Jacobi迭代算法的Python实现详解

Python  /  管理员 发布于 7年前   181

import numpy as npimport 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-1m2=2-mfor 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=1break X=np.copy(x)end=time.perf_counter()#结束计时print('时间：',end-start)print('迭代',i)return X,find,i,end-start

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

def Creat_A(m,n):#生成张量Asize=np.full(m, n)X=np.ones(n)while 1:#随机生成给定形状的张量AA=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]=0for 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]*10print('A:')print(A)return A#由A和给定的X根据Ax**(m-1)=b生成向量bdef 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):#生成对称张量Bsize=np.full(m, n)S=np.zeros(size)print('S',S)for i in range(4):#生成n为向量aa=np.random.random(n)*np.random.randint(-5,6)b=np.copy(a)#对a进行m-1次外积，得到秩1对称张量bfor j in range(m-1):b=outer(b,a)#将不同的b叠加得到低秩对称张量SS=S+bprint('S:')print(S)return Sdef 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)

以上就是本文的全部内容，希望对大家的学习有所帮助，也希望大家多多支持。