【sympy&方程式】pythonで方程式とか微積分を解く

+ クリックでオープン

# 確認用のグラフ
def graph(equation1, equation2, save_name: str):
    def scatter(x, y, name, color: str):
        d = go.Scatter(
            mode='lines+markers',
            x=x, y=y, name=name,
            marker=dict(symbol='circle'),
            line=dict(dash='solid', color=color),
            hovertemplate=f"{name}<br>"
            + 'x: %{x} <br>'
            + 'y: %{y} <br>'
            + "<extra></extra>",
        )
        return d

    # 横軸の値
    yoko = np.arange(-10, 10)
    plot = []
    zipper = zip((equation1, equation2), ('red', 'blue'))
    for num, (equation, color) in enumerate(zipper, 1):
        # sympy.subsだと配列を代入できないので新規で定義
        func = sympy.lambdify([x, y], equation)
        # xには横軸の値yokoを、yは「y=」として計算したいので0を代入
        tate = func(yoko, 0)
        d = scatter(x=yoko, y=tate, name=f"eq{num}", color=color)
        plot.append(d)
    # eq1でのtateの値
    # [166 131 100  73  50  31  16   5  -2  -5  -4   1  10  23  40  61  86 115
    #  148 185]

    # eq2でのtateの値
    # [-50 -45 -40 -35 -30 -25 -20 -15 -10  -5   0   5  10  15  20  25  30  35
    #   40  45]

    layout = go.Layout(template=template.plotly_layout(),)

    fig = go.Figure(data=plot, layout=layout)
    fig.show(config=template.plotly_config())

    pio.write_html(
        fig, f"{save_name}.html",
        auto_open=False, config=template.plotly_config()
    )
    pio.write_image(fig, f"{save_name}.png")

graph(equation1=eq1, equation2=eq2, save_name='simultaneous')

微分系（diff, dsolve）

ここから急に年代が上がって高校数学。微分方程式については特殊解は求めずに一般解のみとする。

微分

# 微分
eq = 2 * x ** 2 + 3 * x - 4
diff = sympy.diff(eq)
print(diff)
# 4*x + 3
print(type(diff))
# <class 'sympy.core.add.Add'>

微分はsympy.diffで行う。ここではxで微分した。もちろん三角関数も微分可能。

# 三角関数の微分も対応
sin = 2 * sympy.sin(x)
diff = sympy.diff(sin)
print(diff)
# 2*cos(x)

2変数の場合は変数を指定しないとエラーになる。

# 2変数の場合は変数を選択
eq = x ** 3 + x + 3 * y ** 2 + 10
# どの変数かを指定しないとエラー
diff = sympy.diff(eq)
# Since there is more than one variable in the expression, the
# variable(s) of differentiation must be supplied to differentiate x**3
# + x + 3*y**2 + 10

diff = sympy.diff(eq, x)
print(diff)
# 3*x**2 + 1
diff = sympy.diff(eq, y)
print(diff)
# 6*y

微分方程式

# fを関数として定義
f = sympy.Function('f')
print(f)
# f
print(type(f))
# <class 'sympy.core.function.UndefinedFunction'>

# 定数としてaを定義
a = sympy.Symbol('a')

# 微分方程式を解く
# Object of type "None" cannot be calledPylancere(portOptionalCall)のエラーあり
eq = sympy.diff(f(x)) - a * f(x)
print(eq)
# -a*f(x) + Derivative(f(x), x)
print(type(eq))
# <class 'sympy.core.add.Add'>

dsolve = sympy.dsolve(eq)
print(dsolve)
# Eq(f(x), C1*exp(a*x))
print(type(dsolve))
# <class 'sympy.core.relational.Equality'>

微分方程式を解くにはdiff solveという意味だろう、sympy.dsolveを使用する。微分方程式と計算結果のEq(f(x), C1*exp(a*x))はJupyter notebookなどでprintを使用せずに出力すると以下のようにキレイな数式となる。

$\begin{aligned} \cfrac{d}{dx}f(x)\ &=\ af(x)\\ f(x)\ &=\ C_1e^{ax} \end{aligned}$

解が2定数を含む式の演算なら自動的に$C_2$のように定数が決められる。

eq = sympy.diff(f(x), x, 2) - a ** 2 * f(x)
print(eq)
# -a**2*f(x) + Derivative(f(x), (x, 2))
dsolve = sympy.dsolve(eq)
print(dsolve)
# Eq(f(x), C1*exp(-a*x) + C2*exp(a*x))

この時の微分方程式とその解は以下。

$\begin{aligned} \cfrac{d^2}{dx^2}f(x)\ &=\ a^2f(x)\\ f(x)\ &=\ C_1e^{-ax}\ +\ C_2e^{ax} \end{aligned}$

また、解けない（と自分は思っている）微分方程式についてはエラーが吐かれる。

# 解けるかどうかは知らんけど、エラーが出る微分方程式もある
eq = sympy.diff(f(x), x, 2) - x**2 + f(x) ** 4
print(eq)
# -x**2 + f(x)**4 + Derivative(f(x), (x, 2))
dsolve = sympy.dsolve(eq)
#     raise NotImplementedError(dummy + "solve" + ": Cannot solve " + str(eq))
# NotImplementedError: solve: Cannot solve -x**2 + f(x)**4 + Derivative(f(x), (x, 2))

この時の微分方程式は以下。

$\cfrac{d^2}{dx^2}f(x)\ =\ -x^2\ +\ f^4(x)$

積分系（integrate）

微分ときたら積分と思うのは私だけだろうか。簡単にざっくりイメージは微分の逆。

不定積分

eq = x ** 2 + 6 * x + 9
integrate = sympy.integrate(eq)
# 積分定数は表記されない
print(integrate)
# x**3/3 + 3*x**2 + 9*x
print(type(integrate))
# <class 'sympy.core.add.Add'>

積分についてはsympy.integrateで可能。不定積分であるが、積分定数（通常は$C$が使われる）は書かれないのでその点には注意。今回の積分は以下。

$\displaystyle\int dx\ x^2\ +\ 6x\ +\ 9$

もちろん三角関数も可能だが、積分結果が$\tan(x)$ではないのが気になる。

eq = 1 / (sympy.cos(x)**2)
integrate = sympy.integrate(eq)
# tan(x)にはなってくれない
print(integrate)
# sin(x)/cos(x)

定積分

eq = x ** 2 + 6 * x + 9
integrate = sympy.integrate(eq, (x, 0, 3))
print(integrate)
# 63
print(type(integrate))
# <class 'sympy.core.numbers.Integer'>

定積分を計算したい場合はsympy.integrateの引数に初めの値と終わりの値を入れればいい。今回は0から3まで積分してみた。式は以下。

$\displaystyle\int_{0}^{3} dx\ x^2\ +\ 6x\ +\ 9$

計算結果は普通に足し算などが可能。

# 普通に計算可能
print(integrate + 1)
# 64

単位を付与（sympy.physics.units）

最後は単位について。というのも執筆者自身、普段からastropyで数値に単位をつけている。単位をつけることでメートルからミリメートルのような単位換算を自動計算してくれて楽。

astropyと単位については以下参照。

astropyの単位だとエラー

# 変数に単位を付与することはできない
m = sympy.Symbol('m') * u.s
#     m = sympy.Symbol('m') * u.s
# TypeError: unsupported operand type(s) for *: 'Symbol' and 'IrreducibleUnit'

今まで通りastropyを使用して単位を付与してみる。まずは変数にそのまま単位を付与してみるがエラー。代入する値に単位をつけて代入しようとしてもエラー。

# 代入するときに単位を付与してもエラー
eq = 2 * x + 5
val = 3 * u.s
subs = eq.subs(x, val)
# RecursionError: maximum recursion depth exceeded in comparison

しかも変数ではない項に単位を付与してもエラー。まじですか。

# 変数ではない項に単位を付与してもエラー
eq = 2 * x + 5 * u.s
# RecursionError: maximum recursion depth exceeded in comparison

と思ったらそもそも単位ありと単位なしを同時に式に書くこと自体無理だった。

eq = 2 + 5 * u.s
# astropy.units.core.UnitConversionError: Can only apply 'add' function to dimensionless quantities when other argument is not a quantity (unless the latter is all zero/infinity/nan)

代入完了後に単位を付与することは可能。

# 代入後の単位の付与は可能
eq = 2 * x + 5
val = 3
subs = eq.subs(x, val) * u.s
print(subs)
# 11.0 s

しかし、これだと単位換算の恩恵を享受することができない。一応、sympyでの単位換算系はsymdimというモジュールで可能らしいがanacondaになかったので執筆者は入れていない。以下に参考URLを書いておく。

sympyの単位だと大丈夫

# sympyの単位を使用することが可能
val = 12 * units.m
print(val)
# 12*meter
print(type(val))
# <class 'sympy.core.mul.Mul'>

じゃあ単位をつけられないのかというとそういうわけではない。astropyと使わなければいい。調べているとちょうどsympyで単位を付与できる機能があるっぽいのでこれを使ってみる。

モジュールのimportはfromを使う例が結構あったが、個人的にはimport sympy.physics.units as unitsの方がしっくりきたのでこれでimport。units.mでメートル単位を使用することが可能。

sympyの単位付与の状態で代入を行うとエラーは出ない。

# sympyの単位だとエラーが起きることはない
eq = 2 * x + 5
subs = eq.subs(x, val)
print(subs)
# 24*meter + 5
print(type(subs))
# <class 'sympy.core.add.Add'>

単位換算はunits.convert_to（変換したい単位つき変数, 変換先単位）で可能。ここではメートルをミリメートル、センチメートル、メートルのまま、そしてキロメートルに変換した。小数が使えないっぽいので数字の表記は分数となっている。

# 単位つきの値の単位変換。小数は使用されず、分数表記となる
print(units.convert_to(val, units.mm))
# 12000*millimeter
print(units.convert_to(val, units.cm))
# 1200*centimeter
print(units.convert_to(val, units.m))
# 12*meter
print(units.convert_to(val, units.km))
# 3*kilometer/250

複数の単位を組み合わせて、単位換算をすることも可能。ここではジュールJと秒sからワットWに変換した。

# ジュールJと秒sからワットWを計算
val1 = 10 * units.s
val2 = 100 * units.J
print(val1)
# 10*second
print(val2)
# 100*joule

val = val2 / val1
print(val)
# 10*joule/second
print(type(val))
# <class 'sympy.core.mul.Mul'>

print(units.convert_to(val, units.W))
# 10*watt

pythonに解いてもらう

今回はsympyを使って変数を定義し、方程式などを解くという方法について解説した。執筆者自身、pythonを学んだのは研究室配属された大学4年生からでこの頃にはもうほとんどの単位を取得してあまり方程式を使用する機会がなかった。

また、大学院の講義でも数学系というよりは宇宙の概論などを主に学んだので今回のsympyを使用する機会はあまりなかった。

しかし、早い時期にこの記事を見ていただいた方は、sympyを使って方程式などを解いてもらうことで検算が可能だ（くれぐれもpythonの結果をそのまま解答にして力をつけないままなのは避けるように）。

今回のような便利な機能をまた見つけたら記事にしていく。

関連記事