常微分方程式を解く手順
sympy.solvers.ode モジュールをインポートすると、次のような種類の常微分方程式を解けるようになります。
・1階変数分離型微分方程式
・1階斉次微分方程式
・1階完全微分方程式
・1階線形微分方程式
・n階定数係数微分方程式
・ベルヌイ型微分方程式
・スツルム=リウヴィル型微分方程式
いくつかの例を使って、SymPy による常微分方程式の基本的な解法手順を説明します。コードは Jupyter Notebook環境を前提に記述されています。
変数分離型微分方程式
最初に
# SYMPY_ODE
# In[1]
# 必要な関数やクラスをインポート
from sympy import var, init_printing, Function
from sympy.solvers import ode
# 数式を MathJax で表示するように設定
init_printing()
数式を TeX で表示する必要がない場合は、init_printing() の一行を削除してください。
sympy.solvers.ode.classify_ode() に微分方程式を渡すと、どの型に分類される方程式なのか表示されます。
# In[2]
# 変数xを定義
var('x')
# yを関数として定義
var('y', cls=Function)
# y'-x(1-y)=0
eq_1 = y(x).diff(x) - x * (1 - y(x))
# 微分方程式を分類
ode.classify_ode(eq_1)
'''
('separable',
'1st_exact',
'1st_linear',
'Bernoulli',
'almost_linear',
'1st_power_series',
'lie_group',
'separable_Integral',
'1st_exact_Integral',
'1st_linear_Integral',
'Bernoulli_Integral',
'almost_linear_Integral')
'''
出力された分類一覧を見ると、
・変数分離型 (separable)
・1階完全 (1st_exact)
・1階線形 (1st_linear)
・ベルヌーイ型 (Bernoulli)
であることがわかります。この分類キーワードを sympy.solvers.ode.dsolve() の hint に渡すと、それぞれのアルゴリズムにしたがって微分方程式を解いてくれますが、一覧の上にあるほど、Sympy が推奨する(つまり効率的に速く解くことができる)キーワードとなっています。したがって、今の場合は hint に ‘separable’ を渡せばよいことになります。
# In[3]
# 微分方程式y'-x(1-y)=0を解く
ode.dsolve(eq_1, hint = "separable")
ics に初期条件を渡して解を得ることもできます。
# In[4]
# 微分方程式y'-x(1-y)=0を変数分離型で解く
ode.dsolve(eq_1, hint="separable", ics={y(0):2})
hint に ‘1st_power_series’ を渡すと級数解を得ます。
# In[5]
# y'-x(1-y)=0を解いて級数解を取得
ode.dsolve(eq_1, hint="1st_power_series", ics={y(0):2})
2階斉次方程式
次は 2階斉次方程式
# 2ND_HOMOGENEOUS_DIFFERENTIAL_EQUATION
# In[1]
from sympy import var, init_printing, Function
from sympy.solvers import ode
# MathJaxで数式を表示
init_printing()
# 記号x,p,qを定義
var('x p q')
# yを関数として定義
var('y', cls=Function)
# 微分方程式y'' + qy = 0
eq_2 = y(x).diff(x, 2) + p * y(x).diff(x) + q * y(x)
# 微分方程式を分類
ode.classify_ode(eq_2)
# ('nth_linear_constant_coeff_homogeneous',
# '2nd_power_series_ordinary')
nth_linear_constant_coeff_homogeneous (定数係数n階斉次方程式) に分類されていることがわかったので、このキーワードを sympy.solvers.ode.dsolve() の hint に渡して微分方程式を解きます。一般解には 2 つの定数
# In[2]
# 二階斉次微分方程式y''+qy=0を変数分離型で解く
ode.dsolve(eq_2,
hint="nth_linear_constant_coeff_homogeneous")
初期条件
# In[3]
# y''+qy=0を初期条件y(0)=1,y'(0)=0で解く
ode.dsolve(eq_2,
hint = "nth_linear_constant_coeff_homogeneous",
ics = {y(0):1, y(x).diff(x, 1).subs(x, 0):0})
ベルヌーイ型微分方程式
次の形の方程式をベルヌーイの微分方程式とよびます。
という方程式を解いてみましょう。すでにベルヌーイ型とわかっていますが、念のために分類を確認しておきます。
# BERNOULLI'S_DIFFERENTIAL_EQUATION
# In[1]
# SymPyから必要な関数やクラスをインポート
from sympy import var, exp, sin, init_printing, Function, simplify
from sympy.solvers import ode
# MathJaxで数式を表示
init_printing()
# 記号x,p,qを定義
var('x p q')
# yを関数として定義
var('y', cls=Function)
# ベルヌーイの微分方程式
eq_3 = y(x).diff(x, 1) + x * y(x) + exp(x) * y(x) ** 2
# 微分方程式を分類
ode.classify_ode(eq_3)
# ('Bernoulli', '1st_power_series', 'lie_group', 'Bernoulli_Integral')
分類リストの最上位に ‘Bernoulli’ が表示されています。
これを dsolve() の hint に与えて微分方程式を解きます。
# In[2]
# 微分方程式y''+xy+exp(x)y**2=0の解
f = ode.dsolve(eq_3, hint = "Bernoulli")
# 式を簡略化
simplify(f)
解の分母には積分が含まれています。
この部分は sympy.integrate() を使って計算することができます。
# In[3]
# 積分関数を読み込む
from sympy import integrate
# 被積分関数
func = exp(x - x**2/2)
# funcをxで積分
z = integrate(func, x)
# 式を整理
simplify(z)
コメント
関数のシンボルは
from sympy import Function
sympy.var( ‘y’, cls=Function))
あるいは
y=sympy.Function(‘y’)
で定義しないとエラーが発生しますので、ご確認をお願いいたします。
申し訳ないです。
ODE_01-1 と ODE_01-2 を修正しておきました。
ODE_01-1 で Function クラスをインポートし、ODE_01-2 で ‘y’ を関数として定義しています。
ODE_02-1 と ODE_03-1 も修正しました。
下記は誤植と思われますので、ご確認ください。
「ベルヌーイ型微分方程式」の説明で、
y” +p(x)y + q(x)y^k = 0 → y’ + p(x)y = q(x)y^k
y” + xy + e^xy^2 = 0 → y’ + xy = e^xy^2
ODE_03-1 プログラムで、
eq_3 = y(x).diff(x, 1) + x * y(x) + exp(x) * y(x)**2 → eq_3 = y(x).diff(x, 1) + x * y(x) – exp(x) * y(x)**2
ありがとうございます。1 階微分のところが 2 階微分になっていましたね。修正しておきました。教科書によってベルヌーイ方程式の表記が異なりますが、q(x) が任意関数なので、y’ + p(x)y + q(x)y^k = 0 と y’ + p(x)y = q(x)y^k は同じ方程式です。
理由はわかりませんが、ODE_03-1 プログラムの実行結果は、
(‘1st_power_series’, ‘lie_group’)
となり、ODE_03-2 プログラムの実行結果は下記のようなエラーとなりました。
ODE x*y(x) + y(x)**2*exp(x) + Derivative(y(x), x) does not match hint Bernoulli
exp(x) * y(x)**2 項の符号を – にすると、両者ともに記事と同じ実行結果が得られました。
※ODE_03-2プログラムの実行結果の積分の中の式は、少し整理が不完全なe^xe^(-x^2/2)となっていました。
Google Colab にインストールされている SymPy のバージョンを確認すると 1.7.1 となっていました。確かに、この環境で ODE_03 を実行すると不具合が生じます。以下のコマンドで SymPy をアップグレードすると正常に動作することが確認できました。
pip install sympy==1.11
アップグレードした後はランタイムを再起動してください。現在の SymPy の公式ドキュメントも 1.11 に準拠しているので、このバージョンを使ったほうがいいと思います。
sympy を 1.11 にアップグレードして正常に動作することを確認いたしました。
ありがとうございました。