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

%定义zuoye7.m文件 function x=zuoye7(a,b,c,d) a1=[0;a];

n=length(b);q=zeros(n,1);p=zeros(n,1); %LU分解

q(1)=b(1);for k=2:n,p(k)=a1(k)/q(k-1); q(k)=b(k)-p(k)*c(k-1); end %解Ly=d

y=zeros(n,1);y(1)=d(1);for k=2:n, y(k)=d(k)-p(k)*y(k-1);end %解Ux=y

x=zeros(n,1); x(n)=y(n)/q(n);for k=n-1:-1:1,x(k)=(y(k)-c(k)*x(k+1))/q(k);end x

在Matlab窗口中执行： a=[-1 -1 -1 -1]'; b=[2 2 2 2 2]'; c=[-1 -1 -1 -1]'; d=[1 0 0 0 0]'; x=zuoye7(a,b,c,d) 运行结果如下： x = 0.83333333333333 0.66666666666667 0.50000000000000 0.33333333333333 0.16666666666667

9

9. 分别用Jacobi迭代法和Gauss-Seidel迭代法求解程组（编写程序）

?10x1+3x2?x3?14????2x1-10x2?3x3?-5? ?x?3x?10x?14?23?1?T取初值X（0）=（0,0,0），精确到小数后面四位。

解：

（1） Jacobi迭代法的Matlab程序：

% x0为初始向量，ep为精度，N为最大次数，x是近似解向量 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); n=length(b);x0=zeros(n,1);x=zeros(n,1);k=0 while K

x(i)=(b(i)-A(I,[1:i-1,i+1:n])*x0(1:i-1,i+1:n))/A(i,i); end

if norm(x-x0,inf)

if k==N,Warning(‘已达到迭代次数上限’);end disp([‘k=’,num2str(k)]),x 运行结果如下： k=15 x=

1.000000327906423

10

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