読者です 読者をやめる 読者になる 読者になる

Gauss法

はじめに

Gauss法は連立1次方程式を解くためのよく知られた方法となります。
条件によっては不整合が発生することもありますが、アルゴリズムとしてはとてもシンプルなものになります。

C言語による最新アルゴリズム事典 (ソフトウェアテクノロジー)

C言語による最新アルゴリズム事典 (ソフトウェアテクノロジー)

式変換

Gauss法は連立方程式の式変換を行う方法のようなものです。
例として、 \eqref{eq:q} を考えます。 $$ \begin{equation} \begin{cases} 2x_1 + 5x_2 + 7x_3 &= 23 \\ 4x_1 + 13x_2 + 20x_3 &= 58 \\ 8x_1 + 29x_2 + 50x_3 &= 132 \end{cases}
\label{eq:q} \end{equation} $$

Gauss法では、$ i $番目の式から$ x_{i} $より番号の若い変数を消していくという処理をしていきます。

$ \eqref{eq:q} $を例にとってみると、 まず2, 3番目の式から$ x_1 $を消します。 $$ \begin{equation} \begin{cases} 2x_1 + 5x_2 + 7x_3 &= 23 \\ 3x_2 + 6x_3 &= 12 \\ 9x_2 + 22x_3 &= 40 \end{cases}
\end{equation} $$

その後、3番目の式から $ x_2 $を消します。 $$ \begin{equation} \begin{cases} 2x_1 + 5x_2 + 7x_3 &= 23 \\ 3x_2 + 6x_3 &= 12 \\ 4x_3 &= 4 \end{cases}
\end{equation} $$

この変換から、まず$ x_3 $が自明になります。
そして、その結果を使って順番に$ x_2, x_1 $と求めていくことができます。

実行できる形にしてみる

これをコードに落としてみます。

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np

def gauss(a, b):
  def getAnswer(a, b):
    res = np.copy(b)
    for i in range(0, len(b))[::-1]:
      for j in range(i + 1, len(b)):
        res[i] -= a[i][j] * res[j]
      res[i] /= a[i][i]
    print(res)

  if not len(a) == len(b):
    raise Exception('matrix and values must be same length')
  for i in range(0, len(a) - 1):
    for j in range(i + 1, len(a)):
      t = a[j][i] / a[i][i]
      for k in range(i, len(a)):
        a[j][k] -= t * a[i][k]
      b[j] -= t * b[i]
  getAnswer(a, b)

こんな感じですかね。

試しに実行してみると、

qa = np.array([[2,5,7],[4,13,20],[8,29,50]])
qb =  np.array([23,58,132])
gauss(qa, qb)
>>> [3 2 1]

これはさきほどの例を入れたものです。
答えも問題なさそうです。

おわりに

Gauss法は係数行列(プログラム中でいう変数a)の対角成分が0に近くなると不正号を起こしてしまいます。
その問題を回避する方法としてピボット選択というのがよく知られています。