Python（SymPy）で微分方程式を解く（力学の問題を題材にして）

物理学の具体的な計算にPython（SymPy）を使って、物理学もPython（SymPy）も同時に学んでしまいましょう。今回はPythonの数式処理ライブラリであるSymPyを使って微分方程式を解きたいと思います。

一般に、未知関数とその導関数を含んだ方程式を微分方程式といいます。一変数関数の場合を常微分方程式、多変数関数の場合を偏微分方程式と呼びます。今回は力学で現れる具体的な常微分方程式をいくつか解いてみます。

Python（SymPy）のごく基本的な使い方については以下の記事を参照してください：

pianofisica.hatenablog.com

不定積分
微分方程式
２階線型常微分方程式：自由落下
１次元調和振動子
Airy（エアリー）方程式

不定 積分

まず微分方程式を解くにあたって関数の積分を復習しておきましょう。

関数 $y=f(x)$ の不定積分

$\qquad\displaystyle{F(x)=\int dx \, f(x) +C}$

（ $C$ は積分定数）は、 $\ f(x)=\exp(-ax)$ を例にとるとPythonでは

import sympy as sp
sp.init_printing()
sp.var('a, x')
sp.integrate(sp.exp(-a*x), x)

などとします（ただし積分定数は省かれてしまいます）。

少しの間、この不定積分を出力する関数を I と定義してこれを用います：

def I(f): return sp.integrate(f, x)

上の例では

I(sp.exp(-a*x))

より

$\qquad\displaystyle{\int dx \,e^{-ax} =\begin{cases} \ -\frac{e^{-ax}}{a}+C &\qquad (a\ne0)\\ \ x +C &\qquad (a=0)\end{cases}}$

です（ただし積分定数は省かれています）。他にも

sp.var('n')
I(x**n)

として $x^n$ の不定積分

$\qquad\displaystyle{\int dx \, x^n =\begin{cases} \ \frac{x^{n+1}}{n+1}+C &\qquad (n\neq-1)\\ \ \log x +C &\qquad (n=-1)\end{cases}}$

が求まりますし、また

I(sp.cos(n*x))

とすれば $\cos(nx)$ の不定積分

$\qquad\displaystyle{\int dx \,\cos(nx) =\begin{cases} \ \frac{\sin(nx)}{n}+C &\qquad (n\neq0)\\ \ x+C &\qquad (n=0)\end{cases}}$

が求まります。さてウォーミングアップは軽く済ませて、早速、微分方程式を解いていきましょう。
　

微分方程式

１階線型常微分方程式

簡単な例として次の１階線型常微分方程式

$\qquad\displaystyle{\frac{dy}{dx}=-ay}$

を考えましょう。ここで $y$ は $x$ の関数、 $a$ は定数とします。

この微分方程式は変数分離法によって解くことができます：

$\qquad\displaystyle{\int dx\, \frac{dy}{dx}\,\frac{1}{y}=-a\int dx}$

より

$\qquad\displaystyle{\log y=-ax + c \quad (c:\text{積分定数})}$

両辺の指数関数をとって

$\qquad\displaystyle{y=C\exp\left(-ax\right) \quad (C:\text{積分定数})}$

と解くことができます。
　

dsolve を使って微分方程式を解く

さて、これをPythonで解かせるには、微分方程式を解く実装 dsolve を使って

import sympy as sp
sp.init_printing()
sp.var('a, x')
y = sp.Function('y')(x)
eq = sp.Eq( sp.diff(y, x),  -a*y )
sp.dsolve(eq)

として

$\qquad\displaystyle{y(x)=Ce^{-ax} \quad (C:\text{積分定数})}$

がわかります（不定積分の場合と違って積分定数も含めて結果を出力します）。

初期値を課して微分方程式を解く

実装 dsolveで初期条件を課した場合に微分方程式の解を求めるには、

たとえば上の例で $y(0)=1$ とするとき

sp.dsolve(eq, ics={y.subs(x,0):1})

から

$\qquad\displaystyle{y(x)=e^{-ax}}$

が求まります。
　

２階線型常微分方程式：自由落下

１次元自由落下は２階常微分方程式（運動方程式）

$\qquad\displaystyle{m\frac{d^2x}{dt^2}=-mg}$

で記述されます。

ここで $x \ [{\rm m}]$ は鉛直方向の座標（上を正の方向にとりました）、

また、 $m \ [{\rm kg}]$ は質点の質量、 $g \ [{\rm m/s^2}]$ は重力加速度です。

この常微分方程式をPythonで解くには

import sympy as sp
sp.init_printing()
sp.var('g, t')
x = sp.Function('x')(t)
eq1 = sp.Eq( sp.diff(x,t,2),  -g )
sp.dsolve(eq1)

とすることで

$\qquad\displaystyle{x(t)=-\frac12 g t^2+C_2 t+C_1}$

がわかります。ここで $C_2$ は初速度、 $C_1$ は初期位置に相当する積分定数です。

さて、初期条件 $\ x(0)=0,\ \frac{dx}{dt}(0)=v$ を課すと

sp.var('v')
sp.dsolve(eq1, ics={x.subs(t,0):0, sp.diff(x,t,1).subs(t,0):v})

によって

$\qquad\displaystyle{x(t)=-\frac12 g t^2+v t}$

が出力されます。
　　

１次元調和振動子

水平右向きを $x$ 軸の正方向にとります。運動方程式は

$\qquad\displaystyle{m\frac{d^2x}{dt^2}=F_x \qquad F_x=-kx}$

です。この状況設定をPythonに入力して解を求めるには

import sympy as sp
sp.init_printing()
sp.var('m,k, t')
x = sp.Function('x')(t)
F = -k*x
eq2 = sp.Eq( m*sp.diff(x,t,2),  F )
sp.dsolve(eq2)

として、一般解が

$\qquad\displaystyle{x(t)=C_{1} e^{-t \sqrt{-\frac{k}{m}}}+C_{2} e^{t \sqrt{-\frac{k}{m}}}}$

とわかります。ここで $C_1, \ C_2$ は積分定数です。

いまバネ定数 $k \ {\rm [ N/m]}$ 、質量 $m\ {\rm [ m]}$ はともに正の値なので

$\qquad\displaystyle{i\omega=\sqrt{-\frac{k}{m}}}$

としてEulerの公式 $\ e^{i\theta}=\cos\theta+i\sin\theta$ を使うと、一般解が

$\qquad\displaystyle{x(t)=A \cos (\omega t)+B\sin(\omega t)}$

で与えられることがわかります。あるいはより具体的に

$\qquad\displaystyle{\frac{d^2x}{dt^2}=-x}$

という場合を考えると

import sympy as sp
sp.init_printing()
sp.var('t')
x = sp.Function('x')(t)
eq3 = sp.Eq( sp.diff(x,t,2),  -x)
sp.dsolve(eq3)

より、その一般解として

$\qquad\displaystyle{x(t)=C_1 \sin (t)+C_2\cos(t)}$

が出力されます。

裳華房の『解析学概論』（矢野・石原）は、常微分方程式（だけでなく、ベクトル解析、複素解析、フーリエ解析）に関してよくまとめられた解説と例題、さまざまな練習問題が載っており、効率よく応用数学が学べる理工系学部生におすすめのテキストです：

<br />

テキストに載っている問題を、手計算で解けるようになることはもとより、Pythonを利用して解いてみるのも、Pythonを有効に活用するための練習として有用です。
　

Airy（エアリー）方程式

もう少し複雑な微分方程式としてAiry方程式

$\qquad\displaystyle{\frac{d^2y}{dx^2}-xy=0}$

を考えてみます：Pythonに

import sympy as sp
sp.init_printing()
sp.var('x')
y = sp.Function('y')(x)
eq4 = sp.Eq( sp.diff(y,x,2), x*y)
sp.dsolve(eq4)

~~と入力するとその結果として~~

~~$\qquad\displaystyle{ y(x)=C_{2}\left(\frac{x^{3}}{6}+1\right)+C_{1} x\left(\frac{x^{3}}{12}+1\right)+O\left(x^{6}\right)}$~~

~~が出力されるだけです。本来であれば~~

~~$\qquad\displaystyle{y(x)=C_{1}{\rm Ai}(x)+C_{2}{\rm Bi}(x)}$~~

~~が出力されるべきですが、特殊関数をともなう計算にはSymPyでは限界があるようです。~~

（追記）SymPyのバージョンが古かったようです。最新版のSymPyにアップデートしたら正しく

$\qquad\displaystyle{y(x)=C_{1}{\rm Ai}(x)+C_{2}{\rm Bi}(x)}$

と出力されました。
　

キーワード：微分方程式、運動方程式、Python、自由落下、調和振動子