[数学] LU分解法をプログラムで書いてみた ref: http://qiita.com/edo_m18/items/1d67532bed4a083cddb3

file0.txt

\begin{bmatrix}
l_{00} & 0 & 0 \\
l_{10} & l_{11} & 0 \\
l_{20} & l_{21} & l_{22}
\end{bmatrix}

file1.txt

\begin{bmatrix}
1 & u_{01} & u_{02} \\
0 & 1 & u_{12} \\
0 & 0 & 1 
\end{bmatrix}

file10.js

(function () {

    'use strict';

    // とりあえず成分数は明示的に書いておく
    let N = 4;

    // 分解したい行列A
    let matA = [
        [8, 16, 24,  32],
        [2,  7, 12,  17],
        [6, 17, 32,  59],
        [7, 22, 46, 105]
    ]; 

    // 分解した下三角行列L用バッファ
    let matL = [
        [0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0],
    ];

    // 分解した上三角行列U用バッファ
    let matU = [
        [1, 0, 0, 0],
        [0, 1, 0, 0],
        [0, 0, 1, 0],
        [0, 0, 0, 1],
    ];

    // 分解後に利用するバッファ用行列
    let lu = [
        [0, 0, 0, 0], 
        [0, 0, 0, 0], 
        [0, 0, 0, 0], 
        [0, 0, 0, 0], 
    ];

    let b = [
        160, 70, 198, 291,
    ];

    for (let i = 0; i < N; i++) {

        let n = N - i - 1;

        // l0_0成分をコピー
        let l0 = matL[i][i] = matA[0][0];

        // l1成分をコピー
        let l1 = [];
        for (let j = 0; j < n; j++) {
            matL[j + i + 1][i] = l1[j] = matA[j + 1][0];
        }

        // →u1^T成分をコピー
        let u1 = [];
        for (let j = 0; j < n; j++) {
            matU[i][j + i + 1] = u1[j] = matA[0][j + 1] / l0;
        }

        // luを求める
        for (let j = 0; j < n; j++) {
            for (let k = 0; k < n; k++) {
                lu[j][k] = l1[j] * u1[k];
            }
        }

        // A1を求める（n-1行列）
        let A1 = [];
        for (let j = 0; j < n; j++) {
            A1[j] = [];
            for (let k = 0; k < n; k++) {
                A1[j][k] = matA[j + 1][k + 1] - lu[j][k];
            }
        }

        // A1を改めてmatAとして再帰的に解く
        matA = A1;
    }

    // 求めたLU行列を使って連立方程式を解く
    let y = new Array(N);
    for (let i = 0; i < N; i++) {
        let sum = 0;
        for (let k = 0; k <= i - 1; k++) {
            sum += matL[i][k] * y[k];
        }
        y[i] = (b[i] - sum) / matL[i][i];
    }

    let x = new Array(N);
    for (let i = N - 1; i >= 0; i--) {
        let sum = 0;
        for (let k = i + 1; k <= N - 1; k++) {
            sum += matU[i][k] * x[k];
        }
        x[i] = y[i] - sum;
    }


    // 以上で求めた$x$の配列が、最終的に求めたい変数$x_0, x_1, x_2, x_3$それぞれの解となります。

}());

file11.txt

x = [4, 3, 2, 1]

file2.txt

\left\{
\begin{array}{ll}
8x_0 + 16x_1 + 24x_2 +  32x_3 = 160 \\
2x_0 +  7x_1 + 12x_2 +  17x_3 = 70 \\
6x_0 + 17x_1 + 32x_2 +  59x_3 = 198 \\
7x_0 + 22x_1 + 46x_2 + 105x_3 = 291
\end{array}
\right.

file3.txt

A = 
\begin{bmatrix}
8 & 16 & 24 &  32 \\
2 &  7 & 12 &  17 \\
6 & 17 & 32 &  59 \\
7 & 22 & 46 & 105
\end{bmatrix}

,\vec{x} = 
\begin{bmatrix}
x_0 \\
x_1 \\
x_2 \\
x_3
\end{bmatrix}

,\vec{b} = 
\begin{bmatrix}
160 \\
70 \\
198 \\
291
\end{bmatrix}

file4.txt

A\vec{x} = \vec{b}

file5.txt

LU\vec{x} = \vec{b}

file6.txt

\vec{y} = U\vec{x}

file7.txt

\left\{
\begin{array}{ll}
L\vec{y} = \vec{b} \\
U\vec{x} = \vec{y}
\end{array}
\right.

file8.txt

A = 
\begin{bmatrix}
8 & 16 & 24 &  32 \\
2 &  7 & 12 &  17 \\
6 & 17 & 32 &  59 \\
7 & 22 & 46 & 105
\end{bmatrix}

,L = 
\begin{bmatrix}
8 & 0 & 0 & 0 \\
2 & 3 & 0 & 0 \\
6 & 5 & 4 & 0 \\
7 & 8 & 9 & 8 \\
\end{bmatrix}

,U = 
\begin{bmatrix}
1 & 2 & 3 & 4 \\
0 & 1 & 2 & 3 \\
0 & 0 & 1 & 5 \\
0 & 0 & 0 & 1 \\
\end{bmatrix}

file9.txt

\left\{
\begin{array}{ll}
L\vec{y} = \vec{b} \\
U\vec{x} = \vec{y}
\end{array}
\right.