多線性方程組迭代演算法——Jacobi迭代演算法的Python實現

多線性方程（張量）組迭代演算法的原理請看這裡：若想看原理部分請留言，不方便公開分享

Gauss-Seidel迭代演算法：多線性方程組迭代演算法——Gauss-Seidel迭代演算法的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

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)