N元連立1次方程式の解法プログラム組んでみた

せっかく数式が使えるようになったしちょっと数学系書いてみる
N元連立1次方程式にはいろいろ解き方あるけど、アルゴリズムが一番楽なクラメルの公式

クラメルの公式とは

$A$を正方行列としたとき $Ax=b$で与えられる連立方程式の$i$番目の解は $$x_i=\frac{det(A_i)}{det(A)}$$

計算の流れ

公式だけでは説明不足で解らいので、例を使って解説

$$\left\{ \begin{array} _2x_1 + 4x_2 = 10 & \\ x_1 + 3x_2 = 6 & \end{array} \right. \tag{1}$$ 式1をクラメルの式を使って解く
行列で$Ax=b$の形に直すと
$$\left( \begin{array}{ccc} 2 & 4 \\ 1 & 3 \end{array} \right) \left( \begin{array}{ccc} x_1 \\ x_2 \end{array} \right) = \left( \begin{array}{ccc} 10 \\ 6 \end{array} \right)$$ よって、 $$A = \left( \begin{array}{ccc} 2 & 4 \\ 1 & 3 \end{array} \right) 、 b = \left( \begin{array}{ccc} 10\\ 6 \end{array} \right)$$

• クラメルの公式の分母
分母は行列$A$の行列式
$$det(A) = 2\times3-4\times1=2$$
• クラメルの公式の分子
$x_i$の場合の行列$A_i$は行列$A$の$i$番目の列と$b$を入れ替えた行列

1. $x_1$の場合 :
$A$の1列目を$b$と入れ替える $$A_1 = \left( \begin{array}{ccc} 10 & 4 \\ 6 & 3 \end{array} \right)$$ 計算すると $$det(A_1) = 10\times3-4\times6=6$$
2. $x_2$の場合 :
$A$の2列目を$b$と入れ替える $$A_2 = \left( \begin{array}{ccc} 2 & 10 \\ 1 & 6 \end{array} \right)$$ 計算すると $$det(A_2) = 2\times6-1\times10=2$$
• 解を求める
$$\left\{ \begin{array} _x_1 = \frac{6}{2}=3 \\ x_2 = \frac{2}{2}=1 \end{array} \right.$$

プログラム

これをN次元でも解けるようにするプログラムを組む
上のアルゴリズムから、

1. 行列式を求めるモジュール
2. 行列の列を入れ替えるモジュール

を作ればいい
数値計算プログラムはCとかFortranで書いた方がいいんだけど、環境がなかったので。。javaを使用
できるだけ、プリミティブ型で書いて、別言語に書き換えやすくしとく

• 行列式を求めるモジュール
余因子展開使って再帰関数で作ってく。N×N行列の余因子展開の公式はコレ
$$\begin{eqnarray} det(A)=\sum^n_{j=0}(-1)^j \cdot a_{0j} \cdot det(A_j) \end{eqnarray}$$ プログラムにするとこんな感じ。


/** 行列式の計算.
 * @param A 行列
 * @param n 行列の長さ
 * @return 計算結果
 */
static double dat(double[][] matrix, int n) {
    if(n == 1) return matrix[0][0];
    if(n == 2) return matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1];
    double sum = 0;
    for(int j = 0; j < n; j++) {
        int sign = ((j % 2) == 0) ? 1 : -1;
        double[][] laplace_marix = createLaplaceMarix(matrix,n,0,j);
        sum += sign * matrix[j][0] * dat(laplace_marix, n-1);
    }
    return sum;
}
  
/** 余因子展開で使う行列を作る.
 * @param A 元の行列
 * @param n 元の行列の長さ
 * @param line 消す行
 * @param col  消す列
 * @return
 */
static double[][] createLaplaceMarix(double[][] matrix, int n, int line, int col) {
    double laplace_marix[][] = new double[n-1][n-1]; // 余因子展開で使う行列
    for(int i = 0, y = 0; i < n; i++) {
        if(i == col)continue;
        for(int j = 0, x = 0; j < n; j++){
            if(j == line)continue;
            laplace_marix[y][x] = matrix[i][j];
            x++;
        }
        y++;
    }
    return laplace_marix;
}

• 行列の列を入れ替えるモジュール
まぁこれは簡単か。

/** 行列の列を入れ替える.
 * @param matrix 元の行列
 * @param n 列の長さ
 * @param col1 入れ替える列
 * @param col2 入れ替わる列 
 * @return 変換後の行列
 */
static double[][] changeColumn(double[][] matrix, int n, int col1, int col2) {
    double[][] temp_matrix = new double[n][n+1];
    // 配列コピー
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n + 1; j++) {
            temp_matrix[i][j] = matrix[i][j];
        }
    }
    // 列の入れ替え
    for(int i = 0; i < n; i++) {
        double temp = temp_matrix[i][col1];
        temp_matrix[i][col1] = temp_matrix[i][col2];
        temp_matrix[i][col2] = temp;
    }
    return temp_matrix;
}

• メイン文
モジュールができたので、組み立てる


public static void main(String[] args) {
    int n = 2; // 連立方程式の元
    double matrix[][]={ // 連立方程式を行列化
            {2, 4, 10}, // 2x + 4y = 10  - 式１
            {1, 3,  6}, //  x + 3y =  6  - 式２
    };
    double denominator = dat(matrix, n); // dat(A) 分母
    double[] result = new double[n]; // 答えの配列
    for(int i = 0; i < n; i++) {
        double[][] molecules_matrix = changeColumn(matrix, n, i, n);// 分子の行列
        result[i] = dat(molecules_matrix, n) / denominator;
    }
    // 解を表示
    for(int i = 0; i < n; i++) {
        System.out.println("result["+ i +"] : " + result[i]);
    }
}

実行してみる

4元連立方程式で試してみる。 $$\left\{ \begin{array} _2x_1 + x_2+3x_3+4x_4 = 2 & \\ 3x_1 + 2x_2+5x_3+2x_4 = 12 & \\ 3x_1 + 4x_2+x_3-x_4 = 4 & \\ -x_1 + -3x_2+x_3+3x_4 = -1 & \\ \end{array} \right. \tag{2}$$ 配列にすると、こう書けばOK

  double matrix[][]={ // 連立方程式を行列化
       // x1   x2   x3   x4    b ※bは移行してるわけではないので、右辺のbを書けばOK
       { 2.0, 1.0, 3.0, 4.0, 2.0}, //式1: 2x +  y + 3z + 4w =  2
       { 3.0, 2.0, 5.0, 2.0, 12.0},//式2: 3x + 2y + 5z + 2w = 12
       { 3.0, 4.0, 1.0,-1.0, 4.0}, //式3: 3x + 4y +  z - 1w =  4
       {-1.0,-3.0, 1.0, 3.0,-1.0}  //式4:-1x - 3y +  z + 3w = -1
  };

回答がこんな感じ

result[0] : 1.0
result[1] : -1.0
result[2] : 3.0
result[3] : -2.0

できたー！