几种常用的矩阵迭代求解算法
矩阵的求解可以分为直接求解方法和迭代求解方法，这里对部分常用的方法进行了列举主要分为迭代法和直接求解法两类，语言基于python，把这些算法做了相应的封装写成了一个类，并给出了一些方程组供测试，很多算法来自互联网，就不一一致谢了。
迭代求解算法，主要有雅克比，高斯-赛德尔，超松弛，共轭梯度，并与numpy库中自带的一些函数做了对比，主要是求解时间和求解结果的。有时间再把原理相关的描述写下来。
"""Created on 2022/1/25 10:45@Author: 123@Email: 985040320@qq.com@Describe: **加入文件描述, 要实现的功能, 注意事项等**"""import numpy as npimport timefrom scipy.sparse import csr_matriximport scipyclass Iteration:
def __init__(self, A, x, b):
self.A = A
self.b = b
self.x = x
def Jacobi(self, n):
# 设Ax= b，其中A=D L U为非奇异矩阵，且对角阵D也非奇异，则当迭代矩阵J的谱半径ρ（J）<1时，雅克比迭代法收敛。
times = 0
# D = np.diag(np.diag(self.A))
L = -np.array(np.tril(self.A, -1))
U = -np.array(np.triu(self.A, 1))
D_inv = np.diag(1/np.diag(self.A))
while times < n:
xnew = self.x
self.x = D_inv.dot(self.b L.dot(self.x) U.dot(self.x))
if abs(self.x-xnew).max() < 0.000000001:
break
times = 1
return self.x.flatten()
def Gauss_Seidel(self, n):
# 需要系数矩阵对称正定或者严格对角占优
times = 0
D = np.diag(np.diag(self.A))
L = -np.array(np.tril(self.A, -1))
U = -np.array(np.triu(self.A, 1))
D_L_inv = np.linalg.inv((D - L))
while times < n:
xnew = self.x
self.x = D_L_inv.dot((b U.dot(self.x)))
if abs(self.x-xnew).max() < 0.000000001:
break
times = 1
return self.x.flatten()
def Successive_Over_Relaxation(self, omega, n):
# 限定条件较少，适用性更普遍
times = 0
D = np.diag(np.diag(self.A))
L = -np.array(np.tril(self.A, -1))
U = -np.array(np.triu(self.A, 1))
D_omega_inv = np.linalg.inv((D - omega * L))
while times < n:
xnew = self.x
self.x = D_omega_inv.dot((omega * b) ((1 - omega) * D).dot(self.x) (omega * U).dot(self.x))
if abs(self.x-xnew).max() < 0.0000000001:
break
times = 1
return self.x.flatten()
def conj_gradient(self, tol, limit):
# 系数矩阵必须对称正定，对称正定可以Cholesky分解
n = self.A.shape[0]
p = np.zeros((n, 1))
r = np.zeros((n, 1))
u = np.zeros((n, 1))
xnew = np.zeros((n, 1))
r = b - np.dot(A, self.x) # 计算r0
p = r.copy() # 计算p0
iters = 0
while True:
iters = iters 1
u = np.dot(self.A, p)
up = np.dot(r.T, r)
alpha = np.dot(r.T, r) / np.dot(p.T, u)
# print("alpha", alpha.flatten())
r = r - u * alpha
# print("r", r.flatten())
xnew = self.x p * alpha
# print("x", xnew.flatten())
beta = np.dot(np.transpose(r), r) / up
p = r p * beta
# print("p", p.flatten())
if abs(xnew - self.x).max() < tol or iters == limit:
break
self.x = xnew
return self.x.flatten()if __name__ == "__main__":
n = 800
emega = 0.5
tol = 0.0000001
itertimemax = 100
# 测试方程组1
A = np.array([[2, -1, 0], [-1, 3, -1], [0, -1, 2]])
x = np.array([[1], [8], [5]])
b = np.array([[1], [8], [5]])
# 测试方程组2
# A = np.array([[1., 1.], [1., -1.]])
# x = np.array([[0.], [0.]])
# b = np.array([[5.], [1.]])
# 测试方程组3
# A = np.array([[16, 4, 8], [4, 5, -4], [8, -4, 22]], dtype=float)
# b = np.array([[4], [2], [5]], dtype=float)
# x = np.array([[1], [1], [1]], dtype=float)
#测试方程组4
# A = np.array([[1, -1, 0, 0], [-0.1, 1, -0.9, 0], [0, -0.9, 1, -0.1], [0, 0, 0, 1]], dtype=float)
# b = np.array([[9], [0], [0], [0]], dtype=float)
# x = np.array([[19], [10], [1],[0]], dtype=float)
# 共轭梯度
t1 = time.time()
Obj = Iteration(A, x, b)
results1 = Obj.conj_gradient(tol, itertimemax)
t2 = time.time()
# 超松弛
t1 = time.time()
Obj = Iteration(A, x, b)
results = Obj.Successive_Over_Relaxation(emega, n)
t2 = time.time()
print(results, f"超松弛耗时{t1-t2}")
# 高斯赛德尔
t1 = time.time()
Obj = Iteration(A, x, b)
results = Obj.Gauss_Seidel(n)
t2 = time.time()
print(results, f"高斯赛德尔耗时{t1-t2}")
# 雅克比
t1 = time.time()
Obj = Iteration(A, x, b)
results = Obj.Jacobi(n)
t2 = time.time()
print(Obj.Jacobi(n), f"雅克比耗时{t1-t2}")
# 测试函数
t1 = time.time()
results = np.linalg.solve(A, b)
t2 = time.time()
print(results.flatten(), f"测试函数耗时{t1 - t2}")除了上面所用到的迭代算法，这里还介绍一种处理CFD这种会遇到的三对角或者五对角矩阵的迭代求解算法，三对角算法，也算迭代算法只不过这种矩阵刚好容易出现在网格离散之后的方程组里面。
import numpy as np class IntSlve: def __init__(self, A, b): self.A = A.copy() self.b = b.copy() self.nf = len(b) def TDMASolver(self): a_1 = np.zeros(self.nf-1) b_1 = np.zeros(self.nf) c_1 = np.zeros(self.nf) for i in range(self.nf): # 矩阵分解为三条对角线上的元素 if i < self.nf-1: c_1[i] = self.A[i, i 1] a_1[i] = self.A[i 1, i] b_1[i] = self.A[i, i] else: b_1[i] = self.A[i, i] ac, bc, cc, dc = map(np.array, (a_1, b_1, c_1, self.b)) for it in range(1, self.nf): mc = ac[it - 1] / bc[it - 1] bc[it] = bc[it] - mc * cc[it - 1] dc[it] = dc[it] - mc * dc[it - 1] xc = bc xc[-1] = dc[-1] / bc[-1] for il in range(self.nf - 2, -1, -1): xc[il] = (dc[il] - cc[il] * xc[il 1]) / bc[il] return xc if __name__ == "__main__": A = np.array([[10, 2, 0, 0],[3,10,4,0],[0,1,7,5],[0,0,3,4]],dtype=float) d = np.array([3, 4, 5, 6.]) sol = IntSlve(A, d) print(sol.TDMASolver()) # 数据对比 print(np.linalg.solve(A, d))
下期写直接求解的相关算法
赞同 6
添加评论
分享
收藏一键生成视频收起
有限体积方法处理非结构网格-多边形和多面体的几何的处理
无语
上一篇文章中已经讲了通过离散求解的流程，目的就是先构造出每个网格的代数方程，代数方程的具体形式如下，碰到的问题就是需要求解三个系数项bc，af，ac
但是这三个系数项就是采用有限体积方法离散获取的，每个都有自己的特点， 有的需要当前网格的值或者这个值的梯度，有的需要与当前网格相邻的单元的值。而且前面也列出了这三个公式的具体表达式，如下
![[公式]](https://public.fangzhenxiu.com/fixComment/commentContent/imgs/1645200537668_abw0gp.jpg?imageView2/0)
前面我们已经知道求解这三项只需要知道gDiff和Tf还有▽ϕ三个就可以了这三个就可以了，而且这三相都是和网格相关的，而且都需要长篇大论的写，所以我们就一个一个的处理，首先是这个gDiff吧，先给出一个网格的图片，图片来自互联网哈。
一个典型的非结构网格
这里把公式先写出来
Ef是表面向量的类正交扩散分量，dCF是两个单元的质心距离，e是两个网格单元质心向量的单元向量，Sf是面的法向向量,已经知道的是网格中所有的点的坐标都。
首先来完成单元的质心计算吧，当然质心的计算其实也是一个比较漫长的步骤，首先是把每个面的面心计算出来，对于三角形来说每个面的面心直接加权就可以，但是对于多边形来说还需要转化为三角形再处理。
1.大致流程是这样的，通过面上的顶点坐标加和然后除以顶点的数目可以计算到几何中心
2.将几何中心与边上的两个点相连得到n边形对应的n个三角形
3.将这n个三角形的面矢量和面积分别计算出来，还有三角形的几何中心
4.将这n个三角形的几何中心和对应的面积加权就可以得到面的面心，三角形面积矢量加和除以2就是这个面的面积矢量Sf。
为了省去你的麻烦，我把函数已经定义好了，当然你也可以自己编写，最好以类的形式写出来，计算结果也存为类里面的属性，可以减少计算量。
"""Created on 2022/1/25 10:45@Author: 123@Email: 985040320@qq.com@Describe: **加入文件描述, 要实现的功能, 注意事项等**""" def face_center_and_area(vertices):#传构成这面的点的列表 # 计算k个点组成的多边形的几何中心 size = len(vertices) geo_center = np.sum(vertices, axis=0) / size face_total_area = 0.0 face_normal_vector = np.zeros((3,)) face_centroid = np.zeros((3,)) # 形成三角子面. # 每一个面都是以边作为基础，以面的几何中心作为顶点然后重建 for i in range(len(vertices)): # 计算子面的点 p1 = vertices[i] p2 = vertices[(i 1) % size] p3 = geo_center # 计算子面的几何中心 subface_geo_center = np.sum([p1, p2, p3], axis=0) / 3.0 # 计算子面的面积和面向量，三角形面积等于面矢量的模长除以2 sf = np.cross((p2 - p1), (p3 - p1)) area = np.linalg.norm(sf) / 2.0 # 计算多边形面的面积和面矢量 face_normal_vector = sf face_total_area = area face_centroid = face_centroid (area * subface_geo_center) face_centroid /= face_total_area face_normal_vector /= 2.0 return face_centroid, face_total_area, face_normal_vector
可见通过上面的四个步骤，我们由点的坐标计算出了面的面积，面心，面心矢量。
面解决了，接下来就是计算体心了，类似于多边形转化为三角形再进行向量求面积，求面心，求面矢量，多面体需要转化为像金字塔的多面体（底面的面积面矢量，面心都已知）之后再进行质心的计算
1.首先由上面的面心计算多面体几何中心，直接加取平均即可
2.通过这个几何中心和各个面的面心（上一步计算给出的）可以构造一个多面体并计算出他的质心（在基面和几何中心的0.25处），这里通过面矢量和这两个点也可以顺便把这个有点像金字塔的几何的体积同时求出。
3.体积加和就是整个网格的体积，第二部求得的质心进行对应的体积加权就可以得到网格的质心。
"""Created on 2022/1/25 10:45@Author: 123@Email: 985040320@qq.com@Describe: **加入文件描述, 要实现的功能, 注意事项等**""" def cell_center_and_volume(faces):#传入面的列表 # 计算单元的几何中心，就是简单的采用所有点的进行平均值计算 geo_centeroid = np.zeros((3,)) for face in faces: geo_centeroid = face_center_and_area(face)[0] geo_centeroid /= len(faces) # 根据每个网格面和几何中心创建类似金字塔的多面体 element_volume = 0.0 element_centroid = np.zeros((3,)) for face in faces: face_centroid, face_area, _ = face_center_and_area(face) # 四面体的中心落在基面和几何中心距离的0.25处 pyramid_centroid = (0.75 * face_centroid) (0.25 * geo_centeroid) pyramid_volume = (1.0 / 3.0) * face_area * np.linalg.norm(face_centroid - geo_centeroid) # 借助四面体的中心计算几何体的体积 element_volume = pyramid_volume # 所有的面计算完成之后用用对应的体积去加一个权 element_centroid = (pyramid_volume * pyramid_centroid) # 最后计算体积加权网格中心 element_centroid = element_centroid / element_volume return element_centroid, element_volume
看现在网格的质心已经求出来了，dCF是两个单元的质心距离就是这两个点构成的向量的模嘛，e就是这两个点构成的向量再除以这个模嘛，sf就是这个面对应的面矢量嘛,所以就可以计算gDiff了。
现在就已经完成了gDiff和Tf还有▽ϕ三个中一个系数的获取，是不是就是通过网格的拓扑关系获取的。后面再依次来处理剩余的两个，同样还是通过网格的拓扑结构。
登录后免费查看全文
著作权归作者所有，欢迎分享，未经许可，不得转载
首次发布时间：2022-02-21
最近编辑：2年前
硕士
|
工程师
签名征集中
课程培训服务行家
课程
培训
服务
行家
