矩阵迭代求解方法

矩阵的求解可以分为直接求解方法和迭代求解方法，这里对部分常用的方法进行了列举主要分为迭代法和直接求解法两类，语言基于，把这些算法做了相应的封装写成了一个类，并给出了一些方程组供测试，很多算法来自互联网，就不一一致谢了。

迭代求解算法，主要有雅克比，高斯-赛德尔，超松弛，共轭梯度，并与numpy库中自带的一些函数做了对比，主要是求解时间和求解结果的。有时间再把原理相关的描述写下来。

"""
Created on 2022/1/25 10:45
@Author: 123
@Email: 985040320@qq.com
@Describe: **加入文件描述, 要实现的功能, 注意事项等**
"""
import numpy as np
import time
from scipy.sparse import csr_matrix
import scipy
class 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))

———END———
限 时 特 惠： 本站每日持续更新海量各大内部创业教程，永久会员只需109元，全站资源免费下载 点击查看详情
站 长 微 信： nanadh666