連立一次方程式の数値計算法

数値計算
  1. ガウス・ジョルダン法
    1. アルゴリズム
    2. ピボット数
    3. python スクリプト
  2. 参考

この記事では，連立1次方程式の数値計算法について説明する．
また，python のサンプルコードで理解を深める．

ガウス・ジョルダン法

ガウス・ジョルダン法 (Gauss-Jordan) は，方程式の各項を順番に削除していくアルゴリズムで，連立1次方程式の数値計算法として用いられる手法である．
この手法は，掃き出し法とも呼ばれる．

アルゴリズム

次のような $n$ 元方程式で説明する．

$$
\begin{aligned}
a_{11} x_{1} + a_{12} x_{2} + a_{13} x_{3} + \cdots + a_{1n} x_{n} &= b_1 \\
a_{21} x_{1} + a_{22} x_{2} + a_{23} x_{3} + \cdots + a_{2n} x_{n} &= b_2 \\
a_{31} x_{1} + a_{32} x_{2} + a_{33} x_{3} + \cdots + a_{3n} x_{n} &= b_3 \\
\vdots & \\
a_{n1} x_{1} + a_{n2} x_{2} + a_{n3} x_{3} + \cdots + a_{nn} x_{n} &= b_n \\
\end{aligned}
$$

まず，この方程式の1式目を $a_{11}$ で除算した後，第1項の係数を消去すると，次の形の式になる．
$a’_{ij}$ は計算後の係数を表している．

$$ \begin{aligned} x_{1} + a’_{12} x_{2} + a’_{13} x_{3} + \cdots + a’_{1n} x_{n} &= b_1 \\\ a’_{22} x_{2} + a’_{23} x_{3} + \cdots + a’_{2n} x_{n} &= b_2 \\\ a’_{32} x_{2} + a’_{33} x_{3} + \cdots + a’_{3n} x_{n} &= b_3 \\\ \vdots & \\\ a’_{n2} x_{2} + a’_{n3} x_{3} + \cdots + a’_{nn} x_{n} &= b_n \end{aligned} $$

次に，2式目を $a’_{22}$ で除算した後，第2項の係数を消去する．

$$ \begin{aligned} x_{1} \quad \qquad + a”_{13} x_{3} + \cdots + a”_{1n} x_{n} &= b_1 \\\ x_{2} + a”_{23} x_{3} + \cdots + a”_{2n} x_{n} &= b_2 \\\ a”_{33} x_{3} + \cdots + a”_{3n} x_{n} &= b_3 \\\ \vdots & \\\ a”_{n3} x_{3} + \cdots + a”_{nn} x_{n} &= b_n \end{aligned} $$

これを繰り返すと，

$$
\begin{aligned}
x_{1} \quad \qquad \qquad \qquad \qquad &= b_1^{(n)} \\
x_{2} \quad \qquad \qquad \qquad &= b_2^{(n)} \\
x_{3} \quad \qquad \qquad &= b_3^{(n)} \\
\vdots & \\
x_{n} &= b_n^{(n)}
\end{aligned}
$$

のように解が求まる．

ピボット数

上記の操作で第1項，第2項と係数を消していく際に $a_{11}$や$a’_{22}$など消去したい項の係数で方程式を除算した．
これらの数はピボット数と呼ばれる．

方程式をピボット数で割るので，ピボット数が $0$ または $0$ に近い値の時に計算が行えないもしくは正しい結果が得られないことがある．

このようなときは，ピボット数が対角成分 $a_{ij} (i=j)$ であることを意識すると，方程式の順序を入れ替え，ピボット数を変えることで対応する，

python スクリプト

次にガウス・ジョルダン法の python スクリプトの例を示す．
スクリプトでは，方程式の係数は2次元の配列として扱う

import numpy as np

## ------ params ------ ##
A_mat = np.array([
    [2.0, 1.0, 3.0, 13.0],
    [1.0, 3.0, 2.0, 13.0],
    [3.0, 2.0, 1.0, 10.0]
    ])

eps = 1e-4
dim = A_mat.shape[0]

## ------ main ------ ##
Ai_mat = A_mat
print(Ai_mat)

for i in range(dim):
    pivot = Ai_mat[i, i]

    ## judge pivot ~ 0
    if np.abs(pivot) < eps:
        print('pivot is near 0, exit')
        exit()

    else:
        ## divide by pivot
        Ai_mat[i, :] = Ai_mat[i, :] / pivot

        ## row reduction
        index = np.arange(dim)
        for j in index[index != i]:
            Ai_mat[j, :] = Ai_mat[j, :] - Ai_mat[j, i]*Ai_mat[i, :]
    print(Ai_mat)

result = Ai_mat[:, -1]
print('result')
print(result)

出力結果

出力結果は次のようになる．
係数が順番に1になっていくのが分かる．

[[ 2.  1.  3. 13.]
 [ 1.  3.  2. 13.]
 [ 3.  2.  1. 10.]]
[[ 1.   0.5  1.5  6.5]
 [ 0.   2.5  0.5  6.5]
 [ 0.   0.5 -3.5 -9.5]]
[[  1.    0.    1.4   5.2]
 [  0.    1.    0.2   2.6]
 [  0.    0.   -3.6 -10.8]]
[[ 1.  0.  0.  1.]
 [ 0.  1.  0.  2.]
 [-0. -0.  1.  3.]]
result
[1. 2. 3.]

参考

[1] 三井田 惇郎，須田 宇宙，数値計算法[第2版]，森北出版，2014