计算方法 课内实验 解线性方程组的直接方法和迭代法

《计算方法》课内实验报告

学生姓名: 及 学 号： 学 院: 班 级: 课程名称： 实验题目: 指导教师 姓名及职称：

张学阳 理学院 数学101 计算方法

解线性方程组的直接方法和迭代法

宋云飞 讲 师 朱秀丽 讲 师 尚宝欣 讲 师

1009300132

2012年12月10日

目 录

一、实验题目........................................................................................ - 1 - 二、实验目的........................................................................................ - 1 - 三、实验内容........................................................................................ - 1 - 四、实验结果........................................................................................ - 2 - 五、实验体会或遇到问题 ................................................................... - 8 -

一、实验题目

解线性方程组的直接方法和迭代法

二、实验目的

1．熟悉Matlab编写及运行数值计算程序的方法。

2．进一步理解求解线性方程组的直接方法和迭代法基础理论。 3．进一步掌握应用不同的方法求解线性方程组的收敛速度及误差分析。

三、实验内容

1．用LU分解及列主元高斯消去法解线性方程组

?701??x1??8?10???32.09999692??x2??5.9000?01??????. ??15?1??x3??5?5???????102??x4??1?2?输出Ax?b中系数A?LU分解法的矩阵L及U，解向量x及detA；列主元法的行交换次序，解向量x及detA；比较两种方法所得的结果. 2．线性方程组Ax?b的A及b为

?10?7 A???8??7787??32??23?565??,b???,

6109??33????5910??31?则解x?(1,1,1,1)T.用Matlab内部函数求detA及A的所有特征值和cond(A)2.若令

78.17.2??10?7.085.0465??， A??A??5.989.899??8??99.98??6.995求解(A??A)(x??x)?b，输出向量?x和?x2，从理论结果和实际计算两方面分析线性方程组Ax?b解的相对误差?x2x2及A的相对误差?A2A2的关系。

- 1 -

3．给出线性方程组Hnx?b，其中系数矩阵Hn为希尔伯特矩阵： Hn??hij??Rn?n，hij?T1,i,j?1,2,?,n.

i?j?1假设x*??1,1,?,1??Rn,b?Hnx*，若取n?6,8，分别用雅可比迭代法及SOR迭代法(??1,1.25,1.5)求解，比较计算结果。

四、实验结果

1.文件程序：

function hl=zhjLU(A,b)

[n n] =size(A); RA=rank(A); if RA~=n

disp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下:'), RA,hl=det(A); return end

if RA==n for p=1:n

h(p)=det(A(1:p, 1:p)); end

hl=h(1:n); for i=1:n if h(1,i)==0

disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：'), hl;RA return end end

if h(1,i)~=0

disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：') for j=1:n

U(1,j)=A(1,j); end

for k=1:n for i=2:n

for j=2:n

L(1,1)=1;L(i,i)=1; if i>j

L(1,1)=1;L(2,1)=A(2,1)/U(1,1); L(i,1)=A(i,1)/U(1,1);

- 2 -

L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k); else

U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j); end end end end

hl;RA,U,L end end

文件程序：

function [RA,RB,n,X]=liezhu(A,b) B=[A b]; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0,

disp('请注意：因为RA~=RB，所以此方程组无解.') return end

if RA==RB if RA==n

disp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1

[Y,j]=max(abs(B(p:n,p))); C=B(p,:); B(p,:)= B(j+p-1,:); B(j+p-1,:)=C; for k=p+1:n

m= B(k,p)/ B(p,p);

B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1); end end

b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1

X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end else

disp('请注意：因为RA=RB

窗口输入输出

>> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2]; >> b=[8 5.9000001 5 1]'; >> zhjLU(A,b)

请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：

- 3 -

RA =4 %U矩阵 U =

10.0000 -7.0000 0 1.0000 0 2.1000 6.0000 2.3000 0 0 -2.1429 -4.2381 0 -0.0000 0 12.7333 %L矩阵 L =

1.0000 0 0 0 -0.3000 1.0000 0 0 0.5000 1.1905 1.0000 -0.0000

0.2000 1.1429 3.2000 1.0000 %矩阵顺序主子试 ans =

10.0000 -0.0000 -150.0001 -762.0001 >> y=inv(L)*b >> x=inv(U)*y %解向量 x =

-0.2749 -1.3299 1.2969

1.4398 %行列式

>> det(L)*det(U) ans = -573.0100

>> [RA,RB,n,X]=liezhu(A,b)

请注意：因为RA=RB=n，所以此方程组有唯一解. A =

10.0000 -7.0000 0 1.0000 0 2.5000 5.0000 -1.5000 0 0 6.0000 2.3000 0 0 0 5.0800 RA =4 RB = 4 n =4

%解向量 X =

0.0000 -1.0000 1.0000

- 4 -

1.0000 %行列式 >> det(A) ans = -762.0001 2．

>> A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; >> b=[32 23 33 31]';

>> cond(A,2)%求cond(A)2. ans =2.9841e+003 >> det(A)%求行列式 ans =1.0000

>> eig(A)%求A的特征值 ans =0.0102 0.8431 3.8581 30.2887

>> A1=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98]; >> x=[1 1 1 1]';

>> I=[1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1];

>> x1=-inv(I+inv(A)*(A1-A))*inv(A)*(A1-A)*x%求?x x1 = -10.5863 17.3741 -4.2258 2.5240

>> norm(x1,2)%求?x2 ans = 20.9322

>> norm(x1,2)/norm(x,2)%相对误差?x2x2 ans =

5.2931

>> norm(A1-A,2)/norm(A,2)%A的相对误差?A2A2

ans =0.0076 3.

雅可比迭代文件一

function X=jacdd(A,b,X0,P,wucha,max1) [n m]=size(A); for j=1:m

a(j)=sum(abs(A(:,j)))-2*(abs(A(j,j))); end

for i=1:n if a(i)>=0

disp('请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛') %return end end

- 5 -

if a(i)<0

disp('请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛 ') end

for k=1:max1 k

for j=1:m

X(j)=(b(j)-A(j,[1:j-1,j+1:m])*X0([1: j-1,j+1:m]))/A(j,j); end

X,djwcX=norm(X'-X0,P); xdwcX=djwcX/(norm(X',P)+eps); X0=X';X1=A\\b; if (djwcX

disp('请注意：雅可比迭代收敛,此方程组的精确解jX和近似解X如下：') return end end

if (djwcX>wucha)&(xdwcX>wucha)

disp('请注意：雅可比迭代次数已经超过最大迭代次数max1 ') end

a,X=X;jX=X1', 文件二 n=8; x0=[]; for i=1:n for j=i:n

A(i,j)=1/(i+j-1); A(j,i)=A(i,j); end

x0=[x0 1]; end x0=x0'; b=A*x0;

jacdd(A,b,x0,2,0.001,100) 运行结果 n=6

ans = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 n=8

ans = Columns 1 through 6

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 Columns 7 through 8

1.0000 1.0000 超松弛迭代文件一

function X=jacdd(A,b,X0,P,wucha,max1,w) [n m]=size(A); for j=1:m

a(j)=sum(abs(A(:,j)))-2*(abs(A(j,j)));

- 6 -

end

for i=1:n if a(i)>=0

disp('请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛') %return end end if a(i)<0

disp('请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛 ') end

for k=1:max1 k

for j=1:m

X(j)=(b(j)-A(j,[1:j-1,j+1:m])*X0([1: j-1,j+1:m]))/A(j,j); end

X,djwcX=norm(X'-X0,P); xdwcX=djwcX/(norm(X',P)+eps); X0=X';X1=A\\b; if (djwcX

disp('请注意：雅可比迭代收敛,此方程组的精确解jX和近似解X如下：') return end end

if (djwcX>wucha)&(xdwcX>wucha)

disp('请注意：雅可比迭代次数已经超过最大迭代次数max1 ') end

a,X=X;jX=X1', 文件二 n=8; w=1; x0=[]; for i=1:n for j=i:n

A(i,j)=1/(i+j-1); A(j,i)=A(i,j); end

x0=[x0 1]; end x0=x0'; b=A*x0;

jacdd(A,b,x0,2,0.001,100,w) 运行结果 n=6,w=1

ans =1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 w=1.25

- 7 -

ans = 1.0e+072 *

0.8986 1.9393 2.5849 3.0372 3.3749 3.6378 (不收敛) w=1.5

ans =1.0e+080 *

1.2736 2.7485 3.6635 4.3046 4.7832 5.1558 (不收敛) n=8,w=1

ans = Columns 1 through 6

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 Columns 7 through 8

1.0000 1.0000 w=1.25

ans = 1.0e+087 *

Columns 1 through 6

0.3650 0.8261 1.1329 1.3580 1.5322 1.6716 Columns 7 through 8 1.7860 1.8818 (不收敛) w=1.5

ans =1.0e+095 *

Columns 1 through 6

0.5139 1.1629 1.5948 1.9118 2.1569 2.3531 Columns 7 through 8

2.5142 2.6490 （不收敛）

五、实验体会或遇到问题

通过这次课内试验熟悉了Matlab编写及运行数值计算程序的方法，进一步理解了求解线性方程组的直接方法和迭代法基础理论，掌握了应用不同的方法求解线性方程组的收敛速度及误差分析，感觉收获很大。

- 8 -

本文来源：https://www.bwwdw.com/article/mipv.html