lu分解

鎖定

lu分解計算

將係數矩陣A轉變成等價兩個矩陣L和U的乘積 ，其中L和U分別是單位下三角矩陣和上三角矩陣。當A的所有順序主子式都不為0時，矩陣A可以唯一地分解為A=LU。其中L是下三角矩陣，U是上三角矩陣。

import java.util.Arrays;
/**
 * 矩陣的直接三角分解 ，調用示例：
 * 
 * DirectDecomposition dd = new DirectDecomposition(data);//data為一個二維double數組，代替一個矩陣
 * 
 * double[][] l = dd.getL();//獲取L
 * 
 * double[][] u = dd.getU();//獲取U
 * 
 * @author jungle
 */
public class DoolittleDecomposition {
    private double[][] data;
    private double[][] l;
    private double[][] u;
    private int n;
    /**
     * 創建一個n階的矩陣
     * 
     * @param n
     */
    public DoolittleDecomposition(double[][] data) {
        if (data == null || data.length == 0 || data.length != data[0].length) {
            throw new RuntimeException("不是一個方陣");
        }
        this.data = data;
        n = data.length;
        l = new double[n][n];
        u = new double[n][n];
        countLU();
    }
    protected void countLU() {
        for (int i = 0; i < n; i++) {// 第一步，計算L的第一列和U的第一行：U1i=A1i,Li1=Ai1/U11
            u[0][i] = data[0][i];
            l[i][0] = data[i][0] / u[0][0];
        }
        //計算U的第r行，L的第r列元素
        //uri=ari-Σ(k=1->r-1)lrkuki    (i=r,r+1,...,n)
        //lir=air-Σ(k=1->r-1)likukr    (i=r+1,r+2,...,n且r≠n)
        for (int r = 1; r < n; r++) {
            for (int i = r; i < n; i++) {
                u[r][i] = data[r][i] - sumLrkUki(r, i);
                if(i==r) l[r][r]=1;
                else if(r==n) l[n][n]=1;
                else l[i][r] = (data[i][r] - sumLikUkr(r, i)) / u[r][r];
            }
        }
    }
    /**
     * 求和：Lrk*Uki 對k求和：1<=k<=r-1
     * 
     * @param r
     * @param i
     * @return
     */
    private double sumLrkUki(int r, int i) {
        double re = 0.0;
        for (int k = 0; k < r; k++) {
            re += l[r][k] * u[k][i];
        }
        return re;
    }
    private double sumLikUkr(int r, int i) {
        double re = 0.0;
        for (int k = 0; k < r; k++) {
            re += l[i][k] * u[k][r];
        }
        return re;
    }
    public double[][] getData() {
        return data;
    }
    public double[][] getL() {
        return l;
    }
    public double[][] getU() {
        return u;
    }
    public static void main(String[] args) {
        double[][] data= {
                {1,2,6},
                {2,5,15},
                {6,15,46},
        };
        DoolittleDecomposition dd = new DoolittleDecomposition(data);
        double[][] l = dd.getL();
        double[][] u = dd.getU();
        int n = l.length;
        System.out.println("L陣：");
        for (int i = 0; i < n; i++) {
            System.out.println(Arrays.toString(u[i]));
        }
        System.out.println("---------------------");
        System.out.println("U陣：");
        for (int i = 0; i < n; i++) {
            System.out.println(Arrays.toString(l[i]));
        }
    }
}

%LU分解法求解Ax=b，假定A矩陣可進行LU分解以及對角線元素均不為0
function [x] = Dool(A,b)
n=length(A);A(2:n,1)=A(2:n,1)/A(1,1);
for t=2:n-1  %進行LU分解
    A(t,t:n)=A(t,t:n)-A(t,1:t-1)*A(1:t-1,t:n);
    A(t+1:n,t)=(A(t+1:n,t)-A(t+1:n,1:t-1)*A(1:t-1,t))/A(t,t);
end
A(n,n)=A(n,n)-A(n,1:n-1)*A(1:n-1,n);
L=tril(A,-1)+eye(n);U=triu(A); %矩陣A分解出L和U
for t=2:n   %解Lx=b
    b(t)=b(t)-L(t,1:t-1)*b(1:t-1);
end
b(n)=b(n)/U(n,n);   %解Ux=b
for t=1:n-1;
    k=n-t;b(k)=(b(k)-U(k,k+1:n)*b(k+1:n))/U(k,k);
end
x=b;  %方程Ax=b的解即為x

參考資料

• 1.    LU分解   ．csdn[引用日期2012-12-01]