SymPyの使い方23 ～ ODE
今回はSymPyを用いて微分方程式を解くODEツールの詳細について紹介したい。 SymPyを用いた微分方程式を解く方法については以前紹介したが、今回はもう少し詳しい使い方になる。
1. 微分方程式ソルバーdsolve
SymPyのdsolveを用いることで微分方程式を解くことができる。例えば
を解くには変数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_linearでGoogle 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をロードして使用する。
例えば
は以下のように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で解ける。
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))]]
