python实现LU分解

引言

这是小编第二篇文章啦！其实纠结了很久作者该用什么自称比较好，于是度娘了一番，发现还真是有些讲究。小编或笔者是最常用的了吧。emmm纠结了一番（其实没纠结）还是自称小编吧，自称笔者感觉很客观很严肃，这里还是活泼一点阿哈哈哈～好，进入正题。

正文

因为矩阵分析的老师要求我们用程序实现LU分解，于是小编认真研究了一下LU分解的步骤，因为涉及矩阵运算，首先想到了用python的numpy实现。然后小编就开始码代码了！
经过小编反反复复的纠结与研究，LU分解的原理其实非常简单。首先高斯消元可以得到U，在高斯消元的过程中把每一次消元过程中的乘数（小编自己想的名词，大家领会精神）存储下来并记录该乘数的位置，就是L下三角矩阵中非主对角元素的值，再通过一个循环将主对角元素赋1就能得到L了。
难点就在于如何设置循环！小编为了探索循环的设置，手动拆解了一个三维矩阵：

A = np.array([[2., 2., 3.],
                  [4., 7., 7.],
                  [-2.,4., 5.]])

首先初始化LU为零矩阵

L = np.zeros(shape=(n,n))
U = np.zeros(shape=(n,n))

然后小编就开始拆解了

L[1,0]=A[1,0]/A[0,0]
A[1]=A[1]-L[1,0]*A[0]

L[2,0]=A[2,0]/A[0,0]
A[2]=A[2]-L[2,0]*A[0]  

L[2,1]=A[2,1]/A[1,1]
A[2]=A[2]-L[2,1]*A[1]

因为小编一开始就想到每次消元都是以一行为基准，后几行来消，所以就想到用base行来定义每次作为基准的行，于是，整理出来了base循环：

for base in range(n-1):
    L[1,base]=A[1,base]/A[base,base]
    A[1]=A[1]-L[1,base]*A[base]

一个循环解决了，之后就很容易发现规律了，在一个确定的基准行下进行消元时，是从第base+1行到最后一行，因此就得到另一个循环了：

for base in range(n-1):
    for i in range(base+1,n):
        L[i,base]=A[i,base]/A[base,base]
        A[i]=A[i]-L[i,base]*A[base]

循环问题解决，LU分解函数诞生！！

import numpy as np
def LU_decompose(A):
    n = len(A)
    L = np.zeros(shape=(n,n))
    U = np.zeros(shape=(n,n))
    
    for base in range(n-1):
        for i in range(base+1,n):
            L[i,base]=A[i,base]/A[base,base]
            A[i]=A[i]-L[i,base]*A[base]
    for i in range(n):#range(n) 范围：[0，n-1]
        L[i,i]=1
    U=np.array(A)
    print("L:")
    print(L)
    print("U:")
    print(U)

A = np.array([[2., 2., 3.],
                  [4., 7., 7.],
                  [-2.,4., 5.]])
LU_decompose(A)

其中还遇到一个问题就是range()函数的范围问题，一定要注意！>o<