数值分析实验Word文档下载推荐.docx

1、 代码如下： sum0=0;for n=1:50000 y=（-1）（n+1）/（2*n-1）; sum1=y+sum1; pi_1=4*sum0; pi_2=4*sum1; error=pi_2-pi_1; sum0=sum1;%d,pi_2=%0.8f, The error:, n,pi_2,error）; if abs（error）10（-4）求得结果是：3画图y=sinx 与其泰勒展开 x=0:pi/100:2*pi;y=sin（x）;y1=0;y2=0;y3=0; for i=0:2 y1=y1+（-1）（i）*x.（2*i+1）/factorial（2*i+1）;5 y2=y2+（

2、-1）（i）*x.（2*i+1）/factorial（2*i+1）;10 y3=y3+（-1）（i）*x.（2*i+1）/factorial（2*i+1）; plot（x,y,*r,x,y1,b,x,y2,-g,x,y3,k） axis（0 2*pi -1.5 1.5） 看出n=2发散，n=10几乎与y=sinx重合。首先用mesh函数.2; k=0.2,0.1,0.05; xi, yi=meshgrid（-10:k（1,i+1）:10）; z=exp（-abs（xi）+cos（xi+yi）+1./（xi.2+yi.2+1）; subplot（2,2,i+1） title（mesh mesh（

3、xi,yi,z）再用surf函数.第二章 实验1.列主元三角分解A. clc A=1 1 1 1 1;1 2 3 4 5;1 3 6 10 15;1 4 10 20 35;1 5 15 35 70E=eye（5） m,n=size（A）; if m=n error（不是方阵） return end if det（A）=0不能分解 u=A;p=eye（m）;l=eye（m）; for i=1:m for j=i: t（j）=u（j,i）; for k=1:i-1 t（j）=t（j）-u（j,k）*u（k,i）; a=i;b=abs（t（i）; for j=i+1: if babs（t（j） b=

4、abs（t（j）; a=j; if a=i for j=1:m c=u（i,j）; u（i,j）=u（a,j）; u（a,j）=c; c=p（i,j）; p（i,j）=p（a,j）; p（a,j）=c; c=t（a）; t（a）=t（i）; t（i）=c; u（i,i）=t（i）; u（j,i）=t（j）/t（i）;i-1 u（i,j）=u（i,j）-u（i,k）*u（k,j）; l=tril（u,-1）+eye（m） u=triu（u,0） B=A/E2.追赶法求方程，电路的电流A=2 -2 0 0 0 0 0 0;-2 5 -2 0 0 0 0 0;0 -2 5 -2 0 0 0 0;0

5、0 -2 5 -2 0 0 0; 0 0 0 -2 5 2 0 0;0 0 0 0 -2 5 -2 0;0 0 0 0 0 -2 5 -2;0 0 0 0 0 0 -2 5bb=220/27 0 0 0 0 0 0 0n=size（A,1）;s=zeros（n,1）;%-b=diag（A）;a=diag（A,-1）;c=diag（A,1）;d=zeros（n,1）;u=zeros（n-1,1）;for i=1:n-1 d（1）=b（1）; u（i）=c（i）/d（i）; d（i+1）=b（i+1）-a（i）*u（i）;%-y=zeros（n,1）;y（1）=bb（1）/d（1）;for i=2

6、:n y（i）=（bb（i）-a（i-1）*y（i-1）/d（i）;%-s（n）=y（n）;for i=n-1:-1:1 s（i）=y（i）-u（i）*s（i+1）;s3.方程组的性态与矩阵条件数的实验clcn=5;A=zeros（n）;C=zeros（n）;b=zeros（1,n）; x（i）=1+0.1*i; A（i,j）=x（i）（j-1）; C（i,j）=1/（i+j-1）; b（i）=b（i）+A（i,j）;ACbd=cond（A,2）DD=Ab（1） 当n=5,10,20 时，cond（A）= 5.3615e+05, 8.6823e+11, 1.5428e+22。看出越来越病态了。

7、（2） 当n=5,10,20 时，解分别如下：说明条件数越大，越病态，发散了。当n=10时就一个值有一定误差。（3） n=10,解方程（A+delta）x=b. clcn=10; % C（i,j）=1/（i+j-1）;A（2,2）=A（2,2）+10（-8）A（10,10）=A（10,10）+10（-8）x=Ab那么结果如下：可以看出很小的扰动导致解误差变得较大了。误差（1）行列式，条件数，特征值（2）（A+A0）（x+x0）=b,求此时的x0与|x0|2.代码如下：A=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10A1=10 7.2 8.1 6.9;7.08 5.07

8、6.02 5;8.2 5.89 9.96 9.01;6.98 5.04 8.97 9.98b=32 23 33 31A0=A1-Ax1=A1bx0=x1-xfanshux0=sqrt（x0（1,1）2+x0（2,1）2+x0（3,1）2+x0（4,1）2）fanshux=sqrt（1+1+1+1）eigAA=eig（A*A）;fanshuA=sqrt（max（eigAA）eigA0A0=eig（A0*A0）;fanshuA0=sqrt（max（eigA0A0）x0_x=fanshux0/fanshuxA0_A=fanshuA0/fanshuAcondA=cond（A）x0_xx=（condA*

9、A0_A）/（1-condA*A0_A）结果如下：然后比较：然后利用公式计算:左边=0.8225，右边=1.0358则有0.82251.0358。第三章 实验1.分别用Jacobi, Gauss-Seild,共轭梯度法解方程（1）Jacobi迭代法A=10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15b=12 -27 14 -17 12n=length（A）; k=0;ep=10（-7）it_max=100;x=zeros（n,1）; y=zeros（n,1）;fprintf（ x1 x2 x3 x4 x5nwhile 1n

10、 y（i）=b（i）; if j=i y（i）=y（i）-A（i,j）*x（j）; if abs（A（i,i）1e-10 | k=it_max return; y（i）=y（i）/A（i,i）; if norm（y-x,inf）ep x=y; k=k+1;%d, %f %f %f %f %fn, k,x（1,1）,x（2,1）,x（3,1）,x（4,1）,x（5,1）得的结果为（2）Gauss-Seild迭代法x0=0 0 0 0 0;ep=10（-7）;N=100;k=1;n*start*n,k）;for k=1:N%dt %calculate x（1） x（1）=b（1）/A（1,1）;

11、for j=2: x（1）=x（1）-A（1,j）*x0（j）/A（1,1）;%ft,x（1）; %calculate x（i） for i=2: x（i）=b（i）/A（i,i）; x（i）=x（i）-A（i,j）*x（j）/A（i,i）; x（i）=x（i）-A（i,j）*x0（j）/A（i,i）;,x（i）; % calculate x（n） x（n）=b（n）/A（n,n）; x（n）=x（n）-A（n,j）*x（j）/A（n,n）;,x（n）; if norm（x-x0,inf） else x0=x;n）;n*end*nx0运行结果如下：可以看出Gauss-Seild只要41次就可以

12、计算结果，而Jacobi需要76次。（3）共轭梯度迭代法b=12 -27 14 -17 12N=length（A）;eps=10（-7）;x=zeros（N,1）;%.normr=b-A*x; %fai=1/2xAX-bxd=r;for k=0:N-1%d,k+1）; a=（norm（r）2）/（d*A*d） x=x+a*d rr=b-A*x; %rr=r（k+1） if （norm（rr）=eps）|（k=N-1） B=（norm（rr）2）/（norm（r）2）; d=rr+B*d; r=rr;计算结果如下：理论上只需要N次得到精确值2.利用共轭梯度法计算矩阵-105阶n=105; if

13、i=j A（i,j）=3; elseif abs（i-j）=1 A（i,j）=-1; elseif （i+j=n+1） & （i=n/2） & （i=n/2+1） A（i,j）=1/2; A（i,j）=0;b（1,1）=2.5;b（1,2:105/2-1）=1.5;b（1,105/2）=1.0;b（1,105/2+1）=1.0;b（1,105/2+2:105-1）=1.5;b（1,105）=2.5;b=bfai=1/2*x A=（norm（r）2）/（d*A*d）; x=x+A*d; x3.利用cgs,bicg,bicgstab,等计算矩阵解利用cgs（）函数A = gallery（wilk,

14、21）b = sum（A,2）;tol = 1e-12; maxit = 15;M1 = diag（10:1 1 1:10）;x = cgs（A,b,tol,maxit,M1）利用bicg（）函数x = bicg（A,b,tol,maxit,M1）利用bicgstab（）函数 maxit = 10;x = bicgstab（A,b,tol,maxit,M1）第五章 实验1.Newton插值f（x）=1/（1+4x2）.图形，误差 （1）输出插值多项式。程序：function p = NewtonChazhi（ x,y,n ）%UNTITLED3 此处显示有关此函数的摘要% 此处显示详细说明f=

15、zeros（n,n）;f（:,1）=yfor j=2: for i=j: f（i,j）=（f（i,j-1）-f（i-1,j-1）/（x（i）-x（i-j+1）;p=f（n,n）;for k=n: p=conv（p,poly（x（k-1）; d=length（p）; p（d）=p（d）+f（k-1,k-1）;pdisp（牛顿插值多项式：polynomial=poly2sym（p,xx=-5:5;y=1./（1+4*x.2）;n=length（x）;p = NewtonChazhi（ x,y,n ）结果：polynomial =- （7319042784910035*x10）/1475739525

16、89676412928 + （256*x8）/93425 + （5*x7）/1152921504606846976 - （7410620163505401*x6）/144115188075855872 + x5/9007199254740992 + （36624*x4）/93425 + x3/9007199254740992 - （5148893614132309*x2）/4503599627370496 + （19*x）/36028797018963968 + 1p = Columns 1 through 6 -0.0000 0 0.0027 0.0000 -0.0514 0.0000 Co

17、lumns 7 through 111.3920 0.0000 -1.1433 0.0000 1.0000（2）在-5,5均匀插入99个节点，计算这些节点上函数f（x）的近似值，在同一张图上画出原函数与插值函数的图形。 x=-5:X=linspace（-5,5,101）;y1=1./（1+4*X.2）;101 N（i）=polyval（p,X（i）; E（i）=N（i）-y1（i）;plot（X,y1,X,N）（4） 误差图形2. f（x）=1/（1+4x2）插值 3.飞机的外形轮廓 y=0,-30,-36,-35,-28.44,-9.4,0;x2=520:-8.67:0; x=520,280

18、,156.6,78,39.62,3.1,0;y2=spline（x,y,x2）;%用样条插值计算近似值x1=0,3.1,39.62,78,156.6,280,520;y1=0,9.4,28.44,35,36,30,0;x3=0:8.67:520;y3=spline（x1,y1,x3）;plot（x2,y2,r-,x3,y3,g-title（机翼外形曲线第九章 四阶龙格库塔方法（自选）四阶龙格-库塔（Runge-Kutta）方法，因为在本人研究方向上，要计算姿态结算方程组，而经典的方法就是四阶RK方法，故这章自己选题做这个方面的实验。1. 解算微分方程组导航姿态方程如下：角度与角速度关系2. Matlab的四阶Rk方法3. 作图比较，一定误差黑线表示实验测试计算结果，红线绿线表示matlab计算结果。两者误差较大。由于计算中全是三角函数计算，误差大。