Jacobi迭代算法的Python实现详解_Python

Jacobi迭代算法的Python实现详解

2021-07-27 00:11Fengqiao_x Python

这篇文章主要介绍了Jacobi迭代算法的Python实现详解,文中通过示例代码介绍的非常详细，对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下

				?

									import numpy as np

									import time

1.1 Jacobi迭代算法

				?

									def Jacobi_tensor_V2(A,b,Delta,m,n,M):

									start=time.perf_counter()#开始计时

									find=0#用于标记是否在规定步数内收敛

									X=np.ones(n)#迭代起始点

									x=np.ones(n)#用于存储迭代的中间结果

									d=np.ones(n)#用于存储Ax**(m-2)的对角线部分

									m1=m-1

									m2=2-m

									for i in range(M):

									print('X',X)

									a=np.copy(A)

									#得Ax**(m-2)

									for j in range(m-2):

									a=np.dot(a,X)

									#得d 和 (2-m)Dx**(m-2)+(L'+U')x**(m-2)

									for j in range(n):

									d[j]=a[j,j]

									a[j,j]=m2*a[j,j]

									#迭代更新

									for j in range(n):

									x[j]=(b[j]-np.dot(a[j],X))/(m1*d[j])

									#判断是否满足精度要求

									if np.max(np.fabs(X-x))<Delta:

									find=1

									break

									X=np.copy(x)

									end=time.perf_counter()#结束计时

									print('时间：',end-start)

									print('迭代',i)

									return X,find,i,end-start<br>

1.2 张量A的生成函数和向量b的生成函数：

				?

									def Creat_A(m,n):#生成张量A

									size=np.full(m, n)

									X=np.ones(n)

									while 1:

									#随机生成给定形状的张量A

									A=np.random.randint(-49,50,size=size)

									#判断Dx**(m-2)是否非奇异，如果是，则满足要求，跳出循环

									D=np.copy(A)

									for i1 in range(n):

									for i2 in range(n):

									if i1!=i2:

									D[i1,i2]=0

									for i in range(m-2):

									D=np.dot(D,X)

									det=np.linalg.det(D)

									if det!=0:

									break

									#将A的对角面张量扩大十倍，使对角面占优

									for i1 in range(n):

									for i2 in range(n):

									if i1==i2:

									A[i1,i2]=A[i1,i2]*10

									print('A:')

									print(A)

									return A

									#由A和给定的X根据Ax**(m-1)=b生成向量b

									def Creat_b(A,X,m):

									a=np.copy(A)

									for i in range(m-1):

									a=np.dot(a,X)

									print('b:')

									print(a)

									return a

1.3 对称张量S的生成函数：

				?

									def Creat_S(m,n):#生成对称张量B

									size=np.full(m, n)

									S=np.zeros(size)

									print('S',S)

									for i in range(4):

									#生成n为向量a

									a=np.random.random(n)*np.random.randint(-5,6)

									b=np.copy(a)

									#对a进行m-1次外积，得到秩1对称张量b

									for j in range(m-1):

									b=outer(b,a)

									#将不同的b叠加得到低秩对称张量S

									S=S+b

									print('S:')

									print(S)

									return S

									def outer(a,b):

									c=[]

									for i in b:

									c.append(i*a)

									return np.array(c)

									return a

1.4 实验一

				?

									def test_1():

									Delta=0.01#精度

									m=3#A的阶数

									n=3#A的维数

									M=200#最大迭代步数

									X_real=np.array( [2,3,4])

									A=Creat_A(m,n) 

									b=Creat_b(A,X_real,m)

									Jacobi_tensor_V2(A,b,Delta,m,n)