數(shù)值分析上機(jī)作業(yè)(總)

..數(shù)值分析上機(jī)實(shí)驗(yàn)一、解線性方程組直接法〔教材49頁14題追趕法程序如下：functionx=followup<A,b>n=rank<A>;for<i=1:n>if<A<i,i>==0>disp<'Error:對角有元素為0'>;return;endend;d=ones<n,1>;a=ones<n-1,1>;c=ones<n-1>;for<i=1:n-1>a<i,1>=A<i+1,i>;c<i,1>=A<i,i+1>;d<i,1>=A<i,i>;endd<n,1>=A<n,n>;for<i=2:n>d<i,1>=d<i,1>-<a<i-1,1>/d<i-1,1>>*c<i-1,1>;b<i,1>=b<i,1>-<a<i-1,1>/d<i-1,1>>*b<i-1,1>;endx<n,1>=b<n,1>/d<n,1>;for<i=<n-1>:-1:1>x<i,1>=<b<i,1>-c<i,1>*x<i+1,1>>/d<i,1>;end主程序如下：functionzhunganfaA=[2-2000000;-25-200000;0-25-20000;00-25-2000;000-25-200;0000-25-20;00000-25-2;000000-25];b=[220/27;0;0;0;0;0;0;0];x=followup<A,b>計算結(jié)果：x=8.14784.07372.03651.01750.50730.25060.11940.0477二、解線性方程組直接法〔教材49頁15題程序如下：functiontiaojianshu<n>A=zeros<n>;forj=1:1:nfori=1:1:nA<i,j>=<1+0.1*i>^<j-1>;endendc=cond<A>d=rcond<A>當(dāng)n=5時c=5.3615e+005d=9.4327e-007當(dāng)n=10時c=8.6823e+011d=5.0894e-013當(dāng)n=20時c=3.4205e+022d=8.1226e-024備注：對于病態(tài)矩陣A來說,d為接近0的數(shù)；對于非病態(tài)矩陣A來說,d為接近1的數(shù)。三、解線性方程組的迭代法〔教材74頁14題〔1用Jacobi迭代法求：Jacobi迭代法程序如下：function[x,n]=jacobi<A,b,x0,eps,varargin>ifnargin==3eps=1.0e-6;M=200;elseifnargin<3errorreturnelseifnargin==5M=varargin{1};endD=diag<diag<A>>;L=-tril<A,-1>;U=-triu<A,1>;B=D\<L+U>;f=D\b;x=B*x0+f;n=1;whilenorm<x-x0>>=epsx0=x;x=B*x0+f;n=n+1;if<n>=M>disp<'Warning:迭代次數(shù)太多,可能不收斂'>;return;endend本題主程序如下：functionyakebidiedaiA=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];b=[12;-27;14;-17;12];x0=[0;0;0;0;0];[x,n]=jacobi<A,b,x0>計算結(jié)果：x=1.0000-2.00003.0000-2.00001.0000n=67經(jīng)過67次迭代,得到最終結(jié)果〔2用Gauss-Seidel迭代法求：Gauss-Seidel迭代法程序如下：function[x,n]=gauseidel<A,b,x0,eps,M>ifnargin==3eps=1.0e-6;M=200;elseifnargin==4M=200;elseifnargin<3errorreturn;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;whilenorm<x-x0>>=epsx0=x;x=G*x0+f;n=n+1;if<n>=M>disp<'Warning:迭代次數(shù)太多,可能不收斂'>;return;endend本題主程序如下：functiongaosidiedaiA=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];b=[12;-27;14;-17;12];x0=[0;0;0;0;0];[x,n]=gauseidel<A,b,x0>計算結(jié)果：x=1.0000-2.00003.0000-2.00001.0000n=38經(jīng)過38次迭代,得到最終結(jié)果。四、矩陣特征值與特征向量的計算〔教材100頁13題冪法求最大特征值的程序：function[l,v,s]=pmethod<A,x0,eps>ifnargin==2eps=1.0e-6;endv=x0;M=5000;m=0;l=0;for<k=1:M>y=A*v;m=max<y>;v=y/m;if<abs<m-l><eps>l=m;s=k;return;elseif<k==M>disp<'收斂速度過慢'>;l=m;s=M;elsel=m;endendend求解本題程序如下：functionmifaA=[19066-8430;6630342-36;336-168147-112;30-3628291];x0=[0001]';[l,v,s]=pmethod<A,x0>求解結(jié)果：l=343.0000v=-0.6667-2.000001.0000s=114結(jié)論：經(jīng)過114次迭代,求得此矩陣的最大特征值為343.0000,及其對應(yīng)特征向量為[-0.6667;-2.0000;0;1.0000]五、函數(shù)逼近〔教材164頁16題本題采用最小二乘法進(jìn)行擬合：線性最小二乘法程序如下：function[a,b]=LZXEC<x,y>if<length<x>==length<y>>n=length<x>;elsedisp<'x和y的維數(shù)不相等'>;return;endA=zeros<2,2>;A<2,2>=n;B=zeros<2,1>;fori=1:nA<1,1>=A<1,1>+x<i>*x<i>;A<1,2>=A<1,2>+x<i>;B<1,1>=B<1,1>+x<i>*y<i>;B<2,1>=B<2,1>+y<i>;endA<2,1>=A<1,2>;s=A\B;a=s<1>;b=s<2>;首先利用1/y代替y,1/x代替x并采用線性最小二乘法求出a與b：functionzuixiaoerchengx1=[2356791011121416171920];y1=[106.42108.26109.58109.50109.86110.00109.93110.59110.60110.72110.90110.76111.10111.30];x2=1./x1;y2=1./y1;[b,a]=LZXEC<x2,y2>計算結(jié)果：b=8.4169e-004a=0.0090繪制圖形：fplot<'x/<0.0090*x+8.4169e-004>',[2,20]>gridontitle<'最小二乘擬合'>六、數(shù)值微分與數(shù)值積分〔教材207頁26題本題采用高斯—勒讓德求積公式求解：高斯—勒讓德求積公式程序如下：functionI=IntGauss<f,a,b,n,AK,XK>if<n<5&&nargin==4>AK=0;XK=0;elseXK1=<<b-a>/2>*XK+<<a+b>/2>;I=<<b-a>/2>*sum<AK.*subs<sym<f>,findsym<f>,XK1>>;Endta=<b-a>/2;tb=<a+b>/2;switchncase0,I=2*ta*subs<sym<f>,findsym<sym<f>>,tb>;case1,I=ta*<subs<sym<f>,findsym<sym<f>>,ta*0.5773503+tb>+...subs<sym<f>,findsym<sym<f>>,-ta*0.5773503+tb>>;case2,I=ta*<0.55555556*subs<sym<f>,findsym<sym<f>>,ta*0.7745967+tb>+...0.55555556*subs<sym<f>,findsym<sym<f>>,-ta*0.7745967+tb>+...0.88888889*subs<sym<f>,findsym<sym<f>>,tb>>;case3,I=ta*<0.3478548*subs<sym<f>,findsym<sym<f>>,ta*0.8611363+tb>+...0.3478548*subs<sym<f>,findsym<sym<f>>,-ta*0.8611363+tb>+...0.6521452*subs<sym<f>,findsym<sym<f>>,ta*0.3398810+tb>...+0.6521452*subs<sym<f>,findsym<sym<f>>,-ta*0.3398810+tb>>;case4,I=ta*<0.2369269*subs<sym<f>,findsym<sym<f>>,ta*0.9061793+tb>+...0.2369269*subs<sym<f>,findsym<sym<f>>,-ta*0.9061793+tb>+...0.4786287*subs<sym<f>,findsym<sym<f>>,ta*0.5384693+tb>...+0.4786287*subs<sym<f>,findsym<sym<f>>,-ta*0.5384693+tb>+...0.5688889*subs<sym<f>,findsym<sym<f>>,tb>>;end本題計算程序如下〔采用4個節(jié)點(diǎn)：functionshuzhijifena=1;b=1;fors=-5:0.05:5q1<1,a>=IntGauss<'cos<1/2*x^2>',0,s,4>;q2<1,b>=IntGauss<'sin<1/2*x^2>',0,s,4>;a=a+1;b=b+1;endplot<q1,q2>;繪制圖形：七、非線性方程及非線性方程組的解法本題采用牛頓法進(jìn)行求解：牛頓法解非線性方程程序如下：function[r,n]=mulNewton<F,x0,eps>ifnargin==2eps=1.0e-4;endx0=transpose<x0>;Fx=subs<F,findsym<F>,x0>;var=findsym<F>;dF=Jacobian<F,var>;dFx=subs<dF,findsym<dF>,x0>;r=x0-inv<dFx>*Fx;n=1;tol=1;whiletol>epsx0=r;Fx=subs<F,findsym<F>,x0>;dFx=subs<dF,findsym<dF>,x0>;r=x0-inv<dFx>*Fx;tol=norm<r-x0>;n=n+1;if<n>100000>disp<'迭代步數(shù)太多,可能不收斂'>;return;endend本題解決方案如下：首先,繪制此方程的圖形,大概確定其與X軸的交點(diǎn)位置。由于,可以得出因此繪制程序如下：ezplot<'log<<513+0.6651*x>/<513-0.6651*x>>-x/<1400*0.0918>',[-772,772,-10,10]>;gridon得到圖形如下圖所示：經(jīng)過放大后,發(fā)現(xiàn)圖形與x軸的交點(diǎn)接近處。計算非零根：令為牛頓法接非線性方程的初值。程序如下：symsxf=log<<513+0.6651*x>/<513-0.6651*x>>-x/<1400*0.0918>;x0=-765[r,n]=mulNewton<f,x0>解得：r=-767.3861x0=765[r,n]=mulNewton<f,x0>解得：r=767.3861結(jié)論：此方程的兩個非零根分別為：八、常微分方程數(shù)值解法〔教材266頁19題本題分別采用四階ADAMS預(yù)測—校正算法和經(jīng)典RK法進(jìn)行求解：四階ADAMS預(yù)測—校正算法如下：functiony=DEYCJZ_yds<f,h,a,b,y0,varvec>formatlong;N=<b-a>/h;y=zeros<N+1,1>;x=a:h:b;y<1>=y0;y<2>=y0+h*Funval<f,varvec,[x<1>y<1>]>;y<3>=y<2>+h*Funval<f,varvec,[x<2>y<2>]>;y<4>=y<3>+h*Funval<f,varvec,[x<3>y<3>]>;fori=5:N+1v1=Funval<f,varvec,[x<i-4>y<i-4>]>;v2=Funval<f,varvec,[x<i-3>y<i-3>]>;v3=Funval<f,varvec,[x<i-2>y<i-2>]>;v4=Funval<f,varvec,[x<i-1>y<i-1>]>;t=y<i-1>+h*<55*v4-59*v3+37*v2-9*v1>/24;ft=Funval<f,varvec,[x<i>t]>;y<i>=y<i-1>+h*<9*ft+19*v4-5*v3+v2>/24;end經(jīng)典RK算法程序如下：functiony=DELGKT4_lungkuta<f,h,a,b,y0,varvec>formatlong;N=<b-a>/h;y=zeros<N+1,1>;y<1>=y0;x=a:h:b;var=findsym<f>;fori=2:N+1K1=Funval<f,varvec,[x<i-1>y<i-1>]>;K2=Funval<f,varvec,[x<i-1>+h/2y<i-1>+K1*h/2]>;K3=Funval<f,varvec,[x<i-1>+h/2y<i-1>+K2*h/2]>;K4=Funval<f,varvec,[x<i-1>+hy<i-1>+h*K3]>;y<i>=y<i-1>+h*<K1+2*K2+2*K3+K4>/6;end其中FUNVAL函數(shù)程序如下：functionfv=Funval<f,varvec,varval>var=findsym<f>;iflength<var><4ifvar<1>==varvec<1>fv=subs<f,varvec<1>,varval<1>>;elsefv=subs<f,varvec<2>,varval<2>>;endelsefv=subs<f,varvec,varval>;end本題解決方案如下〔步長h=0.1：程序：functionchangweifensymsxy;z=-y+2*cos<x>;yy1=DELGKT4_lungkuta<z,0.1,0,pi,1,[xy]>;yy2=DEYCJZ_yds<z,0.1,0,pi,1,[xy]>;a=0:0.1:pi;yy3=cos<a>+sin<a>;plot<a,yy1,'r*'>;gridon;holdon;plot<a,yy2,'b+'>;plot<a,yy3,'-go'>;legend<'RK','Adams','準(zhǔn)確解'>整體圖形：局部放大圖形：繼續(xù)放大：數(shù)據(jù)結(jié)果：yy1<RK>yy2<ADAMS>yy3<準(zhǔn)確解>yy1-yy3yy2-yy3111001.0948374641.11.094837582-1.18389E-070.0051624181.1787356781.1890008331.178735909-2.30293E-070.0102649241.2508563611.2661140651.250856696-3.351E-070.015257371.3104789041.3242156421.310479336-4.32217E-070.0137363051.3570075791.3694466421.3570081-5.21089E-070.0124385421.3899774871.4012422871.389978088-6.012E-070.0112641991.4090592021.4192520421.409059875-6.72089E-070.0101921671.4140620671.4232851431.4140628-7.33353E-070.0092223431.4049360931.413281631.404936878-7.84658E-070.0083447531.3817724651.3893238971.381773291-8.25739E-070.0075506071.3448026251.3516354611.344803481-8.56415E-070.006831981.2943959641.3005785271.29439684-8.76584E-070.0061816861.2310561281.2366502381.231057014-8.86229E-070.0055932241.1554159871.1604775841.155416873-8.85423E-070.0050607111.0682313141.0728110131.068232188-8.74325E-070.0045788250.9703732280.9745168330.970374081-8.53183E-070.0041427520.8628194940.8665684520.862820316-8.2233