自定義求解器之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.50.5],
                [0.5,   10.5],
                [ 0.50.51] ])

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


print(x)

登錄后免費查看全文
立即登錄
App下載
技術鄰APP
工程師必備
  • 項目客服
  • 培訓客服
  • 平臺客服

TOP