4解线性方程组的迭代法
4.1 迭代法和敛散性及其MATLAB 程序
4.1.2 迭代法敛散性的判别及其MATLAB 程序
用谱半径判别迭代法产生的迭代序列的敛散性的MATLAB 主程序
function H=ddpbj(B)
H=eig(B);mH=norm(H,inf); if mH>=1
disp(' 请注意：因为谱半径不小于1，所以迭代序列发散, 谱半径mH 和B 的
所有的特征值H 如下：' )
else
disp(' 请注意：因为谱半径小于1，所以迭代序列收敛, 谱半径mH 和B 的所
有的特征值H 如下：' )
end mH
4.2 雅可比（Jacobi ）迭代及其MATLAB 程序
4.2.2 雅可比迭代的收敛性及其MATLAB 程序
判别雅可比迭代收敛性的MATLAB 主程序
function a=jspb(A) [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
if a(i)
disp(' 请注意：系数矩阵A 是严格对角占优的，此方程组有唯一解，且雅
可比迭代收敛 ') end
例4.2.2 用判别雅可比迭代收敛性的MATLAB 主程序，判别由下列方程组的雅可比迭代产生的序列是否收敛？
⎧10x 1-x 2-2x 3=7. 2, ⎧10x 1-x 2-2x 3=7. 2,
⎪⎪
-x +10x -2x =8. 3, ⎨-x 1+10x 2-2x 3=8. 3, ⎨123
⎪-x -x +0. 5x =4. 2. ⎪-x -x +5x =4. 2;
23123⎩（1）（2）⎩1
解 （1）首先保存名为jspb.m 的M 文件，然后在MATLAB 工作窗口输入程序
>> A=[10 -1 -2;-1 10 -2;-1 -1 5];a=jspb(A)
运行后输出结果
请注意：系数矩阵A 是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛 a =
-8 -8 -1
（2）在MATLAB 工作窗口输入程序
>> A=[10 -1 -2;-1 10 -2;-1 -1 0.5];a=jspb(A)
运行后输出结果
请注意：系数矩阵A 不是严格对角占优的，此雅可比迭代不一定收敛 a =
-8.0000e+000 -8.0000e+000 3.5000e+000
4.2.3 雅可比迭代的两种MATLAB 程序
（一） 雅可比迭代公式的MATLAB 程序
用雅可比迭代解线性方程组AX =b 的MATLAB 主程序
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
if a(i)
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',
例4.2.3 用∞范数和判别雅可比迭代的MATLAB 主程序解例4.2.2 中的方程组，解的精度为0.001，分别取最大迭代次数max1=100，5，初始向量X 0=（0 0 0）T ，并比较它们的收敛速度.
解 （1）取最大迭代次数max1=100时.
①首先保存名为jacdd.m 的M 文件，然后在MATLAB 工作窗口输入程序
>> A=[10 -1 -2;-1 10 -2;-1 -1 5]; b=[7.2;8.3;4.2];
X0=[0 0 0]'; X=jacdd(A,b,X0,inf,0.001,100)
运行后输出结果
请注意：系数矩阵A 是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛
请注意：雅可比迭代收敛, 此方程组的精确解jX 和近似解X 如下： a =
-8 -8 -1 jX =
1.1000 1.2000 1.3000 X =
1.0994 1.1994 1.2993
②在MATLAB 工作窗口输入程序
>> A=[10 -1 -2;-1 10 -2;-1 -1 0.5]; b=[7.2;8.3;4.2]; X0=[0 0
0]';
X=jacdd(A,b,X0,inf, 0.001,100)
运行后输出结果
请注意：系数矩阵A 不是严格对角占优的，此雅可比迭代不一定收敛 请注意：雅可比迭代收敛, 此方程组的精确解jX 和近似解X 如下： a =
-8.0000 -8.0000 3.5000 jX =
24.5000 24.6000 106.6000 X =
24.0738 24.1738 104.7974
（2）取最大迭代次数max1=5时, ①在MATLAB 工作窗口输入程序
>> A=[10 -1 -2;-1 10 -2;-1 -1 5];
b=[7.2;8.3;4.2]; X0=[0 0 0]'; X=jacdd(A,b,X0,inf,0.001,5)
运行后输出结果
请注意：系数矩阵A 是严格对角占优的，此方程组有唯一解，雅可比迭代收敛 请注意：雅可比迭代次数已经超过最大迭代次数max1 a =
-8 -8 -1 jX =
1.1000 1.2000 1.3000 X =
1.0951 1.1951 1.2941
②在MATLAB 工作窗口输入程序
>> A=[10 -1 -2;-1 10 -2;-1 -1 0.5]; b=[7.2;8.3;4.2]; X0=[0 0 0]'; X=jacdd(A,b,X0,inf, 0.001,5)
运行后输出结果
请注意：系数矩阵A 不是严格对角占优的，此雅可比迭代不一定收敛 请注意：雅可比迭代次数已经超过最大迭代次数max1 a =
-8.0000 -8.0000 3.5000 jX =
24.5000 24.6000 106.6000 X =
5.5490 5.6490 27.6553
由（1）和（2）可见，如果系数矩阵A 是严格对角占优的，则雅可比迭代收敛的速度快；如果系数矩阵A 不是严格对角占优的，则雅可比迭代收敛的速度慢. 因此，
a =∑a kj -a kk
j =1i ≠k
n
(k =1, 2, , n ) 的值越小，雅可比迭代收敛的速度越快.
（二）利用雅可比迭代定义编写的解线性方程组的MATLAB 程序
利用雅可比迭代定义编写解线性方程组（4.5）的MATLAB 程序的一般步骤 步骤1 在MATLAB 工作窗口输入程序
>> A=[a11 a12 …a1n; a21 a22 …a2n;…; an1 an2 …ann;]; D=diag(A) U=triu(A,1), L=tril(A,-1)
L ，U ； 运行后即可输出A 的D ，
步骤2 在MATLAB 工作窗口输入程序
>>dD=det(D); if dD==0
disp(' 请注意：因为对角矩阵D 奇异，所以此方程组无解.' ) else
disp(' 请注意：因为对角矩阵D 非奇异，所以此方程组有解.' ) iD=inv(D); B1=iD*(L+U);f1=iD*b; for k=1:max1
X= B1*X0+ f1; X0=X; djwcX=norm(X-X0,P); xdwcX=djwcX/(norm(X,P)+eps); X1=A\b; if (djwcX
' )
return end end
if (djwcX>wucha)|(xdwcX>wucha)
disp(' 请注意：雅可比迭代次数已经超过最大迭代次数max1 ') end end
a,X=X;jX=X1',
4.3 高斯-塞德尔（Gauss-Seidel ）迭代及其MATLAB 程序
4.3.3 高斯-塞德尔迭代两种MATLAB 程序
（一） 高斯-塞德尔迭代定义的MATLAB 程序1
用高斯-塞德尔迭代定义解线性方程组AX b 的MATLAB 主程序1
function X=gsdddy(A,b,X0,P,wucha,max1)
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); B2=iD*U;f2=iD*b;jX=A\b; X=X0; [n m]=size(A); for k=1:max1
X1= B2*X+f2; djwcX=norm(X1-X,P); xdwcX=djwcX/(norm(X,P)+eps);
if (djwcX
return else
k,X1',k=k+1;X=X1;
end end
if (djwcX
disp(' 请注意：高斯-塞德尔迭代收敛, 此A 的分解矩阵D,U,L 和方程
组的精确解jX 和近似解X 如下： ') else
disp(' 请注意：高斯-塞德尔迭代的结果没有达到给定的精度，并且迭代次数已经超过最大迭代次数max1, 方程组的精确解jX 和迭代向量X 如下： ')
X=X';jX=jX' end end
X=X';D,U,L,jX=jX'
例4.3.3 用高斯-塞德尔迭代定义的MATLAB 主程序解下列线性方程组, 取初始值
(k +1)
(0) (0) (0) -x (k ) (x 1, x 2, x 3) =(0, 0, 0) ，要求当x
时，迭代终止.
⎧3x 1+4x 2-5x 3+7x 4=5, ⎪2x -8x +3x -2x =2,
⎧10x 1-x 2-2x 3=7. 2, ⎪1234
⎨⎪
⎨-x 1+10x 2-2x 3=8. 3, ⎪4x 1+51x 2-13x 3+16x 4=-1, ⎪-x -x +0. 5x =4. 2. ⎪7x -2x 2+21x 3+3x 4=21. 23
（1）⎩1（2）⎩1
∞
解 （1）首先保存名为gsdddy.m 的M 文件，然后在MATLAB 工作窗口输入程序
>> A=[10 -1 -2;-1 10 -2;-1 -1 0.5]; b=[7.2;8.3;4.2]; X0=[0 0 0]';
X=gsdddy(A,b,X0,inf, 0.001,100)
运行后输出结果
D =
10.0000 0 0 0 10.0000 0 0 0 0.5000 U =
0 1 2 0 0 2 0 0 0
L =
0 0 0 1 0 0 1 1 0 jX =
24.5000 24.6000 106.6000 X =
24.4996 24.5996 106.5984
请注意：因为对角矩阵D 非奇异，所以此方程组有解.
请注意：高斯-塞德尔迭代收敛, 此A 的分解矩阵D,U,L 和方程组的精确解jX 和
近似解X 如下：
此近似解与例4.2.3中的（1）中②的解X =（24.073 8, 24.173 8, 104.797 4）T 比较，在相同的条件下, 高斯-塞德尔迭代比雅可比迭代得到的近似解的精度更高.
（2）在MATLAB 工作窗口输入程序
>> A=[3 4 -5 7;2 -8 3 -2;4 51 -13 16;7 -2 21 3];b=[5;2;-1;21]; X0=[0 0 0 0]';X=gsdddy(A,b,X0,inf,0.001,100)
运行后输出结果
请注意：因为对角矩阵D 非奇异，所以此方程组有解.
请注意：高斯-塞德尔迭代的记过没有达到给定的精度，并且迭代次数已经超过最
大迭代次数max1, 方程组的精确解jX 和迭代向量X 如下：
jX =
0.1821 -0.2571 0.7286 1.3036
X = 1.0e+142 *
0.2883 0.1062 0.3622 -3.1374
（二） 高斯-塞德尔迭代公式的MATLAB 程序2
用高斯-塞德尔迭代解线性方程组AX b 的MATLAB 主程序2
function X=gsdd(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
if a(i)
disp(' 请注意：系数矩阵A 是严格对角占优的，此方程组有唯一解，且高
斯-塞德尔迭代收敛 ')
end
for k=1:max1 for j=1:m if j==1
X(1)=(b(1)-A(1,2:m)*X0(2:m))/A(1,1) end
if j==m
X(m)=(b(m)-A(m,1:M1)*X(1:M1)')/A(m,m); end
for j=2:M1
X(j)=(b(j)-A(j,1:j-1)*X(1:j-1) -A(j,j+1:m)*X(j+1:m))/A(j,j);
end end
djwcX=norm(X'-X0,P);
xdwcX=djwcX/(norm(X',P)+eps); X0=X';X1=A\b; if (djwcX
X 如下： ')
return end end
if (djwcX>wucha)&(xdwcX>wucha)
disp(' 请注意：高斯-塞德尔迭代次数已经超过最大迭代次数max1 ') end
a,X=X;jX=X1',
4.4 解方程组的超松弛迭代法及其MATLAB 程序
4.4.2 超松弛迭代法收敛性及其MATLAB 程序
用谱半径判别超松弛迭代法产生的迭代序列的敛散性的MATLAB 主程序
function H=ddpbj(A,om)
D=diag(diag(A));U=-triu(A,1); L=-tril(A,-1); iD=inv(D-om*L); B2=iD*(om*U+(1-om)*D); H=eig(B2);mH=norm(H,inf); if mH>=1
disp(' 请注意：因为谱半径不小于1，所以超松弛迭代序列发散, 谱半径
mH 和B 的所有的特征值H 如下：' ) else
disp(' 请注意：因为谱半径小于1，所以超松弛迭代序列收敛, 谱半径mH
和B 的所有的特征值H 如下：' ) end mH
例4.4.1 当取ω=1.15，5时，判别用超松弛迭代法解下列方程组产生的迭代序列是否收敛？
⎧5x 1+x 2-x 3-2x 4=4⎪2x +8x +x +3x =1⎪1234⎨
⎪x 1-2x 2-4x 3-x 4=6⎪⎩-x 1+3x 2+2x 3+7x 4=-3
解 （1）当取ω=1.15时，首先保存名为ddpbj.m 的M 文件，然后在MATLAB 工作窗口输入程序
>> A=[5 1 -1 -2;2 8 1 3;1 -2 -4 -1;-1 3 2 7]; H=ddpbj(A,1.15)
运行后输出结果
请注意：因为谱半径小于1，所以超松弛迭代序列收敛, 谱半径mH 和B 的所有的
特征值H 如下：
mH =
0.1596 H =
0.1049 + 0.1203i 0.1049 - 0.1203i -0.1295 + 0.0556i -0.1295 -
0.0556i
（2）当取ω=5时，然后在MATLAB 工作窗口输入程序
>> H=ddpbj(A, 5)
运行后输出结果
请注意：因为谱半径不小于1，所以超松弛迭代序列发散, 谱半径mH 和B 的所有
的特征值H 如下：
mH =
14.1082 H =
-14.1082 -2.5107 0.5996 + 2.6206i 0.5996 - 2.6206i
4.4.3 超松弛迭代法的MATLAB 程序
用超松弛迭代法解线性方程组AX =b 的MATLAB 主程序
function X=cscdd (A,b,X,om,wucha,max1) D=diag(diag(A));U=-triu(A,1);
L=-tril(A,-1); jX=A\b;[n m]=size(A); iD=inv(D-om*L); B2=iD*(om*U+(1-om)*D); H=eig(B2);mH=norm(H,inf); for k=1:max1
iD=inv(D-om*L); B2=iD*(om*U+(1-om)*D); f2= om*iD*b; X1= B2*X+f2; X=X1; djwcX=norm(X1-jX,inf);
xdwcX=djwcX/(norm(X,inf)+eps);
if (djwcX
disp(' 谱半径mH ，A 的分解矩阵D,U,L 和方程组的精确解jX, 迭代次数
i 如下： ')
mH,D,U,L,jX=jX', i=k-1, return if i> max1
disp(' 迭代次数已经超过最大迭代次数max1, 谱半径mH ，方程组的精
确解jX, 迭代次数i 如下： ')
mH,D,U,L,jX=jX', i=k-1,
end end end
if mH>=1
disp(' 请注意：因为谱半径不小于1，所以超松弛迭代序列发散.' )
disp(' 谱半径mH ，A 的分解矩阵D,U,L 和方程组的精确解jX, 迭代次数
i 和迭代序列X 如下：' ) i=k-1,mH,D,U,L,jX, else
disp(' 因为谱半径小于1，所以超松弛迭代序列收敛，近似解X 如下： ') end 或
function X=cscdd1 (A,b,X,om,wucha,max1)
D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1); jX=A\b;[n
m]=size(A);
iD=inv(D-om*L); B2=iD*(om*U+(1-om)*D); H=eig(B2);mH=norm(H,inf); for k=1:max1
iD=inv(D-om*L); B2=iD*(om*U+(1-om)*D);
f2= om*iD*b; X1= B2*X+f2; X=X1; djwcX=norm(X1-jX,inf); xdwcX=djwcX/(norm(X,inf)+eps); end
if mH>=1
disp(' 请注意：因为谱半径不小于1，所以超松弛迭代序列发散. 谱半径mH ，A 的分解矩阵D,U,L 和方程组的精确解jX 和近似解X 如下：' ) else
disp(' 请注意：因为谱半径小于1，所以超松弛迭代序列收敛.' ) if (djwcX
disp(' 谱半径mH ，A 的分解矩阵D,U,L 和方程组的精确解jX 和近似
解X 如下： ')
mH,D,U,L,jX=jX',
else
disp(' 迭代次数已经超过最大迭代次数max1, 谱半径mH ，方程组的
精确解jX 和迭代向量X 如下： ')
mH,D,U,L,X=X1';jX=jX' return end end
例4.4.3 用超松弛迭代法（取ω=1.15和5）解例4.4.1中的线性方程组.
解 （1）当取ω=1.15时，首先保存名为cscdd.m 的M 文件，然后在MATLAB 工作窗口输入程序
>> A=[5 1 -1 -2;2 8 1 3;1 -2 -4 -1;-1 3 2 7];b=[4;1;6;-3]; X=[0 0 0 0]';X=cscdd (A,b,X,1.15,0.001,100),
运行后输出结果
谱半径mH ，A 的分解矩阵D,U,L 和方程组的精确解jX, 迭代次数i 如下： mH =
0.1596 D =
5 0 0 0 0 8 0 0 0 0 -4 0 0 0 0 7 U =
0 -1 1 2
0 0 -1 -3 0 0 0 1 0 0 0 0 L =
0 0 0 0 -2 0 0 0 -1 2 0 0 1 -3 -2 0 jX =
0.4491 0.2096 -1.4850 -0.0299 i = 3
因为谱半径小于1，所以超松弛迭代序列收敛，近似解X 如下： X =
0.4484 0.2100 -1.4858 -0.0303
（2）当取 =5时，保存名为cscdd.m 的M 文件，然后在MATLAB 工作窗口输入程序
>> A=[5 1 -1 -2;2 8 1 3;1 -2 -4 -1;-1 3 2 7];b=[4;1;6;-3]; X=[0 0 0 0]';X=cscdd (A,b,X,5,0.001,100),
运行后输出结果如下：
请注意：因为谱半径不小于1，所以超松弛迭代序列发散. 谱半径mH ，A 的分解矩阵D,U,L 和方程组的精确解jX, 迭代次数i 和迭代序列X
如下：
i = mH =
99 14.1082 D =
5 0 0 0 0 8 0 0 0 0 -4 0 0 0 0 7 U =
0 -1 1 2 0 0 -1 -3 0 0 0 1 0 0 0 0 L =
0 0 0 0 -2 0 0 0 -1 2 0 0 1 -3 -2 0
jX = X =1.0e+114 * 0.4491 -0.3122 0.2096 1.0497 -1.4850 -3.7174 -0.0299 3.9615