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

  • Post author:
  • Post category:python


参考书籍:数值分析 第五版 李庆杨 王能超 易大义编 第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,一模一样。

在这里插入图片描述

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



版权声明:本文为zhishichengjin原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。