matlab gauss seidel迭代_解线性方程组的经典迭代法(2)-实战

上一节我们介绍了经典迭代法的理论，这节用MATLAB实现它.由于最近在学习微分方程数值解，会用迭代法应用到解方程组中.

目录

1. Jacobi迭代法代码
2. 数值实例与比较
3. 附录（其他迭代法代码）

1.Jacobi迭代法

代码实现基于以下观察：四种迭代法的迭代矩阵和迭代向量都基于最开始的矩阵分解

，然后直接套入公式循环求解即可.

由于四种迭代法的代码差别很小，故正文只写了一个，其余放在附录里.

使用的时候可以根据实际情况调节精度，迭代次数上限等参数.

function

2.数值实例与比较

用有限差分法逼近偏微分方程得到了一个3969阶的线性方程组，并用各个迭代法求解了一下，比较一下各迭代法使用的次数，并和系统自带的线性方程组求解语句比较了一下求解时间.

各算法所用次数和时间

由上图可见，Jacobi用时很长，由于题中的矩阵A最大和最小特征值模相同且操作更复杂些，JOR方法的时间最长，没有起到很好的效果.G-S方法相比Jacobi有很大的提升，SOR相比G-S方法也有很大提升，ADI和共轭梯度法不属于本章内容，但效果也都不错。

MATLAB表示：花里胡哨，我自带的都能秒杀你们，只需1s.

评论区好像有人问到，Matlab内部用的是啥算法，这我也不知道啊，不过Matlab之所以快肯定有一个原因就是它使用了底层的加速，而我们直接写的话肯定没法超过它的。

3.附录

3.1 Gauss-Seidel迭代代码

function [t x1]=Gauss_Seidel(A,b,x0,ep,kmax)
%功能：使用Gauss-Seidel迭代解线性方程组
%输入：参数Ax=b中的A，b,x0为初始值，ep为精度，kmax为最大迭代次数
%输出：参数t x1为迭代次数和解（列向量）
n=length(b);
%默认参数
if nargin<5,kmax=10000;end
if nargin<4,ep=1e-5;end
if nargin<3,x0=zeros(n,1);end
%处理得到U，L，D
U=zeros(n);L=zeros(n);D=zeros(n);
for i=1:n
    for j=1:n
        if i==j D(i,i)=A(i,i);end
        if i>j L(i,j)=A(i,j);end
        if i<j U(i,j)=A(i,j);end
    end
end       
%迭代
x=zeros(n,1);k=1;
G=-((D+L)U);F=(D+L)b;
while k<=kmax
    x=G*x0+F;
    if norm(x-x0,inf)<ep x1=x;break;end
    x0=x;k=k+1;
end
t=k;
x1=x;
if k==kmax,warning('已达迭代次数上限 ');end

3.2 JOR迭代代码

function [t x1]=JOR(A,b,x0,ep,kmax)
%功能：使用JOR迭代解线性方程组
%输入：参数Ax=b中的A，b,x0为初始值，ep为精度，kmax为最大迭代次数
%输出：参数t x1为迭代次数和解（列向量）
n=length(b);
%默认参数
if nargin<5,kmax=10000;end
if nargin<4,ep=1e-5;end
if nargin<3,x0=zeros(n,1);end
%处理得到U，L，D
U=zeros(n);L=zeros(n);D=zeros(n);
for i=1:n
    for j=1:n
        if i==j D(i,i)=A(i,i);end
        if i>j L(i,j)=A(i,j);end
        if i<j U(i,j)=A(i,j);end
    end
end
%求最佳松弛因子w
lambda=eig(-D(L+U));
w=2/(2-max(lambda)-min(lambda));
%迭代
x=zeros(n,1);k=1;
G=(eye(n)-w*DA);F=w*(Db);
while k<=kmax
    x=G*x0+F;
    if norm(x-x0,inf)<ep x1=x;break;end
    x0=x;k=k+1;
end
t=k;
x1=x;
if k==kmax,warning('已达迭代次数上限 ');end
end

3.3 SOR迭代代码

function [t x1]=SOR(A,b,x0,ep,kmax)
%功能：使用SOR迭代解线性方程组
%输入：参数Ax=b中的A，b,x0为初始值，ep为精度，kmax为最大迭代次数
%输出：参数t x1为迭代次数和解（列向量）
n=length(b);
%默认参数
if nargin<5,kmax=10000;end
if nargin<4,ep=1e-5;end
if nargin<3,x0=zeros(n,1);end
%处理得到U，L，D
U=zeros(n);L=zeros(n);D=zeros(n);
for i=1:n
    for j=1:n
        if i==j D(i,i)=A(i,i);end
        if i>j L(i,j)=A(i,j);end
        if i<j U(i,j)=A(i,j);end
    end
end
%求最佳松弛因子w
lambda=eig(-D(L+U));
w=max(abs(lambda));
w=2/(1+sqrt(1-w^2));
%迭代
G=(D+w*L)((1-w)*D-w*U);F=(D+w*L)(w*b);
x=zeros(n,1);k=1;
while k<=kmax
    x=G*x0+F;
    if norm(x-x0,inf)<ep x1=x;break;end
    x0=x;k=k+1;
end
t=k;
x1=x;
if k==kmax,warning('已达迭代次数上限 ');end

完结.