返回 主页

线性方程组求解 - 直接解法

在有限元分析中,经过单元推导—>组装—>边界条件之后,最终求解的方程一般有形如【一般[A]的大小在1w+】:

[A]X=b\tag{1}
Img

对于规模线性方程组的求解,直接法(如高斯消元法,Choleski 分解法等)比迭代法(如牛顿法,SOR法等)会更加受欢迎.本文介绍高斯消元法和Choleski法的数值求解方法.

高斯消元法

公式1展开后如下,上标(0)表示初始值.

Img

首先,消元 x_1,有:

Img

将 x_1 带入剩下的公式,有:

Img

接下,消元 x_2 和其余 x_p(p=3...n) .可以发现,当消元 x_k 时:

Img

经过n-1次消元后,有:

Img

这样就可以以 x_n,x_{n-1},...,x_{1} 的顺序完成求解.

如果每一步的消元过程中,如果对角系数 a_{11}^{(0)},a_{22}^{(1)},a_{33}^{(2)}.... 出现0值,我们需要进行行变换来得到非0值.如果行变换行不通,那就说明[A]是奇异矩阵

Choleski法

该方法基于一个事实:任何方阵[a]都可以表示为上下三角矩阵的乘积.

LU分解

系统方程为 [A]X=b ,[A]可以写作:

Img
Img
Img
Img

满足A=LU的唯一因式分解的[L],[U]可以逐个元素求得,求解顺序是:

Img

l_{ij},u_{ij} 满足:

Img

至此,完成LU分解.

求解方程组

将A=LU带入到[A]X=b可得:LUX=b.方程求解步骤为:

令:

[U]\vec{X}=\vec{Z}

其中Z满足:

[L]\vec{Z}=\vec{b}

展开后:

Img

以上,第一个式子可以求出z1,之后就可以依次求出z2,z3,....,zn【只要 l_{ii} 中没有0值】,一旦 z_i 求得后,则有:

Img

类似高斯消元,以 x_n,x_{n-1},...,x_{1} 的顺序可以完成求解.

对称矩阵的CHOLESKI分解

在大多数的有限元分析中,[A]一般是对称正定矩阵.此时,A有唯一的分解式:

Img

矩阵[U]由下式确定:

Img

对称矩阵的逆

首先进行 A=U^TU 分解,然后:

Img

U^{-1} 满足 U^{-1}U=[I] , U^{-1} 的元素 \lambda_{ij} 满足:

Img

可以看出 U^{-1} 也是上三角阵, U^{-1} 的逆为:

Img

最后,A的逆为:

Img

Crout分解追赶法的编程实现

涉及Crout分解、追赶法的线性方程组求解方法的Python实现.

Python代码

def CroutLU(A:np.ndarray)->Tuple[np.ndarray,np.ndarray]:
    """Crout LU分解算法,A=LU
    input: 
        A (n,n) np.ndarray,方阵
    
    output:
        L: (n,n) np.ndarray,下三角矩阵
        U: (n,n) np.ndarray,上三角矩阵,对角线元素为1.0
    usage:
        A=np.array([[1,2,3,4],
                [1,4,9,16],
                [1,8,27,64],
                [1,16,81,256]])
        L,U=CroutLU(A)
        print("L矩阵:\n", L)
        print("U矩阵:\n", U)
        # 验证分解是否正确
        print("验证LU是否等于A:\n", np.dot(L, U))

    输出:
        L矩阵:
        [[ 1.  0.  0.  0.]
        [ 1.  2.  0.  0.]
        [ 1.  6.  6.  0.]
        [ 1. 14. 36. 24.]]
        U矩阵:
        [[1. 2. 3. 4.]
        [0. 1. 3. 6.]
        [0. 0. 1. 4.]
        [0. 0. 0. 1.]]
        验证LU是否等于A:
        [[  1.   2.   3.   4.]
        [  1.   4.   9.  16.]
        [  1.   8.  27.  64.]
        [  1.  16.  81. 256.]]
    """
    row,col=A.shape
    n=row
    L=np.zeros((n,n))
    U=np.zeros((n,n))
    for k in range(n):
        # 循环列,从k+1列到n列,i=k,...n-1
        for i in range(k,n):
            L[i,k]=A[i,k]-sum([L[i,s]*U[s,k] for s in range(k)])
        
        for j in range(k-1,n):
            U[k,j]=(A[k,j]-sum([L[k,s]*U[s,j] for s in range(k)]))/L[k,k]
    return L,U

def LUChaseMethod(A:np.ndarray,d:np.ndarray)->np.ndarray:
    """LU分解法,追赶法求解线性方程组Ax=d
    求解三对角矩阵A,d的线性方程组Ax=d,其中A为三对角矩阵,d为右端常数
    """
    n=A.shape[0]
    # x: x1,x2...xn
    x=np.zeros(n)
    
    a=np.zeros(n)
    # a:a1,a2...an
    # b:b1...bn
    # c:c1...cn-1
    a[1:],b,c=np.diag(A,k=-1).copy(),np.diag(A,k=0).copy(),np.diag(A,k=1).copy()
    
    L,U=CroutLU(A)
    # size: (n-1,) , u0,u1...u(n-1),u0=0
    us=np.zeros(n)
    us[1:]=np.diag(U,k=1).copy()
    # size: (n,) ,ls : l1...ln
    ls=np.diag(L,k=0).copy()
    # size: (n-1,) , v2...vn-1
    vs=np.diag(L,k=-1).copy()
    # y: y0,y1...yn , y0=0
    y=np.zeros(n+1)
    # i=1....n-1
    for i in range(1,n):
        # print(f"第{i}次迭代")
        y_i_1=y[i-1]
        a_i=a[i-1]
        b_i=b[i-1]
        c_i=c[i-1]
        d_i=d[i-1]
        u_i_1=us[i-1]
        
        l_i=b_i-a_i*u_i_1
        u_i=c_i/l_i
        y_i=(d_i-a_i*y_i_1)/l_i
        
        y[i]=y_i
        ls[i-1]=l_i
        us[i]=u_i
    ls[-1]=b[n-1]-a[n-1]*us[n-1]
    y[n]=(d[n-1]-y[n-1]*a[n-1])/ls[n-1]
    
    x[n-1]=y[n]
    for i in range(n-2,-1,-1):
        # xi=yi-ui*x(i+1),i=n-2...1
        x[i]=y[i+1]-us[i+1]*x[i+1]
    return x

验证

from formu_lib import *
import numpy as np
import matplotlib.pyplot as plt
A=np.array([[1,2,3,4],
            [1,4,9,16],
            [1,8,27,64],
            [1,16,81,256]])

L,U=CroutLU(A)

print("L矩阵:\n", L)
print("U矩阵:\n", U)

# 验证分解是否正确
print("验证LU是否等于A:\n", np.dot(L, U))
L矩阵:
 [[ 1.  0.  0.  0.]
 [ 1.  2.  0.  0.]
 [ 1.  6.  6.  0.]
 [ 1. 14. 36. 24.]]
U矩阵:
 [[1. 2. 3. 4.]
 [0. 1. 3. 6.]
 [0. 0. 1. 4.]
 [0. 0. 0. 1.]]
验证LU是否等于A:
 [[  1.   2.   3.   4.]
 [  1.   4.   9.  16.]
 [  1.   8.  27.  64.]
 [  1.  16.  81. 256.]]
img
from formu_lib import *
import numpy as np
import matplotlib.pyplot as plt
a=np.array([[2,-1,0,0],[-1,3,-2,0],[0,-2,4,-3],[0,0,-3,5]])
d=np.array([6,1,-2,1])
x=LUChaseMethod(a,d)
print(f"x={x}")
# x=[5. 4. 3. 2.]

答案:[5. 4. 3. 2.],ok

返回 主页