北京科技大學(xué)應(yīng)用計(jì)算方法作業(yè)與答案

1、一、第一次作業(yè)（一）2-6計(jì)算下列向量的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 = 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 計(jì)算下列矩陣的行范數(shù)、列范數(shù)、譜范數(shù)、范數(shù)。（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è)用牛頓迭代法求方程在附近的根。要求：給成程序和運(yùn)行結(jié)果.1、 牛頓法的基本原理在求解非線性方程時(shí)，它的困難在于是非線性函數(shù)，為克服這一困難，考慮它的線性展開。設(shè)當(dāng)前點(diǎn)為，在處的展開式為 令，解其方程得到 式為牛頓迭代公式，用牛頓迭代公式求方程根的方法稱為牛頓迭代法。此即牛頓迭代法的設(shè)計(jì)原理。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('兩端點(diǎn)函數(shù)值乘積大于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、運(yùn)行結(jié)果截圖結(jié)論：通過計(jì)算可以看出在附近的根為1.8794.三、第三次作業(yè)編寫高斯順序消元法求解下面方程組的程序并計(jì)算結(jié)果。 1、 高斯順序消元法的設(shè)計(jì)原理高斯順序消元法的基本思想是

6、將線性方程組 通過消元，逐步轉(zhuǎn)化為等價(jià)的上（或下）三角形方程組，然后用回代法求解。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('對(duì)角元素為0！'); %防止對(duì)角元素為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); %通用的求上三角系數(shù)矩陣線性方程組的函數(shù)XA = A; %消元后的系數(shù)矩陣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ùn)行結(jié)果截圖結(jié)論：高斯順序消元法求解出的結(jié)果為。四、第四次作業(yè)編寫迭代法和迭代法求解方程組的程序，并計(jì)算出結(jié)果。 精度要求：（一） 用求解題設(shè)方程組1、迭代原理設(shè)有一個(gè)元線性方程組 它的矩陣形式為，如果非奇異，且。由

8、上式方程組可以得到 而其相應(yīng)的迭代公式 此迭代公式即為迭代。2、迭代法程序代碼functionx,n=jacobi(A,b,x0,eps,varargin)% 求解線性方程組的迭代法其中% A為方程組的系數(shù)矩陣% b為方程組的右端項(xiàng)% eps為精度要求，默認(rèn)值為1e-5% varargin為最大迭代次數(shù)，值100% x為方程組的解% n為迭代次數(shù)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:迭代次數(shù)太多，可能不收斂！'); return; endend3、運(yùn)行結(jié)果截圖結(jié)論：式方程組的解為，需要迭代12步。（二） 用求解題設(shè)方程組1、迭代原理在迭代計(jì)算過程中可看出，方法計(jì)算時(shí)，并未用到已算出的，這時(shí)想到，如果迭代收斂，比更接近方程組的解，若能在迭代過程中盡快用新的信息去替換，則可望收斂更快。由此，可將迭代公式改寫為 此式即為迭代法。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:迭代次數(shù)太多，可能不收斂！'); return;

11、 end end3、運(yùn)行結(jié)果截圖結(jié)論：式方程組的解為，需要迭代5步。五、第五次作業(yè)用歸一化算法（歸一化冪法6-1）求矩陣A的最大模特征值和特征向量，其中1、 冪法的基本思想任取一個(gè)非零初始向量且，由矩陣的乘冪構(gòu)造一迭代序列為 假設(shè)矩陣有個(gè)線性無關(guān)的特征向量，于是給定的初始向量可以用這組特征向量線性表示，即并設(shè)。把代入迭代序列的第一條式子，得同理可得2、程序代碼function l,v,s=pmethod(A,x0,eps)if nargin=2 eps = 1.0e-6;endv = x0; %v為主特征向量M = 5000; %迭代步數(shù)限制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為迭代步數(shù) return; else if(k=M) disp('迭代步數(shù)太多，收斂速度太慢！'); l = m; s = M; else l = m; end endend3、運(yùn)行結(jié)果截圖結(jié)論：由輸出結(jié)果可知，經(jīng)過15步迭代，求得矩陣特征值為9.3597，對(duì)應(yīng)的迭代向量為（1.0000,0.9785，0.6222）T六、第六次作業(yè)已知的函數(shù)值和導(dǎo)數(shù)值如下：求次數(shù)小于等于4的多項(xiàng)式，使并給出余

13、項(xiàng)公式。解：首先構(gòu)造拉格朗日插值函數(shù)化簡得構(gòu)造函數(shù)對(duì)上式求導(dǎo)得根據(jù)導(dǎo)數(shù)信息可以解出所以余項(xiàng)構(gòu)造輔助函數(shù)其中然后反復(fù)利用定理最好可以得到七、第七次作業(yè)已知某地區(qū)在不同月份的平均日照時(shí)間的觀測(cè)數(shù)據(jù)如下表，試分析日照時(shí)間的變化規(guī)律。表1：日照時(shí)間月份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通過圖像可以看出日照隨時(shí)間的大概走勢(shì)，我們也可以推測(cè)出以后每個(gè)月的大概日照長度的大概走勢(shì)。八、第八次作業(yè)在某化學(xué)反應(yīng)中，由實(shí)驗(yàn)得分解物濃度與時(shí)間關(guān)系如下：表2：分解物濃度與時(shí)間關(guān)系時(shí)間0510152025濃度01.272.162.863.443.87時(shí)間303540455055濃度4.154.374.514.584.624.64用最小二乘法求。解：利用matlab數(shù)值逼近可以得到：即最小二乘擬合二次多項(xiàng)式為其中擬合值

15、和實(shí)際值的對(duì)比圖見下:圖2：圖中圓圈為已知數(shù)據(jù)的點(diǎn)，曲線為最小二乘擬合二次多項(xiàng)式九、第九次作業(yè)用龍貝格求積算法計(jì)算積分。計(jì)算過程中數(shù)值保留6位有效數(shù)字。1、 龍貝格算法原理龍貝格求積法的算法為其中逼近的階達(dá)到次。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*(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);%取絕對(duì)值。endquad=R(J+1,J+1);再定義函數(shù)function y=f(x);y=exp(x);%exp(x

17、)以e為底的指數(shù)。命令窗口中輸入clcclearromber('f', 0, 1, 7, 10(-6) 3、 運(yùn)行結(jié)果截圖結(jié)論：求解得最終結(jié)果為1.7183。十、第十次作業(yè)用經(jīng)典R-K方法計(jì)算初值問題，取步長，計(jì)算的近似值 1、 經(jīng)典R-K方法設(shè)計(jì)原理在區(qū)間上取在個(gè)點(diǎn)處的數(shù)值作線性組合，即得到P級(jí)的R-K方法，在P級(jí)R-K方法中，最常用的是標(biāo)準(zhǔn)的（或經(jīng)典的）四階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:一階常微分方程的一般表達(dá)式的右端函數(shù)%h:積分步長%a:自變量取值下限%b:自變量取值上限%y0:函數(shù)初值%varvee：常微分方程的變量組format long;N=(b-a)/h;y=zeros(N+1,1);y