引言
这是小编第二篇文章啦!其实纠结了很久作者该用什么自称比较好,于是度娘了一番,发现还真是有些讲究。小编或笔者是最常用的了吧。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<