掃き出し法による連立方程式の解法プログラム

以前、クラメルの公式を使った連立方程式の解法を紹介したが、13元くらいからめちゃくちゃ遅くなることがわかった
そこで今回は一番よくつかわれてる掃き出し法による解法で作ってみる

・N元連立1次方程式の解法プログラム(マルチスレッド版)
・N元連立1次方程式の解法プログラム組んでみた

掃き出し法とは

$N$元一次方程式$Ax=b$の解法の一つ
係数行列$A$に定数項の列ベクトル$b$を加えた行列 $[A|b]$ に「行列の基本変形」を行い $[E|u]$ に変形したとき、方程式の解$x_i=u_i$となる

行列の基本変形

①	行の入れ替え	:	n行目とm行目を入れ替える
②	行の定数倍	:	n行目を定数倍する　( 定数≠0 )
③	行の加減算	:	n行目を定数倍した値をm行目に加減算する

上記「行列の基本変形」を使うことで、行列を変形させていき、解を導いていく方法が掃き出し法だ
言葉だけじゃよくわからないので例で計算してみる

プログラム

上記の処理をプログラムで自動化させる
まず、行列の基本変形の３つのモジュールを作る

①行の入れ替え


/**
 * 行列の行(line1とline2)を入れ替える
 * @param line1 入れ替え対象の行1
 * @param line2 入れ替え対象の行2
 * @param matrix_n 行列の行数
 * @param matrix 拡大係数行列
 */
void ExchangeLine(int line1, int line2, int matrix_n, double** matrix) {
    double* tmp_line = *(matrix + line2);
    *(matrix + line2) = *(matrix + line1) ;
    *(matrix + line1) = tmp_line;
}

②行の定数倍


/**
 * line行目をc倍にする
 * @param line1 定数倍する行
 * @param c 定数倍にする数
 * @param matrix_n 行列の行数
 * @param matrix 拡大係数行列
 */
void TimesLine(int line, double c, int matrix_n, double** matrix) {
    double* tmp_line = *(matrix + line);
    for(int i = 0; i < matrix_n + 1; i++) {
        *(tmp_line + i) = *(tmp_line + i) * c;
    }
}

③行の加減算


/**
 * line2 行目をc倍したものをline1 行目に加算
 * @param line1  足される行
 * @param line2  足す行
 * @param c      倍数(引き算がしたい場合はマイナス倍にする)
 * @param matrix_n 行列の行数
 * @param matrix 拡大係数行列
 */
void CalcLine(int line1, int line2, double c, int matrix_n, double** matrix) {
    double* tmp1 = *(matrix + line1); // 足される行
    double* tmp2 = *(matrix + line2); // 足す行

    // 一要素ずつ足していく
    for(int i = 0; i < matrix_n + 1; i++) {
        *(tmp1 + i) = *(tmp1 + i) + c * *(tmp2 + i);
    }
}

あとはこれらの組み合わせで実現できる
フローはこんな感じ

一列目から順に単位行列に変形していく形になる

#include<iostream>

#define N 3  // 元数 (3元一次連立方程式だとN = 3)

// 行列の基本変形の３つ
void ExchangeLine(int line1, int line2, int n, double** matrix);
void TimesLine(int line, double c, int n, double** matrix);
void CalcLine(int line1, int line2, double c, int n, double** matrix);

// 「①行の入れ替え」の応用。指定行以下の行をすべて1行ずつスワイプさせる
void RotateMatrix(int line, int matrix_n, double** matrix);

int main() {

    // 拡大係数行列の設定
    double matrix[][N + 1] = {
      {0, 1, 2, 8},
      {0, 2, 1,7},
      {1, 2, 3, 14}
    };

    // 演算しやすいようにポインタに変換
    double* p_matrix[N];
    for(int i = 0;i < N;i++) p_matrix[i] = matrix[i];


    // 行列の変形
    for(int line = 0; line < N; line++) {

        // ①行の入れ替え
        int loop_cout = 0, end_num = N - line;
        while( *(p_matrix[line] + line) == 0) {
            // 対象要素が0のとき、演算できないため、行をスワイプさせる
            RotateMatrix(line, N, p_matrix);
            if( loop_cout++ > end_num ) break;
        }

        // ②行の定数倍 (対象の行を正規化する)
        // 対象の要素が「0」の場合は演算不可能なため終了させる
        double base = *( *(p_matrix + line) + line);
        if(base == 0) break;
        TimesLine(line, 1 / base, N, p_matrix);

        // ③行の加減算
        // 対象行の列は「１」(②で実行済)それ以外の対象列は「０」になるように加減算を行う
        for(int i = 0; i < N; i++) {
            if(i == line) continue;
            double c = *(p_matrix[i] + line);
            CalcLine(i, line, -1 * c, N, p_matrix);
        }
    }


    // 解の表示
    for(int i = 0; i < N; i++) {
        std::cout << "x" << (i + 1) << " = " << *(p_matrix[i] + N) << std::endl;
    }

    return 0;
}


/**
 * line行目以下の行をすべて１つずつスワイプさせる
 * @param line この行以下の行をすべてスワイプさせる
 * @param matrix_n 行列の元数
 * @param matrix スワイプ対象行列
 */
void RotateMatrix(int line, int matrix_n, double** matrix) {
    int line_num = matrix_n - line;
    for(int i = 0; i < line_num - 1; i++) {
        int line1 = i + line;
        int line2 = (i + 1) % line_num + line;
        ExchangeLine(line1, line2, matrix_n, matrix);
    }
}

// 上記してるので、処理内容省略
void ExchangeLine(int line1, int line2, int matrix_n, double** matrix) {上記により省略}
void TimesLine(int line, double c, int matrix_n, double** matrix) {上記により省略}
void CalcLine(int line1, int line2, double c, int matrix_n, double** matrix) {上記により省略}

実行すると。。。

x1 = 1
x2 = 2
x3 = 3

完成！
他の連立方程式を実行する際は、ハイライトがついている箇所(下記の内容)を変更する

・拡大係数行列
・defineで定義している元数N

感想 

クラメルの公式からすると信じられないくらい早いｗｗ
しかも仕組みが簡単で作りやすい
実は連立方程式の解には、「一意の解」、「任意の解」、「解なし」の3種類あり、今回のは「一意の解」に限定している
これらの種類の違いを調べるには行列のランクを調べる必要がある。ランクも掃き出し法の応用に過ぎないので、本記事のロジックが理解できれば余裕
時間があればランクを加味してすべての連立方程式に対応したプログラムもつくってみよかな