北京科技大學應用計算方法作業(yè)與答案

1、一、第一次作業(yè)（一）2-6計算下列向量的1-范數(shù)、-范數(shù)、2-范數(shù)。（1）x=(12,-4,-6,2)T A=12,-4,-6,2A = 12 -4 -6 2 norm(A,1)ans = 24 norm(A,inf)ans = 12 norm(A,2)ans = (2) x=(1,3,-4)T、 A=1,3,-4A = 1 3 -4 norm(A,1)ans = 8 norm(A,inf)ans = 4 norm(A,2)ans = （二）2-9 計算下列矩陣的行范數(shù)、列范數(shù)、譜范數(shù)、范數(shù)。（1） A=3,-1,1;1,1,1;2,1,-1A = 3 -1 1 1 1 1 2 1 -1 no

2、rm(A,1)ans = 6 norm(A,inf)ans = 5 norm(A,2)ans = norm(A,fro)ans = （2）. A=0,1;-1,0A = 0 1 -1 0 norm(A,1)ans = 1 norm(A,inf)ans = 1 norm(A,2)ans = 1 norm(A,fro)ans = 二、第二次作業(yè)用牛頓迭代法求方程在附近的根。要求：給成程序和運行結果.1、 牛頓法的基本原理在求解非線性方程時，它的困難在于是非線性函數(shù)，為克服這一困難，考慮它的線性展開。設當前點為，在處的展開式為 令，解其方程得到 式為牛頓迭代公式，用牛頓迭代公式求方程根的方法稱為牛頓

3、迭代法。此即牛頓迭代法的設計原理。2、程序代碼function root=NewtonRoot(f,a,b,eps)if(nargin=3) eps=;endf1=subs(sym(f),findsym(sym(f),a);f2=subs(sym(f),findsym(sym(f),b);if(f1=0) root=a;endif(f2=0) roor=b;endif(f1*f20) disp(兩端點函數(shù)值乘積大于0！) return;else tol=1; fun=diff(sym(f); fa=subs(sym(f),findsym(sym(f),a); fb=subs(sym(f),fi

4、ndsym(sym(f),b); dfa=subs(sym(fun),findsym(sym(fun),a); dfb=subs(sym(fun),findsym(sym(fun),b); if(dfadfb) root=a-fa/dfa; else root=b-fb/dfb; end while(toleps) r1=root; fx=subs(sym(f),findsym(sym(f),r1); dfx=subs(sym(fun),findsym(sym(fun),r1); root=r1-fx/dfx; tol=abs(root-r1); endend3、運行結果截圖結論：通過計算可以

5、看出在附近的根為.三、第三次作業(yè)編寫高斯順序消元法求解下面方程組的程序并計算結果。 1、 高斯順序消元法的設計原理高斯順序消元法的基本思想是將線性方程組 通過消元，逐步轉化為等價的上（或下）三角形方程組，然后用回代法求解。2、 程序代碼function x,XA=GaussXQByOrder(A,b)%高斯順序消元法N = size(A);n = N(1);for i=1:(n-1) for j=(i+1):n if(A(i,i)=0) disp(對角元素為0！); %防止對角元素為0 return; end l = A(j,i); m = A(i,i); A(j,1:n)=A(j,1:n)-

6、l*A(i,1:n)/m; %消元方程 b(j)=b(j)-l*b(i)/m; endend x=SolveUpTriangle(A,b); %通用的求上三角系數(shù)矩陣線性方程組的函數(shù)XA = A; %消元后的系數(shù)矩陣function x=SolveUpTriangle(A,b)N=size(A);n=N(1);for i=n:-1:1 if(in) s=A(i,(i+1):n)*x(i+1):n,1); else s=0; end x(i,1)=(b(i)-s)/A(i,i);end 3、 運行結果截圖結論：高斯順序消元法求解出的結果為。四、第四次作業(yè)編寫迭代法和迭代法求解方程組的程序，并計算

7、出結果。 精度要求：（一） 用求解題設方程組1、迭代原理設有一個元線性方程組 它的矩陣形式為，如果非奇異，且。由上式方程組可以得到 而其相應的迭代公式 此迭代公式即為迭代。2、迭代法程序代碼functionx,n=jacobi(A,b,x0,eps,varargin)% 求解線性方程組的迭代法其中% A為方程組的系數(shù)矩陣% b為方程組的右端項% eps為精度要求，默認值為1e-5% varargin為最大迭代次數(shù)，值100% x為方程組的解% n為迭代次數(shù)if nargin=3 eps=; M=200;elseif nargin=eps x0=x; x=B*x0+f; n=n+1; if(n=

8、M) disp(Warning:迭代次數(shù)太多，可能不收斂！); return; endend3、運行結果截圖結論：式方程組的解為，需要迭代12步。（二） 用求解題設方程組1、迭代原理在迭代計算過程中可看出，方法計算時，并未用到已算出的，這時想到，如果迭代收斂，比更接近方程組的解，若能在迭代過程中盡快用新的信息去替換，則可望收斂更快。由此，可將迭代公式改寫為 此式即為迭代法。2、迭代法程序代碼functionx,n=gauseidel(A,b,x0,eps,M)if nargin=3 eps=; M=200;elseif nargin=4 M=200;elseif nargin=eps x0=x

9、; x=G*x0+f; n=n+1; if(n=M) disp(Warning:迭代次數(shù)太多，可能不收斂！); return; end end3、運行結果截圖結論：式方程組的解為，需要迭代5步。五、第五次作業(yè)用歸一化算法（歸一化冪法6-1）求矩陣A的最大模特征值和特征向量，其中1、 冪法的基本思想任取一個非零初始向量且，由矩陣的乘冪構造一迭代序列為 假設矩陣有個線性無關的特征向量，于是給定的初始向量可以用這組特征向量線性表示，即并設。把代入迭代序列的第一條式子，得同理可得2、程序代碼function l,v,s=pmethod(A,x0,eps)if nargin=2 eps = ;endv

10、= x0; %v為主特征向量M = 5000; %迭代步數(shù)限制m = 0; l = 0;for(k=1:M) y = A*v; m = max(y); %m為按模最大的分量 v = y/m; if(abs(m - l) x=1:12; y= ; xx=:; yy=spline(x,y,xx); plot(x,y,o,xx,yy)圖1通過圖像可以看出日照隨時間的大概走勢，我們也可以推測出以后每個月的大概日照長度的大概走勢。八、第八次作業(yè)在某化學反應中，由實驗得分解物濃度與時間關系如下：表2：分解物濃度與時間關系時間0510152025濃度0時間303540455055濃度用最小二乘法求。解：利用

11、matlab數(shù)值逼近可以得到：即最小二乘擬合二次多項式為其中擬合值和實際值的對比圖見下:圖2：圖中圓圈為已知數(shù)據(jù)的點，曲線為最小二乘擬合二次多項式九、第九次作業(yè)用龍貝格求積算法計算積分。計算過程中數(shù)值保留6位有效數(shù)字。1、 龍貝格算法原理龍貝格求積法的算法為其中逼近的階達到次。2、 程序代碼functionR,quad,err,h=romber(f,a,b,n,tol)% f是被積函數(shù)。% a, b分別為是積分的上下限。% n1是T數(shù)表的列數(shù)。% tol是允許誤差。% R是T數(shù)表。% quad是所求積分值。M=1;h=b-a;err=1;J=0;R=zeros(4,4);R(1,1)=h*(f

12、eval(f,a)+feval(f,b)/2;while(errtol)&(Jn)|(J4) J=J+1; h=h/2; s=0; for p=1:M x=a+h*(2*p-1); s=s+feval(f,x); end R(J+1,1)=R(J,1)/2+h*s; M=2*M; for K=1:J R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K)/(4K-1); end err=abs(R(J,J)-R(J+1,K+1);%取絕對值。endquad=R(J+1,J+1);再定義函數(shù)function y=f(x);y=exp(x);%exp(x)以e為底的指數(shù)。命令窗口

13、中輸入clcclearromber(f, 0, 1, 7, 10(-6) 3、 運行結果截圖結論：求解得最終結果為。十、第十次作業(yè)用經(jīng)典R-K方法計算初值問題，取步長，計算的近似值 1、 經(jīng)典R-K方法設計原理在區(qū)間上取在個點處的數(shù)值作線性組合，即得到P級的R-K方法，在P級R-K方法中，最常用的是標準的（或經(jīng)典的）四階R-K方法。 2、 程序代碼function fv = Funval(f,varvec,varval)var = findsym(f);varc = findsym(varvec);s1 = length(var);s2 = length(varc);m =floor(s1-1

14、)/3+1);varv = zeros(1,m);if s1 = s2for i=0: (s1-1)/3)k = findstr(varc,var(3*i+1);index = (k-1)/3;varv(i+1) = varval(index+1);endfv = subs(f,var,varv);elsefv = subs(f,varvec,varval);endfunction y=DELGK4_lungkuta(f, h,a,b,y0,varvec)%f:一階常微分方程的一般表達式的右端函數(shù)%h:積分步長%a:自變量取值下限%b:自變量取值上限%y0:函數(shù)初值%varvee：常微分方程的變量組format long;N=(b-a)/h;y=zeros(N+1,1);y(1)=y0;x=a:h:b;var=findsym(f);for i=2:N+1 K1=Funval(f,va