数值分析 - 编程作业(2)

sum=sum+A(i,j)*x0(j); end;

x1(i)=(x1(i)-sum)/A(i,i);

end; sum=0;

for j=1:4

sum=sum+A(5,j)*x0(j); end;

x1(5)=(b(5)-sum)/A(5,5); k=k+1;

end; x1 k

输出结果： x1 =

1.000000024731627 -2.000000022663166 -1.999999993335823 0.999999970325168

k =89

迭代了89次。

2.999999958446139

第四章矩阵特征值与特征向量的计算

1、已知矩阵

?190?66??336??3066303?168?36?84421472830??36?? ?112??291?试用幂法求按模最大的特征值与特征向量； 试用反幂法求按模最小的特征值与特征向量；

解：（1）幂法 程序：

A=[190,66,-84,30;66,303,42,-36;336,-168,147,-112;30,-36,28,291]; x=[0;0;0;1]; deg=1e-5; N=100;

k=0; eigenvalue0=0; [xr,r]=max(abs(x));

a=xr;

y=x/a; x=A*y; eigenvalue=xr;

while abs(eigenvalue-eigenvalue0)>=deg & k

[xr,r]=max(abs(x));

a=xr;

y=x/a; x=A*y; eigenvalue=xr; end; eigenvalue x k 输出结果： eigenvalue =

3.430000571059348e+02 x =

-1.143330070136651e+02 -3.430000489479517e+02 0.000000000000058e+02 1.715010523810607e+02 k = 87

（2）反幂法

A=[190,66,-84,30;66,303,42,-36;336,-168,147,-112;30,-36,28,291]; B=inv(A);

x=[0;0;0;1]; deg=1e-5; N=1000;

k=0; eigenvalue0=0; [xr,r]=max(abs(x));

a=xr;

y=x/a; x=B*y; eigenvalue=xr;

while abs(eigenvalue-eigenvalue0)>=deg & k

[xr,r]=max(abs(x));

a=xr;

y=x/a; x=B*y; eigenvalue=xr; end; eigenvalue x k 输出结果： eigenvalue =

0.002637740424370 x =

-0.001704294528000 0.000852172462123 0.008922284278451 0.000568148573659 k = 68

第五章插值法

3、已知直升机旋转机翼外形轮廓线上某些型值点的数据： kxkyk 1 520 0 2 280 -30 3 156.6 -36 4 78 -35 kxkyk 8 3.1 9.4 9 39.62 28.44 10 78 35 11 156.6 36 5 39.62 -28.44 12 280 30 6 3.1 -9.4 13 520 0 7 0 0 由表中13个节点，用样条插值的方法，增加平面点为101各点，绘制出较平滑的机翼外形曲线图。

解：了解三次样条插值方法的使用。 程序如下：

x=[-520,-280,-156.6,-78,-39.62,-3.1,0,3.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]'; n=13; %求解M for i=1:1:n-1

h(i)=x(i+1)-x(i); end

for i=2:1:n-1

a(i)=h(i-1)/(h(i-1)+h(i)); b(i)=1-a(i);

c(i)=6*((y(i+1)-y(i))/h(i)-(y(i)-y(i-1))/h(i-1))/(h(i-1)+h(i)); end

a(n)=h(n-1)/(h(1)+h(n-1)); b(n)=h(1)/(h(1)+h(n-1));

c(n)=6/(h(1)+h(n-1))*((y(2)-y(1))/h(1)-(y(n)-y(n-1))/h(n-1)); A(1,1)=2; A(1,2)=b(2); A(1,n-1)=a(2); A(n-1,n-2)=a(n); A(n-1,n-1)=2; A(n-1,1)=b(n); for i=2:1:n-2 A(i,i)=2;

A(i,i+1)=b(i+1); A(i,i-1)=a(i+1); end

C=c(2:n); C=C'; m=A\\C;

M(1)=m(n-1); M(2:n)=m;

xx=-520:10.4:520; for i=1:51

for j=1:1:n-1

if x(j)<=xx(i) && xx(i)

end

yy(i)=M(j+1)*(xx(i)-x(j))^3/(6*h(j))-M(j)*(xx(i)-x(j+1))^3/(6*h(j))+(y(j+1)-M(j+1)*h(j)^2/6)*(xx(i)-x(j))/h(j)-(y(j)-M(j)*h(j)^2/6)*(xx(i)-x(j+1))/h(j); end;

for i=52:101

yy(i)=-yy(102-i); end;

for i=1:50

xx(i)=-xx(i); end; plot(xx,yy); hold on; for i=1:1:n/2 x(i)=-x(i); end

plot(x,y,'bd'); title('图形'); 输出结果：

第六章函数逼近

1、钢包问题。炼钢厂出钢时所用的盛钢水的钢包，在使用过程中由于钢液及炉渣对包衬耐火材料的侵蚀，使其容积不断增大。经试验，钢包的容积与相应的使用次数的数据链表如下：

使用次数x容积y 使用次数x容积y 2 106.42 11 3 108.26 110.59