列主元LU直接分解法求n阶线性方程组的解（python）

（一）目的

通过设计、编制、调试2~3个求n阶线性方程组数值解的程序，加深对其数值计算方法及有关的基础理论知识的理解。

（二）要求

    用编程语言实现用高斯(Gauss)消元法求n阶线性方程组的解、用列主元高斯(Gauss)消元法求n阶线性方程组的解、用库郎(Courant)列主直接分解法求n阶线性方程组的解的程序。

二、示例

1、问题

用高斯(Gauss)消元法求n阶线性方程组的解。

2、算法描述

（略）

3、程序中变量说明

    （略）

4、源程序清单及运行结果

（略）

5、按以上4点要求编写上机实验报告。

列主元LU分解法基本原理：

 

 

 

 python实现：

from numpy import *

# 列主元LU分解法  求解Ax=b
def LU(A,b):

    # 获取A的行列
    row = A.shape[0]
    column = A.shape[1]
    # 非方阵情况
    if row != column:
        print('A不是方阵！')
        return

    n = row
    Ip = list(range(n))    # 标记主元交换次序

    for k in range(0,n):

        # 寻找最大主元的行数
        h = argmax(abs(A[:,k]))

        # 只能选择未交换的行进行交换
        if Ip[h] == h:
             temp = copy(A[k, :])
             A[k, :] = A[h, :]
             A[h, :] = temp
             Ip[k] = h

        # 对角线减去对应行列乘积和
        A[k,k] = A[k,k] - dot(A[k,0:k], A[0:k,k])
        if A[k,k] == 0:
            print('第{}个主元为0！'.format(k+1))
            return

        # 直接在A上保存L和U
        for j in range(k+1,n):    # j>k
            A[k,j] = A[k,j] - dot(A[k,0:k], A[0:k,j])
            A[j,k] = (A[j,k] - dot(A[j,0:k], A[0:k,k])) / A[k,k]

    # 单位下三角L和上三角U
    L = tril(A,-1) + identity((n))
    U = triu(A,0)
    L.dtype = double
    U.dtype = double
    print(L,"\n",U)
    # 根据Ip[]的标记，对b进行对应行交换
    for i in range(0,n):
        t = Ip[i]
        if t != i:
            temp = copy(b[i])
            b[i] = b[t]
            b[t] = temp

    # 求解Ly=b
    y = transpose(array([arange(n)], dtype = double))
    y[0] = b[0]
    for i in range(1,n):
        sum = 0
        for t in range(0,i):
            sum += L[i,t]*y[t]
        y[i] = b[i] - sum

    # 求解Ux=y
    x = transpose(array([arange(n)],dtype = double))
    x[n-1] = y[n-1]/U[n-1,n-1]
    for i in range(0,n-1)[::-1]:
        sum = 0
        for t in range(i+1,n):
            sum +=  U[i,t]*x[t]
        x[i] = (y[i]-sum)/U[i,i]

    return x

if __name__ == '__main__':
    A = array([[2.0, -1.0, 3.0, 2.0], [3.0, -3.0, 3.0, 2.0], [3.0, -1.0, -1.0, 2.0], [3.0, -1.0, 3.0, -1.0]], dtype = double)
    b = transpose(array([[6.0, 5.0, 3.0, 4.0]], dtype = double))
    print(LU(A,b))

结果如下：

列主元LU分解法的输出结果：

E:\shuzhifenxishiyan\liezhuyuanLU\venv\Scripts\python.exe E:/shuzhifenxishiyan/liezhuyuanLU/main.py

[[1.         0.         0.         0.        ]

 [0.66666667 1.         0.         0.        ]

 [1.         2.         1.         0.        ]

 [1.         2.         0.33333333 1.        ]]

 [[ 3.         -3.          3.          2.        ]

 [ 0.          1.          1.          0.66666667]

 [ 0.          0.         -6.         -1.33333333]

 [ 0.          0.          0.         -3.88888889]]

[[1.]

 [1.]

 [1.]

 [1.]]

Process finished with exit code 0