【SymPy】常微分方程式

常微分方程式を解く手順

sympy.solvers.ode モジュールをインポートすると、次のような種類の常微分方程式を解けるようになります。
 
・1階変数分離型微分方程式
・1階斉次微分方程式
・1階完全微分方程式
・1階線形微分方程式
・n階定数係数微分方程式
・ベルヌイ型微分方程式
・スツルム＝リウヴィル型微分方程式　
 
いくつかの例を使って、SymPy による常微分方程式の基本的な解法手順を説明します。コードは Jupyter Notebook環境を前提に記述されています。

変数分離型微分方程式

最初に $y’=x(1-y)$ を解いてみます。まず必要なモジュールをまとめてインポートして、実行結果を TeX 形式で表示するように設定しておきます、

# 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')
'''

出力された分類一覧を見ると、$y’=x(1-y)$ は
 
・変数分離型 (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})

$O\left(x^{6}\right)$ は $x^6$ のオーダーの項であることを意味しています。

2階斉次方程式

次は 2階斉次方程式 $y”+py’+qy=0$ を解いてみます。最初にSympy における分類を確認しておきます。

# 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 つの定数 $C_1,\ C_2$ が含まれます。

# In[2]

# 二階斉次微分方程式y''+qy=0を変数分離型で解く
ode.dsolve(eq_2,
           hint="nth_linear_constant_coeff_homogeneous")

初期条件 $y(0)=1,\ y'(0)=0$ によって $C_1$ と $C_2$ を決定すると、次のような解を得ます。

# 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})

ベルヌーイ型微分方程式

次の形の方程式をベルヌーイの微分方程式とよびます。
 \[y’+p(x)y+q(x)y^k=0\]
$p(x)=x,\ q(x)=e^x,\ k=2$ として
 \[y’+xy+e^xy^2=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)

$\operatorname{erf}(x)$ は次式で定義される誤差関数です。
 \[\mathrm{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}\exp (-t^2)dt\]