LU分解解方程程序(python,数值积分)

第二课 因式分解解方程组

如果有看过第一课的朋友，可以注意到，如果在系数矩阵相同，而右边向量不同的情况下，我们需要每次都进行消元去求解，这样会很麻烦。
所以针对系数相同，右边向量b不同的时候，可以使用因式分解的方法更加快捷。
在矩阵中，矩阵A总是能被写成一个下三角矩阵与一个上三角矩阵乘积的形式。也就是

其中[L]为下三角矩阵，[U]为上三角矩阵，表达为下面的形式

系数主对角项的乘积为

所以一般情况下，我们都选择将上三角或者下三角的主对角线值设置为1，也就是

正常一个线性方程的求解形式为

而把分解成上下三角形式的A矩阵代入为

如果让

那么上边方程将变为

下三角[L]和{b}是已知的，这时可以采取’向前迭代法‘(不懂可以看第一课)，求出{y}向量值。当{y}和[U]是已知的，可以采取’向后迭代法‘(不懂可以看第一课)，这样就求出未知量{x}。所以这个过程分为三步，第一步，因式分解；第二步，向前迭代；第三步，向后迭代。
详细的因式分解过程为

第一行： U11=A11,U12=A12,U13=A13
第二行： L21U11=A21,因此L21=A21/U11
L21U12+U22=A22, 因此U22=A22-L21U12
L21U13+U23=A23,因此U23=A23-L21U13
第三行：
L31U11=A31,因此L31=A31/U11
L31U12+L32U22=A32,
因此L32=(A32-L31U12)/U22
L31U13+L32U23+U33=A33,
因此U33=A33-L31U13-L32U23

下面以一个具体算例来详解计算过程：
把上面的系数矩阵分解成下三角和上三角乘积的形式，也就是

求解下三角和上三角，得到

‘向前迭代’[L]{y}={b}，得到

‘向后迭代’[U]{x}={y}，得到

此部分把代码分为一个主程序和三个子程序，好处是为了在之后的程序中可以直接调用子程序，三个子程序分别为因式分解的程序lufac，向后迭代法的程序subbac，向前迭代法的程序subfor。所有子程序全部放入python文件b中，所以以后调用子程序，只需采用import B.子程序名。
下面给出python代码：
主程序：

#使用LU分解的高斯消元法
import numpy as np
import math
import B
n=3
upper=np.zeros((n,n))
lower=np.zeros((n,n))
a=np.array([[2,-3,1],[1,-1<