有两种方法，一个是通过高斯消去法的方式去求一个下三角矩阵，还有一个是直接分解。
显然直接分解好啊233。
然后这是个啥杜尔利特分解公式。
是由矩阵乘法的原理弄出来的。
书上P99页俩公式(5.34 5.35)。
轮流求u和l
三角分解部分的代码：
1 import numpy as np 2 3 mat = [[2,2,3],[4,7,7],[-2,4,5]] 4 5 height = len(mat) 6 width = len(mat[0]) 7 8 # 1 直接分解矩阵 9 # 初始化LU 10 U = [[0 for i in range(width)] for j in range(height)] 11 L = [[0 if i!=j else 1 for i in range(width) ] for j in range(height)] 12 13 print(L,U) 14 15 # 从U第一行开始，对于L来说是从第一列开始 16 for i in range(height): 17 # 依次计算LU的行、列 18 for j in range(i, width): 19 _sum = 0 20 for k in range(i): 21 _sum += L[i][k] * U[k][j] 22 print("U{}{}_sum += L[{}][{}] * U[{}][{}] = {} * {} = {}".format(i,j,i,k,k,j,L[i][k],U[k][j],L[i][k] * U[k][j])) 23 U[i][j] = mat[i][j] - _sum 24 # L的i.j是反着的 25 for j in range(i, width): 26 _sum = 0 27 if j > i: 28 for k in range(i): 29 _sum += L[j][k] * U[k][i] 30 print("L{}{}_sum += L[{}][{}] * U[{}][{}] = {} * {} = {}".format(j,i,j,k,k,i,L[j][k],U[k][i],L[j][k] * U[k][i])) 31 L[j][i] = (mat[j][i] - _sum) / U[i][i] 32 print("L = {}".format(np.array(L))) 33 print("U = {}".format(np.array(U))) 34 35 # 分解完毕
解Ly =b, 回带求解，再解Ux=y, 直接丢完整代码了：
1 import numpy as np 2 # mat = eval(input("请输入系数矩阵：")) 3 # 测试数据[[2,2,3],[4,7,7],[-2,4,5]] 4 mat = [[2,2,3],[4,7,7],[-2,4,5]] 5 # b = eval(input("请输入结果矩阵:")) 6 # 测试数据[[3],[1],[-7]] 7 b = [[3],[1],[-7]] 8 height = len(mat) 9 width = len(mat[0]) 10 11 # 1 直接分解矩阵 12 # 初始化LU 13 U = [[0 for i in range(width)] for j in range(height)] 14 L = [[0 if i!=j else 1 for i in range(width) ] for j in range(height)] 15 16 print(L,U) 17 18 # 从U第一行开始，对于L来说是从第一列开始 19 for i in range(height): 20 # 依次计算LU的行、列 21 for j in range(i, width): 22 _sum = 0 23 for k in range(i): 24 _sum += L[i][k] * U[k][j] 25 print("U{}{}_sum += L[{}][{}] * U[{}][{}] = {} * {} = {}".format(i,j,i,k,k,j,L[i][k],U[k][j],L[i][k] * U[k][j])) 26 U[i][j] = mat[i][j] - _sum 27 # L的i.j是反着的 28 for j in range(i, width): 29 _sum = 0 30 if j > i: 31 for k in range(i): 32 _sum += L[j][k] * U[k][i] 33 print("L{}{}_sum += L[{}][{}] * U[{}][{}] = {} * {} = {}".format(j,i,j,k,k,i,L[j][k],U[k][i],L[j][k] * U[k][i])) 34 L[j][i] = (mat[j][i] - _sum) / U[i][i] 35 print("L = {}".format(np.array(L))) 36 print("U = {}".format(np.array(U))) 37 38 print("分解完毕，开始计算Ly=b") 39 # 分解完毕 解Ly =b 40 y = [[0] for i in range(height)] 41 # print(y) 42 for i in range(height): 43 _sum = 0 44 for k in range(i): 45 _sum += L[i][k] * y[k][0] 46 y[i][0] = (b[i][0] - _sum) # 这里不用除以系数因为对角线为1 47 print(y) 48 print("计算完毕，开始计算Ux=y") 49 x = [[0] for i in range(height)] 50 for i in range(height-1,-1,-1): 51 _sum = 0 52 for k in range(height-1,i,-1): 53 _sum += U[i][k] * x[k][0] 54 x[i][0] = (y[i][0] - _sum) / U[i][i] 55 print(x) 56 print("最终结果x: {}".format(np.array(x)))
over。