matlab 解线性方程组的迭代法

function x=ak(a,b)%a为系数矩阵，b为初始向量（默认为零向量）
%e为精度（默认为1e-6），N为最大迭代次数（默认为100），x为返回解向量
n=length(b);
N=100;
e=1e-6;
x0=zeros(n,1)
%生成一n*1阶零矩阵
x=x0;
x0=x+2*e;
k=0;
d=diag(diag(a));
%生成一个除对角线上元素不为零外其他元素皆为零的矩阵d,且d对角线上的元素为矩阵a 对角线上的元素
l=-tril(a,-1);
%生成一个下三角矩阵
u=-triu(a,1);
%生成一个上三角矩阵
while norm(x0-x,inf)>e & k<N %norm(x0-x,inf)为矩阵（x0-x）的无穷范数
k=k+1;
x0=x;
x=inv(d)*(l+u)*x+inv(d)*b;%雅可比迭代公式k
disp(x')
end
if k==N warning('已达最大迭代次数'); end
function X=BDD(f,x0,TOL)
%X用来存储迭代过程所有的根；
%f是符合不动点迭代要求的迭代方程；%x0设定的迭代初值；
%TOL允许的误差值；
x=feval(f,x0);
n=1;
X(:,n)=x;
while abs(x-x0)>TOL
x0=x;
x=feval(f,x0);
n=n+1;
X(:,n)=x;
end。

jacobi迭代法matlabJacobi迭代法是一种常用的线性方程组求解方法，它是一种迭代法，通过不断迭代来逼近线性方程组的解。

Jacobi迭代法的基本思想是将线性方程组的系数矩阵分解为一个对角矩阵和一个非对角矩阵的和，然后通过迭代求解对角矩阵和非对角矩阵的乘积，最终得到线性方程组的解。

Jacobi迭代法的具体步骤如下：1. 将线性方程组的系数矩阵A分解为一个对角矩阵D和一个非对角矩阵R的和，即A=D+R。

2. 将线性方程组的右端向量b分解为一个对角矩阵D和一个非对角矩阵N的乘积，即b=Dx。

3. 对于任意的初始解向量x0，计算下一次迭代的解向量x1，即x1=D^(-1)(b-Rx0)。

4. 重复步骤3，直到达到预定的精度或迭代次数。

Jacobi迭代法的优点是简单易懂，易于实现，收敛速度较快。

但是，它的缺点也很明显，即收敛速度较慢，需要进行大量的迭代才能达到较高的精度。

在Matlab中，可以使用以下代码实现Jacobi迭代法：function [x,k]=jacobi(A,b,x0,tol,maxit)% Jacobi迭代法求解线性方程组Ax=b% 输入：系数矩阵A，右端向量b，初始解向量x0，精度tol，最大迭代次数maxit% 输出：解向量x，迭代次数kn=length(b); % 系数矩阵A的阶数D=diag(diag(A)); % 对角矩阵DR=A-D; % 非对角矩阵Rx=x0; % 初始解向量for k=1:maxitx1=D\(b-R*x); % 计算下一次迭代的解向量if norm(x1-x)<tol % 判断是否达到精度要求break;endx=x1; % 更新解向量end输出结果可以使用以下代码实现：A=[4 -1 0; -1 4 -1; 0 -1 4]; % 系数矩阵b=[15; 10; 10]; % 右端向量x0=[0; 0; 0]; % 初始解向量tol=1e-6; % 精度要求maxit=1000; % 最大迭代次数[x,k]=jacobi(A,b,x0,tol,maxit); % Jacobi迭代法求解线性方程组fprintf('解向量x=[%f; %f; %f]\n',x(1),x(2),x(3)); % 输出解向量fprintf('迭代次数k=%d\n',k); % 输出迭代次数以上就是Jacobi迭代法的主要内容，通过Matlab实现Jacobi迭代法可以更好地理解其基本思想和具体步骤。

线性方程组求解的迭代法及Matlab的应用2 定常迭代法 32.1 雅可比迭代法 42.2 高斯－赛德尔迭代法 42.3 超松弛迭代法 52.4 迭代的收敛性分析 62.5 实例 73 不定常迭代 93.1 最速下降法 93.2 共轭梯度法 103.3 实例 114 Matlab在定常迭代与不定常迭代中的应用 124.1 雅可比迭代法的程序 124.2 高斯-赛德尔迭代的程序 134.3 超松弛迭代的程序 134.4 最速下降法的程序 144.5 共轭梯度法的程序 144.6 Matlab实现的实例 154.6.1 定常迭代的收敛速度的比较 154.6.2 超松弛迭代法松弛因子的选择 164.6.3 不定常迭代的收敛速度的比较 18参考文献 20致谢 211 引言1.1 课题的目的和意义数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。

为计算数学的主体部分。

它的主要内容有插值法，函数逼近，曲线拟和，数值积分，数值微分，解线性方程组的直接方法，解线性方程组的迭代法，非线性方程求根，常微分方程的数值解法。

数值分析这门学科有如下特点：面向计算机；有可靠的理论分析；要有好的计算复杂性；要有数值实验；要对算法进行误差分析。

:在许多工程实际应用中，超大规模的线性方程组的数值解法是时常要遇到的问题。

由于线性方程组的维数巨大，给具体的计算带来很大的问题——算法对计算机的内存需求大，算法的收敛速度慢以及计算舍人误差的累积扩张。

这些往往使理论上较好的算法无法真正的应用到工程实际中，因此寻求一种真正能实际应用的数值算法一直是人们关注的问题。

通常求解线性方程组一般可以分为直接解法和迭代解法。

现在流行的算法一般采用迭代的算法来求解线性方程组，这主要是为了加快求解的速度。

另外由于计算机的发展，在许多领域里涌现了一些新型的算法如神经网络，遗传算法，粒子群算法，模拟退火算法以及蚁群算法等。

解线性方程组b AX =的迭代法是从初始解出发，根据设计好的步骤用逐次求出的近似解逼近精确解.在第三章中介绍的解线性方程组的直接方法一般适合于A 为低阶稠密矩阵（指n 不大且元多为非零）的情况，而在工程技术和科学计算中常会遇到大型稀疏矩阵（指n 很大且零元较多）的方程组，迭代法在计算和存贮两方面都适合后一种情况.由于迭代法是通过逐次迭代来逼近方程组的解，所以收敛性和收敛速度是构造迭代法时应该注意的问题.另外，因为不同的系数矩阵具有不同的性态，所以大多数迭代方法都具有一定的适用范围.有时，某种方法对于一类方程组迭代收敛，而对另一类方程组迭代时就发散.因此，我们应该学会针对具有不同性质的线性方程组构造不同的迭代.4.1 迭代法和敛散性及其MATLAB 程序4.1.2 迭代法敛散性的判别及其MATLAB 程序根据定理4.1和谱半径定义，现提供一个名为pddpb.m 的M 文件，用于判别迭代公用谱半径判别迭代法产生的迭代序列的敛散性的MATLAB 主程序 输入的量：线性方程组b AX =的迭代公式（4.7）中的B ； 输出的量：矩阵B 的所有特征值和谱半径mH )(B ρ=及其有关迭代法产生的迭代序列的敛散性的相关信息.function H=ddpbj(B)H=eig(B);mH=norm(H,inf); if mH>=1disp('请注意：因为谱半径不小于1，所以迭代序列发散,谱半径mH 和B 的所有的特征值H 如下：')elsedisp('请注意：因为谱半径小于1，所以迭代序列收敛,谱半径mH 和B 的所有的特征值H 如下：')end mH4.1.3 与迭代法有关的MATLAB 命令（一） 提取（产生）对角矩阵和特征值可以用表4–1的MATLAB 命令提取对角矩阵和特征值.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的MATLAB 命令提取矩阵的上三角形矩阵和下三角形矩阵.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 .(三)稀疏矩阵的处理对稀疏矩阵在存贮和运算上的特殊处理，是MA TLAB 进行大规模科学计算时的特点和优势之一.可以用表4–3的MATLAB 命令，输入稀疏矩阵的非零元（零元不必输入），即可进行运算.MATLAB 命令 功 能ZA=sparse(r,c,v,m,n)表示在第r 行、第c 列输入数值v ，矩阵共m 行n 列，输出ZA ，给出 (r , c ) 及v 为一稀疏矩阵.MA=full(ZA) 输入稀疏矩阵ZA ，输出为满矩阵MA （包含零元）4.2 雅可比（Jacobi ）迭代及其MATLAB 程序4.2.2 雅可比迭代的收敛性及其MATLAB 程序判别雅可比迭代收敛性的MATLAB 主程序输入的量：线性方程组b AX =的系数矩阵A ； 输出的量：系数矩阵=A ()nn ija ⨯的kk nki j kja a a -=∑≠=1 ),,2,1(n k =的值和有关雅可比迭代收敛性的相关信息.[n m]=size(A); for j=1:ma(j)=sum(abs(A(:,j)))-2*(abs(A(j,j))); end for i=1:n if a(i)>=0disp('请注意：系数矩阵A 不是严格对角占优的，此雅可比迭代不一定收敛')return end end if a(i)<0disp('请注意：系数矩阵A 是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛 ') end例4.2.2 用判别雅可比迭代收敛性的MATLAB 主程序，判别由下列方程组的雅可比迭代产生的序列是否收敛？（1）⎪⎩⎪⎨⎧=+--=-+-=--;2.45,3.8210,2.7210321321321x x x x x x x x x （2）⎪⎩⎪⎨⎧=+--=-+-=--.2.45.0,3.8210,2.7210321321321x x x x x x x x x 解 （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+0004.2.3 雅可比迭代的两种MATLAB 程序利用MATLAB 程序和雅可比迭代解线性方程组b AX =的常用的方法有两种，一种方法是根据雅可比迭代公式（4.11）、（4.12）式、定理4.3和公式（4.14）编写一个名为jacdd.m 的M 文件并保存，然后在MATLAB 工作窗口输入对应的命令，执行此文件；另一种方法是根据雅可比迭代的定义，利用提取对角矩阵和上、下三角矩阵的MATLAB 程序解线性方程组b AX =.下面我们分别介绍这两种方法.用雅可比迭代解线性方程组b AX =的MATLAB 主程序输入的量：线性方程组b AX =的系数矩阵A 和b , 初始向量X 0,范数的名称P = 1, 2, inf 或 'fro .'，近似解X 的误差（精度）wucha 和迭代的最大次数max1；输出的量：系数矩阵=A ()nn ija ⨯的kk nki j kja a a -=∑≠=1 ),,2,1(n k =的值和有关雅可比迭代收敛性的相关信息及其b AX =的精确解jX 和近似解X .的M 文件如下：function X=jacdd(A,b,X0,P,wucha,max1) [n m]=size(A); for j=1:ma(j)=sum(abs(A(:,j)))-2*(abs(A(j,j))); end for i=1:n if a(i)>=0disp('请注意：系数矩阵A 不是严格对角占优的，此雅可比迭代不一定收敛')return end end if a(i)<0disp('请注意：系数矩阵A 是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛 ') endfor k=1:max1kfor j=1:mX(j)=(b(j)-A(j,[1:j-1,j+1:m])*X0([1: j-1,j+1:m]))/A(j,j);endX,djwcX=norm(X'-X0,P); xdwcX=djwcX/(norm(X',P)+eps); X0=X';X1=A\b;if (djwcX<wucha)&(xdwcX<wucha)disp('请注意：雅可比迭代收敛,此方程组的精确解jX和近似解X如下：')returnendendif (djwcX>wucha)&(xdwcX>wucha)disp('请注意：雅可比迭代次数已经超过最大迭代次数max1 ')enda,X=X;jX=X1',例4.2.3用 范数和判别雅可比迭代的MATLAB主程序解例4.2.2 中的方程组，解的精度为0.001，分别取最大迭代次数max1=100，5，初始向量X0=（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 -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]; 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 0 0]'; X=jacdd(A,b,X0,inf,0.001,5)运行后输出结果请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，雅可比迭代收敛请注意：雅可比迭代次数已经超过最大迭代次数max1a =-8 -8 -1jX =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 不是严格对角占优的，则雅可比迭代收敛的速度慢.因此，kk nki j kj a a a -=∑≠=1 ),,2,1(n k =的值越小，雅可比迭代收敛的速度越快.（二）利用雅可比迭代定义编写的解线性方程组的MATLAB 程序利用雅可比迭代定义编写解线性方程组（4.5）的MATLAB 程序的一般步骤 步骤1 将线性方程组（4.5）的系数矩阵A 分解为U L D A --=，其中=D ⎪⎪⎪⎪⎪⎭⎫⎝⎛=nn nn a a a a a a diag22112211),,(, -=L⎪⎪⎪⎪⎪⎭⎫⎝⎛-0001,2121n n n n a a aa -=U ,⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-000,1112nn n a a a . 在MATLAB 工作窗口输入程序>> A=[a11 a12 …a1n; a21 a22 …a2n;…; an1 an2 …ann;]; D=diag(A) U=triu(A,1), L=tril(A,-1)运行后即可输出U L D A ，，的； 步骤2 若对角矩阵D 非奇异（即),,1,0n i a ii =≠，则（4.5）化为 b D X U L D X 11)(--++=.若记b D f U L D B 1111),(--=+=.则方程组的迭代形式可写作1)(1)1(f X B Xk k +=+ )2,1,0( =k 可以利用下面的MATLAB 程序完成>>dD=det(D); if dD==0disp('请注意：因为对角矩阵D 奇异，所以此方程组无解.') 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=A\b;if (djwcX<wucha)&(xdwcX<wucha)disp('请注意：雅可比迭代收敛,此方程组的精确解jX和近似解X如下：')returnendendif (djwcX>wucha)|(xdwcX>wucha)disp('请注意：雅可比迭代次数已经超过最大迭代次数max1 ')endenda,X=X;jX=X1',4.3 高斯-塞德尔（Gauss-Seidel）迭代及其MATLAB程序4.3.3 高斯-塞德尔迭代两种MATLAB程序AX=的常用方法有两种，一利用MATLAB程序和高斯-塞德尔迭代解线性方程组b种方法是根据高斯-塞德尔迭代公式（4.16）、（4.17）、定理4.3和公式（4.14）编写一个名为gsdd.m的M文件并保存，然后在MATLAB工作窗口输入对应的命令，执行此文件；另一种方法是根据高斯-塞德尔迭代的定义，利用提取对角矩阵和上、下三角矩阵的AX=.下面我们分别介绍这两种方法.MATLAB程序解线性方程组b（一）高斯-塞德尔迭代定义的MATLAB程序1AX=的MATLAB主程序1 用高斯-塞德尔迭代定义解线性方程组bAX=的系数矩阵A和b, 初始向量X0,范数的名称P = 1, 2, 输入的量：线性方程组binf, 或'fro.'，近似解X的误差（精度）wucha和迭代的最大次数max1.A()n n ij a⨯的对角元构成的对角矩阵D、A的上三角形矩阵输出的量：以系数矩阵=U，但对角元为0、A的下三角形矩阵L，但对角元为0和有关高斯-塞德尔迭代收敛性的AX=的精确解jX和近似解X.相关信息及其bD=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=A\b; 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 (djwcX<wucha)|(xdwcX<wucha)returnelsek,X1',k=k+1;X=X1;endendif (djwcX<wucha)|(xdwcX<wucha)disp('请注意：高斯-塞德尔迭代收敛,此A的分解矩阵D,U,L和方程组的精确解jX和近似解X如下：')elsedisp('请注意：高斯-塞德尔迭代发散，并且迭代次数已经超过最大迭代次数max1,方程组的精确解jX和迭代向量X如下：')X=X';jX=jX' end endX=X';D,U,L,jX=jX'例4.3.3 用高斯-塞德尔迭代定义的MATLAB 主程序解下列线性方程组,取初始值)0,0,0(),,()0(3)0(2)0(1=x x x ，要求当3)()1(10-∞+<-k k x x 时，迭代终止.（1）⎪⎩⎪⎨⎧=+--=-+-=--.2.45.0,3.8210,2.7210321321321x x x x x x x x x （2）⎪⎪⎩⎪⎪⎨⎧=++--=+-+=-+-=+-+.2132127,11613514,22382,575434321432143214321x x x x x x x x x x x x x x x x解 （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 非奇异，所以此方程组有解.请注意：高斯-塞德尔迭代收敛,此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.3036X = 1.0e+142 *0.2883 0.1062 0.3622 -3.1374（二） 高斯-塞德尔迭代公式的MATLAB 程序2 根据高斯-塞德尔迭代公式（4.16）、（4.17）、定理4.3和公式（4.14），现提供名为用高斯-塞德尔迭代解线性方程组b AX =的MATLAB 主程序2 输入的量：线性方程组b AX =的系数矩阵A 和b , 初始向量X 0,范数的名称P = 1, 2, inf, 或 'fro.'，近似解X 的误差（精度）wucha 和迭代的最大次数max1.输出的量：系数矩阵=A ()nn ija ⨯的kk nki j kja a a -=∑≠=1 ),,2,1(n k =的值和有关D =10.0000 0 0 0 10.0000 0 0 0 0.5000 U =0 1 2 0 0 2 0 0 0L =0 0 0 1 0 0 1 1 0 jX =24.5000 24.6000 106.6000 X =24.4996 24.5996 106.5984AX=的精确解jX和近似解X.高斯-塞德尔迭代收敛性的相关信息及其b[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)<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-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=A\b;if (djwcX<wucha)|(xdwcX<wucha)disp('请注意：高斯-塞德尔迭代收敛,此方程组的精确解jX和近似解X 如下：')returnendendif (djwcX>wucha)&(xdwcX>wucha)disp('请注意：高斯-塞德尔迭代次数已经超过最大迭代次数max1 ') enda,X=X;jX=X1',4.4 解方程组的超松弛迭代法及其MATLAB程序用雅可比迭代法和高斯-塞德尔迭代法解线性方程组时，有时收敛速度很慢，为了提高收敛速度，我们介绍超松弛迭代法，它是对高斯-塞德尔迭代的一种加速方法，适用于大型稀疏矩阵方程组的求解.4.4.2 超松弛迭代法收敛性及其MATLAB程序根据定理4.5和谱半径定义，现提供名为ddpbj.m的M文件，用于判别超松弛迭代用谱半径判别超松弛迭代法产生的迭代序列的敛散性的MATLAB主程序AX=的系数矩阵A和松弛因子om；输入的量：线性方程组b输出的量：矩阵])1([)(1D U L D B ωωωω-+-=-的所有特征值和谱半径mH)(ωρB =及其有关超松弛迭代法产生的迭代序列的敛散性的相关信息.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 的所有的特征值H 如下：') elsedisp('请注意：因为谱半径小于1，所以超松弛迭代序列收敛,谱半径mH和B 的所有的特征值H 如下：') end mH例4.4.1 当取ω=1.15，5时，判别用超松弛迭代法解下列方程组产生的迭代序列是否收敛？⎪⎪⎩⎪⎪⎨⎧-=+++-=---=+++=--+372364213824254321432143214321x x x x x x x x x x x x x x x x 解 （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.6206i4.4.3 超松弛迭代法的MATLAB 程序用超松弛迭代法解线性方程组b AX =的MATLAB 主程序输入的量：线性方程组b AX =的系数矩阵A 和b , 初始向量X ,范数的名称P = 1, 2, inf, 或 'fro.'，松弛因子om ,近似解X 的误差（精度）wucha 和迭代的最大次数max1.输出的量：谱半径mH ，以系数矩阵A 的对角元构成的对角矩阵D 、A 的上三角形矩阵U ，但对角元为0、A 的下三角形矩阵L ，但对角元为0, 迭代次数i ,有关超松弛迭代收敛性的相关信息及其b AX =的精确解jX 和近似解X .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: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);if (djwcX<wucha)|(xdwcX<wucha)disp('谱半径mH，A的分解矩阵D,U,L和方程组的精确解jX,迭代次数i如下：')mH,D,U,L,jX=jX', i=k-1,returnif i> 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,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: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，所以超松弛迭代序列发散.谱半径mH，A的分解矩阵D,U,L和方程组的精确解jX和近似解X如下：')elsedisp('请注意：因为谱半径小于1，所以超松弛迭代序列收敛.')if (djwcX<wucha)|(xdwcX<wucha)disp('谱半径mH，A的分解矩阵D,U,L和方程组的精确解jX和近似解X 如下：')mH,D,U,L,jX=jX',elsedisp('迭代次数已经超过最大迭代次数max1,谱半径mH，方程组的精确解jX和迭代向量X如下：')mH,D,U,L,X=X1';jX=jX'returnendend例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.1596D =5 0 0 00 8 0 00 0 -4 00 0 0 7U =0 -1 1 20 0 -1 -30 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),运行后输出结果如下：请注意：因为谱半径不小于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 20 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。

基于Matlab的解线性方程组的几种迭代法的实现及比较线性方程组的解法有很多种，其中一类常用的方法是迭代法。

迭代法根据一个初值逐步逼近方程组的解，在每一次迭代中利用现有的信息产生新的近似值，并不断地修正。

下面介绍基于Matlab的三种迭代法：雅可比迭代法、高斯-赛德尔迭代法和超松弛迭代法，并进行比较。

1. 雅可比迭代法雅可比迭代法是迭代法中最简单的一种方法。

对于线性方程组Ax=b，雅可比迭代法的迭代公式为：x_{i+1}(j)=1/a_{jj}(b_j-\\sum_{k=1,k\eq j}^n a_{jk}x_i(k))其中，i表示迭代次数，j表示未知数的下标，x_i表示第i次迭代的近似解，a_{jk}表示系数矩阵A的第j行第k列元素，b_j 表示方程组的常数项第j项。

在Matlab中，可以使用以下代码实现雅可比迭代：function [x,flag]=jacobi(A,b,X0,tol,kmax)n=length(b);x=X0;for k=1:kmaxfor i=1:nx(i)=(b(i)-A(i,:)*x+A(i,i)*x(i))/A(i,i);endif norm(A*x-b)<tolflag=1;returnendendflag=0;return其中，参数A为系数矩阵，b为常数项列向量，X0为初值列向量，tol为迭代误差容许值（默认为1e-6），kmax为最大迭代次数（默认为1000）。

函数返回值x为近似解列向量，flag表示是否满足容许误差要求。

2. 高斯-赛德尔迭代法高斯-赛德尔迭代法是雅可比迭代法的改进。

其基本思想是，每次迭代时，利用已经求出的新解中的信息来更新其他未知数的值。

迭代公式为：x_{i+1}(j)=(1/a_{jj})(b_j-\\sum_{k=1}^{j-1}a_{jk}x_{i+1}(k)-\\sum_{k=j+1}^n a_{jk}x_i(k))与雅可比迭代法相比，高斯-赛德尔迭代法的每一次迭代都利用了前面已求得的近似解，因此可以更快地收敛。

matlab逐次超松弛迭代法
逐次超松弛迭代法（Gauss-Seidel Overrelaxation Method）
是一种用于求解线性方程组的数值方法，常用于解决大型稀疏矩阵
的方程组。

在MATLAB中，可以通过编写逐次超松弛迭代法的代码来
实现该算法。

首先，让我们回顾一下逐次超松弛迭代法的基本原理。

该方法
是基于迭代的思想，通过不断迭代计算得到线性方程组的近似解。

在每一次迭代中，通过更新当前解向量的各个分量来逐步逼近方程
组的精确解。

逐次超松弛迭代法引入了松弛因子，可以加速收敛速度。

在MATLAB中，可以使用以下步骤来实现逐次超松弛迭代法：
1. 首先，编写一个函数来表示线性方程组的系数矩阵和右侧向量。

这个函数可以接受系数矩阵、右侧向量和当前解向量作为输入，并返回更新后的解向量。

2. 接下来，编写主程序来调用这个函数，并设置迭代的终止条件。

可以选择设置最大迭代次数或者设定一个收敛精度作为终止条
件。

3. 在主程序中，使用一个循环来进行迭代计算，直到满足设定的终止条件为止。

在每一次迭代中，调用上述编写的函数来更新解向量。

4. 最后，输出得到的近似解向量作为结果。

需要注意的是，逐次超松弛迭代法的收敛性与松弛因子的选择有关，通常需要根据具体的线性方程组进行调整。

总之，在MATLAB中实现逐次超松弛迭代法需要编写系数矩阵和右侧向量的函数以及主程序来进行迭代计算，并且需要注意收敛性和松弛因子的选择。

希望这个回答能够帮助你更好地理解和实现逐次超松弛迭代法。

MATLAB代码解线性方程组的迭代法解线性方程组的迭代法1.rs里查森迭代法求线性方程组Ax=b的解function[x,n]=rs(A,b,x0,eps,M)if(nargin==3)eps=1.0e-6;%eps表示迭代精度M=10000;%M表示迭代步数的限制值elseif(nargin==4) M=10000;endI=eye(size(A));n=0;x=x0;tol=1;%迭代过程while(tol>eps)x=(I-A)*x0+b;n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x;if(n>=M)disp('Warning:迭代次数太多，可能不收敛！'); return;endend2.crs里查森参数迭代法求线性方程组Ax=b的解function[x,n]=crs(A,b,x0,w,eps,M)if(nargin==4)eps=1.0e-6;%eps表示迭代精度M=10000;%M表示迭代步数的限制值elseif(nargin==5)M=10000;endI=eye(size(A));n=0;x=x0;tol=1;%迭代过程while(tol>eps)x=(I-w*A)*x0+w*b;n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x;if(n>=M)disp('Warning:迭代次数太多，可能不收敛！'); return;endend3.grs里查森迭代法求线性方程组Ax=b的解function[x,n]=grs(A,b,x0,W,eps,M)if(nargin==4)eps=1.0e-6;%eps表示迭代精度M=10000;%M表示迭代步数的限制值elseif(nargin==5)M=10000;endI=eye(size(A));n=0;x=x0;tol=1;%前后两次迭代结果误差%迭代过程while(tol>eps)x=(I-W*A)*x0+W*b;%迭代公式n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x;if(n>=M)disp('Warning:迭代次数太多，可能不收敛！'); return;endend4.jacobi雅可比迭代法求线性方程组Ax=b的解function[x,n]=jacobi(A,b,x0,eps,varargin)if nargin==3eps=1.0e-6;M=200;elseif nargin<3errorreturnelseif nargin==5M=varargin{1};endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B=D\(L+U);f=D\b;x=B*x0+f;n=1;%迭代次数while norm(x-x0)>=epsx0=x;x=B*x0+f;n=n+1;if(n>=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend5.gauseidel高斯-赛德尔迭代法求线性方程组Ax=b的解function[x,n]=gauseidel(A,b,x0,eps,M)if nargin==3eps=1.0e-6;M=200;elseif nargin==4M=200;elseif nargin<3errorreturn;endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵G=(D-L)\U;f=(D-L)\b;x=G*x0+f;n=1;%迭代次数while norm(x-x0)>=epsx0=x;x=G*x0+f;n=n+1;if(n>=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend6.SOR超松弛迭代法求线性方程组Ax=b的解function[x,n]=SOR(A,b,x0,w,eps,M)if nargin==4eps=1.0e-6;M=200;elseif nargin<4errorreturnelseif nargin==5M=200;endif(w<=0||w>=2)error;return;endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B=inv(D-L*w)*((1-w)*D+w*U);f=w*inv((D-L*w))*b;x=B*x0+f;n=1;%迭代次数while norm(x-x0)>=epsx0=x;x=B*x0+f;n=n+1;if(n>=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend7.SSOR对称逐次超松弛迭代法求线性方程组Ax=b的解function[x,n]=SSOR(A,b,x0,w,eps,M)if nargin==4eps=1.0e-6;M=200;elseif nargin<4errorreturnelseif nargin==5M=200;endif(w<=0||w>=2)error;return;endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B1=inv(D-L*w)*((1-w)*D+w*U);B2=inv(D-U*w)*((1-w)*D+w*L);f1=w*inv((D-L*w))*b;f2=w*inv((D-U*w))*b;x12=B1*x0+f1;x=B2*x12+f2;n=1;%迭代次数while norm(x-x0)>=epsx0=x;x12=B1*x0+f1;x=B2*x12+f2;n=n+1;if(n>=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend8.JOR雅可比超松弛迭代法求线性方程组Ax=b的解function[x,n]=JOR(A,b,x0,w,eps,M)if nargin==4eps=1.0e-6;M=10000;elseif nargin==5M=10000;endif(w<=0||w>=2)%收敛条件要求error;return;endD=diag(diag(A));%求A的对角矩阵B=w*inv(D);%迭代过程x=x0;n=0;%迭代次数tol=1;%迭代过程while tol>=epsx=x0-B*(A*x0-b);n=n+1;tol=norm(x-x0);x0=x;if(n>=M)disp('Warning:迭代次数太多，可能不收敛！'); return;endend9.twostep两步迭代法求线性方程组Ax=b的解function[x,n]=twostep(A,b,x0,eps,varargin) if nargin==3eps=1.0e-6;M=200;elseif nargin<3errorreturnelseif nargin==5M=varargin{1};endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B1=(D-L)\U;B2=(D-U)\L;f1=(D-L)\b;f2=(D-U)\b;x12=B1*x0+f1;x=B2*x12+f2;n=1;%迭代次数while norm(x-x0)>=epsx0=x;x12=B1*x0+f1;x=B2*x12+f2;n=n+1;if(n>=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend10.fastdown最速下降法求线性方程组Ax=b的解function[x,n]=fastdown(A,b,x0,eps)if(nargin==3)eps=1.0e-6;endx=x0;n=0;tol=1;while(tol>eps)%以下过程可参考算法流程r=b-A*x0;d=dot(r,r)/dot(A*r,r);x=x0+d*r;tol=norm(x-x0);x0=x;n=n+1;end11.conjgrad共轭梯度法求线性方程组Ax=b的解function[x,n]=conjgrad(A,b,x0)r1=b-A*x0;p=r1;n=0;for i=1:rank(A)%以下过程可参考算法流程if(dot(p,A*p)< 1.0e-50)%循环结束条件break;endalpha=dot(r1,r1)/dot(p,A*p);x=x0+alpha*p;r2=r1-alpha*A*p;if(r2< 1.0e-50)%循环结束条件break;endbelta=dot(r2,r2)/dot(r1,r1);p=r2+belta*p;n=n+1;end12.preconjgrad预处理共轭梯度法求线性方程组Ax=b的解function[x,n]=preconjgrad(A,b,x0,M,eps)if nargin==4eps=1.0e-6;endr1=b-A*x0;iM=inv(M);z1=iM*r1;p=z1;n=0;tol=1;while tol>=epsalpha=dot(r1,z1)/dot(p,A*p);x=x0+alpha*p;r2=r1-alpha*A*p;z2=iM*r2;belta=dot(r2,z2)/dot(r1,z1);p=z2+belta*p;n=n+1;tol=norm(x-x0);x0=x;%更新迭代值r1=r2;z1=z2;end13.BJ块雅克比迭代法求线性方程组Ax=b的解function[x,N]=BJ(A,b,x0,d,eps,M)if nargin==4eps=1.0e-6;M=10000;elseif nargin<4errorreturnelseif nargin==5M=10000;%参数的默认值endNS=size(A);n=NS(1,1);if(sum(d)~=n)disp('分块错误！');return;endbnum=length(d);bs=ones(bnum,1);for i=1:(bnum-1)bs(i+1,1)=sum(d(1:i))+1;%获得对角线上每个分块矩阵元素索引的起始值endDB=zeros(n,n);for i=1:bnumendb=bs(i,1)+d(i,1)-1;DB(bs(i,1):endb,bs(i,1):endb)=A(bs(i,1):endb,bs(i,1):endb);%求A的对角分块矩阵endfor i=1:bnumendb=bs(i,1)+d(i,1)-1;invDB(bs(i,1):endb,bs(i,1):endb)=inv(DB(bs(i,1):endb,bs(i,1):e ndb));%求A的对角分块矩阵的逆矩阵endN=0;tol=1;while tol>=epsx=invDB*(DB-A)*x0+invDB*b;%由于LB+DB=DB-AN=N+1;%迭代步数tol=norm(x-x0);%前后两步迭代结果的误差x0=x;if(N>=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend14.BGS块高斯-赛德尔迭代法求线性方程组Ax=b的解function[x,N]=BGS(A,b,x0,d,eps,M)if nargin==4eps=1.0e-6;M=10000;elseif nargin<4errorreturnelseif nargin==5M=10000;endNS=size(A);n=NS(1,1);bnum=length(d);bs=ones(bnum,1);for i=1:(bnum-1)bs(i+1,1)=sum(d(1:i))+1;%获得对角线上每个分块矩阵元素索引的起始值endDB=zeros(n,n);for i=1:bnumendb=bs(i,1)+d(i,1)-1;DB(bs(i,1):endb,bs(i,1):endb)=A(bs(i,1):endb,bs(i,1):endb); %求A的对角分块矩阵endLB=-tril(A-DB);%求A的下三角分块阵UB=-triu(A-DB);%求A的上三角分块阵N=0;tol=1;while tol>=epsinvDL=inv(DB-LB);x=invDL*UB*x0+invDL*b;%块迭代公式N=N+1;tol=norm(x-x0);x0=x;if(N>=M)disp('Warning:迭代次数太多，可能不收敛！');return;end15.BSOR块逐次超松弛迭代法求线性方程组Ax=b的解function[x,N]=BSOR(A,b,x0,d,w,eps,M)if nargin==5eps=1.0e-6;M=10000;elseif nargin<5errorreturnelseif nargin==6M=10000;%参数默认值endNS=size(A);n=NS(1,1);bnum=length(d);bs=ones(bnum,1);for i=1:(bnum-1)bs(i+1,1)=sum(d(1:i))+1;%获得对角线上每个分块矩阵元素索引的起始值endDB=zeros(n,n);for i=1:bnumendb=bs(i,1)+d(i,1)-1;DB(bs(i,1):endb,bs(i,1):endb)=A(bs(i,1):endb,bs(i,1):endb); %求A的对角矩阵endLB=-tril(A-DB);%求A的下三角阵UB=-triu(A-DB);%求A的上三角阵N=0;iw=1-w;while tol>=epsinvDL=inv(DB-w*LB);x=invDL*(iw*DB+w*UB)*x0+w*invDL*b;%块迭代公式N=N+1;tol=norm(x-x0);x0=x;if(N>=M)disp('Warning:迭代次数太多，可能不收敛！'); return;endend。