読者です 読者をやめる 読者になる 読者になる

yukiのブログ

作ったものなど

Sympyで方程式を解いてみる

Sympyで方程式を解いてみる

Sympyがまだ面白そうだったのでもっと色々やりました。

とりあえず前に受けた微分方程式講義ノートを引っ張り出してそこに書いてある練習問題を解いていきました。

微分方程式はすぐ解ける形もあればあれこれ手を焼いてやらないと解けないものもあるようです。すぐに解けないときは小手先の変更(hintを変えるとか)だけではうまくいかないことが多いので、さっさと地道な方法で解いたほうがいいと思います。

この記事はjupyter notebookで書いたのをmarkdownに変換したものです。notebookはここにあります。

# 初期処理
%matplotlib inline
from sympy import *
from sympy.plotting import plot
from IPython.display import display
from  sympy.solvers.ode import *
init_printing()

# デフォルトだと微分方程式を解いたときに誤差関数(error function)を使うのでそれを書き直す関数
def expand_erfs(expr):
    """ 誤差関数を開く """
    return expr.replace(erf, lambda x: 2/sqrt(pi) * Integral(exp(-t**2), (t, 0, x)))\
               .replace(erfi, lambda x: 2/sqrt(pi) * Integral(exp(t**2), (t, 0, x))).simplify()

# 積分定数を潰す関数constantsimp()は積分定数を引数として
# 渡さなければいけないので、C1,C2,...をシンボルとして定義しなければいけない。
# それだとめんどくさいので、積分定数でないシンボルを渡せばすむ関数があると便利。
def constsimp(expr, variable):
    """ variable以外の自由変数をconstantsimp()に渡す """
    return constantsimp(expr, [c for c in ans.free_symbols if c != variable])

x = symbols("x")
a, b, c, t = symbols("a b c t")
y, f, g, P, Q, R = symbols('y f g P Q R', cls=Function)

方程式

expr = Eq(x**2+3*x-4)
display(expr, solveset(expr, x)) # solve()という関数もあるがsolveset()が推奨されている

\(\displaystyle x^{2} + 3 x - 4 = 0\)

\(\displaystyle \left\{-4, 1\right\}\)

expr = x**3+3*x**2+3*x-4
display(Eq(expr), solveset(expr))

\(\displaystyle x^{3} + 3 x^{2} + 3 x - 4 = 0\)

\(\displaystyle \left\{-1 + \sqrt[3]{5}, -1 - \frac{\sqrt[3]{5}}{2} - \frac{\sqrt{3} i}{2} \sqrt[3]{5}, -1 - \frac{\sqrt[3]{5}}{2} + \frac{\sqrt{3} i}{2} \sqrt[3]{5}\right\}\)

a, b, c = symbols("a b c")
solve(a*x**2+b*x+c, x)

\(\displaystyle \left [ \frac{1}{2 a} \left(- b + \sqrt{- 4 a c + b^{2}}\right), \quad - \frac{1}{2 a} \left(b + \sqrt{- 4 a c + b^{2}}\right)\right ]\)

solveset(Gt(x**2, 0), domain=S.Reals)

\(\displaystyle \left(-\infty, 0\right) \cup \left(0, \infty\right)\)

solveset(sin(x) - sqrt(3)/2, domain=S.Reals)

\(\displaystyle \left\{2 n \pi + \frac{2 \pi}{3}\; \middle|\; n \in \mathbb{Z}\right\} \cup \left\{2 n \pi + \frac{\pi}{3}\; \middle|\; n \in \mathbb{Z}\right\}\)

微分方程式

変数分離型

eq = Eq(Derivative(y(x)), g(x)*g(y(x)))
display(eq, dsolve(eq, y(x)))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = g{\left (x \right )} g{\left (y{\left (x \right )} \right )}\)

\(\displaystyle \int^{y{\left (x \right )}} \frac{1}{g{\left (y \right )}}\, dy = C_{1} + \int g{\left (x \right )}\, dx\)

例題: ロジスティック方程式

# シンボルを定義
t, r, K = symbols("t r K")
N = symbols("N", cls=Function)
# ロジスティック方程式を定義
eq = Eq(Derivative(N(t)), r*(1-(N(t)/K))*N(t))
display(eq)
# 一般解を求める
ans = dsolve(eq)
display(ans)

# 特殊化する
args = [(K, 100), (r, 0.75)]
init_val = 1.0
func = ans.rhs.subs(args)
# 初期値を方程式の解として求める
c = solve(Eq(func.subs([(t, 0)]), init_val), domain=S.Reals)
# プロットする
C1 = symbols("C1") # 積分定数用のシンボルを定義しなきゃいけないのがめんどくさい
plot(func.subs([(C1, c[0])]), (t, -10, 15))

\(\displaystyle \frac{d}{d t} N{\left (t \right )} = r \left(1 - \frac{1}{K} N{\left (t \right )}\right) N{\left (t \right )}\)

\(\displaystyle N{\left (t \right )} = \frac{K e^{C_{1} K + r t}}{e^{C_{1} K + r t} - 1}\)

png

<sympy.plotting.plot.Plot at 0x23613929b70>

LiouVille(リウビル)型

LiouVille = Eq(y(x).diff(x,2) + P(x)*y(x).diff(x) + Q(y(x))*y(x).diff(x)**2)
display(LiouVille, dsolve(LiouVille, y(x), hint="Liouville"))

\(\displaystyle P{\left (x \right )} \frac{d}{d x} y{\left (x \right )} + Q{\left (y{\left (x \right )} \right )} \frac{d}{d x} y{\left (x \right )}^{2} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = 0\)

\(\displaystyle C_{1} + C_{2} \int e^{- \int P{\left (x \right )}\, dx}\, dx + \int^{y{\left (x \right )}} e^{\int Q{\left (y \right )}\, dy}\, dy = 0\)

例題

eq = Eq(y(x).diff(x,2) + 2*x*y(x).diff(x) + 2*y(x)*y(x).diff(x) ** 2)
ans = dsolve(eq, hint="Liouville")
display(eq, ans, constsimp(expand_erfs(ans), x))

\(\displaystyle 2 x \frac{d}{d x} y{\left (x \right )} + 2 y{\left (x \right )} \frac{d}{d x} y{\left (x \right )}^{2} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = 0\)

\(\displaystyle C_{1} + C_{2} \operatorname{erf}{\left (x \right )} + \frac{\sqrt{\pi}}{2} \operatorname{erfi}{\left (y{\left (x \right )} \right )} = 0\)

\(\displaystyle C_{1} + C_{2} \int_{0}^{x} e^{- t^{2}}\, dt + \int_{0}^{y{\left (x \right )}} e^{t^{2}}\, dt = 0\)

eq = Eq(y(x).diff(x), sin(x)*tan(y(x)))
display(eq, dsolve(eq, y(x)))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = \sin{\left (x \right )} \tan{\left (y{\left (x \right )} \right )}\)

\(\displaystyle \left [ y{\left (x \right )} = - \operatorname{a{}c{}o{}s}{\left (- \sqrt{C_{1} e^{- 2 \cos{\left (x \right )}} + 1} \right )} + 2 \pi, \quad y{\left (x \right )} = - \operatorname{a{}c{}o{}s}{\left (\sqrt{C_{1} e^{- 2 \cos{\left (x \right )}} + 1} \right )} + 2 \pi, \quad y{\left (x \right )} = \operatorname{a{}c{}o{}s}{\left (- \sqrt{C_{1} e^{- 2 \cos{\left (x \right )}} + 1} \right )}, \quad y{\left (x \right )} = \operatorname{a{}c{}o{}s}{\left (\sqrt{C_{1} e^{- 2 \cos{\left (x \right )}} + 1} \right )}\right ]\)

同次型

eq = Eq(y(x).diff(x), f(y(x)/x))
display(eq, dsolve(eq, y(x)))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = f{\left (\frac{1}{x} y{\left (x \right )} \right )}\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- \int^{\frac{x}{y{\left (x \right )}}} \frac{f{\left (\frac{1}{u_{2}} \right )}}{u_{2} f{\left (\frac{1}{u_{2}} \right )} - 1}\, du_{2}}\)

eq = Eq(y(x).diff(x) , exp(y(x)/x) + y(x)/x)
display(eq, dsolve(eq))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = e^{\frac{1}{x} y{\left (x \right )}} + \frac{1}{x} y{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = \log{\left (\left(\frac{1}{C_{1} - \log{\left (x \right )}}\right)^{x} \right )}\)

分数型

eq = Eq(y(x).diff(x), (x-2*y(x)+3)/(2*x+y(x)-4))
display(eq, dsolve(eq))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = \frac{x - 2 y{\left (x \right )} + 3}{2 x + y{\left (x \right )} - 4}\)

\(\displaystyle \left [ y{\left (x \right )} = - 2 x - \sqrt{C_{1} + 5 x^{2} - 10 x} + 4, \quad y{\left (x \right )} = - 2 x + \sqrt{C_{1} + 5 x^{2} - 10 x} + 4\right ]\)

eq = Eq(y(x).diff(x), x*y(x)/(x**2+y(x)))
display(eq, dsolve(eq, hint="lie_group"))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = \frac{x y{\left (x \right )}}{x^{2} + y{\left (x \right )}}\)

\(\displaystyle \left [ y{\left (x \right )} = C_{1} \left(- \sqrt{- 2 C_{1} x^{2} + 1} - 1\right), \quad y{\left (x \right )} = C_{1} \left(\sqrt{- 2 C_{1} x^{2} + 1} - 1\right)\right ]\)

# 自動で解けないのもあるので、そのときは答を代入する
eq = Eq(y(x).diff(x), x**2*y(x)/(x**3+y(x)))
# この答えは手で解いたもの
ans_implicit = Eq(y(x)**3/(2*x**3+3*y(x)), 6*C1)
# 実数解だけ取る(虚数解でも成り立つだろうけどとりあえず無視)
ans = solve(ans_implicit, y(x))[0]
# これがTrueなのでans_implicitが正しいと分かる
display(eq, ans_implicit, Eq(y(x), ans), eq.subs([(y(x), ans)]).doit().simplify())

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = \frac{x^{2} y{\left (x \right )}}{x^{3} + y{\left (x \right )}}\)

\(\displaystyle \frac{y^{3}{\left (x \right )}}{2 x^{3} + 3 y{\left (x \right )}} = 6 C_{1}\)

\(\displaystyle y{\left (x \right )} = - \frac{18 C_{1}}{\sqrt[3]{- 162 C_{1} x^{3} + \frac{1}{2} \sqrt{- 629856 C_{1}^{3} + 104976 C_{1}^{2} x^{6}}}} - \frac{1}{3} \sqrt[3]{- 162 C_{1} x^{3} + \frac{1}{2} \sqrt{- 629856 C_{1}^{3} + 104976 C_{1}^{2} x^{6}}}\)

\(\displaystyle \mathrm{True}\)

定数変化法

eq = Eq(y(x).diff(x) + P(x)*y(x), Q(x))
display(eq, dsolve(eq, y(x), hint="almost_linear"))

\(\displaystyle P{\left (x \right )} y{\left (x \right )} + \frac{d}{d x} y{\left (x \right )} = Q{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + \int Q{\left (x \right )} e^{\int P{\left (x \right )}\, dx}\, dx\right) e^{- \int P{\left (x \right )}\, dx}\)

eq = Eq(x*y(x).diff(x) - 2*y(x), x**3*cos(x))
display(eq, dsolve(eq, y(x)))

\(\displaystyle x \frac{d}{d x} y{\left (x \right )} - 2 y{\left (x \right )} = x^{3} \cos{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = x^{2} \left(C_{1} + \sin{\left (x \right )}\right)\)

# 一般解から解いてみる
eq = Eq(y(x).diff(x) + P(x)*y(x), Q(x))
display(eq.subs([(P(x), -2/x), (Q(x), x**2*cos(x))]),
        dsolve(eq, hint="almost_linear").subs([(P(x), -2/x), (Q(x), x**2*cos(x))]).doit())

\(\displaystyle \frac{d}{d x} y{\left (x \right )} - \frac{2}{x} y{\left (x \right )} = x^{2} \cos{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = x^{2} \left(C_{1} + \sin{\left (x \right )}\right)\)

Bernoulli(ベルヌーイ)型

n = symbols("n")
eq = Eq(y(x).diff(x) + P(x)*y(x), Q(x)*y(x)**n)
display(eq, dsolve(eq, hint="Bernoulli"))

\(\displaystyle P{\left (x \right )} y{\left (x \right )} + \frac{d}{d x} y{\left (x \right )} = Q{\left (x \right )} y^{n}{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = \left(\left(C_{1} - \left(n - 1\right) \int Q{\left (x \right )} \left(e^{- n \int P{\left (x \right )}\, dx}\right) e^{\int P{\left (x \right )}\, dx}\, dx\right) e^{- \left(- n + 1\right) \int P{\left (x \right )}\, dx}\right)^{\frac{1}{- n + 1}}\)

eq = Eq(y(x).diff(x)+x*y(x), x*y(x)**3)
display(eq, dsolve(eq, hint="Bernoulli"))

\(\displaystyle x y{\left (x \right )} + \frac{d}{d x} y{\left (x \right )} = x y^{3}{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = \frac{1}{\sqrt{\left(C_{1} + e^{- x^{2}}\right) e^{x^{2}}}}\)

Riccati(リッカチ)型

# 自動で解くのは無理っぽい
eq = Eq(y(x).diff(x) - y(x)**2 - (2-x)*y(x) + 2*x - 1)
spec_ans = x
display(eq.subs([(y(x), spec_ans)]).doit().simplify())
bern = eq.subs([(y(x), P(x)+spec_ans)]).doit().simplify()
ans = dsolve(bern, hint="Bernoulli")
display(eq, bern, ans, ans.subs([(P(x), y(x)-x)]))

\(\displaystyle \mathrm{True}\)

\(\displaystyle 2 x - \left(- x + 2\right) y{\left (x \right )} - y^{2}{\left (x \right )} + \frac{d}{d x} y{\left (x \right )} - 1 = 0\)

\(\displaystyle - x P{\left (x \right )} - P^{2}{\left (x \right )} - 2 P{\left (x \right )} + \frac{d}{d x} P{\left (x \right )} = 0\)

\(\displaystyle P{\left (x \right )} = \frac{e^{x \left(\frac{x}{2} + 2\right)}}{C_{1} - \frac{\sqrt{2} \sqrt{\pi}}{2 e^{2}} \operatorname{erfi}{\left (\sqrt{2} \left(\frac{x}{2} + 1\right) \right )}}\)

\(\displaystyle - x + y{\left (x \right )} = \frac{e^{x \left(\frac{x}{2} + 2\right)}}{C_{1} - \frac{\sqrt{2} \sqrt{\pi}}{2 e^{2}} \operatorname{erfi}{\left (\sqrt{2} \left(\frac{x}{2} + 1\right) \right )}}\)

Clairaur(クレロー)型

eq = Eq(y(x), x*y(x).diff(x) - y(x).diff(x)**2)
dsolve(eq)

\(\displaystyle y{\left (x \right )} = \frac{x^{2}}{4} - \frac{1}{4} \left(C_{1} + x\right)^{2}\)

減衰振動

m, k, u = symbols("m k u")
eq = Eq(m*y(x).diff(x,2), -k*y(x) - u * y(x).diff(x))
dsolve(eq)

\(\displaystyle y{\left (x \right )} = C_{1} e^{\frac{x}{2 m} \left(- u - \sqrt{- 4 k m + u^{2}}\right)} + C_{2} e^{\frac{x}{2 m} \left(- u + \sqrt{- 4 k m + u^{2}}\right)}\)

定数係数斉次線形

# 一番簡単なやつ
eq = Eq(y(x).diff(x,2) - y(x).diff(x) - 6*y(x))
display(eq, dsolve(eq))

\(\displaystyle - 6 y{\left (x \right )} - \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = 0\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- 2 x} + C_{2} e^{3 x}\)

# 重解のパターン
eq = Eq(y(x).diff(x, 2) - 6*y(x).diff(x) + 9*y(x))
display(eq, dsolve(eq))

\(\displaystyle 9 y{\left (x \right )} - 6 \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = 0\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + C_{2} x\right) e^{3 x}\)

eq = Eq(y(x).diff(x,3) - 5*y(x).diff(x,2) + 8*y(x).diff(x) - 4*y(x))
display(eq, dsolve(eq))

\(\displaystyle - 4 y{\left (x \right )} + 8 \frac{d}{d x} y{\left (x \right )} - 5 \frac{d^{2}}{d x^{2}} y{\left (x \right )} + \frac{d^{3}}{d x^{3}} y{\left (x \right )} = 0\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + \left(C_{2} + C_{3} x\right) e^{x}\right) e^{x}\)

eq = Eq(y(x).diff(x,4) - 8*y(x).diff(x,2) + 16*y(x))
display(eq, dsolve(eq))

\(\displaystyle 16 y{\left (x \right )} - 8 \frac{d^{2}}{d x^{2}} y{\left (x \right )} + \frac{d^{4}}{d x^{4}} y{\left (x \right )} = 0\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + C_{2} x\right) e^{- 2 x} + \left(C_{3} + C_{4} x\right) e^{2 x}\)

定数係数非斉次線形

eq = Eq(y(x).diff(x, 2) + 3*y(x).diff(x) + 2*y(x), exp(3*x))
display(eq, dsolve(eq))

\(\displaystyle 2 y{\left (x \right )} + 3 \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = e^{3 x}\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- 2 x} + C_{2} e^{- x} + \frac{e^{3 x}}{20}\)

eq = Eq(y(x).diff(x, 2) + 2*y(x).diff(x) + y(x), exp(-x))
display(eq, dsolve(eq))

\(\displaystyle y{\left (x \right )} + 2 \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = e^{- x}\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + C_{2} x + \frac{x^{2}}{2}\right) e^{- x}\)

eq = Eq(y(x).diff(x, 2) + 3*y(x).diff(x) + 2*y(x), x)
display(eq, dsolve(eq))

\(\displaystyle 2 y{\left (x \right )} + 3 \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = x\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- 2 x} + C_{2} e^{- x} + \frac{x}{2} - \frac{3}{4}\)

eq = Eq(y(x).diff(x, 2) - y(x).diff(x) - 2*y(x), x**2)
display(eq, dsolve(eq))

\(\displaystyle - 2 y{\left (x \right )} - \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = x^{2}\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- x} + C_{2} e^{2 x} - \frac{x^{2}}{2} + \frac{x}{2} - \frac{3}{4}\)

eq = Eq(y(x).diff(x, 2) + 2*a*y(x).diff(x) + a**2*y(x), x)
display(eq, dsolve(eq))

\(\displaystyle a^{2} y{\left (x \right )} + 2 a \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = x\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + C_{2} x\right) e^{- a x} + \frac{x}{a^{2}} - \frac{2}{a^{3}}\)