python实现LU分解

数学原理

对线性方程组

$$A_{n\times n}{X}_{n\times 1}=b_{n\times 1}$$

求解思路为进行$n$次迭代，在第$k$次迭代中计算U的第$k$行，L的第$k$列。
于是编程思路为两个for循环嵌套，外层循环控制迭代次数，内层循环计算L，U矩阵行列具体值。

$$U_{k,i}=A_{k,i}-(L_{k,0}{U}_{0,i}+\cdots +L_{k,k-1}{U}_{k-1,i}),i=k,k+1,\cdots ,n.$$
$$L_{j,k}=(A_{j,k}-L_{j,0}{U}_{0,k}-\cdots -L_{j,k-1}{U}_{k-1,k})/U_{k,k},j=k+1,\cdots n.$$

注：为了便于编程，公式中矩阵第一个元素索引为0。

Python代码

/'''
功能：使用LU分解求解线性方程组                                      |符号
输入：线性方程组的系数矩阵和常数阵                                   |A,b
输出：系数矩阵A经分解后得到的L,U矩阵，以及线性方程组的解向量x          |L,U,x
'''
import numpy as np

A = np.array([[4, 2, 1, 5],
              [8, 7, 2, 10],
              [4, 8, 3, 6],
              [6, 8, 4, 9]])
b = np.array([[7, 4, 6, 8]])
n = np.shape(A)[0]
L = np.zeros((n, n)) + np.eye(n)
U = np.zeros((n, n))
for k in range(n):
    for i in range(k, n):
        U[k, i] = A[k, i] - L[k, 0:k].dot(U[0:k, i])
    for j in range(k + 1, n):
        L[j, k] = (A[j, k] - L[j, 0:k].dot(U[0:k, k])) / U[k, k]
 A_inv = np.linalg.inv(U).dot(np.linalg.inv(L))
x = A_inv.dot(b.T)
print("L:\n{}".format(L))
print("U:\n{}".format(U))
print("x:{}".format(x))