【数值方法-Python实现】Crout分解+追赶法实现

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

原文链接：https://www.cnblogs.com/aksoam/p/18366119

Codes

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:
    ```python
    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

Vaildation

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.]]

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