Pythonで線形方程式を解く

How much gold :: CheckiOを解く方法として、何か簡単な方法はあるんだろうなと思いながら、線形方程式を解く方法を採用したわけです。

線形方程式を簡単に解く方法はないかな？と探したところ、 SymPyを使って式を与えるだけで解くことが出来ました。テストその1をベタ書きしています。

ソースコードおよびセットアップ方法は、 Pythonを使って一瞬で連立方程式を解く – Qiita を参考にしています。要するに、右辺がゼロとなる方程式を、連立方程式を解くのに必要な分渡しているわけです。

from fractions import Fraction
from sympy import *

g, t, i, c = symbols('g t i c')
out = solve([
    g * Fraction(1, 1) + t * Fraction(1, 1) - Fraction(1, 2),
    g * Fraction(1, 1) + i * Fraction(1, 1) - Fraction(1, 3),
    g * Fraction(1, 1) + c * Fraction(1, 1) - Fraction(1, 4),
    g * Fraction(1, 1) + t * Fraction(1, 1) + i * Fraction(1, 1) + c * Fraction(1, 1) - Fraction(1, 1),
], [g, t, i, c])

out[g]  # -> Fraction(1, 24)

ただ、CheckiOでは外部ライブラリの利用はできません。また、式を動的に生成する方法があるかも分かりませんので、複数ケースに対応する方法が分かりませんでした。

結局、ガウスの消去法 – Wikipedia を使って、自前実装で解くことにしました。

def gauss_elimination(m):
    '''mは[g, t, i, c, value]というベクトルを4本持つ行列'''
    n = len(m)
    # 前進消去
    for i in range(n):
        v = m[i][i]
        if v == 0:
            # ゼロ以外の値を持つ行と交換する
            for k in range(i + 1, n):
                if m[k][i] != 0:
                    m[i], m[k] = m[k], m[i]
                    v = m[i][i]
                    break
        for j in range(n + 1):
            m[i][j] /= v
        for k in range(i + 1, n):
            v = -m[k][i]
            for j in range(n + 1):
                m[k][j] += m[i][j] * v

    # 後退代入
    for i in range(n-1, 0, -1):
        v = m[i][i]
        m[i][i] /= v
        m[i][n] /= v
        for k in range(i - 1, -1, -1):
            v2 = m[k][i] * m[i][n]
            m[k][n] -= v2
            m[k][i] = Fraction(0, 1)
    return m

これでgauss_elimination()が解けたときにはすべての変数の値が求められた状態になります。良かった良かった。

…と思ったら、最短1行で書くことができるらしい。 G+T=a, G+I=b, G+C=c, G+T+I+C=1から、G = (a + b + c - 1) / 2が導出できて、 a, b, cを、それぞれGが含まれている場合の式に変換することで計算できると (テスト2の方は、1 – 式とすればGが含まれる方になる)。なるほど、勉強になりました。

技術情報棚卸し(平日限定)

todoa2cの技術情報棚卸しです。平日限定ってことはアレだ。言わせんな恥ずかしい。

Comments