Doolittle分解法（LU分解法）的Python实现

在解一般的非奇异矩阵线性方程组的时候，或者在迭代改善算法中，需要使用LU分解法。
对于一个一般的非奇异矩阵A=(a11, a12,…,a1n,a21,…ann)，可分解为一个下三角矩阵L和一个上三角矩阵U。其中L的主对角线元素都是1.
希望得到一个M，最后在需要的时候将M拆分为L和U

M=[[u11 u12 u13 ...... u1n
    l21 u22 u23 ...... u2n
    l31 l32 u33 ...... u3n
    ......                 ]]
L=[[1   0   0......
	l21 1   0......
	l31 l32 1......
	......         ]]
U=[[u11 u12 u13......
	0   u22 u23......
	0   0   u33......
	......         ]]

doolittle分解法分两步进行
首先，根据已知A直接做出一个基矩阵M，其中，矩阵的第一行就是A的第一行，矩阵的第一列为A的第一列除以矩阵M的第一列的第一个元。即
u1j=a1j j=1,2,…,n
li1=ai1/u11, i=1,2,…,n
接着，从第二行开始，从第二列到最后一列一次求解M的元素，列数大于等于行数则属于U，否则属于L，有公式
ukj=akj-sum(lkmumj)(m∈(1, k-1)), j=1,2,…,n
lik=(aik-sum(limumk)(m∈(1, k-1)))/ukk, i=1,2,…,n
可以发现lik与ukj公式的区别只有lik多除以了一个ukk即这一列的对角线元
下面给出一个例子

A=[[1  2  -3
    2  -1  3
    3  -2  2]]

将A进行doolittle分解

代码如下

# -*- coding: utf-8 -*-
"""
Created on Tue Feb 12 11:27:19 2019

@author: 鹰皇
"""

#Doolittle分解法
import numpy as np
A=[[1.0,2.0,-3.0],[2.0,-1.0,3.0],[3.0,-2.0,2.0]]

def get_base(A):
    base=list(np.zeros((len(A),len(A))))
    D=[]
    for i in base:
        D.append(list(i))
    return D

def U_initiate(A):
    u=get_base(A)
    for i in range(0,len(A)):
        u[0][i]=A[0][i]
    return u

def L_initiate(A, u):
    l=get_base(A)
    for i in range(0,len(A)):
        l[i][i]=1.0
        if i==0:
            pass
        else:
            l[i][0]=A[i][0]/u[0][0]
    return l

def get_u(A,L,U,k,j):
    sum1=0.0
    for m in range(0,k):
        sum1=sum1+L[k][m]*U[m][j]
    U[k][j]=A[k][j]-sum1
    return U

def get_l(A,L,U,i,k):
    sum1=0.0
    for m in range(0,k):
        sum1=sum1+L[i][m]*U[m][k]
    L[i][k]=(A[i][k]-sum1)/U[k][k]
    return L

def main(A):
    U=U_initiate(A)
    L=L_initiate(A,U)
    n=1
    while n<len(A):
        t=1
        while t<len(A):
            if t>=n:
                U=get_u(A,L,U,n,t)
            else:
                L=get_l(A,L,U,n,t)
            t=t+1
            print t
        n=n+1
    return U,L

print main(A)

算得L，U，解答完毕