计算方法 | 方程组求解的三角分解

有两种方法，一个是通过高斯消去法的方式去求一个下三角矩阵，还有一个是直接分解。

显然直接分解好啊233。

然后这是个啥杜尔利特分解公式。

是由矩阵乘法的原理弄出来的。

书上P99页俩公式(5.34 5.35)。

轮流求u和l

三角分解部分的代码：

import numpy as np

mat = [[2,2,3],[4,7,7],[-2,4,5]]

height = len(mat)
width = len(mat[0])

# 1 直接分解矩阵
# 初始化LU
U = [[0 for i in range(width)] for j in range(height)]
L = [[0 if i!=j else 1 for i in range(width) ] for j in range(height)]

print(L,U)

# 从U第一行开始，对于L来说是从第一列开始
for i in range(height):
    # 依次计算LU的行、列
    for j in range(i, width):
        _sum = 0
        for k in range(i):
            _sum += L[i][k] * U[k][j]
            print("U{}{}_sum += L[{}][{}] * U[{}][{}] = {} * {} = {}".format(i,j,i,k,k,j,L[i][k],U[k][j],L[i][k] * U[k][j]))
        U[i][j] = mat[i][j] - _sum
    # L的i.j是反着的
    for j in range(i, width):
        _sum = 0
        if j > i:
            for k in range(i):
                _sum += L[j][k] * U[k][i]
                print("L{}{}_sum += L[{}][{}] * U[{}][{}] = {} * {} = {}".format(j,i,j,k,k,i,L[j][k],U[k][i],L[j][k] * U[k][i]))
            L[j][i] = (mat[j][i] - _sum) / U[i][i]
    print("L = 
{}".format(np.array(L)))
    print("U = 
{}".format(np.array(U)))
    
# 分解完毕

解Ly =b, 回带求解，再解Ux=y, 直接丢完整代码了：

import numpy as np
# mat = eval(input("请输入系数矩阵："))
# 测试数据[[2,2,3],[4,7,7],[-2,4,5]]
mat = [[2,2,3],[4,7,7],[-2,4,5]]
# b = eval(input("请输入结果矩阵:"))
# 测试数据[[3],[1],[-7]]
b = [[3],[1],[-7]]
height = len(mat)
width = len(mat[0])

# 1 直接分解矩阵
# 初始化LU
U = [[0 for i in range(width)] for j in range(height)]
L = [[0 if i!=j else 1 for i in range(width) ] for j in range(height)]

print(L,U)

# 从U第一行开始，对于L来说是从第一列开始
for i in range(height):
    # 依次计算LU的行、列
    for j in range(i, width):
        _sum = 0
        for k in range(i):
            _sum += L[i][k] * U[k][j]
            print("U{}{}_sum += L[{}][{}] * U[{}][{}] = {} * {} = {}".format(i,j,i,k,k,j,L[i][k],U[k][j],L[i][k] * U[k][j]))
        U[i][j] = mat[i][j] - _sum
    # L的i.j是反着的
    for j in range(i, width):
        _sum = 0
        if j > i:
            for k in range(i):
                _sum += L[j][k] * U[k][i]
                print("L{}{}_sum += L[{}][{}] * U[{}][{}] = {} * {} = {}".format(j,i,j,k,k,i,L[j][k],U[k][i],L[j][k] * U[k][i]))
            L[j][i] = (mat[j][i] - _sum) / U[i][i]
    print("L = 
{}".format(np.array(L)))
    print("U = 
{}".format(np.array(U)))
    
print("分解完毕，开始计算Ly=b")
# 分解完毕 解Ly =b
y = [[0] for i in range(height)]
# print(y)
for i in range(height):
    _sum = 0
    for k in range(i):
        _sum += L[i][k] * y[k][0]
    y[i][0] = (b[i][0] - _sum)          # 这里不用除以系数因为对角线为1
    print(y)
print("计算完毕，开始计算Ux=y")
x = [[0] for i in range(height)]
for i in range(height-1,-1,-1):
    _sum = 0
    for k in range(height-1,i,-1):
        _sum += U[i][k] * x[k][0]
    x[i][0] = (y[i][0] - _sum) / U[i][i]
    print(x)
print("最终结果x: 
{}".format(np.array(x)))

over。