数值分析上机作业(总)

..数值分析上机实验一、解线性方程组直接法〔教材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>计算结果：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>当n=5时c=5.3615e+005d=9.4327e-007当n=10时c=8.6823e+011d=5.0894e-013当n=20时c=3.4205e+022d=8.1226e-024备注：对于病态矩阵A来说,d为接近0的数；对于非病态矩阵A来说,d为接近1的数。三、解线性方程组的迭代法〔教材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:迭代次数太多,可能不收敛'>;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>计算结果：x=1.0000-2.00003.0000-2.00001.0000n=67经过67次迭代,得到最终结果〔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:迭代次数太多,可能不收敛'>;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>计算结果：x=1.0000-2.00003.0000-2.00001.0000n=38经过38次迭代,得到最终结果。四、矩阵特征值与特征向量的计算〔教材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>求解结果：l=343.0000v=-0.6667-2.000001.0000s=114结论：经过114次迭代,求得此矩阵的最大特征值为343.0000,及其对应特征向量为[-0.6667;-2.0000;0;1.0000]五、函数逼近〔教材164页16题本题采用最小二乘法进行拟合：线性最小二乘法程序如下：function[a,b]=LZXEC<x,y>if<length<x>==length<y>>n=length<x>;elsedisp<'x和y的维数不相等'>;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>计算结果：b=8.4169e-004a=0.0090绘制图形：fplot<'x/<0.0090*x+8.4169e-004>',[2,20]>gridontitle<'最小二乘拟合'>六、数值微分与数值积分〔教材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个节点：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>;绘制图形：七、非线性方程及非线性方程组的解法本题采用牛顿法进行求解：牛顿法解非线性方程程序如下：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<'迭代步数太多,可能不收敛'>;return;endend本题解决方案如下：首先,绘制此方程的图形,大概确定其与X轴的交点位置。由于,可以得出因此绘制程序如下：ezplot<'log<<513+0.6651*x>/<513-0.6651*x>>-x/<1400*0.0918>',[-772,772,-10,10]>;gridon得到图形如下图所示：经过放大后,发现图形与x轴的交点接近处。计算非零根：令为牛顿法接非线性方程的初值。程序如下：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结论：此方程的两个非零根分别为：八、常微分方程数值解法〔教材266页19题本题分别采用四阶ADAMS预测—校正算法和经典RK法进行求解：四阶ADAMS预测—校正算法如下：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经典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函数程序如下：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','准确解'>整体图形：局部放大图形：继续放大：数据结果：yy1<RK>yy2<ADAMS>yy3<准确解>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