迭代解法的matlab实现

1、逐次求出的A为低阶稠解线性方程组AX =b逐次求出的A为低阶稠适合于算中常会遇至吠型稀近似解逼近精适合于算中常会遇至吠型稀密矩阵（指n不大且元疏矩阵（指n很大且零兀较多）的方程组，迭代法在计算和存贮两方面都适合后一种情 况.因为迭代法是通过逐次迭代来逼近方程组的解，所以收敛性和收敛速度是构造迭代法 时应该注意的问题.另外，因为不同的系数矩阵具有不同的性态，所以大多数迭代方法都 具有一定的适用范围.有时，某种方法对于一类方程组迭代收敛，而对另一类方程组迭代 时就发散.因此，我们应该学会针对具有不同性质的线性方程组构造不同的迭代4.1迭代法和敛散性及其MATLAB程序4.1.2迭代法敛散性的判别及

2、其MATLAB程序根据定理4.1和谱半径定义，现提供一个名为pddpb.m的M文件，用于判别迭代公 式（4.7）产生的迭代序列的敛散性.用谱半径判别迭代法产生的迭代序列的敛散性的MATLAB主程序输入的量：线性方程组AX = b的迭代公式（4.7）中的B ；输出的量：矩阵B的所有特征值和谱半径mH = P （B）及其有关迭代法产生的迭代序 列的敛散性的相关信息.function H=ddpbj（B）H=eig（B）;mH=norm（H,inf）;if mH=1disp（请注意：因为谱半径不小于1,所以迭代序列发散,谱半径mH和B的所 有的特征值H如下：）elsedisp（请注意：因为谱半径小于

3、1，所以迭代序列收敛,谱半径mH和B的所有 的特征值H如下：）endmH4.1.3与迭代法有关的MATLAB命令（一）提取（产生）对角矩阵和特征值可以用表4 - 1的MATLAB命令提取对角矩阵和特征值.表4-1提取（产生）对角矩阵和特征值MATLAB命令功能DX=diag(X)若输入向量X，则输出DX是以X为对角元的对角矩阵； 若输入矩阵X，则输出DX是以X的对角元构成的向量；DX=diag(diag(X)输入矩阵X，输出DX是以X的对角元构成的对角矩阵， 可用于迭代法中从A中提取D.lm=eig(A)输入矩阵A，输出lm是A的所有特征值.（二）提取（产生）上（下）三角形矩阵可以用表4-2的

4、MATLAB命令提取矩阵的上三角形矩阵和下三角形矩阵.表4-2 提取矩阵的上三角形矩阵和下三角形矩阵MATLAB命令功能U=triu(A)输入矩阵A，输出U是A的上三角形矩阵.L=tril(A)输入矩阵A，输出L是A的下三角形矩阵.U=triu(A,1)输入矩阵A，输出U是A的上三角形矩阵，但对角元为0，可 用于迭代法中从A中提取U .L=tril(A,-1)输入矩阵A，输出L是A的下三角形矩阵，但对角元为0，可 用于迭代法中从A中提取L.（三）稀疏矩阵的处理对稀疏矩阵在存贮和运算上的特殊处理，是MATLAB进行大规模科学计算时的特点和优势之一.可以用表4-3的MATLAB命令，输入稀疏矩阵的

5、非零元（零元不必输入），即可进行运算.表4-3 稀疏矩阵的MATLAB命令MATLAB命令功能ZA=sparse(r,c,v,m,n)表示在第r仃、第c列输入数值力 矩阵共m仃n 列，输出ZA，给出Sc）及V为一稀疏矩阵.MA=full(ZA)输入稀疏矩阵Z4,输出为满矩阵M4 （包含零元）4.2雅可比（Jacobi）迭代及其MATLAB程序4.2.2雅可比迭代的收敛性及其MATLAB程序根据定理4.3和公式（4.14），现提供一个名为jspb.m的M文件如下:判别雅可比迭代收敛性的MATLAB主程序输入的量：线性方程组AX=b的系数矩阵A；输出的量：系数矩阵A =(.)万nxnj输出的量：系

6、数矩阵A =(.)万nxnj=ii丰kkjkk(k = 1, 2,n）的值和有关雅可比迭代收敛性的相关信息.function a=jspb(A)n m=size(A);or j=1:ma(j)=sum(abs(A(:,j)-2*(abs(A(j,j);end for 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）在MAT

7、LAB工作窗口输入程序 A=10 -1 -2;-1 10 -2;-1 -1 0.5;a=jspb（A）运行后输出结果请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛a =-8.0000e+000 -8.0000e+000 3.5000e+0004.2.3雅可比迭代的两种MATLAB程序利用MATLAB程序和雅可比迭代解线性方程组AX = b的常用的方法有两种，一种方法是根据雅可比迭代公式(4.11)、(4.12)式、定理4.3和公式(4.14)编写一个名为jacdd.m的M文件并保存，然后在MATLAB工作窗口输入对应的命令，执行此文件；另一 种方法是根据雅可比迭代的定义，利用提取

8、对角矩阵和上、下三角矩阵的MATLAB程序解 线性方程组AX = b .下面我们分别介绍这两种方法.(一)雅可比迭代公式的MATLAB程序用雅可比迭代解线性方程组AX = b的MATLAB主程序输入的量：线性方程组AX = b的系数矩阵A和b,初始向量X0范数的名称P= 1,才aakk(k = 1, 2, , n)的值和有关雅万nxn2, inf或fro.，近似解X的误差(精度)wucha和迭代的最大次数maxi； (k = 1, 2, , n)的值和有关雅万nxnj=ii丰k可比迭代收敛性的相关信息及其AX = b的精确解jX和近似解X.根据雅可比迭代(4.11)、(4.12)式、定理4.3

9、和公式(4.14)，现提供一个名为jacdd.m 的M文件如下：function X=jacdd(A,b,X0,P,wucha,max1) n m=size(A);or j=1:ma(j)=sum(abs(A(:,j)-2*(abs(A(j,j);end for i=1:nif a(i)=0disp(请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛)returnendendif a(i)0disp(请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅 可比迭代收敛)endfor k=1:max1kfor j=1:mX(j)=(b(j)-A(j,1:j-1,j+1:m)*X0(

10、1: j-1,j+1:m)/A(j,j);endX,djwcX=norm(X-X0,P); xdwcX=djwcX/(norm(X,P)+eps); X0=X;X1=Ab;if (djwcXwucha)&(xdwcXwucha)&(xdwcXwucha)disp(请注意：雅可比迭代次数已经超过最大迭代次数max1 )enda,X=X;jX=X1,例4.2.3用8范数和判别雅可比迭代的MATLAB主程序解例4.2.2中的方程组，解 的精度为0.001，分别取最大迭代次数max1=100, 5,初始向量X0=(0 0 0)T，并比较 它们的收敛速度.解 (1)取最大迭代次数max1=100时.首先

11、保存名为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-1jX =1.10001.20001.3000X =1.09941.19941.2993在MATLAB工作窗口输入程序 A=10 -1 -2;-1 10 -2;-1 -1 0.5； b=7.2;8.3；4

12、.2; X0=0 0 0;X=jacdd（A,b,X0,inf, 0.001,100）运行后输出结果请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛请注意：雅可比迭代收敛，此方程组的精确解jX和近似解X如下： a =-8.0000-8.00003.5000jX =24.500024.6000 106.6000X = 24.073824.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,

13、5） 运行后输出结果请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，雅可比迭代收敛 请注意：雅可比迭代次数已经超过最大迭代次数max1 a =-8-8-1jX = 1.10001.20001.3000X =1.09511.19511.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不是严格对角占优的，此雅可比迭代不一定收敛 请注意：雅可比迭代次数已经超过最大迭代次数max1a =-8.0000-8.

14、00003.5000jX =24.500024.6000106.6000X =5.54905.649027.6553a = fj=1 i壬a = fj=1 i壬k-|akJ （k = 1,2, , n）的值越小，雅可比迭代收敛的速度越快.（二）利用雅可比迭代定义编写的解线性方程组的MATLAB程序利用雅可比迭代定义编写解线性方程组（4.5）的MATLAB程序的一般步骤步骤1将线性方程组（4.5）的系数矩阵A分解为A = D - L-U，其中(0)(0 a a )121na 00:21,U =: I. I. an-1,na a -a 00 A=a11 a12 -a1n; a21 a22 a2n;

15、an1 D=diag(A) U=triu(A,1), L=tril(A,-1)运行后即可输出A的D,L, U ；步骤2若对角矩阵D非奇异(即a,。0, i = 1,n)，则(4.5)化为X = D -1(L + U) X + D-ib .若记B = D -1 (L + U), f = D -1b .则方程组的迭代形式可写作11X (k+1) = B X (k) + f(k = 0,1,2 )可以利用下面的MATLAB程序完成1dD=det(D);ifdD=0disp(请注意：因为对角矩阵D奇异，所以此方程组无解.)elsedisp(请注意：因为对角矩阵D非奇异，所以此方程组有解.)iD=inv

16、(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)&(xdwcXwucha)|(xdwcXwucha)disp(请注意：雅可比迭代次数已经超过最大迭代次数max1 )endenda,X=X;jX=X1,4.3高斯-塞德尔(Gauss-Seidel)迭代及其MATLAB程序4.3.3高斯-塞德尔迭代两种MATLAB程序利用MATLAB程序和高斯-塞德尔迭代解线性方程组AX = b的常用方法有两种

17、，一 种方法是根据高斯-塞德尔迭代公式(4.16)、(4.17)、定理4.3和公式(4.14)编写一 个名为gsdd.m的M文件并保存，然后在MATLAB工作窗口输入对应的命令，执行此文件; 另一种方法是根据高斯-塞德尔迭代的定义，利用提取对角矩阵和上、下三角矩阵的 MATLAB程序解线性方程组AX = b.下面我们分别介绍这两种方法.(一)高斯-塞德尔迭代定义的MATLAB程序1根据高斯-塞德尔迭代定义，现提供名为gsdddy.m的M文件如下：用高斯-塞德尔迭代定义解线性方程组AX = b的MATLAB主程序1输入的量:线性方程组AX = b的系数矩阵A和b,初始向量X0,范数的名称P =

18、1, 2,inf,或fro,近似解X的误差(精呷 wucha和迭代的最大次数maxi.输出的量：以系数矩阵A =L. 的对角元构成的对角矩阵D. A的上三角形矩阵 nxnU,但对角元为0、A的下三角形矩阵L,但对角元为0和有关高斯-塞德尔迭代收敛性的 相关信息及其AX = b的精确解jX和近似解X.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=0disp(请注意：因为对角矩阵D奇异，所以此方程组无解.)elsedisp(请注意：因为对角矩阵D非奇异

19、，所以此方程组有解.) 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/(norm(X,P)+eps);if (djwcXwucha)|(xdwcXwucha)returnelsek,X1,k=k+1;X=X1;endendif (djwcXwucha)|(xdwcX 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.0

20、01,100)运行后输出结果请注意：因为对角矩阵D非奇异，所以此方程组有解.请注意：高斯-塞德尔迭代收敛，此A的分解矩阵D,U,L和方程组的精确解jX和 近似解X如下：D =L =10.000000000010.00000100000.5000110U =jX =01224.500024.6000 106.6000002X =00024.499624.5996 106.5984此近似解与例4.2.3中的（1）中的解X = （24.0738, 24.1738, 104.7974） t比 较，在相同的条件下，高斯-塞德尔迭代比雅可比迭代得到的近似解的精度更高.（2）在MATLAB工作窗口输入程序

21、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.25710.72861.3036X = 1.0e+142 *0.28830.10620.3622-3.1374（二）高斯-塞德尔迭代公式的MATLAB程序2根据高斯-塞德尔迭代公式（4.16）、（4.17）、定

22、理4.3和公式（4.14），现提供名为gsdd.m的M文件如下用高斯-塞德尔迭代解线性方程组AX = b的MATLAB主程序2输入的量:线性方程组AX = b的系数矩阵A和b,初始向量X0,范数的名称P = 1,2, inf,或fro.，近似解X的误差（精度）wucha和迭代的最大次数maxi.才aakjkkj=ii丰k（k = 1, 2, , n）j=ii丰k（k = 1, 2, , n）的值和有关高万nxn斯-塞德尔迭代收敛性的相关信息及其AX = b的精确解/X和近似解X.function X=gsdd(A,b,X0,P,wucha,max1)n m=size(A);for j=1:ma

23、(j)=sum(abs(A(:,j)-2*(abs(A(j,j);endfor i=1:nif a(i)=0disp(请注意：系数矩阵A不是严格对角占优的，此高斯塞德尔迭代不定收敛)return endend if a(i)0disp(请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且高 斯-塞德尔迭代收敛)endfor k=1:max1for j=1:mif j=1X(1)=(b (1)-A(1,2:m)*X0(2:m)/A(1,1)endif j=mX(m)=(b(m)-A(m,1:M1)*X(1:M1)/A(m,m);endfor j=2:M1X(j)=(b(j)-A(j,1:j-

24、1)*X(1:j-1) -A(j,j+1:m)*X(j+1:m)/A(j,j);endenddjwcX=norm(X-X0,P);xdwcX=djwcX/(norm(X,P)+eps); X0=X;X1=Ab;if (djwcXwucha)|(xdwcXwucha)&(xdwcXwucha)disp(请注意：高斯-塞德尔迭代次数已经超过最大迭代次数maxi )enda,X=X;jX=X1,4.4解方程组的超松弛迭代法及其MATLAB程序用雅可比迭代法和高斯-塞德尔迭代法解线性方程组时，有时收敛速度很慢，为了提 高收敛速度，我们介绍超松弛迭代法，它是对高斯-塞德尔迭代的一种加速方法，适用于 大型

25、稀疏矩阵方程组的求解.4.4.2超松弛迭代法收敛性及其MATLAB程序根据定理4.5和谱半径定义，现提供名为ddpbj.m的M文件，用于判别超松弛迭代 公式(4.23)产生的迭代序列的敛散性.用谱半径判别超松弛迭代法产生的迭代序列的敛散性的MATLAB主程序输入的量：线性方程组AX = A的系数矩阵A和松弛因子om；输出的量：矩阵% = (-诚)-1U + (1 -)D 的所有特征值和谱半径 mH = P (B“)及其有关超松弛迭代法产生的迭代序列的敛散性的相关信息function H=ddpbj(A,om)D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1);

26、iD=inv(D-om*L);B2=iD*(om*U+(1-om)*D);H=eig(B2);mH=norm(H,inf);if mH=idisp(请注意：因为谱半径不小于1,所以超松弛迭代序列发散,谱半径mH 和B的所有的特征值H如下：)elsedisp(请注意：因为谱半径小于1,所以超松弛迭代序列收敛,谱半径mH 和B的所有的特征值H如下：)end mH例4.4.1当取=1.15, 5时，判别用超松弛迭代法解下列方程组产生的迭代序列是 否收敛？5 % + x 2 - x 3 一 2 x 4 = 42 x + 8 x + x + 3 x = 1 A=5 1 -1 -2;2 8 1 3;1 -

27、2 -4 -1;-1 3 2 7; H=ddpbj(A,1.15)运行后输出结果请注意：因为谱半径小于1，所以超松弛迭代序列收敛,谱半径mH和B的所有 的特征值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.6

28、206i 0.5996 - 2.6206i4.4.3超松弛迭代法的MATLAB程序根据超松弛迭代公式(4.23)和定理4.5,现提供名为cscdd.m的M文件如下：用超松弛迭代法解线性方程组AX = b的MATLAB主程序输入的量：线性方程组AX = b的系数矩阵A和b,初始向量X,范数的名称P = 1, 2, inf,或fro.，松弛因子om,近似解X的误差(精度)wucha和迭代的最大次数maxi.输出的量：谱半径mH，以系数矩阵A的对角元构成的对角矩阵D、A的上三角形矩 阵，但对角元为0、A的下三角形矩阵乙，但对角元为0,迭代次数z，有关超松弛迭代收 敛性的相关信息及其AX = b的精确

29、解jX和近似解X.function 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=djwcX/(norm(X,inf)+

30、eps);if (djwcXwucha)|(xdwcX 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,所以超松弛迭代序列收敛，近似解乂如下：)end或function X=cscdd1 (A,b,X,om,wucha,max1)D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1); jX=Ab;n m=size(A