自定义求解器之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)