自定义求解器之LDLT分解法
分解法解方程组是一些有限元软件的主流求解器常用的方法,比如PKPM软件就采用这个方法。对称正定矩阵 可以分解为 ,这里 为下三角矩阵且主对角元素皆为1。令 ,则 。这种分解也是 分解的特殊形式。对于方程组 ,可以写成 ,令 ,依次解三个简单方程组即可。考察一个3X3矩阵: 之后的矩阵分解算法同 分解。对于多工况分析时形成的荷载矩阵 ,若令 ,则 ,亦即 ,求解更简单。import numpy as npdef 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)来源:数值分析与有限元编程