数值分析 (2)

1、数值计算方法编程作业化工与环境学院 2120111381 王雪梅实验题一第二章 14解：（1）编程：function B=ZGF(A,d)m,n=size(A);b(1)=A(1,1);c(1)=A(1,2);a(n)=A(m,m-1);b(n)=A(m,m);for i=2:(m-1) a(i)=A(i,i-1); b(i)=A(i,i); c(i)=A(i,i+1);endu(1)=b(1);y(1)=d(1);for k=2:m l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*c(k-1); y(k)=d(k)-l(k)*y(k-1);endI(m)=y(m)/u(m)

2、;for k=m-1:-1:1 I(k)=(y(k)-c(k)*I(k+1)/u(k);end（2）计算：>> 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 5A=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

3、 0 0 0 -2 5 -2 0 0 0 0 0 0 -2 5 -2 0 0 0 0 0 0 -2 5>> d=220/27 0 0 0 0 0 0 0d =8.1481 0 0 0 0 0 0 0>> B=ZGF(A,d)I =8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477实验题二第三章 14 解：（1）用Jacobi迭代法解方程。编程：function x,k=Jacobi(A,b,x0,espi,N)m,n=size(A);k=1;x=x0;for i=1:m sum=0; for j=1:m if

4、j=i sum=sum+A(i,j)*x0(j); end end x(i)=(b(i)-sum)/A(i,i);endwhile abs(norm(x-x0)>espi&k<N x0=x; for i=1:m sum=0; for j=1:m if j=i sum=sum+A(i,j)*x0(j); end end x(i)=(b(i)-sum)/A(i,i); end k=k+1;end计算：>> 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,15 A = 10 1 2 3 4

5、 1 9 -1 2 -3 2 -1 7 3 -5 3 2 3 12 -1 4 -3 -5 -1 15>> b=12,-27,14,-17,12' b = 12 -27 14 -17 12>> x0=0,0,0,0,0' x0 = 0 0 0 0 0>> N=1000 N = 1000>> espi=1e-8 espi = 1.0000e-008>> x,k=Jacobi (A,b,x0,espi,N)x = 1.0000 -2.0000 3.0000 -2.0000 1.0000k = 91 （2）用Gauss_Sei

6、del迭代法解方程。编程：function x,k=GaussSeidel (A,b,x0,espi,N)m,n=size(A);k=1;x=x0;for i=1:m sum=0; for j=1:m if j<i sum=sum+A(i,j)*x(j); end if j>i sum=sum+A(i,j)*x0(j); end end x(i)=(b(i)-sum)/A(i,i);endwhile abs(norm(x-x0)>espi&k<N x0=x; for i=1:m sum=0; for j=1:m if j<i sum=sum+A(i,j)*

7、x(j); end if j>i sum=sum+A(i,j)*x0(j); end end x(i)=(b(i)-sum)/A(i,i); end k=k+1;end计算：>> 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,15 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 15>> b=12,-27,14,-17,12' b = 12 -27 14 -17 12>> x0=0,0,0

8、,0,0' x0 = 0 0 0 0 0>> N=1000 N = 1000>> espi=1e-8 espi = 1.0000e-008>> x,k= GaussSeidel (A,b,x0,espi,N)x = 1.0000 -2.0000 3.0000 -2.0000 1.0000k = 50 （3）用共轭梯度法解方程。编程：function x,k=GEF (A,b,x0,espi,N)d0=b-A*x0;r0=d0;lmd0=r0'*d0/(d0'*(A*d0);x=x0+lmd0*d0;k=1;while abs(norm

9、(x-x0)>espi & k<N r0=b-A*x; beta=-(r0'*(A*d0)/(d0'*(A*d0); d0=r0+beta*d0; lmd0= r0'*d0/(d0'*(A*d0); x0=x; x=x+lmd0*d0; k=k+1;end计算：>> 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,15 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 15>

10、;> b=12,-27,14,-17,12' b = 12 -27 14 -17 12>> x0=0,0,0,0,0' x0 = 0 0 0 0 0>> N=1000 N = 1000>> espi=1e-8 espi = 1.0000e-008>> x,k= GEF (A,b,x0,espi,N)x = 1.0000 -2.0000 3.0000 -2.0000 1.0000k = 6实验题三第五章 23.（1）用MatLab计算三弯矩方程参数 x=-520,-280,-156.6,-78,-39.62,-3.1,0,3.

11、1,39.62,78,156.6,280,520; y=0,-30,-36,-35,-28.44,-9.4,0,9.4,28.44,35,36,30,0;h=240,123.4,78.6,38.38,36.52,3.1,3.1,36.52,38.38,78.6,123.4,240; a=zeros(12,1);c=a;b=a;for i=1:11 a(i+1)=h(i)/(h(i)+h(i+1); b(i)=1-a(i); c(i)=6*(y(i+2)-y(i+1)/h(i+1)-(y(i+1)-y(i)/h(i)/(h(i)+h(i+1); end; a(1)=h(12)/(h(1)+h(1

12、2);b(12)=1-a(1);c(12)= 6/(h(1)+h(12)*(y(2)-y(1)/h(1)-(y(13)-y(12)/h(12);（2）解出M for i=1:12 for j=1:12 if j=i A(i,i)=2; elseif j=i+1 A(i,j)=b(i); elseif i=j+1 A(i,j)=a(i); end; end; end; A(1,12)=a(1);A(12,1)=b(12); X=Ac; M=zeros(13,1); for i=2:13 M(i)=X(i-1); end; M(1)=M(13);（3）将M代入插值方程式，计算方程的值 xx=-52

13、0:10.4:520;xx=xx; yy=zeros(101,1); for i=1:51 if x(1)<=xx(i)<x(2) yy(i)=M(2)*(xx(i)-x(1)3/(6*h(1)-M(1)*(xx(i)-x(2)3/(6*h(1)+ (y(2)-M(2)*h(1)2/6)*(xx(i)-x(1)/h(1)- (y(1)-M(1)*h(1)2/6)*(xx(i)-x(2)/h(1); elseif x(2)<=xx(i)<x(3) yy(i)=M(3)*(xx(i)-x(2)3/(6*h(2)-M(2)*(xx(i)-x(3)3/(6*h(2)+ (y(3

14、)-M(3)*h(2)2/6)*(xx(i)-x(2)/h(2)- (y(2)-M(2)*h(2)2/6)*(xx(i)-x(3)/h(2); elseif x(3)<=xx(i)<x(4) yy(i)=M(4)*(xx(i)-x(3)3/(6*h(3)-M(3)*(xx(i)-x(4)3/(6*h(3)+ (y(4)-M(4)*h(3)2/6)*(xx(i)-x(3)/h(3)- (y(3)-M(3)*h(3)2/6)*(xx(i)-x(4)/h(3); elseif x(4)<=xx(i)<x(5) yy(i)=M(5)*(xx(i)-x(4)3/(6*h(4)-M

15、(4)*(xx(i)-x(5)3/(6*h(4)+ (y(5)-M(5)*h(4)2/6)*(xx(i)-x(4)/h(4)- (y(4)-M(4)*h(4)2/6)*(xx(i)-x(5)/h(4); elseif x(5)<=xx(i)<x(6) yy(i)=M(6)*(xx(i)-x(5)3/(6*h(5)-M(5)*(xx(i)-x(6)3/(6*h(5)+ (y(6)-M(6)*h(5)2/6)*(xx(i)-x(5)/h(5)- (y(5)-M(5)*h(5)2/6)*(xx(i)-x(6)/h(5); elseif x(6)<=xx(i)<x(7) yy(

16、i)=M(7)*(xx(i)-x(6)3/(6*h(6)-M(6)*(xx(i)-x(7)3/(6*h(6)+ (y(7)-M(7)*h(6)2/6)*(xx(i)-x(6)/h(6)- (y(6)-M(6)*h(6)2/6)*(xx(i)-x(7)/h(6); elseif x(7)<=xx(i)<=x(8) yy(i)=M(8)*(xx(i)-x(7)3/(6*h(7)-M(7)*(xx(i)-x(8)3/(6*h(7)+ (y(8)-M(8)*h(7)2/6)*(xx(i)-x(7)/h(7)- (y(7)-M(7)*h(7)2/6)*(xx(i)-x(8)/h(7); en

17、d end（4）作图并画出另一半for i=52:101 yy(i)=-yy(102-i); %将y另一半求出； end; for i=1:50 xx(i)=-xx(i); %将x负值改为原值； end plot(xx,yy);所得的机翼外形图如下图所示：实验题四第六章 16.x0=2,3,5,6,7,9,10,11,12,14,16,17,19,20; y0=106.42,108.26,109.58,109.50,109.86,110.00,109.93,110.59,110.60,110.72,110.90,110.76,111.10,111.30;x=x0;y=y0;x=1./x;y=1

18、./y;n=size(x,2);tolx1=sum(x(:)；toly1=sum(y(:)；y=y.*x;toly2=sum(y(:)；x=x.2;tolx2=sum(x(:)；if n>=tolx1if n>=tolx2a1=(toly2-tolx1*toly1/n)/(tolx2-tolx1*tolx1/n);a0=(toly1-a1*tolx1)/n;elsea0=(toly1-toly2*tolx1/tolx2)/(n-tolx1*tolx1/tolx2);a1=(toly2-a0*tolx1)/tolx2；endelseif tolx1>tolx2a1=(toly1

19、-n*toly2/tolx1)/(tolx1-n*tolx2/tolx1)；a0=(toly2-tolx2*a1)/tolx1；elsea0=(toly1-toly2*tolx1/tolx2)/(n-tolx1*tolx1/tolx2);a1=(toly2-a0*tolx1)/tolx2；endendx=x0;x=1./x;y=a0+a1.*x;y=1./y;x=1./x;plot(x,y);hold on;plot(x0,y0,'bd');输出图形为：实验题五第七章 26.epsx=1e-3;epsy=1e-7;k=50;det=5/k;a=0:10/k:5;x(1)=0;y

20、(1)=0;fx2max=0;fy2max=0;for i = 2:1:k+1 ix(i)=0;y(i)=0;a=0;b=0+(i-1)*det;num=20;h=(b-a)/num;xk=a:h:b;for j =1:1:num+1fx2(j)=-sin(1/2*xk(j)*xk(j)-xk(j)*xk(j)*cos(1/2*xk(j)*xk(j);fy2(j)=cos(1/2*xk(j)*xk(j)-xk(j)*xk(j)*sin(1/2*xk(j)*xk(j);endfx2=abs(fx2); temp=max(fx2);if fx2max<temp fx2max=temp; en

21、drx=abs(b-a)*h*h*fx2max);fy2=abs(fy2); temp=max(fy2);if fy2max<temp fy2max=temp; endry=abs(b-a)*h*h*fy2max);while rx>epsx | ry>epsynum=num*1.5;h=(b-a)/num;xk=a:h:b;for j =1:1:num+1 fx2(j)=-sin(1/2*xk(j)*xk(j)-xk(j)*xk(j)*cos(1/2*xk(j)*xk(j);fy2(j)=cos(1/2*xk(j)*xk(j)-xk(j)*xk(j)*sin(1/2*xk(

22、j)*xk(j);endfx2=abs(fx2);temp=max(fx2); if fx2max<temp fx2max=temp; endrx=abs(b-a)*h*h*fx2max);fy2=abs(fy2); temp=max(fy2);if fy2max<temp fy2max=temp; endry=abs(b-a)*h*h*fy2max);end;x(i)=x(i)+cos(1/2*a*a)+cos(1/2*b*b);y(i)=y(i)+sin(1/2*a*a)+sin(1/2*b*b); for j =1:1:num+1x(i)=x(i)+2*cos(1/2*xk(

23、j)*xk(j);y(i)=y(i)+2*sin(1/2*xk(j)*xk(j); endx(i)=x(i)*h/2;y(i)=y(i)*h/2;endxx=x;yy=y;for i=1:1:k yy(k+i)=y(i); xx(k+i)=x(i); yy(i)=y(k-i+1); xx(i)=-x(k-i+1);end;plot(xx,yy);title('7.26考纽螺线图');输出结果：实验题六第八章 20.（1）简单迭代法x(1)=0.5;i=1;x(i+1)=exp(-x(i);while abs(x(i+1)-x(i)>1e-8; x(i+2)=exp(-x(i+1) i=i+1;end运行结果如下：0.50.6065306597126330.5452392118926050.5797030948780680.5600646279389020.5711721489772150.5648629469803230.5684380475700660.5664094527469210.5675596342622420.5669072129354710.5672771959707790.5670673518537280.5671863600876380.5671188642569860.5