線性方程組迭代算法——Gauss-Seidel迭代算法的python實現

原理：

請看本人博客：線性方程組的迭代求解算法——原理

代碼：

 

import numpy as np
max=100#迭代次數上限
Delta=0.01
m=2#階數：矩陣為2階
n=3#維數：3X3的矩陣

shape=np.full(m, n)
shape=tuple(shape)

def read_tensor(f,shape):#讀取張量
    data=[]
    for i in range(n**(m-1)):
        line = f.readline()
        data.append(list(map(eval, line.split(","))))
    return np.array(data).reshape(shape)
    
def read_vector(f):#讀取向量
    line = f.readline()
    line = line.replace("\n","")
    line=list(map(eval, line.split(",")))
    return np.array(line)
    
#讀取數據    
f = open("jacobi_data.txt")
A=read_tensor(f,shape)#讀取矩陣A
b=read_vector(f)#讀取b
f.close()
print('A:')
print(A)
print('b:',b)

U=np.copy(A)#求U
DL=np.copy(A)#求D-L
for i in range(n):
    for j in range(n):
        if j<=i:
            U[i,j]=0
        else:
            DL[i,j]=0
U=0-U

#迭代求解
x=np.ones(n)#用於存儲迭代過程中x的值
y=np.ones(n)#用於存儲中間結果
DLU=np.dot(np.linalg.inv(DL),U)#對DL求逆，然后和U相乘
DLb=np.dot(np.linalg.inv(DL),b)#對DL求逆，然后和b相乘
print('x:',x)
for iteration in range(max):
    #迭代計算
    y=np.dot(DLU,x)+DLb
        
    #判斷是否達到精度要求    
    if np.max(np.fabs(x-y))<Delta:
        print('iteration:',iteration)
        break
    #將y幅值到x，開始下一輪迭代    
    x=np.copy(y)  
    print('x:',x)



    
    
    
    
    
    

 

數據：

組織形式：前3行是A的數據，最后1行是b的數據。

結果：