数值分析第六章节 用Python实现解线性方程组的迭代法

参考书籍：数值分析 第五版 李庆杨 王能超 易大义编 第5章 解线性方程组的迭代法
文章声明：如有发现错误，欢迎批评指正

迭代法的基本概念

6.1.1引言：定义：(1)对于给定的线性方程组$x=Bx+f$，用公式$x^{(k+1)}=B{x}^{(k)}+f$逐步带入求近似解的方法称为迭代法(或称为一阶定常迭代法，这里$B$与$k$无关)(2)如果$\underset{k\to \infty }{\lim }{x}^{(k)}$存在(记为$x^{\ast }$)，称此迭代法收敛，显然$x^{\ast }$就是此方程组的解，否则称此迭代法发散。6.1.2：向量序列与矩阵序列的极限：给定线性方程组$x=Bx+f$及一阶定常迭代法$x^{(k+1)}=B{x}^{(k)}+f$式，对任意选取初始向量$x^{(0)}$，迭代法$x^{(k+1)}=B{x}^{(k)}+f$式收敛的充要条件是矩阵$B$的谱半径$\rho (B)<1$。其他跳过。

雅可比迭代法与高斯-塞格尔迭代法

雅可比迭代法

$$\{\begin{array}{c}{x}^{(0)}\\ x^{(k+1)}=B{x}^{(k)}+f,k=0,1,\dots ,\end{array}{x}^{(0)}为初始向量，B=D^{-1}(L+U),f=D^{-1}b$$
我感觉我写得挺好，可以算作通用代码，前提必须保证收敛。输入：输入系数矩阵行数，系数矩阵，初始向量，迭代次数。输出：解的向量。命名十分规范，懂了理论不难看懂。

def func1(B,x):#不通用的矩阵乘法global nlt=[]for i in range(n):cnt=0for j in range(n):cnt+=B[i][j]*x[j]lt.append(cnt)return lt
def func2(Bx,f):#不通用的矩阵加法global nlt=[]for i in range(n):lt.append(Bx[i]+f[i])return lt
n=int(input())
lt=[]
for _ in range(n):lt.append([eval(_) for _ in input().strip().split()])
D_inv=[[0 for _ in range(n)] for _ in range(n)]
for i in range(n):D_inv[i][i]=1/lt[i][i]
L_sum_U=[[0 for _ in range(n)] for _ in range(n)]
for i in range(1,n):for j in range(i):L_sum_U[i][j]=-lt[i][j]
for i in range(n-1):for j in range(i+1,n):L_sum_U[i][j]=-lt[i][j]
B=[[0 for _ in range(n)] for _ in range(n)]
for i in range(n):for j in range(n):B[i][j]=L_sum_U[i][j]*D_inv[i][i]
f=[0 for _ in range(n)]
for i in range(n):f[i]=D_inv[i][i]*lt[i][-1]
x=[eval(_) for _ in input().strip().split()]
num=int(input())
for _ in range(1,num+1):x=func2(func1(B,x),f)
print(x)

用的例1，一模一样。

高斯-塞格尔迭代法

$$\{\begin{array}{c}{x}^{(0)}\\ x^{(k+1)}=B{x}^{(k)}+f,k=0,1,\dots ,\end{array}{x}^{(0)}为初始向量，B=(D-L{)}^{-1}U,f=(D-L{)}^{-1}b$$
我感觉我写得挺好，可以算作通用代码，前提必须保证收敛。输入：输入系数矩阵行数，系数矩阵，初始向量，迭代次数。输出：解的向量。命名十分规范，懂了理论不难看懂。

def func1(lt1,lt2):#矩阵乘法a,b=len(lt1),len(lt2[0])lt=[[0 for _ in range(b)] for _ in range(a)]for i in range(a):for j in range(b):for p in range(len(lt1[0])):lt[i][j]+=lt1[i][p]*lt2[p][j]return lt
def func2(lt1,lt2):#不通用的矩阵加法global nlt=[]for i in range(n):lt.append([lt1[i][0]+lt2[i][0]])return lt
n=int(input())
lt=[]
for _ in range(n):lt.append([eval(_) for _ in input().strip().split()])
D=[[0 for _ in range(n)] for _ in range(n)]
for i in range(n):D[i][i]=lt[i][i]
L=[[0 for _ in range(n)] for _ in range(n)]
for i in range(1,n):for j in range(i):L[i][j]=-lt[i][j]
U=[[0 for _ in range(n)] for _ in range(n)]
for i in range(n-1):for j in range(i+1,n):U[i][j]=-lt[i][j]
D_minus_L=[[0 for _ in range(n)] for _ in range(n)]
for i in range(n):for j in range(n):D_minus_L[i][j]=D[i][j]-L[i][j]
#这里涉及一个求解下三角阵的逆矩阵
D_minus_L_inv=[[0 for _ in range(n)] for _ in range(n)]
for i in range(n):for j in range(i):cnt=0for k in range(i):cnt-=D_minus_L[i][k]*D_minus_L_inv[k][j]D_minus_L_inv[i][j]=cnt/D_minus_L[i][i]D_minus_L_inv[i][i]=1/D_minus_L[i][i]
B=func1(D_minus_L_inv,U)
f=func1(D_minus_L_inv,[[lt[_][-1]] for _ in range(n)])
x=[[eval(_)] for _ in input().strip().split()]
num=int(input())
for _ in range(1,num+1):x=func2(func1(B,x),f)
print(x)

用的例1，一模一样。

就这样吧，剩下方法，自己研究。