【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")

y (x) = C_{1} e^{- \frac{x^{2}}{2}} + 1

ics に初期条件を渡して解を得ることもできます。

# In[4]

# 微分方程式y'-x(1-y)=0を変数分離型で解く
ode.dsolve(eq_1, hint="separable", ics={y(0):2})

y (x) = 1 + e^{- \frac{x^{2}}{2}}

hint に ‘1st_power_series’ を渡すと級数解を得ます。

# In[5]

# y'-x(1-y)=0を解いて級数解を取得
ode.dsolve(eq_1, hint="1st_power_series", ics={y(0):2})

y (x) = 2 - \frac{x^{2}}{2} + \frac{x^{4}}{8} + O (x^{6})

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

2階斉次方程式

次は 2階斉次方程式 $y ” + p y^{'} + q y = 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 (x) = C_{1} e^{\frac{x (- p - \sqrt{p^{2} - 4 q})}{2}} + C_{2} e^{\frac{x (- p + \sqrt{p^{2} - 4 q})}{2}}

初期条件 $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 (x) = (- \frac{p}{2 \sqrt{p^{2} - 4 q}} + \frac{1}{2}) e^{\frac{x (- p - \sqrt{p^{2} - 4 q})}{2}} + (\frac{p}{2 \sqrt{p^{2} - 4 q}} + \frac{1}{2}) e^{\frac{x (- p + \sqrt{p^{2} - 4 q})}{2}}

ベルヌーイ型微分方程式

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

y (x) = \frac{e^{- \frac{x^{2}}{2}}}{C_{1} + \int e^{- \frac{x^{2}}{2} + x} d x}

解の分母には積分が含まれています。
この部分は sympy.integrate() を使って計算することができます。

# In[3]

# 積分関数を読み込む
from sympy import integrate

# 被積分関数
func = exp(x - x**2/2)

# funcをxで積分
z = integrate(func, x)

# 式を整理
simplify(z)

\frac{\sqrt{2} \sqrt{π} e^{\frac{1}{2}} \erf (\frac{\sqrt{2} (x - 1)}{2})}{2}

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

HNaito より:

2023年3月5日 10:27 PM

関数のシンボルは
from sympy import Function
sympy.var( ‘y’, cls=Function))
あるいは
y=sympy.Function(‘y’)
で定義しないとエラーが発生しますので、ご確認をお願いいたします。

返信
- あとりえこばとより:
  
  2023年3月6日 8:59 AM
  
  申し訳ないです。
  ODE_01-1 と ODE_01-2 を修正しておきました。
  ODE_01-1 で Function クラスをインポートし、ODE_01-2 で ‘y’ を関数として定義しています。
  ODE_02-1 と ODE_03-1 も修正しました。
  
  返信
HNaito より:

2023年3月6日 1:38 AM

下記は誤植と思われますので、ご確認ください。
「ベルヌーイ型微分方程式」の説明で、
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

返信
- あとりえこばとより:
  
  2023年3月6日 9:22 AM
  
  ありがとうございます。1 階微分のところが 2 階微分になっていましたね。修正しておきました。教科書によってベルヌーイ方程式の表記が異なりますが、q(x) が任意関数なので、y’ + p(x)y + q(x)y^k = 0 と y’ + p(x)y = q(x)y^k は同じ方程式です。
  
  返信
  - HNaito より:
    
    2023年3月6日 9:51 AM
    
    理由はわかりませんが、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)となっていました。
    
    返信
    - あとりえこばとより:
      
      2023年3月6日 11:18 AM
      
      Google Colab にインストールされている SymPy のバージョンを確認すると 1.7.1 となっていました。確かに、この環境で ODE_03 を実行すると不具合が生じます。以下のコマンドで SymPy をアップグレードすると正常に動作することが確認できました。
      pip install sympy==1.11
      アップグレードした後はランタイムを再起動してください。現在の SymPy の公式ドキュメントも 1.11 に準拠しているので、このバージョンを使ったほうがいいと思います。
      
      返信
      - HNaito より:
        
        2023年3月6日 5:47 PM
        
        sympy を 1.11 にアップグレードして正常に動作することを確認いたしました。
        ありがとうございました。
        
        返信