つれづれなる備忘録

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

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ツールを用いて、非線形微分方程式と連立微分方程式の解き方について紹介した。

「凄麺 冷し中華 海藻サラダ風」

近くのドラッグストアで売っていた冷やし中華カップ麺。「凄麺 冷し中華 海藻サラダ風」

https://www.newtouch.co.jp/ec/upload/save_image/0407101807_606d084f289ae.jpg

www.newtouch.co.jp

焼きそばと同じ要領で、お湯で戻して湯切りするが最後に冷水を入れて水切りすることで冷やし中華になる。フライ麺でつくる冷やし中華は袋麺では中華三昧シリーズが有名だが、カップ麺では珍しい。

かやくもカップ麺としては珍しい白きくらげ、わかめ、赤つのまた。

食べた感想としては、麺は冷水でしめることで、もちもち感が強く、かやくもインスタントとは思えないほど弾力が強い。特に何も言われなければカップ麺とは思えないほど、冷やし中華を再現している。

カップヌードル旨辛カルビ味焼きそば

カップヌードルシリーズの新製品:カップヌードル旨辛カルビ味焼きそばを食べてみた。

www.nissin.com

https://cdn.nissin.com/image?id=10590&s=720

甘辛のソースをカップヌードルの麺に和えた味は想像通りのBBQ風味のカップ焼きそば。思ったよりも辛くはない。

麺の量に対してソースが足りず、麺全体を均等に和えることができなかった。

日清といえばUFO焼きそばがあるが、カップヌードルサイズの方がカバンに入れて持ち歩くには便利。