1、 第四章 解线性方程组的迭代法的 MATLAB 程序48. 高等教育出版社 教育电子音像出版社 作者：任玉杰 解线性方程组 的迭代法是从初始解出发，根据设计好的步骤用逐次求出的bAX近似解逼近精确解.在第三章中介绍的解线性方程组的直接方法一般适合于 为低阶稠A密矩阵（指 n 不大且元多为非零）的情况，而在工程技术和科学计算中常会遇到大型稀疏矩阵（指 n 很大且零元较多）的方程组，迭代法在计算和存贮两方面都适合后一种情况.由于迭代法是通过逐次迭代来逼近方程组的解，所以收敛性和收敛速度是构造迭代法时应该注意的问题.另外，因为不同的系数矩阵具有不同的性态，所以大多数迭代方法都具有一定的适用范围.有时
2、，某种方法对于一类方程组迭代收敛，而对另一类方程组迭代时就发散.因此，我们应该学会针对具有不同性质的线性方程组构造不同的迭代.4.1 迭代法和敛散性及其 MATLAB 程序4.1.2 迭代法敛散性的判别及其 MATLAB 程序根据定理 4.1 和谱半径定义，现提供一个名为 pddpb.m 的 M 文件，用于判别迭代公式（4.7）产生的迭代序列的敛散性.用谱半径判别迭代法产生的迭代序列的敛散性的 MATLAB 主程序输入的量：线性方程组 的迭代公式（4.7）中的 ；bAXB输出的量：矩阵 的所有特征值和谱半径 mH 及其有关迭代法产生的迭代B)(序列的敛散性的相关信息.function H=dd
3、pbj(B)H=eig(B);mH=norm(H,inf);if mH=1disp(请注意：因为谱半径不小于1，所以迭代序列发散, 谱半径mH和B的所有的特征值H如下： )elsedisp(请注意：因为谱半径小于1，所以迭代序列收敛, 谱半径mH和B的所有的特征值H如下： )endmH4.1.3 与迭代法有关的 MATLAB 命令（一） 提取（产生）对角矩阵和特征值可以用表 41 的 MATLAB 命令提取对角矩阵和特征值. 表 41 提取（产生）对角矩阵和特征值MATLAB 命令 功 能DX=diag(X) 若输入向量 X，则输出 DX 是以 X 为对角元的对角矩阵；若输入矩阵 X，则输出
4、DX 是以 X 的对角元构成的向量；DX=diag(diag(X) 输入矩阵 X，输出 DX 是以 X 的对角元构成的对角矩阵，可用于迭代法中从 A 中提取 D.lm=eig(A) 输入矩阵 A，输出 lm 是 A 的所有特征值.第四章 解线性方程组的迭代法第四章 解线性方程组的迭代法的 MATLAB 程序49. 高等教育出版社 教育电子音像出版社 作者：任玉杰 (二) 提取（产生）上（下）三角形矩阵 可以用表 42 的 MATLAB 命令提取矩阵的上三角形矩阵和下三角形矩阵. 表 42 提取矩阵的上三角形矩阵和下三角形矩阵MATLAB 命令 功 能U=triu(A) 输入矩阵 ，输出 是 的
5、上三角形矩阵.AUL=tril(A) 输入矩阵 ，输出 是 的下三角形矩阵.LU=triu(A,1) 输入矩阵 ，输出 是 的上三角形矩阵，但对角元为 0，可用于迭代法中从 中提取 .L=tril(A,-1) 输入矩阵 ，输出 是 的下三角形矩阵，但对角元为 0，可用于迭代法中从 中提取 .(三)稀疏矩阵的处理对稀疏矩阵在存贮和运算上的特殊处理，是 MATLAB 进行大规模科学计算时的特点和优势之一.可以用表 43 的 MATLAB 命令，输入稀疏矩阵的非零元（零元不必输入） ，即可进行运算. 表 43 稀疏矩阵的 MATLAB 命令MATLAB 命令 功 能ZA=sparse(r,c,v,m
6、,n) 表示在第 r 行、第 c 列输入数值 v，矩阵共 m 行 n列，输出 ZA，给出 (r, c) 及 v 为一稀疏矩阵.MA=full(ZA) 输入稀疏矩阵 ZA，输出为满矩阵 MA（包含零元）4.2 雅可比（Jacobi）迭代及其 MATLAB 程序4.2.2 雅可比迭代的收敛性及其 MATLAB 程序根据定理 4.3 和公式（4.14）,现提供一个名为 jspb.m 的 M 文件如下：判别雅可比迭代收敛性的 MATLAB 主程序输入的量：线性方程组 的系数矩阵 ；bAXA输出的量：系数矩阵 的 的值和有关nijakkijja1 ),21(n雅可比迭代收敛性的相关信息.function
7、 a=jspb(A)n m=size(A);for j=1:ma(j)=sum(abs(A(:,j)-2*(abs(A(j,j);endfor i=1:nif a(i)=0disp(请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛 )returnendendif a(i) 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)运行后输
8、出结果请注意：系数矩阵 A 不是严格对角占优的，此雅可比迭代不一定收敛a = -8.0000e+000 -8.0000e+000 3.5000e+0004.2.3 雅可比迭代的两种 MATLAB 程序利用 MATLAB 程序和雅可比迭代解线性方程组 的常用的方法有两种，一种bAX方法是根据雅可比迭代公式（4.11） 、 （4.12）式、定理 4.3 和公式（4.14）编写一个名为 jacdd.m 的 M 文件并保存，然后在 MATLAB 工作窗口输入对应的命令，执行此文件；另一种方法是根据雅可比迭代的定义，利用提取对角矩阵和上、下三角矩阵的 MATLAB程序解线性方程组 .下面我们分别介绍这两
9、种方法.bAX（一） 雅可比迭代公式的 MATLAB 程序用雅可比迭代解线性方程组 的 MATLAB 主程序输入的量：线性方程组 的系数矩阵 和 b, 初始向量 X0,范数的名称 P = A1, 2, inf 或 fro.，近似解 X 的误差（精度）wucha 和迭代的最大次数 max1；输出的量：系数矩阵 的 的值和有关Anijakkijja1 ),21(n雅可比迭代收敛性的相关信息及其 的精确解 jX 和近似解 X.b根据雅可比迭代（4.11） 、 （4.12）式、定理4.3和公式（4.14） ，现提供一个名为jacdd.m的M文件如下：function X=jacdd(A,b,X0,P,
10、wucha,max1)n m=size(A);for j=1:ma(j)=sum(abs(A(:,j)-2*(abs(A(j,j);endfor i=1:nif a(i)=0disp(请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛)returnendendif a(i)wucha)jX=X1,例 4.2.3 用 范数和判别雅可比迭代的 MATLAB 主程序解例 4.2.2 中的方程组，解的精度为 0.001，分别取最大迭代次数 max1=100，5， 初始向量 X0=（0 0 0） T，并比较它们的收敛速度.解 （1）取最大迭代次数 max1=100 时.首先保存名为 jacdd
11、.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 -1jX = 1.1000 1.2000 1.3000X = 1.0994 1.1994 1.2993在 MATLAB 工作窗口输入程序 A=10 -1 -2;-1 10 -2;-1 -1 0.5;
12、 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.5000jX = 24.5000 24.6000 106.6000X = 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
13、0 0; X=jacdd(A,b,X0,inf,0.001,5)运行后输出结果请注意：系数矩阵 A 是严格对角占优的，此方程组有唯一解，雅可比迭代收敛 请注意：雅可比迭代次数已经超过最大迭代次数 max1 第四章 解线性方程组的迭代法的 MATLAB 程序52. 高等教育出版社 教育电子音像出版社 作者：任玉杰 a = -8 -8 -1jX = 1.1000 1.2000 1.3000X = 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
14、,b,X0,inf, 0.001,5)运行后输出结果请注意：系数矩阵 A 不是严格对角占优的，此雅可比迭代不一定收敛请注意：雅可比迭代次数已经超过最大迭代次数 max1 a = -8.0000 -8.0000 3.5000jX = 24.5000 24.6000 106.6000X = 5.5490 5.6490 27.6553由（1）和（2）可见，如果系数矩阵 是严格对角占优的，则雅可比迭代收敛的A速度快；如果系数矩阵 不是严格对角占优的，则雅可比迭代收敛的速度慢.因此， 的值越小，雅可比迭代收敛的速度越快. knkijjaa1 ),21(n（二）利用雅可比迭代定义编写的解线性方程组的 MA
15、TLAB 程序利用雅可比迭代定义编写解线性方程组（4.5）的 MATLAB 程序的一般步骤步骤 1 将线性方程组（4.5）的系数矩阵 分解为 ，其中AULD,Dnnaadig 2121),(.L001,21nna U0,112n 在 MATLAB 工作窗口输入程序 A=a11 a12 a1n; a21 a22 a2n; an1 an2 ann;D=diag(A) U=triu(A,1), L=tril(A,-1)运行后即可输出 ；LDA，的步骤 2 若对角矩阵 非奇异（即 ，则（4.5）化为),1,0niai.bDXU1)(若记 .fB11,则方程组的迭代形式可写作1)()(kk )2,0(k
16、可以利用下面的 MATLAB 程序完成dD=det(D);if dD=0disp(请注意：因为对角矩阵 D 奇异，所以此方程组无解. )第四章 解线性方程组的迭代法的 MATLAB 程序53. 高等教育出版社 教育电子音像出版社 作者：任玉杰 elsedisp(请注意：因为对角矩阵D 非奇异，所以此方程组有解. )iD=inv(D); B1=iD*(L+U);f1=iD*b;for k=1:max1X= B1*X0+ f1; X0=X; djwcX=norm(X-X0,P);xdwcX=djwcX/(norm(X,P)+eps); X1=Ab;if (djwcXwucha)|(xdwcXwuc
17、ha)disp(请注意：雅可比迭代次数已经超过最大迭代次数max1 )endenda,X=X;jX=X1,4.3 高斯-塞德尔（Gauss-Seidel）迭代及其 MATLAB 程序4.3.3 高斯-塞德尔迭代两种 MATLAB 程序利用 MATLAB 程序和高斯-塞德尔迭代解线性方程组 的常用方法有两种，一bAX种方法是根据高斯-塞德尔迭代公式（4.16） 、 （4.17） 、定理 4.3 和公式（4.14）编写一个名为 gsdd.m 的 M 文件并保存，然后在 MATLAB 工作窗口输入对应的命令，执行此文件；另一种方法是根据高斯-塞德尔迭代的定义，利用提取对角矩阵和上、下三角矩阵的MAT
18、LAB 程序解线性方程组 .下面我们分别介绍这两种方法.bAX（一） 高斯-塞德尔迭代定义的 MATLAB 程序 1根据高斯-塞德尔迭代定义，现提供名为 gsdddy.m 的 M 文件如下：用高斯-塞德尔迭代定义解线性方程组 的 MATLAB 主程序 1b输入的量：线性方程组 的系数矩阵 和 b, 初始向量 X0,范数的名称 P = 1, A2, inf, 或 fro.，近似解 X 的误差（精度）wucha 和迭代的最大次数 max1.输出的量：以系数矩阵 的对角元构成的对角矩阵 D、A 的上三角形矩阵AnijaU，但对角元为 0、A 的下三角形矩阵 L，但对角元为 0 和有关高斯-塞德尔迭代
19、收敛性的相关信息及其 的精确解 jX 和近似解 X.bfunction 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=0disp(请注意：因为对角矩阵D 奇异，所以此方程组无解.)elsedisp(请注意：因为对角矩阵D 非奇异，所以此方程组有解. )iD=inv(D-L); B2=iD*U;f2=iD*b;jX=Ab; X=X0; n m=size(A);for k=1:max1X1= B2*X+f2; djwcX=norm(X1-X,P);xdwcX=djwcX
20、/(norm(X,P)+eps);if (djwcX 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 非奇异，所以此方程组有解 .请注意：高斯-塞德尔迭代收敛,此 A 的分解矩阵 D,U,L 和方程组的精确解 jX和近似解 X 如下： 此近似解与例 4.2.3 中的（1）中的解 =（24.073 8, 24.173 8, 104.797 4）XT比较，在相同的条件下, 高斯-塞德尔迭代比雅可比迭代得到的近似解的精度更高.（2）在 M
21、ATLAB 工作窗口输入程序 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.3036X = 1.0e+142 *0.2883 0.1062 0.3622 -3.1374（二） 高斯-塞德尔迭代公式的 MATL
22、AB 程序 2根据高斯-塞德尔迭代公式（4.16） 、 （4.17） 、定理 4.3 和公式（4.14） ，现提供名为 gsdd.m 的 M 文件如下用高斯-塞德尔迭代解线性方程组 的 MATLAB 主程序 2bAX输入的量：线性方程组 的系数矩阵 和 b, 初始向量 X0,范数的名称 P = 1, D =10.0000 0 00 10.0000 00 0 0.5000U =0 1 20 0 20 0 0L =0 0 01 0 01 1 0jX =24.5000 24.6000 106.6000X =24.4996 24.5996 106.5984第四章 解线性方程组的迭代法的 MATLAB
23、程序55. 高等教育出版社 教育电子音像出版社 作者：任玉杰 2, inf, 或 fro.，近似解 X 的误差（精度）wucha 和迭代的最大次数 max1.输出的量：系数矩阵 的 的值和有关Anijakkijja1 ),21(n高斯-塞德尔迭代收敛性的相关信息及其 的精确解 jX 和近似解 X.bXfunction X=gsdd(A,b,X0,P,wucha,max1)n m=size(A);for j=1:ma(j)=sum(abs(A(:,j)-2*(abs(A(j,j);endfor i=1:nif a(i)=0disp(请注意：系数矩阵A不是严格对角占优的，此高斯- 塞德尔迭代不一定
24、收敛)returnendendif a(i)wucha)jX=X1, 4.4 解方程组的超松弛迭代法及其 MATLAB 程序用雅可比迭代法和高斯-塞德尔迭代法解线性方程组时，有时收敛速度很慢，为了提高收敛速度，我们介绍超松弛迭代法，它是对高斯-塞德尔迭代的一种加速方法，适用于大型稀疏矩阵方程组的求解.4.4.2 超松弛迭代法收敛性及其 MATLAB 程序第四章 解线性方程组的迭代法的 MATLAB 程序56. 高等教育出版社 教育电子音像出版社 作者：任玉杰 根据定理 4.5 和谱半径定义，现提供名为 ddpbj.m 的 M 文件，用于 判别超松弛迭代公式（4.23）产生的迭代序列的敛散性.用
25、谱半径判别超松弛迭代法产生的迭代序列的敛散性的 MATLAB 主程序输入的量：线性方程组 的系数矩阵 和松弛因子 om；bAXA输出的量：矩阵 的所有特征值和谱半径 mH )1()(1DULDB及其有关超松弛迭代法产生的迭代序列的敛散性的相关信息.)(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=1disp(请注意：因为谱半径不小于1，所以超松弛迭代序列发散,谱半径mH和B的
26、所有的特征值H如下： )elsedisp(请注意：因为谱半径小于1，所以超松弛迭代序列收敛,谱半径mH和B的所有的特征值H如下： )endmH例 4.4.1 当取 =1.15，5 时，判别用超松弛迭代法解下列方程组产生的迭代序列是否收敛？ 372364181321xx解 （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 的所有的特征值
27、H 如下： mH = 0.1596H =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.1082H =-14.1082 -2.5107 0.5996 + 2.6206i 0.5996 - 2.6206i4.4.3 超松弛迭代法的 MATLAB 程序根据超松弛迭代公式（4.23）和定理 4.
28、5，现提供名为 cscdd.m 的 M 文件如下：用超松弛迭代法解线性方程组 的 MATLAB 主程序bAX输入的量：线性方程组 的系数矩阵 和 b, 初始向量 X,范数的名称 P = 1, 第四章 解线性方程组的迭代法的 MATLAB 程序57. 高等教育出版社 教育电子音像出版社 作者：任玉杰 2, inf, 或 fro.，松弛因子 om,近似解 X 的误差（精度）wucha 和迭代的最大次数 max1.输出的量：谱半径 mH，以系数矩阵 A 的对角元构成的对角矩阵 D、 A 的上三角形矩阵 U，但对角元为 0、A 的下三角形矩阵 L，但对角元为 0, 迭代次数 i,有关超松弛迭代收敛性的
29、相关信息及其 的精确解 jX 和近似解 X.bfunction X=cscdd (A,b,X,om,wucha,max1)D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1); jX=Ab;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:max1iD=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=dj
30、wcX/(norm(X,inf)+eps);if (djwcX max1disp(迭代次数已经超过最大迭代次数max1, 谱半径mH ，方程组的精确解jX, 迭代次数i如下： )mH,D,U,L,jX=jX, i=k-1,endendendif mH=1disp(请注意：因为谱半径不小于1，所以超松弛迭代序列发散.)disp(谱半径mH，A的分解矩阵 D,U,L和方程组的精确解jX,迭代次数i和迭代序列X 如下： )i=k-1,mH,D,U,L,jX,elsedisp(因为谱半径小于1，所以超松弛迭代序列收敛，近似解X如下： )end或function X=cscdd1 (A,b,X,om,w
31、ucha,max1)D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1); jX=Ab;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:max1iD=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);endif mH=1disp(请注意：因为谱半径不小于1
32、，所以超松弛迭代序列发散.谱半径mH，A的分解矩阵D,U,L和方程组的精确解jX和近似解X如下：)elsedisp(请注意：因为谱半径小于1，所以超松弛迭代序列收敛.)if (djwcX 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.1596D =5 0 0 00 8 0 00 0 -4 00 0 0 7U =0 -1 1 20 0 -1 -30
33、 0 0 10 0 0 0L =0 0 0 0-2 0 0 0-1 2 0 01 -3 -2 0jX =0.4491 0.2096 -1.4850 -0.0299i =3因为谱半径小于 1，所以超松弛迭代序列收敛，近似解 X 如下： X =0.44840.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),运行后输出结果如下：请注意：
34、因为谱半径不小于 1，所以超松弛迭代序列发散.谱半径 mH，A 的分解矩阵 D,U,L 和方程组的精确解 jX,迭代次数 i 和迭代序列X 如下：i = mH =99 14.1082D =5 0 0 00 8 0 00 0 -4 00 0 0 7U =0 -1 1 2第四章 解线性方程组的迭代法的 MATLAB 程序59. 高等教育出版社 教育电子音像出版社 作者：任玉杰 0 0 -1 -30 0 0 10 0 0 0L =0 0 0 0-2 0 0 0-1 2 0 01 -3 -2 0jX = X =1.0e+114 *0.4491 -0.31220.2096 1.0497-1.4850 -3.7174-0.0299 3.9615
