python的史蒂芬加速迭代法_Jacobi迭代法的python实现

设有线性方程组Ax=b，其中n阶矩阵A=(aij)使非奇异矩阵

将A分解为A=M-N，因而Mx=Nx+b，求解Ax=b转化为求解x=M的逆Nx+M的逆b

可构造一阶梯定常迭代法x(k+1)=Bx(k)+f

令A=D-L-U，其中D是有A的主对角元组成的对角矩阵，L是A的下三角矩阵每一行乘以-1，U是A的上三角矩阵每一行乘以-1

Jacobi迭代法中，令M=D，A=D-N，因而B=M的逆N=M的逆(M-A)=D的逆*(L+U)

f=D的逆b

确定初值x0，开始迭代

若x最终收敛于x，则线性方程组有解x*，否则无解

确定需要的绝对误差限或相对误差限e，确定相应的迭代次数

下面给出一个例子

矩阵A

10.0 3.0 1.0

2.0 -10.0 3.0

1.0 3.0 10.0

b=(14.0,-5.0,14.0)’

要求绝对误差限e=0.00001即精确到小数点后5位

代码如下

# -*- coding: utf-8 -*-

"""

Created on Wed Feb 13 08:09:46 2019

@author: 鹰皇

"""

#jacobi迭代法

import numpy as np

A=[[10.0,3.0,1.0],[2.0,-10.0,3.0],[1.0,3.0,10.0]]

b=[[14.0,-5.0,14.0]]

def get_base(A):#获得一个基，在基上修改

base=list(np.zeros((len(A),len(A))))

D=[]

for i in base:

D.append(list(i))

return D

def get_D(A):#获得对角矩阵D

D=get_base(A)

for i in A:

D[A.index(i)][A.index(i)]=A[A.index(i)][A.index(i)]

return D

def get_LU(A,D):#获得L+U

AA=np.array(A)

DD=np.array(D)

LLUU=DD-AA

LU=[]

for i in LLUU:

LU.append(list(i))

return LU

def matrix_to_list(x):#将矩阵转换为列表的函数，numpy中应当有这一函数

d=[]

ans=[]

for i in x:

d.append(i.tolist())

for i in d:

ans.append(i[0])

return ans

def roll(A,b,x0):#主循环结构

D=get_D(A)

LU=np.mat(get_LU(A,D))

D=np.mat(get_D(A))

B=D.I*LU

b1=np.mat(b).T

f=D.I*b1

x=np.mat(x0).T

y=B*x+f

return matrix_to_list(y.T)

x1=[[1,1,2]]

x0= roll(A,b,x1)

def main(A,b,x0,e):#输入系数矩阵，方程组右端b，初始值x0，以及需要的绝对误差限

n=0

ans=[]

ans.append(x1)

ans.append(x0)

ans1=[]

ans1.append(x1[0])

ans1.append(x0[0])

while abs(ans[-1][0][1]-ans[-2][0][1])>e:

n=n+1

x0=roll(A,b,x0)

ans.append(x0)

for i in ans:

ans1.append(i[0])

return ans1,n#输出迭代结果和迭代次数

e=0.00001

print main(A,b,x0,e)