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

1 线性方程组迭代法概述

设线性方程组为

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

其中，

如果序列{

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

设

则等价方程可表示为：

其中，,。 

任取始解向量，代入等价方程组可得迭代公式：

其分量形式可以表示为：

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

如果将上式改为：

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

G-S法和J法的不同点在于：每得到一个新分量时，在计算以下各分量时，用代替旧的,因为新分量比旧分量更接近于真实解。

3 程序实现流程

Jacobi迭代法：

Gauss-Seidel迭代法：

 4 案例及Python代码

已知线性方程组：

 采用 Jacobi迭代法和Gauss-Seidel迭代法求解上述方程组的解，误差要求为，输出线性方程组的解以及迭代次数。

Python代码：

① Jacobi迭代法

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

运行结果：

x1=[[1.1]
 [1. ]
 [1.8]]
[[6.65227873e-10]
 [3.32616823e-10]
 [0.00000000e+00]]
L_J=16

 ② Gauss-Seidel迭代法

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

运行结果如下：

x1=[[1.1]
 [1. ]
 [1.8]]
[[8.81430928e-10]
 [4.40715464e-10]
 [0.00000000e+00]]
L_G=17