アルゴリズム

Tweet

掃き出し法

ヨビノリヨビノリ　【大学数学】線形代数入門⑤(連立方程式：掃き出し法)【線形代数】を参考に吐き出し法を実装してみた。

# coding UTF-8
import numpy as np


def find_one(mat: np.ndarray, t: int) -> [bool, int]:
    """
    t行目以下で係数が1の行を探す
    :param mat: 式全て
    :param t: 対角線行列の位置
    :return: 見つかった場合[True,行番号]を返す。
    見つからなかった場合は[False, -1.]を返す。
    """
    row_num = mat.shape[0]
    for r in range(t, row_num, 1):
        if mat[r][t] == 1.:
            return [True, r]

    return [False, 0]


def swap(mat: np.ndarray, sr: int, dr: int) -> None:
    """
    行を入れ替える
    :param mat: データ
    :param sr: 入れ替え前の行番号
    :param dr: 入れ替え先の行番号
    :return: なし(matの配列を書き換える)
    """
    tmp = np.copy(mat[sr])
    mat[sr] = mat[dr]
    mat[dr] = tmp


def normarize(mat: np.ndarray, t: int) -> None:
    """
    ターゲットの行を体格行列の成分で割る
    :param mat:　データ
    :param t: 対角行列の位置
    :return: なし。解のない式の場合は例外を返す
    """
    origin = mat[t][t]
    if origin == 0.:
        raise Exception("解がありません")

    mat[t] = mat[t] / origin


def sub_mat(mat: np.ndarray, t: int) -> None:
    """
    ターゲット以外の行を引く
    :param mat: データ
    :param t: 引き算のもとになる行番号
    :return: なし
    """
    row_num = mat.shape[0]
    for i in range(row_num):
        if i != t:
            origin = np.copy(mat[t])
            compo = mat[i][t]
            origin *= compo
            mat[i] -= origin


def main():
    mat = [
        [3., 1., 2., 2., 5.],
        [1., -1., 2., 2., 7.],
        [2., 1., 3., 2., 6.],
        [1., 1., 1., 1., 2.]
    ]

    mat = np.array(mat)

    for t in range(len(mat)):
        # ターゲットの対角行列以下の行に1があるか？
        ret, r = find_one(mat, t)
        if ret:
            # 1が見つかった
            if t != r:
                # ターゲットより下で1が見つかったので入れ替え
                swap(mat, t, r)
        else:
            # 1が見つからないので、自分自身の行を正規化する
            try:
                normarize(mat, t)
            except Exception as e:
                print("解がありません")
                return

        # 自分以外の行を引く
        sub_mat(mat, t)

    print(mat)


if __name__ == '__main__':
    main()