Pythonでのガウス・ジョルダン法の実装

Python

ガウス・ジョルダン法（Gauss-Jordan method）は、線形方程式を解くための手法の1つです。この方法は、行列を操作して連立方程式を単純な対角線行列に変換し、それから簡単に解を求めることができます。この記事では、Pythonでガウス・ジョルダン法を実装する方法を紹介します。

ガウス・ジョルダン法とは

ガウス・ジョルダン法は、18世紀末から19世紀初めにかけて、数学者のカール・フリードリッヒ・ガウスとフーゴ・ジョルダンによって独立に発明されました。ガウス・ジョルダン法は、線形方程式を解くための手法であり、係数行列と定数項の行列を結合させた拡大係数行列を作り、行の基本変形を繰り返すことで、拡大係数行列を上三角行列に変形させ、最後に後退代入によって解を求める方法です。この手法を使うことで、線形方程式の解を効率的に求めることができます。以下の手順で解を求めます。

  1. 係数行列と定数項の行列を結合させた拡大係数行列を作る。
  2. 拡大係数行列の左上の要素を1にするため、1行目を左から順に、その行の要素を1で割る。
  3. 1列目の要素を0にするため、2行目から最終行まで、以下の手順を繰り返す。
    1. 2行目の1列目の要素を0にするため、2行目を以下のように更新する。
      • 2行目 - (2行目の1列目の要素 ÷ 1行目の1列目の要素) × 1行目
    2. 3行目の1列目の要素を0にするため、3行目を以下のように更新する。
      • 3行目 - (3行目の1列目の要素 ÷ 1行目の1列目の要素) × 1行目
    3. 最終行の1列目の要素を0にするため、最終行を以下のように更新する。
      • 最終行 - (最終行の1列目の要素 ÷ 1行目の1列目の要素) × 1行目
  4. 拡大係数行列の2行目から最終行まで、左から2列目、3列目、…と同様の手順を繰り返す。
  5. 係数行列と定数項の行列に戻し、解を求める。

この手法を使うことで、線形方程式の解を効率的に求めることができます。

ガウス・ジョルダン法で線型方程式を解く具体例

例として、次のような3つの方程式を考えてみます。

2x + y - z = 8
-3x - y + 2z = -11
-2x + y + 2z = -3

この方程式を行列で表すと以下のようになります。

| 2  1 -1 |   | x |   | 8 |
| -3 -1  2 | x | y | = |-11|
| -2  1  2 |   | z |   | -3|

これをガウス・ジョルダン法を用いて解いてみましょう。まず、この行列を拡大係数行列として表します。

| 2  1 -1  8 |
| -3 -1  2 -11|
| -2  1  2 -3 |

次に、最初の列の2行目と3行目の係数を0にするために、1行目を2で割ったものを2行目に足し、1行目を-1でかけて3行目に足します。

| 1 1/2 -1/2 4 |
| 0 -5/2 3/2 -3|
| 0  2   3  -5 |

次に、2行目の2列目を1にするために、2行目を-2/5倍します。

| 1 1/2 -1/2 4 |
| 0   1  -3/5 3/5|
| 0  2   3  -5 |

次に、1行目の2列目と3列目を0にするために、2行目を-1/2倍して1行目に足し、3行目を-2倍して1行目に足します。

| 1   0 -1/5 17/5|
| 0   1 -3/5  3/5 |
| 0   0  9/5 -27/5|

最後に、3行目の3列目を1にするために、3行目を5/9倍します。

| 1   0 -1/5 17/5|
| 0   1 -3/5  3/5 |
| 0   0   1  -3  |

これで、拡大係数行列を単純な対角線行列に変換することができました。最後に、この行列から解を求めることができます。

x = 2
y = 3
z = -3

これが、元の方程式の解となります。

ガウス・ジョルダン法のPythonでの実装

それでは、Pythonでこのガウス・ジョルダン法を実装してみましょう。以下は、この例題を解くためのPythonコードです。

def gauss_jordan(A):
    n = len(A)
    for i in range(n):
        pivot = A[i][i]
        for j in range(i, n+1):
            A[i][j] /= pivot
        for j in range(n):
            if i != j:
                pivot = A[j][i]
                for k in range(i, n+1):
                    A[j][k] -= pivot * A[i][k]
    return [A[i][n] for i in range(n)]

A = [[2, 1, -1, 8], [-3, -1, 2, -11], [-2, 1, 2, -3]]
print(gauss_jordan(A))

このコードを実行すると、以下のような出力が得られます。

[2.0, 3.0, -3.0]

これが、先ほど求めた解と一致していることが確認できます。この実装では、拡大係数行列を2重のforループで操作しています。まず、各行を対角成分で割り、その後、各行の対角成分以外の要素を0にするために、他の行の対角成分を利用して行列の各要素を変換しています。

この実装を改良することで、パフォーマンスを向上させることができます。例えば、部分選択を行うことで、係数行列の行の順序を変えることで、計算誤差を減らすことができます。また、部分選択を行わない場合には、ピボット要素が0になる場合には別のピボット要素を選ぶ必要があります。さらに、行列式が0になる場合には、解が存在しないか、無限に存在する可能性があります。

以上が、Pythonでガウス・ジョルダン法を実装する方法と、その具体例とコードです。この方法は、線形方程式を解く場合に非常に有用な手法であり、Pythonを用いて簡単に実装することができます。

問題解決のための「アルゴリズム×数学」が基礎からしっかり身につく本 [ 米田優峻 ]