AstroBlog

宇宙系大学院生の戯言

PythonのSymPyでの3次方程式の解き方

SPONSORED LINK

ちょっと3次方程式を解く必要に迫られて、しかしこれを自分で解くのは非常にめんどくさいなと。 未知数の次数について整理するのすらめんどくさい。 そこで少し調べていたら、PythonのSymPyというモジュールを使えば、方程式を簡単に解くことができるということを知りました。

今回は、SymPyで方程式を解くときのコード例を紹介します。

まず、解きたい方程式はこちらです。 $$\frac{nt}{\left(mt + nt\sigma^{2}\frac{A}{p}\right)^{½}} = 10$$ これをtについて解きたいという状況です。
t以外の文字は本当は具体的な値があるのですが、今回は適当な定数として扱います。

さて、肝心のSymPyでの書き方ですが、これにはどうやら2つのタイプがあって、1つは

from sympy import *

# t以外の文字の値はここでは適当な値とする

def main():
    var("t")      # tを未知数とする
    n = 10
    m = 20
    sigma = 30
    A = 40
    p = 50
    eq = Eq(n*t/(m*t + n*t + sigma*A/p)**(1/2) - 10)    # 方程式を書く
    ans = solve(eq)      # 方程式の解を計算
    print(ans)      # 表示

if __name__ == '__main__':
    main()

というやり方です。
この場合結果はこうなります。

[30.7797338380595]

関数には別にしなくても良さそうですが、なんとなくここでは関数を作りました。

もう1つの書き方は、

import sympy

# t以外の文字の値はここでは適当な値とする
n = 10
m = 20
sigma = 30
A = 40
p = 50

t = sympy.Symbol('t')     # ここでtを未知数と定義
exp = n*t/(m*t + n*t + sigma*A/p)**(1/2) - 10     # 方程式を定義
print(sympy.solve(exp))    # 表示

という書き方です。 importのところが少し違いますね。 この場合は結果はこうなります。

[30.7797338380595]

当然ですが、上のやり方と同じ結果になりますね。

今回はこんな感じです。SymPyは微分方程式でも使えるようなので、少しずつ習得していきたいです(あまり使わないかも)。