线性方程组迭代求解——Gauss-Seidel迭代算法(Python实现)
程序员文章站
2022-10-08 17:55:24
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的数据。
结果:
上一篇: 孝康章皇后13岁入宫为妃,15岁生下康熙,为什么她24岁就去世了?
下一篇: AI绘制逼真的叶子图形