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

1、一、第一次作業(yè)（一）2-6計算下列向量的1-范數、-范數、2-范數。（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 = 14.1421(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 = 5.0990（二）2

2、-9 計算下列矩陣的行范數、列范數、譜范數、范數。（1）>> A=3,-1,1;1,1,1;2,1,-1A = 3 -1 1 1 1 1 2 1 -1>> norm(A,1)ans = 6>> norm(A,inf)ans = 5>> norm(A,2)ans = 3.7888>> norm(A,'fro')ans = 4.4721（2）>> A=0,1;-1,0A = 0 1 -1 0>> norm(A,1)ans = 1>> norm(A,inf)ans = 1>>

3、norm(A,2)ans = 1>> norm(A,'fro')ans = 1.414215二、第二次作業(yè)用牛頓迭代法求方程在附近的根。要求：給成程序和運行結果.1、 牛頓法的基本原理在求解非線性方程時，它的困難在于是非線性函數，為克服這一困難，考慮它的線性展開。設當前點為，在處的展開式為 令，解其方程得到 式為牛頓迭代公式，用牛頓迭代公式求方程根的方法稱為牛頓迭代法。此即牛頓迭代法的設計原理。2、程序代碼function root=NewtonRoot(f,a,b,eps)if(nargin=3) eps=1.0e-4;endf1=subs(sym(f),find

4、sym(sym(f),a);f2=subs(sym(f),findsym(sym(f),b);if(f1=0) root=a;endif(f2=0) roor=b;endif(f1*f2>0) disp('兩端點函數值乘積大于0！') return;else tol=1; fun=diff(sym(f); fa=subs(sym(f),findsym(sym(f),a); fb=subs(sym(f),findsym(sym(f),b); dfa=subs(sym(fun),findsym(sym(fun),a); dfb=subs(sym(fun),findsym(sy

5、m(fun),b); if(dfa>dfb) root=a-fa/dfa; else root=b-fb/dfb; end while(tol>eps) 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、運行結果截圖結論：通過計算可以看出在附近的根為1.8794.三、第三次作業(yè)編寫高斯順序消元法求解下面方程組的程序并計算結果。 1、 高斯順序消元法的設計原理高斯順序消元法的基本思想是

6、將線性方程組 通過消元，逐步轉化為等價的上（或下）三角形方程組，然后用回代法求解。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)-l*A(i,1:n)/m; %消元方程 b(j)=b(j)-l*b(i)/m; endend x=SolveUpTria

7、ngle(A,b); %通用的求上三角系數矩陣線性方程組的函數XA = A; %消元后的系數矩陣function x=SolveUpTriangle(A,b)N=size(A);n=N(1);for i=n:-1:1 if(i<n) 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è)編寫迭代法和迭代法求解方程組的程序，并計算出結果。 精度要求：（一） 用求解題設方程組1、迭代原理設有一個元線性方程組 它的矩陣形式為，如果非奇異，且。由

8、上式方程組可以得到 而其相應的迭代公式 此迭代公式即為迭代。2、迭代法程序代碼functionx,n=jacobi(A,b,x0,eps,varargin)% 求解線性方程組的迭代法其中% A為方程組的系數矩陣% b為方程組的右端項% eps為精度要求，默認值為1e-5% varargin為最大迭代次數，值100% x為方程組的解% n為迭代次數if nargin=3 eps=1.0e-6; M=200;elseif nargin<3 error returnelseif nargin=5 M=varargin1;endD=diag(diag(A);L=-tril(A,-1);U=-tr

9、iu(A,1);B=D(L+U);f=Db;x=B*x0+f;n=1;while norm(x-x0)>=eps x0=x; x=B*x0+f; n=n+1; if(n>=M) disp('Warning:迭代次數太多，可能不收斂！'); return; endend3、運行結果截圖結論：式方程組的解為，需要迭代12步。（二） 用求解題設方程組1、迭代原理在迭代計算過程中可看出，方法計算時，并未用到已算出的，這時想到，如果迭代收斂，比更接近方程組的解，若能在迭代過程中盡快用新的信息去替換，則可望收斂更快。由此，可將迭代公式改寫為 此式即為迭代法。2、迭代法程序代碼f

10、unctionx,n=gauseidel(A,b,x0,eps,M)if nargin=3 eps=1.0e-3; M=200;elseif nargin=4 M=200;elseif nargin<3 error return;endD=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);G=(D-L)U;f=(D-L)b;x=G*x0+f;n=1;while norm(x-x0)>=eps x0=x; x=G*x0+f; n=n+1; if(n>=M) disp('Warning:迭代次數太多，可能不收斂！'); return;

11、 end end3、運行結果截圖結論：式方程組的解為，需要迭代5步。五、第五次作業(yè)用歸一化算法（歸一化冪法6-1）求矩陣A的最大模特征值和特征向量，其中1、 冪法的基本思想任取一個非零初始向量且，由矩陣的乘冪構造一迭代序列為 假設矩陣有個線性無關的特征向量，于是給定的初始向量可以用這組特征向量線性表示，即并設。把代入迭代序列的第一條式子，得同理可得2、程序代碼function l,v,s=pmethod(A,x0,eps)if nargin=2 eps = 1.0e-6;endv = x0; %v為主特征向量M = 5000; %迭代步數限制m = 0; l = 0;for(k=1:M) y

12、= A*v; m = max(y); %m為按模最大的分量 v = y/m; if(abs(m - l)<eps) l = m; %到所需精度，退出，l為主特征值 s = k; %s為迭代步數 return; else if(k=M) disp('迭代步數太多，收斂速度太慢！'); l = m; s = M; else l = m; end endend3、運行結果截圖結論：由輸出結果可知，經過15步迭代，求得矩陣特征值為9.3597，對應的迭代向量為（1.0000,0.9785，0.6222）T六、第六次作業(yè)已知的函數值和導數值如下：求次數小于等于4的多項式，使并給出余

13、項公式。解：首先構造拉格朗日插值函數化簡得構造函數對上式求導得根據導數信息可以解出所以余項構造輔助函數其中然后反復利用定理最好可以得到七、第七次作業(yè)已知某地區(qū)在不同月份的平均日照時間的觀測數據如下表，試分析日照時間的變化規(guī)律。表1：日照時間月份123456日照h/月80.967.267.150.532.033.6月份789101112日照h/月36.646.852.362.064.171.2解：首先在matlab中輸入語句>> x=1:12;>> y=80.9 67.2 67.1 50.5 32.0 33.6 36.6 46.8 52.3 62.0 64.1 71.2;

14、>> xx=0.0:0.1:12.0;>> yy=spline(x,y,xx);>> plot(x,y,'o',xx,yy)圖1通過圖像可以看出日照隨時間的大概走勢，我們也可以推測出以后每個月的大概日照長度的大概走勢。八、第八次作業(yè)在某化學反應中，由實驗得分解物濃度與時間關系如下：表2：分解物濃度與時間關系時間0510152025濃度01.272.162.863.443.87時間303540455055濃度4.154.374.514.584.624.64用最小二乘法求。解：利用matlab數值逼近可以得到：即最小二乘擬合二次多項式為其中擬合值

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

16、#39;,b)/2;while(err>tol)&(J<n)|(J<4) 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);再定義函數function y=f(x);y=exp(x);%exp(x

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

18、 = length(varc);m =floor(s1-1)/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:一階常微分方程的一般表達式的右端函數%h:積分步長%a:自變量取值下限%b:自變量取值上限%y0:函數初值%varvee：常微分方程的變量組format long;N=(b-a)/h;y=zeros(N+1,1);y