自定義求解器之LDLT分解法
瀏覽:2458
分解法解方程組是一些有限元軟件的主流求解器常用的方法,比如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)
技術鄰APP
工程師必備
工程師必備
- 項目客服
- 培訓客服
- 平臺客服
TOP




















