首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >自定义求解器之LDLT分解法

自定义求解器之LDLT分解法

作者头像
fem178
发布2021-04-16 15:39:08
1.1K0
发布2021-04-16 15:39:08
举报
LDL^T

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

对称正定矩阵

A

可以分解为

A = LDL^T

,这里

L

为下三角矩阵且主对角元素皆为1。令

U^{*}=DL^T

,则

A=LU^{*}

。这种分解也是

LU

分解的特殊形式。

对于方程组

Ax=b

,可以写成

LDL^Tx=b

,令

Ly=b,Dz=y,L^Tx=z

,依次解三个简单方程组即可。

考察一个3X3矩阵:

\begin{bmatrix} A_{11} & A_{12} & A_{13}\\ A_{21} & A_{22} & A_{23}\\ A_{31} & A_{32} & A_{33}\\ \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0\\ L_{21} & 1 & 0\\ L_{31} & L_{32} & 1\\ \end{bmatrix} \begin{bmatrix} D_1 & 0 & 0\\ 0 & D_2 & 0\\ 0 & 0 & D_3\\ \end{bmatrix} \begin{bmatrix} 1 & L_{21} & L_{31}\\ 0 & 1 & L_{32}\\ 0 & 0 & 1\\ \end{bmatrix}
U^{*}=DL^T= \begin{bmatrix} D_1 & D_1L_{21} & D_1L_{31}\\ 0 & D_2 & D_2L_{32}\\ 0 & 0 & D_3\\ \end{bmatrix}

之后的矩阵分解算法同

LL^T

分解。

对于多工况分析时形成的荷载矩阵

P

,若令

P=LDQ

,则

LDL^Tu= LDQ

,亦即

L^Tu=Q

,求解更简单。


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)

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2021-04-12,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 数值分析与有限元编程 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档