参考书籍:数值分析 第五版 李庆杨 王能超 易大义编 第5章 解线性方程组的迭代法
文章声明:如有发现错误,欢迎批评指正
迭代法的基本概念
6.1.1引言:定义:(1)对于给定的线性方程组
x
=
B
x
+
f
x=Bx+f
x
=
B
x
+
f
,用公式
x
(
k
+
1
)
=
B
x
(
k
)
+
f
x^{(k+1)}=Bx^{(k)}+f
x
(
k
+
1
)
=
B
x
(
k
)
+
f
逐步带入求近似解的方法称为迭代法(或称为一阶定常迭代法,这里
B
B
B
与
k
k
k
无关)(2)如果
lim
k
→
∞
x
(
k
)
\lim\limits_{k\rightarrow\infty}x^{(k)}
k
→
∞
lim
x
(
k
)
存在(记为
x
∗
x^*
x
∗
),称此迭代法收敛,显然
x
∗
x^{*}
x
∗
就是此方程组的解,否则称此迭代法发散。6.1.2:向量序列与矩阵序列的极限:给定线性方程组
x
=
B
x
+
f
x=Bx+f
x
=
B
x
+
f
及一阶定常迭代法
x
(
k
+
1
)
=
B
x
(
k
)
+
f
x^{(k+1)}=Bx^{(k)}+f
x
(
k
+
1
)
=
B
x
(
k
)
+
f
式,对任意选取初始向量
x
(
0
)
x^{(0)}
x
(
0
)
,迭代法
x
(
k
+
1
)
=
B
x
(
k
)
+
f
x^{(k+1)}=Bx^{(k)}+f
x
(
k
+
1
)
=
B
x
(
k
)
+
f
式收敛的充要条件是矩阵
B
B
B
的谱半径
ρ
(
B
)
<
1
\rho(B)<1
ρ
(
B
)
<
1
。其他跳过。
雅可比迭代法与高斯-塞格尔迭代法
雅可比迭代法
{
x
(
0
)
x
(
k
+
1
)
=
B
x
(
k
)
+
f
,
k
=
0
,
1
,
…
,
x
(
0
)
为初始向量,
B
=
D
−
1
(
L
+
U
)
,
f
=
D
−
1
b
\left\{\begin{matrix}x^{(0)}\\x^{(k+1)}=Bx^{(k)}+f,k=0,1,\dots,\end{matrix}\right.x^{(0)}为初始向量,B=D^{-1}(L+U),f=D^{-1}b
{
x
(
0
)
x
(
k
+
1
)
=
B
x
(
k
)
+
f
,
k
=
0
,
1
,
…
,
x
(
0
)
为初始向量,
B
=
D
−
1
(
L
+
U
)
,
f
=
D
−
1
b
我感觉我写得挺好,可以算作通用代码,前提必须保证收敛。输入:输入系数矩阵行数,系数矩阵,初始向量,迭代次数。输出:解的向量。命名十分规范,懂了理论不难看懂。
def func1(B,x):
#不通用的矩阵乘法
global n
lt=[]
for i in range(n):
cnt=0
for j in range(n):
cnt+=B[i][j]*x[j]
lt.append(cnt)
return lt
def func2(Bx,f):
#不通用的矩阵加法
global n
lt=[]
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,一模一样。
高斯-塞格尔迭代法
{
x
(
0
)
x
(
k
+
1
)
=
B
x
(
k
)
+
f
,
k
=
0
,
1
,
…
,
x
(
0
)
为初始向量,
B
=
(
D
−
L
)
−
1
U
,
f
=
(
D
−
L
)
−
1
b
\left\{\begin{matrix}x^{(0)}\\x^{(k+1)}=Bx^{(k)}+f,k=0,1,\dots,\end{matrix}\right.x^{(0)}为初始向量,B=(D-L)^{-1}U,f=(D-L)^{-1}b
{
x
(
0
)
x
(
k
+
1
)
=
B
x
(
k
)
+
f
,
k
=
0
,
1
,
…
,
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 n
lt=[]
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=0
for 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,一模一样。
就这样吧,剩下方法,自己研究。