北京科技大学计算方法大作业(3)

0.999999801006905 1.000000327906432

（2）Gauss-Seidel迭代法的Matlab程序：

% x0为初始向量，ep为精度，N为最大次数，x是近似解向量 Format long;clear; A=[10 3 1;2 -10 3;1 3 10];

b=[14 -5 14];n=length(b);N=500;ep=1e-6;x0=zero(n,1);P=inf; %以下是Guass-Seidal迭代法程序

D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);dD=det(D); if dD==0

disp(‘请注意：因为对角矩阵D奇异，所以此方程组无解。’) else

disp(‘请注意：因为对角矩阵D非奇异，所以此方程组无解。’) iD=inv(D-L);B_GS=ID*U;d_GS=iD*b; end k=0; while k

if k==N,warning(‘已达到迭代次数上限’)；end disp([‘k=’,num2str(k)]),x,%D,U,L,B_GS,d_GS,%jx=A\\b, 运行结果如下：

11

请注意：因为对角矩阵D非奇异，所以此方程组有解。 k=9

x=1.000000043651052 1.000000031224555 0.999999986267528

10．编写幂法程序求矩阵

10.5??1?A??110.25???2??0.50.25?

?x1(??10?6)1按模最大的特征值和对应的特征向量。 解：幂法程序如下：

% A 为n 阶方阵，u0 为初始向量，

%epsilon 为上限，max1 为循环次数。lambda 返回按模最大的特征值，u 返回

对应的特征向量。

clear;format long;clc;max1=1000; epsilon=1e-6; A=[1 1 0.5;1 1 0.25;0.5 0.25 2]; u0=[1;1;1]; u=u0;v0=u0; lambda=0;k=0;err=1;R=[]; while ((kepsilon))

v=A*u; [m,j]=max(abs(v)); m=v(j); u=v/m; err=abs(lambda-m); lambda=m; k=k+1;

Temp=[k;m;v;u];R=[R Temp];

12

end

k, lambda, u, disp(R), xlswrite('d:\\\\1.xls', R,'sheet1', 'b1') ;

%[D,Lmd]=eig(A) %D 之列为A 的特征向量，Lmd 之对角线上为A 的特征值

运行结果如下： k = 23

lambda = 2.536527202038736 u = 0.748222175139906 0.649662222855908 1.000000000000000

11. 编写用原点位移加速反幂法程序求矩阵

? 7 3 -2??A?? 3 4 -1????-2 -1 3 ??

?10(??10) q?1.9最接近于的特征值和相应的特征向量。

解：

反幂法程序如下：

% A 为n 阶方阵，v0,u0 为初始向量，

%epsilon 为上限，max1 为循环次数，q 为近似特征值。lambda 为提高精度的

特征值，v 给出对应的特征向量。

13

clear;format long;max1=100; epsilon=1e-10;

A=[7 3 -2;3 4 -1;-2 -1 3]; u0=[1;1;1]; u=u0;v=u0;q=1.9; lambda=0;k=0;err=1;mu=0.5;

n=length(u0); z=zeros(n,1);A=A-q*eye(n); %LU 分解

n=length(u0);UU=zeros(n,n);L=eye(n,n); UU(1,:)=A(1,:); L(2:n,1)=A(2:n,1)/UU(1,1); for s=2:n

UU(s,s:n)=A(s,s:n)-L(s,1:s-1)*UU(1:s-1,s:n);

L(s+1:n,s)=(A(s+1:n,s)-L(s+1:n,1:s-1)*UU(1:s-1,s))/UU(s,s); end

while ((kepsilon)) [m,j]=max(abs(v)); m=v(j); u=v/m; %解下三角方程组Lz=u

z(1)=u(1);for j=2:n, z(j)=u(j)-L(j,1:j-1)*z(1:j-1); end %解上三角方程组Uv=z v(n)=z(n)/UU(n,n);

for j=n-1:-1:1, v(j)=(z(j)-UU(j,j+1:n)*v(j+1:n))/UU(j,j); end err=abs(1/m-1/mu); k=k+1; mu=m; end

disp(['迭代次数k= ',num2str(k)])

disp('要求的矩阵A 的特征值lambda='), lambda= q+1/m disp('要求的矩阵A 的特征向量u'), u

14

运行结果如下： 迭代次数k= 17

要求的矩阵A 的特征值lambda= lambda =

2.000000000009091

要求的矩阵A 的特征向量u u =

0.999999999973922 -0.999999999956935 1.000000000000000

12. 已知插值点

(-2.00,17.00), (0.00,1.00), (1.00,2.00), (2.00,17.00),

求三次插值多项式，并计算f(0.6). 解：

其Matlab程序如下： function yy=zuoye10(x,y,xx) m=length(x);n=length(y);

if m~=n, error('向量x与y的长度必须一致');end

15