【ガウスの消去法(Gaussian elimination method)】Pythonによるガウスの消去法アルゴリズム

#**はじめに**
Pythonでガウスの消去法アルゴリズムを実装しましたので、掲載します。

また、Fortran、Rubyでも同様にプログラムを書いてますでので、必要に応じてご参考ください。

#**関連**
**・Fortran**
[【ガウスの消去法(Gaussian elimination method)】Fortranによるガウスの消去法アルゴリズム](https://qiita.com/takuyanin/items/3addb4d6e4ff85673e3b)

**・Ruby**
[【ガウスの消去法(Gaussian elimination method)】Rubyによるガウスの消去法アルゴリズム(構造化プログラミング/ オブジェクト指向プログラミング)](https://qiita.com/takuyanin/items/61fc7189550f4c969306)

#**プログラム**
###**メインプログラム**
計算を行うメインプログラムファイルを準備します。
このプログラムに必要モジュールをインポートします。

```calc.py
import sys
from gauss_elimination import GaussElimination as GE

if __name__ == "__main__":
    try:
        obj = GE()
        # obj = GaussElimination()
        obj.display_matrixA_and_matrixB()
        obj.calculate()
        obj.display_matrixX()
    except Exception as e:
        print("=== エラー内容 ===")
        print('type:\t{0}'.format(type(e)))
        tr = sys.exc_info()[2]
        print("message:{0}".format(e.with_traceback(tr)))
    else:
        print("finish (no error)")
    finally:
        print("all finish")
```

###**モジュール**
**ガウスの消去法クラス**を用意し、その中に必要メソッドを定義します。

```gauss_elimination.py
import pandas as pd
import numpy as np
from math import fabs
from sys import argv

class GaussElimination():
    """
    A. インスタンス変数
        a: 行列A
        b: 行列B
        n_a: 行列Aの次元
        n_b_col: 行列Bの列数

    B. メソッド
        __init__: コンストラクタ
        display_matrixA_and_matrixB: 行列A, Bを出力
        display_matrixX: 行列X(解)を出力
        calculate: ガウスの消去法で計算を行う
        __partial_pivoting: 部分ピボット選択
    """

    def __init__(self):
        # 初期値を設定
        self.a = [
          [-1, 2, 2],
          [2, 3, -1],
          [1, 1, -1]
        ]
        self.b = [
          [1, 2],
          [-3, 2],
          [-2, 2]
        ]

        if len(argv) == 3:
            self.a = np.array(pd.read_csv(argv[1], header=-1, sep=',', dtype='float32'))[:,:]
            self.b = np.array(pd.read_csv(argv[2], header=-1, sep=',', dtype='float32'))[:,:]

        self.n_a = len(self.a)
        self.n_b_col = len(self.b[0])

    def display_matrixA_and_matrixB(self):
        # 行列Aを出力
        print("A =")
        for i in range(0,self.n_a):
            for j in range(0, self.n_a):
                print('{:+.5f}\t'.format(self.a[i][j]), end="")
            print()
        # 行列Bを出力
        print("B =")
        for i in range(0,self.n_a):
            for j in range(0, self.n_b_col):
                print('{:+.5f}\t'.format(self.b[i][j]), end="")
            print()

    def display_matrixX(self):
        # 行列X(解)を出力
        print("X (satisfying AX=B) =")
        for i in range(0, self.n_a):
            for j in range(0, self.n_b_col):
                print('{:+.5f}\t'.format(self.b[i][j]), end="")
            print()

    def calculate(self):
        # 前進消去
        for j in range(0, self.n_a):
            self.__partial_pivoting(j)
            # ピボット列の掃き出し
            for i in range(j + 1, self.n_a):
                r = float(self.a[i][j]) / float(self.a[j][j])
                # 各行減算
                # 1. 行列Aについて
                for k in range(j, self.n_a):
                    self.a[i][k] = float(self.a[i][k]) - float(self.a[j][k]) * r
                # 2. 行列Bについて
                for k in range(0, self.n_b_col):
                    self.b[i][k] = float(self.b[i][k]) - float(self.b[j][k]) * r
        # 後退代入
        for j in range(0, self.n_b_col):
            for i in reversed(range(0, self.n_a)):
                r = float(self.b[i][j])
                # シグマ演算
                for k in range(i + 1, self.n_a):
                    # ここでのrはx
                    r = r - float(self.a[i][k]) * float(self.b[k][j])
                self.b[i][j] = r / float(self.a[i][i])

    def __partial_pivoting(self, j):
        # 最大値を記憶
        amax = fabs(self.a[j][j])
        imax = j

        for i in range(j + 1, self.n_a):
            r = fabs(self.a[i][j])
            if r > amax:
                # 最大値を変更
                amax = r
                imax = i

        if j == imax:
            return

        # ピボット行と最大値行が異なれば交換
        # 1. 行列Aについて
        for k in range(0, self.n_a):
            r = self.a[imax][k]
            self.a[imax][k] = self.a[j][k]
            self.a[j][k] = r
        # 2. 行列Bについて
        for k in range(0, self.n_b_col):
            r = self.b[imax][k]
            self.b[imax][k] = self.b[j][k]
            self.b[j][k] = r
```

###**実行**
```
$ python3 calc.py
```
行列を定義したCSVファイル(例： a.csv, b.csv)を指定することもできます。

```
$ python3 calc.py a.csv b.csv
```

一応、CSVファイルの例としてa.csvとb.csvの中身を下記に載せておきます。

###**CSVファイル**
```a.csv
-1, 2, 2
2, 3, -1
1, 1, -1
```

```b.csv
1, 2
-3, 2
-2, 2
```

###**結果**
```
$ python3 calc.py a.csv b.csv
A =
-1.00000	+2.00000	+2.00000
+2.00000	+3.00000	-1.00000
+1.00000	+1.00000	-1.00000
B =
+1.00000	+2.00000
-3.00000	+2.00000
-2.00000	+2.00000
X (satisfying AX=B) =
+1.00000	-6.00000
-1.00000	+3.00000
+2.00000	-5.00000
finish (no error)
all finish
```

#**おわりに**
Python, Rubyそれぞれで同様の処理をプログラムに組んでみましたが、感じたことは、**Pythonのインデント記法って素晴らしい**、ということです。笑
これがPythonコードの可読性を大きく高めていることは間違いないでしょう。

Python, Rubyどちらのコードも見た目が似通っているだけに、可読性の観点から優劣をつけるのは難しいですが、インデントの恩恵を感じられてよかったです。

こちらの記事が役にたったという方は、いいね、よろしくお願いします(^^)