• 线性方程组求解的迭代方法&Python实现

    目录

    1 线性方程组迭代法概述

    2 Jacobi迭代法和Gauss-Seidel迭代法

    3 程序实现流程

     4 案例及Python代码

    1 线性方程组迭代法概述

    设线性方程组为

    Ax=b

    其中A为非奇异，且b≠0，则线性方程组存在唯一非零解x*。令A=M-N，其中M为非奇异，则上述公式可以改写成等价方程组：

    x=Gx+f

    其中，G=M^{-1}N,f=M^{-1}b。则可以得到迭代公式：

    x_{n+1}=Gx_n+f

    如果序列{{x_n}}是收敛的，则有：

    \lim_{n\rightarrow \infty }x_n=x*.

    2 Jacobi迭代法和Gauss-Seidel迭代法

    设

    M=\begin{bmatrix} a_1_1 & 0&...& 0\\ 0&a_2_2 & ...& 0\\ ...& ...& ... &... \\ 0& 0& ...& a_n_n \end{bmatrix}" role="presentation" style="position: relative;">M=\begin{bmatrix} a_1_1 & 0&...& 0\\ 0&a_2_2 & ...& 0\\ ...& ...& ... &... \\ 0& 0& ...& a_n_n \end{bmatrix},      N=\begin{bmatrix} 0 &a_1_2 &...&a_1_n \\ a_2_1& 0 & ...&a_2_n \\ ... & ...&... & ...\\ a_n_1& a_n_2 & ... & 0 \end{bmatrix}" role="presentation" style="position: relative;">N=\begin{bmatrix} 0 &a_1_2 &...&a_1_n \\ a_2_1& 0 & ...&a_2_n \\ ... & ...&... & ...\\ a_n_1& a_n_2 & ... & 0 \end{bmatrix}

    则等价方程x=Gx+f可表示为：

    其中，G=M^{-1}N,f=M^{-1}b。 

    任取始解向量x^{(0)}=(x^{(0)}_1,x^{(0)}_2,...,x^{(0)}_n)^T，代入等价方程组可得迭代公式：

    x^{(k+1)}=Gx^{(k)}+f,(k=0,1,2,...,n)

    其分量形式可以表示为：

    上述的迭代公式称为 Jacobi迭代法，简称J法。

    如果将上式改为：

    则称为Gauss-Seidel迭代法，简称G-S法。

    G-S法和J法的不同点在于：每得到一个新分量x^{k+1}_i时，在计算以下各分量时，用x^{k+1}_i代替旧的x^k_i,因为新分量比旧分量更接近于真实解x_i*。

    3 程序实现流程

    Jacobi迭代法：

    Gauss-Seidel迭代法：

     

     4 案例及Python代码

    已知线性方程组：

     采用 Jacobi迭代法和Gauss-Seidel迭代法求解上述方程组的解，误差要求为10^-^9" role="presentation" style="position: relative;">10^-^9，输出线性方程组的解以及迭代次数。

    Python代码：

    ① Jacobi迭代法

    1. #-----Jacobi迭代法求解上述方程组的解
    2. import numpy as np
    3. A=np.array([[4,-1,2],[2,-5,1],[-2,1,4]])
    4. b=np.array([[7,-1,6]]).T
    5. e=10**-9 #误差
    6. N=10000 #最大迭代次数
    7. n=len(b)
    8. x0=np.zeros((n,1)) #迭代初值
    9. x1=np.zeros((n,1)) #输出矩阵初始化
    10. L_J=0 #初始化Jacobi迭代法的迭代次数
    11. #-----Jacobi迭代法-----#
    12. for i in range(N):
    13.     for j in range(n):
    14.       index=np.append(np.arange(0,j),np.arange(j+1,n)) #剔除掉j之后的线性序列
    15.       x1[j]=(b[j]-A[j,index]@x0[index])/A[j,j] #迭代公式
    16.     L_J=L_J+1 #累计迭代次数
    17.     if max(abs(A@x1-b))#利用残差判断
    18.         break
    19.     x0 = x1  # 更新初始向量
    20. print(f"x1={x1}") #输出线性方程组的解
    21. print(A@x1-b) #验证解的正确性
    22. print(f"L_J={L_J}") #输出迭代次数

    运行结果：

    1. x1=[[1.1]
    2.  [1. ]
    3.  [1.8]]
    4. [[6.65227873e-10]
    5.  [3.32616823e-10]
    6.  [0.00000000e+00]]
    7. L_J=16

     ② Gauss-Seidel迭代法

    1. #Gauss-Seidel迭代法求解上述方程组的解
    2. import numpy as np
    3. A=np.array([[4,-1,2],[2,-5,1],[-2,1,4]])
    4. b=np.array([[7,-1,6]]).T
    5. e=10**-9 #误差
    6. max_number=10000 #最大迭代次数
    7. n=len(b)
    8. x0=np.zeros((n,1)) #迭代初值
    9. x1=np.zeros((n,1)) #输出矩阵初始化
    10. L_G=0 #初始化迭代次数
    11. for k in range(max_number):
    12.     for j in range(n):
    13.       if j==0:
    14.           x1[0]=(b[0]-A[0,1:n]@x0[1:n])/A[0,0]
    15.       elif j==n:
    16.           x1[n]=(b[n-1]-A[n-1,0:n]@x1[0:n])/A[0,0]
    17.       else:
    18.           x1[j]=(b[j]-A[j,0:j]@x1[0:j])/A[j,j]-A[j,j+1:n]@x0[j+1:n]/A[j,j] #迭代公式
    19.     L_G=L_G+1 #更新迭代次数
    20.     if max(abs(A@x1-b))#利用残差判断
    21.         break
    22.     x0 = x1  # 更新迭代初值
    23. print(f"x1={x1}") #输出线性方程组的解
    24. print(A@x1-b) #验证解的正确性
    25. print(f"L_G={L_G}") #输出迭代次数

    运行结果如下：

    1. x1=[[1.1]
    2.  [1. ]
    3.  [1.8]]
    4. [[8.81430928e-10]
    5.  [4.40715464e-10]
    6.  [0.00000000e+00]]
    7. L_G=17

  • 相关阅读:
    如何通过沉浸式投影技术提升文旅夜游的互动体验？
    【LeetCode热题100】--283.移动零
    如何使用闲置的云服务器搭建一个属于自己的私人云网盘（可道云kodbox）
    Vue实战篇二十八：实现一个手机版的购物车
    【编程题】【Scratch一级】2022.09 踢足球
    龙口联合化学通过注册：年营收5.5亿 李秀梅控制92.5%股权
    【Python】取火柴小游戏（巴什博弈）
    数据湖存储在大模型中的应用
    ROS | ros::NodeHandle
    aspose-slides-22.5-jdk16
  • 原文地址：https://blog.csdn.net/dongke1991/article/details/127641843