つれづれなる備忘録

日々の発見をあるがままに綴る

SymPyの使い方23 ~ ODE

 今回はSymPyを用いて微分方程式を解くODEツールの詳細について紹介したい。 SymPyを用いた微分方程式を解く方法については以前紹介したが、今回はもう少し詳しい使い方になる。

atatat.hatenablog.com

1. 微分方程式ソルバーdsolve

 SymPyのdsolveを用いることで微分方程式を解くことができる。例えば

 \displaystyle \frac{d^{2}f(x)}{dx^{2}} = 9 f(x)

を解くには変数x,関数fを定義して以下のようにdsolveを使用すると解が得られる。

from sympy import *
init_printing(False)
x=symbols('x')
f = Function('f')
 dsolve(Derivative(f(x), x, x) + 9*f(x), f(x))
>Eq(f(x), C1*sin(3*x) + C2*cos(3*x))

2. 非線形微分方程式

線形微分方程式だけでなく、非線形微分方程式にも適用できる。(あくまで解が存在する場合) このときにdsolveのオプションhintをで解法を指定することで、計算が早く収束する場合がある。

以下hintで1st_exactまたはalmolst_linearGoogle Colab上で3秒程度の実行時間だった。

eq = sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x)
dsolve(eq, hint='1st_exact')
>[Eq(f(x), -acos(C1/cos(x)) + 2*pi), Eq(f(x), acos(C1/cos(x)))]
 dsolve(eq, hint='almost_linear')
>[Eq(f(x), -acos(C1/cos(x)) + 2*pi), Eq(f(x), acos(C1/cos(x)))]

以下hint設定なしでも解けるが、実行時間としては25秒だった。

dsolve(eq)
>[Eq(f(x), -acos(C1/cos(x)) + 2*pi), Eq(f(x), acos(C1/cos(x)))]

いずれも解は同じだが実行時間がhintで指定することで短縮できた。

hintで設定できるキーワードは多数あり、以下のようにallhintsにより調べることができる。

from sympy.solvers.ode.ode import allhints
allhints
>('factorable',
 'nth_algebraic',
 'separable',
 '1st_exact',
 '1st_linear',
 'Bernoulli',
 'Riccati_special_minus2',
 '1st_homogeneous_coeff_best',
 '1st_homogeneous_coeff_subs_indep_div_dep',
 '1st_homogeneous_coeff_subs_dep_div_indep',
 'almost_linear',
 'linear_coefficients',
 'separable_reduced',
 '1st_power_series',
 'lie_group',
 'nth_linear_constant_coeff_homogeneous',
 'nth_linear_euler_eq_homogeneous',
 'nth_linear_constant_coeff_undetermined_coefficients',
 'nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients',
 'nth_linear_constant_coeff_variation_of_parameters',
 'nth_linear_euler_eq_nonhomogeneous_variation_of_parameters',
 'Liouville',
 '2nd_linear_airy',
 '2nd_linear_bessel',
 '2nd_hypergeometric',
 '2nd_hypergeometric_Integral',
 'nth_order_reducible',
 '2nd_power_series_ordinary',
 '2nd_power_series_regular',
 'nth_algebraic_Integral',
 'separable_Integral',
 '1st_exact_Integral',
 '1st_linear_Integral',
 'Bernoulli_Integral',
 '1st_homogeneous_coeff_subs_indep_div_dep_Integral',
 '1st_homogeneous_coeff_subs_dep_div_indep_Integral',
 'almost_linear_Integral',
 'linear_coefficients_Integral',
 'separable_reduced_Integral',
 'nth_linear_constant_coeff_variation_of_parameters_Integral',
 'nth_linear_euler_eq_nonhomogeneous_variation_of_parameters_Integral',
 'Liouville_Integral')

hintで1st_power_seriesを指定すると、名前の通り多項式1次近似なので1st_exactなどとは解が異なる。

dsolve(eq,hint='1st_power_series')
>Eq(f(x), -x**2*cos(C1)/(2*sin(C1)) + x**4*(-5 - 3*cos(C1)**2/sin(C1)**2)*cos(C1)/(24*sin(C1)) + C1 + O(x**6))

3. 連立微分方程式

 連立微分方程式の場合はsympy.solvers.ode.systemsからdsolve_systemをロードして使用する。

例えば

\displaystyle \begin{eqnarray}
       \frac{df(x)}{dx} &=& g(x) \\
       \frac{dg(x)}{dx} &=& f(x)   \end{eqnarray}

は以下のように2本の微分方程式を配列eqsとして定義する。

f, g = symbols("f g", cls=Function)
eqs = [Eq(f(x).diff(x), g(x)), Eq(g(x).diff(x), f(x))]

dsolveと同じようにdsolve_system(eqs)とすれば解くことができる。なおdsolveのようにhintオプションはない。

from sympy.solvers.ode.systems import dsolve_system
dsolve_system(eqs)
>[[Eq(f(x), -C1*exp(-x) + C2*exp(x)), Eq(g(x), C1*exp(-x) + C2*exp(x))]]

初期値を指定する場合はdictionary形式でics={f(0): 1, g(0): 0}とする。

dsolve_system(eqs,ics={f(0): 1, g(0): 0})
>[[Eq(f(x), exp(x)/2 + exp(-x)/2), Eq(g(x), exp(x)/2 - exp(-x)/2)]

もう少し複雑な連立微分方程式で、以下解法を紹介したものもdsolve_systemで解ける。

\displaystyle \begin{eqnarray}
       \frac{df(x)}{dx} &=& af(x)+bg(x) \\
       \frac{dg(x)}{dx} &=& cf(x)+dg(x)   \end{eqnarray}

atatat.hatenablog.com

a,b,c,d=symbols('a b c d')
eqs = [Eq(f(x).diff(x),a*f(x)+b*g(x)), Eq(g(x).diff(x), c*f(x)+d*g(x))]
dsolve_system(eqs)
>[[Eq(f(x), 2*C1*b*exp(x*(a + d + sqrt(a**2 - 2*a*d + 4*b*c + d**2))/2)/(-a + d + sqrt(a**2 - 2*a*d + 4*b*c + d**2)) - 2*C2*b*exp(x*(a + d - sqrt(a**2 - 2*a*d + 4*b*c + d**2))/2)/(a - d + sqrt(a**2 - 2*a*d + 4*b*c + d**2))), Eq(g(x), C1*exp(x*(a + d + sqrt(a**2 - 2*a*d + 4*b*c + d**2))/2) + C2*exp(x*(a + d - sqrt(a**2 - 2*a*d + 4*b*c + d**2))/2))]]

4. まとめ

 今回はSymPyのODEツールを用いて、非線形微分方程式と連立微分方程式の解き方について紹介した。