Gauss-Seidel迭代算法的Python實現(xiàn)詳解

import numpy as np
import time

def GaussSeidel_tensor_V2(A,b,Delta,m,n,M):
start=time.perf_counter()
find=0
X=np.ones(n)
d=np.ones(n)
m1=m-1
m2=2-m
for i in range(M):
print('X',X)
x=np.copy(X)
#迭代更新
for j in range(n):
a=np.copy(A)
for k in range(m-2):
a=np.dot(a,x)
for k in range(n):
d[k]=a[k,k]
a[k,k]=m2*a[k,k] 
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

def Creat_A(m,n):#生成張量A
size=np.full(m, n)
X=np.ones(n)
while 1:
#隨機(jī)生成給定形狀的張量A
A=np.random.randint(-49,50,size=size)
#判斷Dx**(m-2)是否非奇異，如果是，則滿足要求，跳出循環(huán)
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的對角面張量擴(kuò)大十倍，使對角面占優(yōu)
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根據(jù)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

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進(jìn)行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

def test_2():
Delta=0.01#精度
m=3#A的階數(shù)
n=3#A的維數(shù)
M=200#最大迭代步數(shù)
X_real=np.array( [2,3,4])
A=Creat_A(m,n) 
b=Creat_b(A,X_real,m)
GaussSeidel_tensor_V2(A,b,Delta,m,n)