雅克比迭代和高斯赛德尔迭代法可以用来求线性方程组的近似解，与雅克比迭代不同的是，高斯赛德尔迭代法在第k+1次迭代求 $x_i^{k+1}$时，会用到 $x_1^k,x_2^k..x_{i-1}^k$ 的值。
 假设有如下的线性方程组：
 $$\begin{array}{rl}{a}_{11}\cdot x_1+a_{12}\cdot x_2+\dots +a_{1n}\cdot x_n& =b_1\\ a_{21}\cdot x_1+a_{22}\cdot x_2+\dots +a_{2n}\cdot x_n& =b_2\\ ⋮& =⋮\\ a_{n1}\cdot x_1+a_{n2}\cdot x_2+\dots +a_{nn}\cdot x_n& =b_n.\end{array}$$
 雅克比迭代求解公式如下：
 $$x_m^{(k+1)}=\frac{1}{a_{mm}}\left(b_m-\sum _{j\ne m}{a}_{mj}{x}_j^{(k)}\right),1\le m\le n$$
 高斯赛德尔迭代求解公式如下：
 $$x_m^{k+1}=\frac{1}{a_{mm}}\left(b_m-\sum _{j=1}^{m-1}{a}_{mj}\cdot x_j^{k+1}-\sum _{j=m+1}^n{a}_{mj}\cdot x_j^k\right),1\le m\le n$$
 
该表达式的 $x_m$由线性方程组的第m行转换而来，同时考虑到了在计算 $x_m^{k+1}$时， $x_1^{k+1},x_2^{k+1}..x_{m-1}^{k+1}$已经迭代产生的情况，利用最新的参数进行迭代更新。
 
可以看到在两种迭代方法的求解公式中，$x_m$ 的计算需要除以其系数 $a_{mm}$ ,所以高斯赛德尔迭代法要求线性方程组的主对角矩阵元素非0。

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()

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

 
 
刚学 Jacobi算法和Gauss_Siedel算法不久，觉的对以后学习会有帮助，所以记下来，希望感兴趣的朋友共勉！
 
 

  
  
  雅克比迭代 
  
 
   

   

    #include
    <
    iostream
    >
    
#include
    "
    math.h
    "
    

    using
     
    namespace
     std;

    #define
     n 3 
    

    double
     a[n][n]
    =
    {{
    5
    ,
    2
    ,
    1
    },{
    -
    1
    ,
    4
    ,
    2
    },{
    2
    ,
    -
    3
    ,
    10
    }}, x[n]
    =
    {
    -
    3
    ,
    2
    ,
    1
    }, b[n]
    =
    {
    -
    12
    ,
    20
    ,
    3
    }, c[n];

    int
     i
    =
    0
    ,j,k
    =
    0
    ,counter
    =
    15
    ; 

    int
     main()
{


    int
     k
    =
    0
    ,i
    =
    0
    ; 

    double
     num1,E; 

    float
     y[n];

    do
    
 {
 E
    =
    0
    ; 

    for
    (j
    =
    0
    ;j
    <
    n;j
    ++
    )
 {
 c[j]
    =
    x[j];
 } 
 num1
    =
    0
    ;

    for
    (j
    =
    0
    ;j
    <
    n;j
    ++
    )
 {

    if
    (j
    !=
    i)
 num1
    =
    num1
    +
    a[i][j]
    *
    x[j]; 
 }

    for
    (j
    =
    0
    ;j
    <
    n;j
    ++
    )
 {
 c[j]
    =
    x[j];
 }

 y[i]
    =
    (b[i]
    -
    num1)
    /
    a[i][i];


    if
    (i
    <
    n)
 { i
    ++
    ; }


    else
    
 {

    for
    (j
    =
    0
    ;j
    <
    n;j
    ++
    )
 { x[j]
    =
    y[j]; }


    for
    (j
    =
    0
    ;j
    <
    n;j
    ++
    )
 {

    if
    (E
    <
    (
    float
    )fabs(c[j]
    -
    y[j])) E
    =
    (
    float
    )fabs(c[j]
    -
    y[j]);
 }

    if
    (k
    >=
    7
    )
 { cout
    <<
    "
    Iterations is 
    "
    <<
    k
    +
    1
    <<
    "
    ,
    "
    <<
    "
    Iterative precision is 
    "
     ;
 printf(
    "
    %.6lf\n
    "
    ,E);

    for
    (j
    =
    0
    ;j
    <
    n;j
    ++
    )
 printf(
    "
    x[%d]=%.5lf\t
    "
    ,j,x[j]); 
 cout
    <<
    endl
    <<
    endl;
 }
 i
    =
    0
    ;
 k
    ++
    ;

    if
    (E
    <
    0.001
    ) 
    break
    ;
 } 
 }
    while
    (k
    <=
    counter);
 cout
    <<
    "
    Iterations is 
    "
    <<
    k
    <<
    endl;
 cout
    <<
    "
    Iterative terminating condition 
    "
    ;
 printf(
    "
    x[%d]-x[%d]oo<%.6lf\n
    "
    ,k,k
    -
    1
    ,E);
 cout
    <<
    "
    \t\t---------Solution-----------
    "
    <<
    endl;

    for
    (i
    =
    0
    ;i
    <
    n;i
    ++
    )
 {
 printf(
    "
    x[%d]=%.5lf\n
    "
    ,i,x[i]);
 }
}
   
 
  
 
 
 
 
 
 
 
 
 
 

  
  
  高斯-赛德尔迭代法 
  
 
   

   

    #include 
    <
    iostream
    >
    
#include 
    <
    math.h
    >
    

    using
     
    namespace
     std;

    int
     main()

{

    double
     a[
    3
    ][
    3
    ]
    =
    {{
    1
    ,
    2
    ,
    -
    2
    },{
    3
    ,
    1
    ,
    -
    1
    },{
    2
    ,
    4
    ,
    1
    }},b[
    3
    ]
    =
    {
    10
    ,
    12
    ,
    -
    20
    };


    double
     x[
    3
    ]
    =
    {
    0
    ,
    0
    ,
    0
    },sum1,sum2;


    int
     i,j,k,n
    =
    3
    ;


    for
     (k
    =
    0
    ;k
    <
    5
    ;k
    ++
    )

{ 
    for
    (i
    =
    0
    ;i
    <
    n;i
    ++
    )

 { sum1
    =
    0
    ;sum2
    =
    0
    ;


    for
    (j
    =
    0
    ;j
    <
    i
    -
    1
    ;j
    ++
    )

 { sum1
    =
    sum1
    +
    a[i][j]
    *
    x[j];

 }


    for
    (j
    =
    i
    +
    1
    ;j
    <
    n;j
    ++
    )

 {sum2
    =
    sum2
    +
    a[i][j]
    *
    x[j];}

 x[i]
    =
    (b[i]
    -
    sum1
    -
    sum2)
    /
    a[i][i];

 }


    for
    (i
    =
    0
    ;i
    <
    n;i
    ++
    )

 { printf(
    "
    x%d=%-15f
    "
    ,i
    +
    1
    ,x[i]);}

 printf(
    "
    \n
    "
    ); 

} 

}

   
 
  
 
 
 
 
 
 
 
转载于:https://www.cnblogs.com/FCWORLD/archive/2010/12/04/1896063.html

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

 
前言
 
今天我们要说的就是数值微积分，赶紧看看他和高等数学中的微积分有什么区别吧。本文是科学计算与MATLAB语言专题六第2小节的学习笔记，如果大家有时间的话，可以去听听课，没有的话，可以看看我的笔记，还是很不错的。
 
1 直接法
 
1.线性方程组的直接解法
 高斯（Gauss）消去法
 列主元消去法
 矩阵的三角分解法
 高斯（Gauss）消去法是一个经典的直接法，由它改进得到的列主元消去法，是目前计算机上求解线性方程组的标准算法，其特点就是通过消元将一般线性方程组的求解问题转化为三角方程组的求解问题。此外，还有矩阵的三角分解法等许多直接求解算法。
 （1）利用左除运算符的直接解法
 MATLAB提供了一个左除运算符“\”用于求解线性方程组，它使用列主元消去法，使用起来十分方便。对于线性方程组$Ax=b$，可以利用左除运算符反斜杠求解，b左除以A可获得线性方程组的数值解x。
 $Ax=b$→$x=A\backslash b$
 注意，如果矩阵A是奇异的或接近奇异的，则MATLAB会给出警告信息。
 例1 用左除运算符求解下列线性方程组。
 $\{\begin{array}{r}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{array}$
 
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\to Ax=b\to LUx=b$→ $\{\begin{array}{r}Ly=b\\ Ux=y\end{array}$
 对于三角方程很容易求解，于是可以首解向量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\ast 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.线性方程组的迭代解法
 迭代法是一种不断用变量的原值推出它的新值的过程，是用计算机解决问题的一种基本方法。
 $\{\begin{array}{r}10x_1-x_2=9\\ -x_1+10x_2-2x_3=7\\ -2x_2+10x_3=6\end{array}$$⟹$ $\{\begin{array}{rl}{x}_2& =10x_1-9\\ x_1& =10x_2-2x_3-7\\ x_3& =(6+2x_2)10\end{array}$$⟹$$\{\begin{array}{rl}{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{array}$
 （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}{ccccc}0& & & & \\ a_{21}& 0& & & \\ a_{31}& a_{32}& 0& & \\ ⋮& ⋮& ⋮& \ddots & \\ a_{n1}& a_{n2}& \cdots & \cdots & 0\end{array}\right]$
 $U=-\left[\begin{array}{ccccc}0& a_{12}& a_{13}& \cdots & a_{1n}\\ & 0& a_{23}& \cdots & a_{2n}\\ & & \ddots & \ddots & ⋮\\ & & & \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$ $⟹(B=D^{-1}(L+U),f=D^{-1}b$) $⟹x^{k+1}=B{x}^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）迭代法
 $$\begin{gathered} D{x}^{k+1}=（L+U）x^k+b\to D{x}^{k+1}=L{x}^{k+1}+U{x}^k+b\to （D-L）x^{k+1}=U{x}^k+b \\ {x}^{k+1}=（D-L{）}^{-1}U\times 0+（U-L{）}^{-1}b\to B=D-L^{-1}U,f=D-L^{-1}b\to x^{k+1}=B{x}^k+f \end{gathered}$$
 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{array}{ccc}4& -2& -1\\ -2& 4& 3\\ -1& -3& 3\end{array}\right]$$\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]$=$\left[\begin{array}{c}1\\ 5\\ 0\end{array}\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{array}{ccc}1& -2& 2\\ 1& 1& 1\\ 2& 2& 1\end{array}\right]$$\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]$=$\left[\begin{array}{c}9\\ 7\\ 6\end{array}\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$ ，为了编辑出优美的公式，大家可以多看看官方支持文档.最后没有一键三连，但欢迎大家
 点赞👍，收藏⭐，转发🚀，
 如有问题、建议，请您在评论区留言💬哦。