矩陣的求解可以分為直接求解方法和迭代求解方法，這里對部分常用的方法進行了列舉主要分為迭代法和直接求解法兩類，語言基于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))

喜歡的朋友可以給個關注或者聯系我