ガウス・ジョルダン消去法は、連立一次方程式の解を求めるための数値計算アルゴリズムです。
基本概念
- ガウス消去法：連立一次方程式を解くための基本的な手法
- ガウス・ジョルダン消去法：ガウス消去法を改良した方法で、コードが簡潔になり、簡約化された行階段形行列を一度に求められる
拡大行列の構築
ガウス消去法では、通常n個のn元連立一次方程式をn行×(n+1)列の「拡大行列」として表現します。具体的な手順は以下の通りです：
- 方程式を整理し、左辺がすべて「係数×変数」の和の形になるように移項します。特定の変数が存在しない場合は、係数0を補います。定数項はすべて右辺に移動します。
- これにより、各方程式の左辺にはn個の一次項（係数0も含む）が、右辺には1つの定数項が存在します。
- 行列のi行j列は、i番目の方程式中のx_jの係数を表します。特に、i行(n+1)列は方程式右辺の定数項を表します。
（実際の問題では、拡大行列が直接与えられることが多いので、手動での処理は不要な場合が多いです。）
例： [\begin{cases} 3 x_1 + 4 x_3 - 2 = 0\ x_1 + 5 x_3 + 3 = 0\ 4x_1 + 2x_2 + 3x_3 - 5 = 0 \end{cases} ]ステップ1：方程式の整理（移項、0係数の補完）
[\begin{cases} 3 x_1 &+& 0 x_2 &+& 4 x_3 &=& 2\ 1 x_1 &+& 0 x_2 &+& 5 x_3 &=& -3\ 4 x_1 &+& 2 x_2 &+& 3 x_3 &=& 5 \end{cases} ]ステップ2：拡大行列の構築（上記の係数から行列を作成）
[\left[ \begin{array}{ccc|c} 3 & 0 & 4 & 2\ 1 & 0 & 5 & -3\ 4 & 2 & 3 & 5 \end{array} \right] ]縦線の左側が係数、右側が定数項です。以降のガウス消去法はこの行列に対して操作を行います。
消去過程
この行列に対して、以下の操作を行うことができます：
- 行の入れ替え
- 特定の行のすべての要素を定数倍または定数で除算
- 特定の行に別の行を加えるまたは引く
これらの操作を任意に行った後の行列は元の行列と等価であり、これら3つの操作のみでガウス消去法の全過程を完了できます。
ガウス・ジョルダン消去法の特徴は、有効な方程式を「単位行列」（縦線右側の拡大部分を除き、左側の主対角線がすべて1）に一度に消去できる点です。
具体的には、行列のx行y列をM(x, y)とした場合、以下の手順で進めます：
- 列ごとに順次処理し、i列目について、M(p, i)が0ではないp ≥ iを探します。
- p行とi行を入れ替えます（後続のコード記述を大幅に簡略化できます）。これにより、この行を使って他のすべての行のi列目を0にできます。
- まずi行のすべての要素をM(i, i)で除算し、M(i, i)をすべて1にします。これにより行列は元の行列と等価のままです。
- その他のすべての行について、各j ≠ iに対し、この行全体からM(j, i)倍のi行を引きます。これによりM(j, i)が消去されます。 （厳密にはM(j, i) ÷ M(i, i)倍すべきですが、M(i, i) = 1なので約分されます）
- これにより、i列目はi行を除いてすべて0になります。各i ∈ [1, n]に対してこの操作を繰り返すと、最終的に「簡約化された行階段形行列」または「正規化された行階段形行列」が得られます。つまり、縦線左側が単位行列となり、その拡大部分が直接方程式の解となります。
解の判定
上記は方程式が解を持つ場合の通常の処理過程ですが、無解をどう判定するのでしょうか。
ここでいう無解は「一意な解を持たない」、つまり「解がない」と「無数の解を持つ」の2つのケースを指します。
これには線形代数の知識（「線形従属」の知識）が必要です。簡単に言えば、特定の列で有効なpが見つからない場合、つまりi行以降がすべて0の場合、ここには他の式から導き出せる式が含まれていることを意味します。これは実質的に有用な式ではなく、n個の未知数に対して(n-1)個以下の方程式しかない状態です。この場合、方程式は一意な解を持たなくなります。
この場合は直接「一意な解を持たない」と判定できます。しかし、問題によっては「解がない」か「無数の解を持つ」かを区別する必要がある場合があり、ロジックを変更する必要があります。
具体的には、「完全な消去」ではなく「可能な限りの消去」を目指します。特定の列で有効なpが見つからない場合は、その列をスキップして次の列でマッチングを続けます。
これにより、最終的に「完全なマッチング」ができない場合、特定の列でp行のマッチングが見つからず、行列の下部にいくつかの行が残ります。これらの行の係数はすべて0であることが保証されます（この行の各列について、この列がマッチング成功していたら0に消去され、マッチング成功しなかった場合は有効なpが見つからなかったのでここも0です）。
これらの式を個別に抽出します：
[\begin{cases} 0 x_1 &+& 0 x_2 &+& \cdots &+& 0 x_n &=& v_k\ 0 x_1 &+& 0 x_2 &+& \cdots &+& 0 x_n &=& v_{k + 1}\ \vdots && \vdots && \ddots && \vdots && \vdots\ 0 x_1 &+& 0 x_2 &+& \cdots &+& 0 x_n &=& v_n\ \end{cases} ]上記のすべてのvが0である場合、前のxはどのような値でも取り得ます—無数の解があります。vの中に0以外の数値が現れる場合、xをどのように取っても0以外の数字を組み合わせることはできません—解がありません。
コード例
コードでは、実数演算が大量に行われるため必ず精度問題が発生するため、0の判定はv = 0ではなく|v| < EPSを使用する必要があります。ここでEPSは極小の数値（例：10^-6）を表します。
以下に「洛谷P2455 [SDOI2006] 線形方程式」問題の参考コードを示します。
#include <iomanip>
#include <iostream>
#include <cmath>
using namespace std;
const int MAXN = 105;
int equation_count;
double coefficient_matrix[MAXN][MAXN];
const double EPS = 1e-8;
int solve_linear_system()
{
int current_row = 1; // 未処理の最も早い行番号
for (int current_col = 1; current_col <= equation_count; current_col++) // 列を順に処理
{
/* ステップ1：非ゼロの要素を探す */
int pivot_row = current_row; // 未処理の行からpを探す
while (pivot_row <= equation_count && abs(coefficient_matrix[pivot_row][current_col]) < EPS) pivot_row++;
if (pivot_row > equation_count) continue; // 有効なpが見つからない場合はこの列をスキップ
/* ステップ2：行の入れ替え */
for (int i = 1; i <= equation_count + 1; i++)
swap(coefficient_matrix[pivot_row][i], coefficient_matrix[current_row][i]);
/* ステップ3：ピボット要素を1にする */
double pivot_value = coefficient_matrix[current_row][current_col]; // 事前に保存
for (int i = 1; i <= equation_count + 1; i++)
coefficient_matrix[current_row][i] /= pivot_value; // coefficient_matrix[current_row][current_col]を1に調整
/* ステップ4：他の行からこの列を消去 */
for (int i = 1; i <= equation_count; i++)
{
if (i == current_row) continue; // 自分自身を消去してはいけない
double factor = coefficient_matrix[i][current_col]; // 事前に保存
for (int j = 1; j <= equation_count + 1; j++)
coefficient_matrix[i][j] -= coefficient_matrix[current_row][j] * factor; // coefficient_matrix[i][current_col]を消去
}
current_row++;
}
if (current_row == equation_count + 1) return 1; // 一意な解が存在
// 残りの行を確認
for (int i = current_row; i <= equation_count; i++)
if (abs(coefficient_matrix[i][equation_count + 1]) >= EPS) return -1; // 解が存在しない
return 0; // 無数の解が存在
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr), cout.tie(nullptr);
cin >> equation_count;
for (int i = 1; i <= equation_count; i++)
for (int j = 1; j <= equation_count + 1; j++)
cin >> coefficient_matrix[i][j];
int result = solve_linear_system();
if (result == 1)
{
for (int i = 1; i <= equation_count; i++)
cout << 'x' << i << '=' << fixed << setprecision(6) << coefficient_matrix[i][equation_count + 1] << '\n';
}
else cout << result << '\n';
return 0;
}
「洛谷P3389 【模板】高斯消元法」の簡略化されたコード例を以下に示します：
const double EPS=1e-8;
bool solve_system()
{
for(int current_row=1; current_row<=equation_count; current_row++)
{
int pivot_row=current_row;
while(abs(coefficient_matrix[pivot_row][current_row])<EPS && pivot_row<=equation_count)
pivot_row++;
if(pivot_row>equation_count) return false; // 一意な解が存在しない
// 行の入れ替え
for(int j=1; j<=equation_count+1; j++)
swap(coefficient_matrix[current_row][j], coefficient_matrix[pivot_row][j]);
// ピボット行を正規化
double pivot_val=coefficient_matrix[current_row][current_row];
for(int j=current_row; j<=equation_count+1; j++)
coefficient_matrix[current_row][j]/=pivot_val;
// 他の行から消去
for(int j=1; j<=equation_count; j++)
{
if(j==current_row) continue;
double factor=coefficient_matrix[j][current_row];
for(int k=current_row; k<=equation_count+1; k++)
coefficient_matrix[j][k]-=factor*coefficient_matrix[current_row][k];
}
}
return true;
}