python矩阵的分解及其应用（LU分解）

目录

1、代码

2、结果

---

1、代码

author:hewang
qq:207962168

import numpy as np 
from scipy import linalg
def getMatrix(epsion):
    A=np.array([[epsion,1],[1,0]],dtype=np.float64)
    return A 
def getRHS():
    b=np.array([[1],[3.0]],dtype=np.float64)
    return b
def LUsolution(A,b,epsion):
    L=np.array([[1,0],[1/epsion,1]],dtype=np.float64)
    U=np.array([[epsion,1],[0,-1/epsion]],dtype=np.float64)
    y=np.linalg.solve(L,b)
    x=np.linalg.solve(U,y)
    return L,U,x
def LUPivotsolution(A,b,epsion):
    P=np.array([[0,1],[1,0]],dtype=np.float64)
    b=np.matmul(P,b)
    L=np.array([[1,0],[epsion,1]],dtype=np.float64)
    U=np.array([[1,0],[0,1]],dtype=np.float64)
    y=np.linalg.solve(L,b)
    x=np.linalg.solve(U,y)
    return L,U,x
def printDetails(L,U,x):
    print("L")
    print(L)
    print("U")
    print(U)
    print("solution of Ax=b is:")
    print(x)


def main():
   
    #epsion =1.0*2**(-52)
    #epsion =1.01*2**(-52)
    epsion =1.1*2**(-57)
    A=getMatrix(epsion)
    b=getRHS()
    print("epsilon:{}".format("%lf",epsion))
    print("=====WITHOUT pivot======")
    L,U,x=LUsolution(A,b,epsion)
    printDetails(L,U,x)
    print("=====WITH PIVOT=====")
    L,U,x=LUPivotsolution(A,b,epsion)
    printDetails(L,U,x)

if __name__ == '__main__':
    main()

2、结果

epsilon:%lf
=====WITHOUT pivot======
L
[[1.00000000e+00 0.00000000e+00]
 [1.31013807e+17 1.00000000e+00]]
U
[[ 7.63278329e-18  1.00000000e+00]
 [ 0.00000000e+00 -1.31013807e+17]]
solution of Ax=b is:
[[0.]
 [1.]]
=====WITH PIVOT=====
L
[[1.00000000e+00 0.00000000e+00]
 [7.63278329e-18 1.00000000e+00]]
U
[[1. 0.]
 [0. 1.]]
solution of Ax=b is:
[[3.]
 [1.]]

Process finished with exit code 0