用LU分解法来解矩阵方程

最新推荐文章于 2022-05-14 20:30:45 发布

LazzyBoi懒惰男孩

最新推荐文章于 2022-05-14 20:30:45 发布

阅读量2.1k

点赞数

分类专栏：数值计算与优化

本文链接：https://blog.csdn.net/m0_37717751/article/details/81050141

版权

数值计算与优化专栏收录该内容

4 篇文章 0 订阅

订阅专栏

运用的是克劳特（Crout）分解，也就是U是单位上三角形矩阵

import numpy as np
#np.set_printoptions(threshold=np.inf)

def LU(A, n):
    L = [[0 for x in range(n)] for y in range(n)]
    U = [[0 for x in range(n)] for y in range(n)]  
    for i in range(n, 0, -1):
	    for tmp1 in range(n - i, n):
		    sum = 0
		    for tmp2 in range(n - i):
			    sum += L[tmp1][tmp2] * U[tmp2][n - i]
		    L[tmp1][n - i] = A[tmp1][n - i] - sum
	    for tmp1 in range(n - i + 1, n):
	    	sum = 0
	    	for tmp2 in range(n - i):
	    		sum += L[n - i][tmp2] * U[tmp2][tmp1]
	    	U[n - i][tmp1] = (A[n - i][tmp1] - sum) / float(L[n - i][n - i])
    for i in range(n):
    	U[i][i] = 1
    #print("L:")
    #printMat(L, 3)	
    #print("U")
    #printMat(U, 3)
    return L, U
    
def LUSolve(A, b, n):
	L = [[0 for x in range(n)] for y in range(n)]
	U = [[0 for x in range(n)] for y in range(n)]
	L, U = LU(A, n)
	Y = [0 for x in range(n)]
	X = [0 for x in range(n)]
	Y[0] = b[0][0] / L[0][0]
	for i in range(1, n):
		sum = 0
		for tmp in range(0, i):
			sum += L[i][tmp] * Y[tmp]
		Y[i] = (b[0][i] - sum) / float(L[i][i])
	X[n - 1] = Y[n - 1]
	for i in range(n - 2, -1, -1):
		sum = 0
		for tmp in range(i + 1, n):
			sum += U[i][tmp] * X[tmp]
		X[i] = Y[i] - sum
	return X

if __name__ == "__main__":
    M = np.random.random(size = (100, 100))
    #print("M:")
    #print(M)
    A = M
    for i in range(0, 100):
        A[i][i] += 1  # A = M + I
    #print("A:")
    #print(A)
    X = np.random.random(size = (100, 1))
    for i in range(0, 100):
        X[i][0] = i + 1
    #print("X:")
    #print(X)
    B = np.dot(A, X)
    b = np.transpose(B)
    #(b)
    result = LUSolve(A, b, 100)
    print(result)