Python实现高斯赛德尔迭代法（Gauss-Seidel）与雅可比迭代法（Jacobi）
数值分析：Python实现高斯赛德尔迭代法（Gauss-Seidel）与雅可比迭代法（Jacobi）
一、基本迭代法

二、雅可比迭代法（Jacobi）



三、高斯—赛德尔迭代法（Gauss-Seidel）
题目：分别使用雅可比迭代法与高斯—赛德尔迭代法求解Ax=b线性方程组，其中A为10阶希尔伯特矩阵输出解向量与迭代次数。
希尔伯特矩阵 Hilbert matrix
希尔伯特矩阵（Hilbert matrix），矩阵的一种，其元素A（i,j）=1/(i+j-1)，i,j分别为其行标和列标。


代码实现：
import numpy as np
n=10
def Jacobi(A,b,x0,xstar):
    k=0
    while True:
        for i in range(0, n):
            sum = 0.0
            for j in range(0, i):
                sum = sum + A[i][j] * x0[j]
            for j in range(i + 1, n):
                sum = sum + A[i][j] * x0[j]
            xstar[i] = (b[i] - sum) / A[i][i]
        temp = np.fabs(xstar[0] - x0[0])
        for j in range(1, n):
            if np.fabs(xstar[j] - x0[j]) > temp:
                temp = np.fabs(xstar[j] - x0[j])
        for j in range(0, n):
            x0[j] = xstar[j]
        k = k + 1
        if (temp < 1.0e-6 or k > 1000) :
            break
    print("Jacobi:",k)

def Gauss(A,b,x0,xstar):
    k = 0
    while True:
        for i in range(0, n):
            sum = 0.0
            for j in range(0, i):
                sum = sum + A[i][j] * xstar[j]
            for j in range(i + 1, n):
                sum = sum + A[i][j] * x0[j]
            xstar[i] = (b[i] - sum) / A[i][i]
        temp = np.fabs(xstar[0] - x0[0])
        for j in range(1, n):
            if np.fabs(xstar[j] - x0[j]) > temp:
                temp = np.fabs(xstar[j] - x0[j])
        for j in range(0, n):
            x0[j] = xstar[j]
        k = k + 1
        if (temp < 1.0e-6 or k > 1000) :
            break
    print("Gauss:",k)

def main():
    A=np.zeros([n,n],float)
    b=np.zeros(n,float)
    x0=np.zeros(n,float)
    Jex = np.zeros(n,float)
    GSex = np.zeros(n,float)
    for i in range(0,n):
        for j in range(0,n):
            A[i][j] = 1./(i+1+j)
    for i in range(0,n):
        A[i][i] = A[i][i] + 0.5
        b[i] = 1
        x0[i] = 0

    Jacobi(A,b,x0,Jex)
    for i in range(0,n):
        x0[i] = 0
    Gauss(A,b,x0,GSex)
    for i in range(0,n):
        print("Jex[",i,"]=",Jex[i])
    print("\n")
    for i in range(0,n):
        print("GSex[",i,"]=",GSex[i])

if __name__ == '__main__':
main()


通过编码计算，对比两种算法的收敛速度与结果。可以看出，高斯赛德

尔迭代法的收敛速度比雅可比迭代法收敛速度要快的多。

文章目录
前言
1 直接法
2 迭代法
小结

前言
今天我们要说的就是数值微积分，赶紧看看他和高等数学中的微积分有什么区别吧。本文是科学计算与MATLAB语言专题六第2小节的学习笔记，如果大家有时间的话，可以去听听课，没有的话，可以看看我的笔记，还是很不错的。
1 直接法
1.线性方程组的直接解法

高斯（Gauss）消去法

列主元消去法

矩阵的三角分解法

高斯（Gauss）消去法是一个经典的直接法，由它改进得到的列主元消去法，是目前计算机上求解线性方程组的标准算法，其特点就是通过消元将一般线性方程组的求解问题转化为三角方程组的求解问题。此外，还有矩阵的三角分解法等许多直接求解算法。

（1）利用左除运算符的直接解法

MATLAB提供了一个左除运算符“\”用于求解线性方程组，它使用列主元消去法，使用起来十分方便。对于线性方程组$Ax=b$，可以利用左除运算符反斜杠求解，b左除以A可获得线性方程组的数值解x。

$Ax=b$→$x=A\backslash b$

注意，如果矩阵A是奇异的或接近奇异的，则MATLAB会给出警告信息。

例1  用左除运算符求解下列线性方程组。

$\left\{
\begin{aligned}
2x_1+x_2-5x_3+x_4=13\\
x_1-5x_2+7x_4=-9\\
2x_2+x_3-x_4=6\\
x_1+6x_2-x_3-4x_4=0
\end{aligned}
\right.$
A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';
x=A\b

（2）利用矩阵分解求解线性方程组矩阵分解是设计算法的重要技巧，是指将一个给定的矩阵分解成若干个特殊类型矩阵的乘积，从而将一个一般的矩阵计算问题转化为几个易求的特殊矩阵的计算问题。通过矩阵分解方法求解线性方程组的优点是运算速度快，可以节省存储空间。

LU分解

QR分解

Cholesky分解

①LU分解的基本思想

矩阵的LU分解就是将一个n阶矩阵表示为一个下三角矩阵和一个上三角矩阵的乘积。线性代数中已经证明，只要方阵是非奇异的，LU分解总是可以进行的。

$A=LU →Ax=b→LUx=b$→ $\left\{
\begin{aligned}Ly=b\\Ux=y
\end{aligned}
\right.$

对于三角方程很容易求解，于是可以首解向量y使Ly=b，再求解Ux=y，从而达到求解线性方程组Ax=b的目的。

②MATLAB的LU分解函数

LU分解函数是根据列主元LU分解算法定义的，具有较好的数据稳定性。lu函数有两种调用格式：

[L，U]=lu（A）：

产生一个上三角阵U和一个变换形式的下三角阵L，使之满足A=LU。注意，这里的矩阵A必须是方阵。

[L，U，P]=lu（A）：

产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PA=LU。同样，矩阵A必须是方阵。

当使用第一种格式时，矩阵L往往不是一个下三角阵，但可以通过行交换成为一个下三角阵。

③用LU分解求解线性方程组

$Ax=b$→$LUx=b$→$x=U \backslash (L\backslash b）$

或

$Ax=b$→$PAx=Pb$→$LUx=Pb$→$x=U \backslash（L*b）$

通过LU分解后可以大大提高运算速度。

例2  用LU分解求解例1中的线性方程组。
A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';
[L,U]=lu(A);
x=U\(L\b)

x =
  -66.5556
   25.6667
  -18.7778
   26.5556

2 迭代法
2.线性方程组的迭代解法

迭代法是一种不断用变量的原值推出它的新值的过程，是用计算机解决问题的一种基本方法。

$\left\{
\begin{aligned}
10x_1-x_2=9\\
-x_1+10x_2-2x_3=7\\
-2x_2+10x_3=6
\end{aligned}
\right.$$\implies$ $\left\{
\begin{aligned}
x_2&=10x_1-9\\
x_1&=10x_2-2x_3-7\\
x_3&=(6+2x_2)10
\end{aligned}
\right.$$\implies$$\left\{
\begin{aligned}
x_1^{k+1}&=10x_2^{k}-2x_3^{k}-7\\
x_2^{k+1}&=10x_1^{k}-9\\
x_3^{k+1}&=0.6+0.2x_2^{k}
\end{aligned}
\right.$

（1）雅可比（Jacobi）迭代法

$Ax=b,
A=D-L-U，
(D-L-U)x=b$

$D=\left [
\begin{array}{cccc}
 a_{11} & & & \\
               & a_{22} &  & \\
                &  & \ddots & \\
                                          & &   & a_{nn} 
\end{array}
\right]$

$L=-\left [
\begin{array}{cccc}
  0 &  &&  &  \\
  a_{21}  & 0 &&  & \\
  a_{31}  &  a_{32}   &0&  &  \\
\vdots &\vdots&\vdots &\ddots&\\
  a_{n1} &  a_{n2}&\cdots& \cdots&0
\end{array}
\right]$

$U=-\left [
\begin{array}{cccc}
  0 & a_{12} &a_{13}& \cdots & a_{1n} \\
               & 0 &a_{23}& \cdots &a_{2n} \\
                &  &\ddots& \ddots & \vdots \\
                          &   & &\ddots&a_{n-1n}\\
                                          &&&   &0
\end{array}
\right]$
求解公式为：

$x=D^{-1}（L+U）x+D^{-1}b$

与之对应的迭代公式为：

$x^{k+1}=D^{-1}（L+U）x^{k}+D^{-1}b$ $\implies(B=D^{-1}(L+U),f=D^{-1}b$) $\implies x^{k+1}=Bx^{k}+f$

雅可比迭代法的函数文件jacobi.m：
function [y,n]=jacobi(A,b,x0,ep)
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=D\(L+U);
f=D\b;
y=B*x0+f;
n=1;
while norm(y-x0)>=ep
    x0=y;
    y=B*x0+f;
    n=n+1;
end

（2）高斯-赛德尔（Gauss-Serdel）迭代法

$Dx^{k+1}=（L+U）x^{k}+b
→Dx^{k+1}=Lx^{k+1}+Ux^k+b→
（D-L）x^{k+1}=Ux^{k}+b \\
x^{k+1}=（D-L）^{-1}U×0+（U-L）^{-1}b
→B=D-L^{-1}U,f=D-L^{-1}b→ x^{k+1}=Bx^{k}+f$

Gauss-Serdel迭代法的函数文件gauseidel.m
function [y,n]=gauseidel(A,b,x0,ep)
D=diag(diag(A));
L=-tril(A,-1); 
U=-triu(A,1);
B=(D-L)\U;
f=(D-L)\b;
y=B*x0+f;
n=1;
while norm(y-x0)>=ep
    x0=y;
    y=B*x0+f;
    n=n+1;
end



例3  分别用雅可比迭代法和高斯-赛德尔迭代法求解线性方程组。设迭代初值为0，迭代精度为10^(-6)。

$\left[\begin{matrix}
4&-2&-1\\
-2&4&3\\
-1&-3&3
\end{matrix}\right]$$\left[\begin{matrix}
x_1\\
x_2\\
x_3
\end{matrix}\right]$=$\left[\begin{matrix}
1\\
5\\
0
\end{matrix}\right]$
A=[4,-2,-1;-2,4,3;-1,-3,3];
b=[1,5,0]';
[x,n]=jacobi(A,b,[0,0,0]',1.0e-6)
[x,n]=gauseidel(A,b,[0,0,0]',1.0e-6)

x =
    0.9706
    0.8529
    1.1765
n =
    35
x =
    0.9706
    0.8529
    1.1765
n =
    16

例4  分别用雅可比迭代法和高斯-赛德尔迭代法求解下列线性方程组，看是否收敛。

$\left[\begin{matrix}
1&-2&2\\
1&1&1\\
2&2&1
\end{matrix}\right]$$\left[\begin{matrix}
x_1\\
x_2\\
x_3
\end{matrix}\right]$=$\left[\begin{matrix}
9\\
7\\
6
\end{matrix}\right]$
A=[1,2,-2;1,1,1;2,2,1];
b=[9;7;6];
[x,n]=jacobi(A,b,[0;0;0],1.0e-6)
[x,n]=gauseidel(A,b,[0;0;0],1.0e-6)

x =
   -27
    26
     8
n =
     4
x =
   NaN
   NaN
   NaN
n =
    1012

小结
直接法：以矩阵初等变换为基础，可以求得方程组的精确解；占用的内存空间大、程序实现较为复杂；一般适合求解低阶稠密线性方程组。

迭代法：从给定初始值逐步逼近精确解的过程，求解过程占用存储空间小，程序设计简单；适用于求解大型稀疏矩阵线性方程组；要考虑算法的收敛性。

大家学会了吗？另外分享给大家一个Markdown的公式神器-$\ Katex$ ，为了编辑出优美的公式，大家可以多看看官方支持文档.最后没有一键三连，但欢迎大家

点赞👍，收藏⭐，转发🚀，

如有问题、建议，请您在评论区留言💬哦。

雅可比迭代法
雅可比迭代法-高斯-赛德尔迭代法求线性方程组

雅可比迭代法-高斯-赛德尔迭代法求线性方程组

求解方程组



matlab程序实现

function [x,k] = Fjacobi(A,b,x0,nm,eps)
    %用雅可比迭代法求解方程组Ax=b
    %输入：A为方程组的系数矩阵，b为方程组右端的列向量，
    %输入：x0为迭代初值构成的列向量，nm为最大迭代次数,eps为误差精度
    %输出：x为求得的方程组的解构成的列向量，k为迭代次数
    D = diag(diag(A));
    L=-tril(A,-1);
    U=-triu(A,1);
    B=D\(L+U);
    d=D\b;
    k=1;
    disp('雅可比迭代法');
    while k<=nm
       x=B*x0+d; %雅可比迭代式
       if norm(x-x0)<eps
           disp('迭代次数为');k
           disp('方程组的解为');x
       return;
       end
       x0=x; 
       k=k+1;
    end
    %下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
    disp('在最大迭代次数内不收敛!');
    disp('最大迭代次数后的结果为');x

function [x,k] = Fgseid(A,b,x0,nm,eps)
    %用高斯-塞德尔迭代法法求解方程组Ax=b
    %输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数,eps为误差精度
    %输出: x为求得的方程组的解构成的列向量，k为迭代次数
    D = diag(diag(A));
    L = -tril(A,-1);
    U = -triu(A,1);
    G = (D-L)\U;
    d = (D-L)\b;
    
    k = 1;
    disp('高斯-塞德尔迭代法');
    while k<=nm
        x  = G * x0 + d;%计算迭代矩阵
        if norm(x-x0)<eps
            disp('迭代次数为');k
            disp('方程组的解为');x
        return;
        %上面：达到精度要求就结束程序，输出迭代次数和方程组的解
        end
        x0 = x;
        k  = k + 1;
    end
    %下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
    disp('在最大迭代次数内不收敛!');
    disp('最大迭代次数后的结果为');x



求解函数；

a = [5 -1 -1 -1; -1 10 -1 -1; -1 -1 5 -1; -1 -1 -1 10];
b = [ 4 13 8 34]';
x0 = [0 0 0 0]';
nm = 100;
eps = 0.0001;
[x,k] = Fjacobi(a,b,x0,nm,eps);
[x,k] = Fgseid(a,b,x0,nm,eps);

结果：

一般而言，因雅可比迭代收敛速度不够快，所以在工程中应用不多。并且在雅可比迭代收敛速
度很慢的情况下，通常高斯-赛德尔方法也不会很快。我们可以对高斯-赛德尔方法做出一定的
修改，来提高收敛速度。本节，我们介绍超松弛迭代法。
下面我们对该算法进行收敛性分析：
这里我们解释一下什么叫严格对角占优矩阵：严格对角占优矩阵就是矩阵主对角线元素的的模
大于与他同行的其他元素的模的总和。
下面我们具体代码实现一下：
%超松弛（SOR）迭代法，计算线性方程组的解
function [x,k] = SOR_iteration(A,b,x0,w,tol)
% tol为输入误差容限，x0为迭代初始值

% 默认最多迭代300次，超出300次会给出警示
max = 300;

if(w<=0||w>=2)
    % MATLAB中error语句用于报错跳出，并可以给出相应提示
    error('错啦！w的值不符合要求');
    
    % 在MATLAB中遇到return的意思就是结束整个函数的执行，
    % break是仅仅跳出循环体
    return;
end

%  取出X矩阵的对角元，然后构建一个以X对角元为对角的对角矩阵
D = diag(diag(A));

% 求A的下三角矩阵,对角线元素为0，再每个矩阵元素取负号
L = -tril(A,-1);

% 求A的上三角矩阵，对角线元素为0，再每个矩阵元素取负号
U = -triu(A,1);

% 在MATLAB中inv是求矩阵的逆矩阵的意思，同具有一样的功能
B = inv(D-L*w)*((1-w)*D+w*U);
f = w*inv((D-L*w))*b;
x = B*x0+f;
k = 1;

while norm(x-x0)>=tol
    x0 = x;
    x = B*x0+f;
    k = k+1;
    if(k>=max)
        disp('迭代次数超过',max1,'次，方程组可能不收敛');
        return;
    end
    
    [k,x']
end
我们用超松弛迭代法求解一下下面的方程组：
在编辑界面输入如下命令：
A = [4,-2,-1;-2,4,-2;-1,-2,3];

% 注意必须加'表示转置
b = [0,-2,3]';
x0 = [1,1,1]';
[x,k] = SOR_iteration(A,b,x0,1.45,1e-6)
最后计算结果如下：