自定义求解器之LDLT分解法_UM-仿真秀干货文章

自定义求解器之LDLT分解法

分解法解方程组是一些有限元软件的主流求解器常用的方法，比如PKPM软件就采用这个方法。

对称正定矩阵可以分解为 ,这里为下三角矩阵且主对角元素皆为1。令 ,则。这种分解也是分解的特殊形式。

对于方程组，可以写成，令，依次解三个简单方程组即可。

考察一个3X3矩阵:

之后的矩阵分解算法同分解。

对于多工况分析时形成的荷载矩阵 ,若令 ,则，亦即 ,求解更简单。


import numpy as np

def LDLTSolver(A, b):
    n = len(A)
    D = np.zeros((n))
    for k in range(n):
        D[k] = A[k, k] - np.dot(A[k, 0:k], A[0:k, k])
        for i in range(k+1, n):           
            A[i, k] = ( A[i, k] - np.dot(A[i, 0:k], A[0:k, k]) ) / D[k]
            A[0:i, i] = D[0:i] * A[i, 0:i]            

    for k in range(1,n): 
        A[0:k,k] = 0.0

    for k in range(n):
        A[k,k] = 1

    # 求解 [L]{y} = {b} 
    y = np.zeros((n))
    for k in range(n):
        h = b[k] - np.dot(A[k,0:k], y[0:k]) 
        y[k] = h

    
    # 求解 [D]{z} = {y}
    b = y/D
    # 求解 [L^T]{x} = {z} 
    for k in range(n-1,-1,-1):
        h = b[k] - np.dot(A[k+1:n,k], b[k+1:n])
        b[k] = h

    return b
    


A = np.array([   [1,  0.5, 0.5],
                [0.5,   1, 0.5],
                [ 0.5, 0.5, 1] ])

b = np.array([1, -2, 3])
x = LDLTSolver(A, b)


print(x)

来源：数值分析与有限元编程

自定义求解器之Cholesky分解法

对称正定矩阵可以分解为，这种分解被称为Cholesky分解，是LU分解的一个重要特例，可以显著降低计算量。在计算机程序中常常用到这种方法解线性代数方程组。它的优点是节省存储量，得到的L矩阵可以覆盖原来的A矩阵。对于方程组，可以写成，令，则。考察一个3X3矩阵: 则分解的算法： import numpy as npimport mathclass LinerSolver: def __init__(self, A, b): self.A = A self.b = b def CholeskiSolver(self): n = len(self.A) # 分解 [A] = [L] [L^T] for k in range(n): self.A[k,k] = math.sqrt(self.A[k,k] - np.dot(self.A[k,0:k], self.A[k,0:k])) for i in range(k+1,n): self.A[i,k] = (self.A[i,k] - np.dot(self.A[i,0:k], self.A[k,0:k])) / self.A[k,k] for k in range(1,n): self.A[0:k,k] = 0.0 # 求解 [L]{y} = {b} for k in range(n): self.b[k] = (self.b[k] - np.dot(self.A[k,0:k], self.b[0:k])) / self.A[k,k] # 求解 [L^T]{x} = {y} for k in range(n-1,-1,-1): self.b[k] = (self.b[k] - np.dot(self.A[k+1:n,k], self.b[k+1:n])) / self.A[k,k] return self.bA = np.array([ [ 4, -1, 1], [-1, 4.25, 2.75], [1, 2.75, 3.5] ])b = np.array([4, 6, 7.25])cls = LinerSolver(A, b) #创建一个求解器的实例clsx = cls.CholeskiSolver() #调用Choleski法求解print(x) 与高斯消去法相比，LL分解的优点在于，一旦A被分解，我们就可以对任意多个常量向量b求解Ax=b。每个附加解决方案的成本相对较小，因为前向和后向替换操作比分解过程花费的时间少得多.来源：数值分析与有限元编程