數(shù)值分析實驗

數(shù)值實驗?zāi)夸洈?shù)值實驗 1第一章實驗 11.根號5的極限 12.求pi的近似值 23．畫圖y=sinx與其泰勒展開 34．三維圖mesh與surf實現(xiàn)。 4第二章實驗 61.列主元三角分解A. 62.追趕法求方程，電路的電流 83.方程組的性態(tài)與矩陣條件數(shù)的實驗 94.Wilson矩陣。求特征值，條件數(shù)。誤差 12第三章實驗 151.分別用Jacobi,Gauss-Seild,共軛梯度法解方程 152.利用共軛梯度法計算矩陣--10^5階 193.利用cgs,bicg,bicgstab,等計算矩陣解 20第五章實驗 221.Newton插值f(x)=1/(1+4x^2).圖形，誤差 222.f(x)=1/(1+4x^2)插值 243.飛機的外形輪廓 24第九章四階龍格庫塔方法（自選） 251. 解算微分方程組 252. Matlab的四階Rk方法 253. 作圖比較，一定誤差 27第一章實驗1.根號5的極限代碼如下:x0=5;fori=0:1:100xi=sqrt(x0);ei=x0-xi;fprintf('NO:%d,xi=%0.8f,Theerroris:%0.8f\n',i,xi,ei);x0=xi;ifei<10^(-8)break;endend可以看出極限是1.2.求pi的近似值代碼如下：sum0=0;forn=1:1:50000y=(-1)^(n+1)/(2*n-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,Theerror:%0.8f\n',n,pi_2,error);ifabs(error)<10^(-4)break;endend求得結(jié)果是：3．畫圖y=sinx與其泰勒展開x=0:pi/100:2*pi;y=sin(x);y1=0;y2=0;y3=0;fori=0:2y1=y1+(-1)^(i)*x.^(2*i+1)/factorial(2*i+1);endfori=0:5y2=y2+(-1)^(i)*x.^(2*i+1)/factorial(2*i+1);endfori=0:10y3=y3+(-1)^(i)*x.^(2*i+1)/factorial(2*i+1);endplot(x,y,'*r',x,y1,'b',x,y2,'-g',x,y3,'k')axis([02*pi-1.51.5])看出n=2發(fā)散，n=10幾乎與y=sinx重合。4．三維圖mesh與surf實現(xiàn)。首先用mesh函數(shù).fori=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再用surf函數(shù).fori=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.列主元三角分解A.clcA=[11111;12345;1361015;14102035;15153570]E=eye(5)[m,n]=size(A);ifm~=nerror('不是方陣')returnendifdet(A)==0error('不能分解')endu=A;p=eye(m);l=eye(m);fori=1:mforj=i:mt(j)=u(j,i);fork=1:i-1t(j)=t(j)-u(j,k)*u(k,i);endenda=i;b=abs(t(i));forj=i+1:mifb<abs(t(j))b=abs(t(j));a=j;endendifa~=iforj=1:mc=u(i,j);u(i,j)=u(a,j);u(a,j)=c;endforj=1:mc=p(i,j);p(i,j)=p(a,j);p(a,j)=c;endc=t(a);t(a)=t(i);t(i)=c;endu(i,i)=t(i);forj=i+1:mu(j,i)=t(j)/t(i);endforj=i+1:mfork=1:i-1u(i,j)=u(i,j)-u(i,k)*u(k,j);endendendl=tril(u,-1)+eye(m)u=triu(u,0)B=A/E2.追趕法求方程，電路的電流A=[2-2000000;-25-200000;0-25-20000;00-25-2000;000-25200;0000-25-20;00000-25-2;000000-25]bb=[220/270000000]n=size(A,1);s=zeros(n,1);%è?3?èy????b=diag(A);a=diag(A,-1);c=diag(A,1);d=zeros(n,1);u=zeros(n-1,1);fori=1:n-1d(1)=b(1);u(i)=c(i)/d(i);d(i+1)=b(i+1)-a(i)*u(i);end%×·μ?1y3ìy=zeros(n,1);y(1)=bb(1)/d(1);fori=2:ny(i)=(bb(i)-a(i-1)*y(i-1))/d(i);end%??μ?1y3ìs(n)=y(n);fori=n-1:-1:1s(i)=y(i)-u(i)*s(i+1);ends3.方程組的性態(tài)與矩陣條件數(shù)的實驗clcn=5;A=zeros(n);C=zeros(n);b=zeros(1,n);fori=1:nx(i)=1+0.1*i;forj=1:nA(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=A\b'當(dāng)n=5,10,20時，cond(A)=5.3615e+05,8.6823e+11,1.5428e+22?？闯鲈絹碓讲B(tài)了。當(dāng)n=5,10,20時，解分別如下：說明條件數(shù)越大，越病態(tài)，發(fā)散了。當(dāng)n=10時就一個值有一定誤差。n=10,解方程（A+delta）x=b.clcn=10;A=zeros(n);C=zeros(n);b=zeros(1,n);fori=1:nx(i)=1+0.1*i;forj=1:nA(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=A\b'那么結(jié)果如下：可以看出很小的擾動導(dǎo)致解誤差變得較大了。4.Wilson矩陣。求特征值，條件數(shù)。誤差（1）行列式，條件數(shù)，特征值（2）(A+A0)(x+x0)=b,求此時的x0與||x0||2.代碼如下：clcA=[10787;7565;86109;75910]A1=[107.28.16.9;7.085.076.025;8.25.899.969.01;6.985.048.979.98]b=[32233331]A0=A1-Ax=A\b'x1=A1\b'x0=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*A0_A)/(1-condA*A0_A)結(jié)果如下：然后比較：然后利用公式計算:左邊=0.8225，右邊=1.0358則有0.8225<1.0358。

第三章實驗1.分別用Jacobi,Gauss-Seild,共軛梯度法解方程（1）Jacobi迭代法clcA=[101234;19-12-3;2-173-5;32312-1;4-3-5-115]b=[12-2714-1712]n=length(A);k=0;ep=10^(-7)it_max=100;x=zeros(n,1);y=zeros(n,1);fprintf('NO:x1x2x3x4x5\n')while1fori=1:ny(i)=b(i);forj=1:nifj~=iy(i)=y(i)-A(i,j)*x(j);endendifabs(A(i,i))<1e-10|k==it_maxreturn;endy(i)=y(i)/A(i,i);endifnorm(y-x,inf)<epbreak;endx=y;k=k+1;fprintf('NO:%d,%f%f%f%f%f\n',k,x(1,1),x(2,1),x(3,1),x(4,1),x(5,1))end得的結(jié)果為（2）Gauss-Seild迭代法clcA=[101234;19-12-3;2-173-5;32312-1;4-3-5-115]b=[12-2714-1712]x0=[00000]';n=length(A);ep=10^(-7);N=100;k=1;x=zeros(n,1);fprintf('\n**********start**************\n',k);fork=1:Nfprintf('%d\t',k);%calculatex(1)x(1)=b(1)/A(1,1);forj=2:nx(1)=x(1)-A(1,j)*x0(j)/A(1,1);endfprintf('%f\t',x(1));%calculatex(i)fori=2:n-1x(i)=b(i)/A(i,i);forj=1:i-1x(i)=x(i)-A(i,j)*x(j)/A(i,i);endforj=i+1:nx(i)=x(i)-A(i,j)*x0(j)/A(i,i);endfprintf('%f\t',x(i));end%calculatex(n)x(n)=b(n)/A(n,n);forj=1:n-1x(n)=x(n)-A(n,j)*x(j)/A(n,n);endfprintf('%f\t',x(n));ifnorm(x-x0,inf)<epbreak;elsek=k+1;x0=x;endfprintf('\n');endfprintf('\n***********end*************\n',k);x0運行結(jié)果如下：可以看出Gauss-Seild只要41次就可以計算結(jié)果，而Jacobi需要76次。（3）共軛梯度迭代法clcA=[101234;19-12-3;2-173-5;32312-1;4-3-5-115]b=[12-2714-1712]'x0=[00000]';N=length(A);eps=10^(-7);x=zeros(N,1);%μü′ú?ü???òá?.normê????ó￡??òá?·?êyr=b-A*x;%??êμ±íê?μ?ê??o???÷·??ò?￡fai=1/2x'AX-bxd=r;fork=0:N-1fprintf('μú%d′?μü′ú￡o',k+1);a=(norm(r)^2)/(d'*A*d)x=x+a*drr=b-A*x;%rr=r(k+1)if(norm(rr)<=eps)||(k==N-1)break;endB=(norm(rr)^2)/(norm(r)^2);d=rr+B*d;r=rr;end計算結(jié)果如下：理論上只需要N次得到精確值2.利用共軛梯度法計算矩陣--10^5階clcn=10^5;fori=1:nforj=1:nifi==jA(i,j)=3;elseifabs(i-j)==1A(i,j)=-1;elseif(i+j==n+1)&&(i~=n/2)&&(i~=n/2+1)A(i,j)=1/2;elseA(i,j)=0;endendendb(1,1)=2.5;b(1,2:10^5/2-1)=1.5;b(1,10^5/2)=1.0;b(1,10^5/2+1)=1.0;b(1,10^5/2+2:10^5-1)=1.5;b(1,10^5)=2.5;b=b';N=length(A);eps=10^(-7);x=zeros(N,1);%μü′ú?ü???òá?.normê????ó￡??òá?·?êyr=b-A*x;%??êμ±íê?μ?ê??o???÷·??ò?￡fai=1/2*x'AX-bxd=r;fork=0:N-1fprintf('μú%d′?μü′ú￡o',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;endB=(norm(rr)^2)/(norm(r)^2);d=rr+B*d;r=rr;end3.利用cgs,bicg,bicgstab,等計算矩陣解利用cgs（）函數(shù)A=gallery('wilk',21)b=sum(A,2);tol=1e-12;maxit=15;M1=diag([10:-1:111:10]);x=cgs(A,b,tol,maxit,M1)利用bicg（）函數(shù)A=gallery('wilk',21)b=sum(A,2);tol=1e-12;maxit=15;M1=diag([10:-1:111:10]);x=bicg(A,b,tol,maxit,M1)利用bicgstab（）函數(shù)A=gallery('wilk',21)b=sum(A,2);tol=1e-12;maxit=10;M1=diag([10:-1:111:10]);x=bicgstab(A,b,tol,maxit,M1)第五章實驗1.Newton插值f(x)=1/(1+4x^2).圖形，誤差（1）輸出插值多項式。程序：function[p]=NewtonChazhi(x,y,n)%UNTITLED3此處顯示有關(guān)此函數(shù)的摘要%此處顯示詳細說明f=zeros(n,n);f(:,1)=y';forj=2:nfori=j:nf(i,j)=(f(i,j-1)-f(i-1,j-1))/(x(i)-x(i-j+1));endendp=f(n,n);fork=n:-1:2p=conv(p,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)結(jié)果：牛頓插值多項式：polynomial=-(7319042784910035*x^10)/1475739525896764129