HFUT计算方法实验：用python实现Doolittle法矩阵三角分解

import numpy as np
def dolittle(d_num,arr_a,arr_m,arr_u):      #dolittle三角分解函数
    arr_u[0]=arr_a[0]
    arr_m[0][0]=1
    for i in range(1,d_num):
        arr_m[i][0]=arr_a[i][0]/arr_u[0][0]

    for i in range(1,d_num):    #i表示当前求的这一行
        for j in range(i,d_num):
            sum=0               #j表示当前求的这一列
            for k in range(i):
                sum+=(arr_m[i][k]*arr_u[k][j])
            arr_u[i][j]=arr_a[i][j]-sum
        for j in range(i,d_num):
            if(i==j):
                arr_m[i][i]=1
            if(i!=j):
                sum=0
                for k in range(i):
                    sum += (arr_m[j][k] * arr_u[k][i])
                arr_m[j][i]=(arr_a[j][i]-sum)/arr_u[i][i]

    return arr_m,arr_u

def fun_m(n,a,b,y):     #l*y=b
    y[0]=b[0]
    for i in range(1,n):
        sum=0
        for j in range(0,i):
            sum+=a[i][j]*y[j]
        y[i]=b[i]-sum
    return y

def fun_u(n, a, x, y):  # u*x=y
    i = n - 1
    while i >= 0:
        sum = 0
        j = n - 1
        while j > i:
            sum += a[i][j] * x[j]
            j -= 1
        x[i] = (y[i] - sum) / a[i][i]
        i -= 1
    return x

def fun(a,x,b,n):       #三角分解解方程组
    l = np.zeros((n,n))
    u = np.zeros((n, n))
    l,u=dolittle(n,a,l,u)
    y=np.zeros(n)
    fun_u(n,u,x,fun_m(n,l,b,y))
    return x

if __name__=='__main__':
    print("计算方程组：ax=b")
    a=np.array([[1,2,3],[1,4,6],[2,10,18]])
    print("系数矩阵 a=\n",a)
    l,u=dolittle(3,a,np.zeros((3,3)),np.zeros((3,3)))
    print("单位下三角矩阵 l=\n",l)
    print("上三角矩阵 u=\n",u)
    b=np.array([10,17,44])
    x=np.zeros(3)
    fun(a,x,b,3)
    print("解向量 x=\n",x)