数值分析实验

1、数值实验目录数值实验1第一章 实验11.根号5的极限12.求pi的近似值23画图y=sinx 与其泰勒展开34三维图mesh 与surf实现。4第二章 实验61.列主元三角分解A.62.追赶法求方程，电路的电流83.方程组的性态与矩阵条件数的实验94.Wilson矩阵。求特征值，条件数。误差12第三章 实验151.分别用Jacobi, Gauss-Seild,共轭梯度法解方程152.利用共轭梯度法计算矩阵-105阶193.利用cgs,bicg,bicgstab,等计算矩阵解20第五章 实验221.Newton插值f(x)=1/(1+4x2).图形，误差222. f(x)=1/(1+4x2)插值

2、243.飞机的外形轮廓24第九章 四阶龙格库塔方法（自选）251.解算微分方程组252.Matlab的四阶Rk方法253.作图比较，一定误差27第一章 实验1.根号5的极限代码如下: x0=5;for i=0:1:100 xi=sqrt(x0); ei=x0-xi; fprintf('NO:%d,xi=%0.8f, The error is: %0.8fn', i,xi,ei); x0=xi; if ei<10(-8) break; endend可以看出极限是1.2.求pi的近似值 代码如下： sum0=0;for n=1:1:50000 y=(-1)(n+1)/(2*n

3、-1); sum1=y+sum1; pi_1=4*sum0; pi_2=4*sum1; error=pi_2-pi_1; sum0=sum1; fprintf('NO:%d,pi_2=%0.8f, The error: %0.8fn', n,pi_2,error); if abs(error)<10(-4) break; endend求得结果是：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); end f

4、or i=0:5 y2=y2+(-1)(i)*x.(2*i+1)/factorial(2*i+1); end for i=0:10 y3=y3+(-1)(i)*x.(2*i+1)/factorial(2*i+1); end 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重合。4三维图mesh 与surf实现。首先用mesh函数.for i=0:1:2; k=0.2,0.1,0.05; xi, yi=mesh

5、grid(-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(xi,yi,z)end再用surf函数.for i=0:1: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(xi,yi,z)end第二章 实验1

6、.列主元三角分解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 error('不能分解') end u=A;p=eye(m);l=eye(m); for i=1:m for j=i:m t(j)=u(j,i); for k=1:i-1 t(j)=t(j)-u(j,k)*u(k,i); end end a=i;b=abs(t(i); for j=i+1

7、:m if b<abs(t(j) b=abs(t(j); a=j; end end if a=i for j=1:m c=u(i,j); u(i,j)=u(a,j); u(a,j)=c; end for j=1:m c=p(i,j); p(i,j)=p(a,j); p(a,j)=c; end c=t(a); t(a)=t(i); t(i)=c; end u(i,i)=t(i); for j=i+1:m u(j,i)=t(j)/t(i); end for j=i+1:m for k=1:i-1 u(i,j)=u(i,j)-u(i,k)*u(k,j); end end end l=tril(

8、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 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=

9、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);end%-×·µÄ¹ý³Ì-y=zeros(n,1);y(1)=bb(1)/d(1);for i=2:n y(i)=(bb(i)-a(i-1)*y(i-1)/d(i);end%-¸ÏµÄ¹ý³Ì

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

11、（2） 当n=5,10,20 时，解分别如下：说明条件数越大，越病态，发散了。当n=10时就一个值有一定误差。（3） n=10,解方程（A+delta）x=b. clcn=10;A=zeros(n);C=zeros(n);b=zeros(1,n);for i=1:n x(i)=1+0.1*i; for j=1:n A(i,j)=x(i)(j-1); % C(i,j)=1/(i+j-1); b(i)=b(i)+A(i,j); endendA(2,2)=A(2,2)+10(-8)A(10,10)=A(10,10)+10(-8)Ad=cond(A,2)x=Ab'那么结果如下：可以看出很小的扰

12、动导致解误差变得较大了。4.Wilson矩阵。求特征值，条件数。误差（1）行列式，条件数，特征值（2）(A+A0)(x+x0)=b,求此时的x0与|x0|2.代码如下：clcA=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 6.02 5;8.2 5.89 9.96 9.01;6.98 5.04 8.97 9.98b=32 23 33 31A0=A1-Ax=Ab'x1=A1b'x0=x1-xfanshux0=sqrt(x0(1,1)2+x0(2,1)2+x0(3,1)2+x0(4,1)2)fanshux

13、=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*A0_A)/(1-condA*A0_A)结果如下：然后比较：然后利用公式计算:左边=0.8225，右边=1.0358则有0.8225<1.0358。第三章 实验1.分别用Jacobi, Gauss-Seild,共轭梯度法解方程（1）Jacobi

14、迭代法clcA=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('NO: x1 x2 x3 x4 x5n')while 1 for i=1:n y(i)=b(i); for j=1:n if j=i y(i)=y(i)-A(i,j)*x(j); end end if abs(A(i,i)<1e-10 | k=it_max r

15、eturn; end y(i)=y(i)/A(i,i); end if norm(y-x,inf)<ep break; end x=y; k=k+1; fprintf('NO:%d, %f %f %f %f %fn', k,x(1,1),x(2,1),x(3,1),x(4,1),x(5,1)end得的结果为（2）Gauss-Seild迭代法 clcA=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 12x0=0 0 0 0 0'n=length(A); ep=1

16、0(-7);N=100;k=1;x=zeros(n,1); fprintf('n*start*n',k);for k=1:N fprintf('%dt',k); %calculate x(1) x(1)=b(1)/A(1,1); for j=2:n x(1)=x(1)-A(1,j)*x0(j)/A(1,1); end fprintf('%ft',x(1); %calculate x(i) for i=2:n-1 x(i)=b(i)/A(i,i); for j=1:i-1 x(i)=x(i)-A(i,j)*x(j)/A(i,i); end for

17、j=i+1:n x(i)=x(i)-A(i,j)*x0(j)/A(i,i); end fprintf('%ft',x(i); end % calculate x(n) x(n)=b(n)/A(n,n); for j=1:n-1 x(n)=x(n)-A(n,j)*x(j)/A(n,n); end fprintf('%ft',x(n); if norm(x-x0,inf)<ep break; else k=k+1; x0=x; end fprintf('n');end fprintf('n*end*n',k);x0运行结果如下：

18、可以看出Gauss-Seild只要41次就可以计算结果，而Jacobi需要76次。（3）共轭梯度迭代法 clcA=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 12'x0=0 0 0 0 0'N=length(A); eps=10(-7);x=zeros(N,1);%µü´ú½üËÆÏòÁ¿.normÊǾØ

19、3;ó£¬ÏòÁ¿·¶Êýr=b-A*x; %Æäʵ±íʾµÄÊǸºËÑË÷·½Ïò¡£fai=1/2x'AX-bxd=r;for k=0:N-1 fprintf('µÚ%d´Îµ

20、;ü´ú£º',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) break; end B=(norm(rr)2)/(norm(r)2); d=rr+B*d; r=rr;end计算结果如下：理论上只需要N次得到精确值2.利用共轭梯度法计算矩阵-105阶clcn=105;for i=1:n for j=1:n if i=j A(i,j)=3; elseif abs(i-j)=1 A(i,j)=-1; else

21、if (i+j=n+1) && (i=n/2) && (i=n/2+1) A(i,j)=1/2; else A(i,j)=0; end endendb(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=b'N=length(A); eps=10(-7);x=zeros(N,1);%µü´ú½üËÆÏòÁ

22、¿.normÊǾØÕó£¬ÏòÁ¿·¶Êýr=b-A*x; %Æäʵ±íʾµÄÊǸºËÑË÷·½Ïò¡£fai=1/2*x'AX-bxd=r;for k=0:N-1 fpri

23、ntf('µÚ%d´Îµü´ú£º',k+1); A=(norm(r)2)/(d'*A*d); x=x+A*d; x' rr=b-A*x; %rr=r(k+1) if (norm(rr)<=eps)|(k=N-1) break; end B=(norm(rr)2)/(norm(r)2); d=rr+B*d; r=rr;end3.利用cgs,bicg,bicgstab,等计算矩阵解利用cgs（）函数A = gallery('wilk',21)b

24、 = sum(A,2);tol = 1e-12; maxit = 15; M1 = diag(10:-1:1 1 1:10);x = cgs(A,b,tol,maxit,M1)利用bicg（）函数A = gallery('wilk',21)b = sum(A,2);tol = 1e-12; maxit = 15; M1 = diag(10:-1:1 1 1:10);x = bicg(A,b,tol,maxit,M1)利用bicgstab（）函数A = gallery('wilk',21)b = sum(A,2);tol = 1e-12; maxit = 10;

25、M1 = diag(10:-1:1 1 1:10);x = bicgstab(A,b,tol,maxit,M1)第五章 实验1.Newton插值f(x)=1/(1+4x2).图形，误差 （1）输出插值多项式。程序：function p = NewtonChazhi( x,y,n )%UNTITLED3 此处显示有关此函数的摘要% 此处显示详细说明f=zeros(n,n);f(:,1)=y'for j=2:n for i=j:n f(i,j)=(f(i,j-1)-f(i-1,j-1)/(x(i)-x(i-j+1); endendp=f(n,n);for k=n:-1:2 p=conv(p

26、,poly(x(k-1); d=length(p); p(d)=p(d)+f(k-1,k-1);endpdisp('牛顿插值多项式：');polynomial=poly2sym(p,'x')endx=-5:1:5;y=1./(1+4*x.2);n=length(x);p = NewtonChazhi( x,y,n )结果：牛顿插值多项式： polynomial = - (7319042784910035*x10)/147573952589676412928 + (256*x8)/93425 + (5*x7)/1152921504606846976 - (7410620163505401*x6)/144115188075855872 + x5/9007199254740992 + (36624*x4)/93425 + x3/9007199254740992 - (5148893614132309*x2)/4503599627370496