欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页  >  IT编程

线性方程组迭代求解——Gauss-Seidel迭代算法(Python实现)

程序员文章站 2022-04-10 21:07:28
Jacobi迭代求解线性方程组:Gauss-Seidel实现 ......

原理:

请看本人博客:线性方程组的迭代求解算法——原理

代码:

 

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的数据。

线性方程组迭代求解——Gauss-Seidel迭代算法(Python实现)

结果:

 

线性方程组迭代求解——Gauss-Seidel迭代算法(Python实现)