一、简介
在线性代数中, LU分解(LU Decomposition)是矩阵分解的一种,可以将一个矩阵分解为一个单位下三角矩阵和一个上三角矩阵的乘积(有时是它们和一个置换矩阵的乘积)。LU分解主要应用在数值分析中,用来解线性方程、求反矩阵或计算行列式。
详解
二、实现
# -*- coding: utf-8 -*-
"""
Created on Tue Dec 13 15:59:27 2016
LU PLU矩阵分解 LU求解方程
@author: Administrator
"""
from numpy import *
import matplotlib.pyplot as plt
def lu_decomp(A):
m,n = A.shape
if m != n:
return
temp_A = A.copy()
L = eye(n)
U = zeros((m,n))
for k in range(0,n-1):
for i in range(k+1,n):
if abs(temp_A[i,k]) > 1.0e-9:
lam = temp_A[i,k] / temp_A[k,k]
temp_A[i,k+1:n] = temp_A[i,k+1:n] - lam * temp_A[k,k+1:n]
temp_A[i,k] = 0
L[i,k] = lam
U[:n,:n] = temp_A[:n,:n]
return L,U
def plu_decomp(A):
m,n = A.shape
if m != n:
return
temp_A = A.copy()
L = eye(n)
U = zeros((m,n))
index = [0 for i in range(n-1)]
for k in range(0,n-1):
index[k] = argmax(abs(temp_A[k:n,k])) + k
if(index[k] != k):
for j in range(n):
temp_A[index[k],j],temp_A[k,j] = temp_A[k,j],temp_A[index[k],j]
for i in range(k+1,n):
if abs(temp_A[i,k]) > 1.0e-9:
if temp_A[k,k] == 0.:
return
lam = temp_A[i,k] / temp_A[k,k]
temp_A[i,k+1:n] = temp_A[i,k+1:n] - lam * temp_A[k,k+1:n]
temp_A[i,k] = 0
L[i,k] = lam
U[:n,:n] = temp_A[:n,:n]
return U
#列主元求解方程 ex=true
def solve_ex(A,b,ex=True):
m,n = A.shape
if m != n:
return
temp_A = A.copy()
x = zeros(n)
index = [0 for i in range(n-1)]
#消去
for k in range(0,n-1):
index[k] = argmax(abs(temp_A[k:n,k])) + k
if(index[k] != k and ex):
for j in range(n):
temp_A[index[k],j],temp_A[k,j] = temp_A[k,j],temp_A[index[k],j]
b[k],b[index[k]] = b[index[k]],b[k]
for i in range(k+1,n):
if abs(temp_A[i,k]) > 1.0e-9:
lam = temp_A[i,k] / temp_A[k,k]
temp_A[i,k+1:n] = temp_A[i,k+1:n] - lam * temp_A[k,k+1:n]
temp_A[i,k] = 0
b[i] = b[i] - lam * b[k]
#回代
x[n-1] = b[n-1] / temp_A[n-1,n-1]
for k in range(n-2,-1,-1):
temp = b[k]
for j in range(k+1,n):
temp = temp - temp_A[k,j] * x[j]
x[k] = temp / temp_A[k,k]
return x
if __name__ == '__main__':
A = array([[1,-6,-1],[1,2,2],[2,-2,1]]) #square matrix
b = array([-2,4,1],dtype=float)
A = array([[2,1,1,0],[4,3,3,1],[8,7,9,5],[6,7,9,8]]) #square matrix
A = matrix(A,dtype=float)
m,n = A.shape
L,U = lu_decomp(A)
e = eye(n,dtype=float)
A_ = []
for i in range(0,n):
A_.append(solve_ex(U,solve_ex(L,e[:,i]),ex=False))
print matrix(A_,dtype=float)