SymPy で1階線形常微分方程式や2階定数係数線形常微分方程式を解析的に解く。
SymPy の import
from sympy.abc import *
from sympy import *
常微分方程式
未定義関数の宣言 Function()
$y$ を（変数ではなく）関数として宣言するには y = Function('y') とします。
y = Function('y')
上記のように関数として宣言すると y のタイプは UndefinedFunction 未定義関数になります。
type(y)
未定義関数として宣言された $y$ を変数 $x$ の関数 $y(x)$ とする場合は，y(x) と書きます。
y(x)
参考：複数の未定義関数の宣言
複数の未定義関数を宣言する際，以下のような書き方もできるようです。
# f, g, h = Function('f g h') などとは書けない。
#
# f = Function('f')
# g = Function('g')
# h = Function('h') と3行書くよりも簡単
var('f g h', cls = Function)
微分の表記 Derivative(), diff()
微分 $\displaystyle \frac{dy}{dx}$ を表すには Derivative(y(x), x) または diff(y(x), x)。未定義関数については diff(y(x), x) は Derivative(y(x), x) と同じです。
$\displaystyle \frac{dy}{dx} = $ Derivative(y(x), x) $=$ diff(y(x), x) $=$ y(x).diff(x)
diff(y(x), x)
# このような書き方もできる
y(x).diff(x)
# 未定義関数 y を t の関数としたいなら...
display(y(t))
display(y(t).diff(t))
微分方程式の定義 Eq()
例として微分方程式 $\displaystyle \frac{dy}{dx} = -2 x\, y$ の定義の仕方についてまとめます。
方程式とは未知変数や未知関数を含む等式のことですから，等式と同じように
Eq(左辺, 右辺) のように定義します。
# y を関数として宣言
y = Function('y')
# 解くべき微分方程式を定義して eq に代入
eq = Eq(diff(y(x), x), -2*x*y(x))
# eq に代入された方程式を表示
eq
微分方程式を解く dsolve()
SymPy で微分方程式を解くには dsolve() 関数を使います。
eq に代入された微分方程式
$\displaystyle \frac{dy}{dx} = -2 x\, y$ を未知関数 y(x) について解くには，以下のように書きます。
# dsolve() で微分方程式 eq を y(x) について解く
dsolve(eq, y(x))
上記の出力結果に現れる $C_1$ は積分定数です。
ちなみに，微分方程式が
$$\frac{dy}{dx} + 2 x y = 0$$
のように右辺がゼロの形になっている場合は，左辺のみを dsolve() の引数にしても同じ答えが出ます。
sahen = diff(y(x), x) + 2*x*y(x)
display(sahen)
dsolve(sahen, y(x))
初期条件 ics={}
初期条件を課して積分定数の値を決めることができます。
微分方程式を解く際の初期条件は dsolve() のオプションに ics={} のように「辞書」形式で設定します。
例として
$$\frac{dy}{dx} + 2 x y = 0$$
を $x = 0$ で $y(0) = 1$ という初期条件のもとで解く場合：
dsolve(sahen, y(x), ics={y(0):1})
この初期条件によって積分定数が $C_1 = 1$ と決まることがわかります。
1階線形常微分方程式
いくつか例題を解いてみます。
マルサスの人口モデル
ある時刻 $t$ における，とある国（地域）の総人口を $N(t)$ とすると，マルサスの人口モデルは以下の微分方程式で表されます。
$$\frac{dN}{dt} = \gamma\, N$$
この微分方程式を初期条件 $t=t_0$ で $N(t_0) = N_0$ のもとで解きます。
# N を関数として宣言
N = Function('N')
# gamma を変数として宣言（ガンマ関数とカブるため）
var('gamma')
# 微分方程式の定義
eqm = Eq(diff(N(t), t), gamma*N(t))
eqm
# 初期値変数の宣言
var('t0, N0')
# 初期条件は ics={} で以下のように...
solm = dsolve(eqm, N(t), ics={N(t0):N0})
solm
ヴェアフルストによる修正人口モデル
マルサスの人口モデルは，$\gamma > 0$ の場合には人口の際限なき増加（まさに指数関数的な増加）を予測する。しかし，実際にはさまざまな要因により，このような無制限な増加は続かない。
人口過密の要因を考慮にいれた修正を取り入れたヴェアフルストのモデルは，以下の微分方程式で表すことができる。
$$\frac{dN}{dt} = \gamma N \left(1 – \frac{N}{N_{\rm max}}\right)$$
この方程式を，同じ初期条件 $t=t_0$ で $N(t_0) = N_0$ のもとで解く。
N = Function('N')
var('gamma N_max')
eqv = Eq(diff(N(t), t), gamma*N(t)*(1-N(t)/N_max))
eqv
var('t0 N0')
ansv = dsolve(eqv, N(t), ics = {N(t0):N0})
ansv.simplify()
積分因子法の例題
（非同次項 $Q(x)$ がある）一般的な1階線形微分方程式
$$\frac{dy}{dx} + P(x) y = Q(x)$$
を人力で解く際には，積分因子法を使います。結果のみまとめると，
$$g(x) = \exp \left\{\int^x P(x’)\,dx’ \right\}$$
で与えられる積分因子 $g(x)$ を使うと，解が
$$y = \frac{1}{g(x)} \left\{\int^x g(x’) Q(x’) \, dx’ + C \right\}$$
と書くことができる，という方法です。
積分因子法を使って解く例題として，以下のような微分方程式を考えます。
$$\frac{dy}{dx} + \frac{y}{x} = \frac{\sin x}{x}$$
SymPy の dsolve() を使えば，特に積分因子法を意識することなく解くことができます。
y = Function('y')
eq = Eq(diff(y(x), x) + y(x)/x, sin(x)/x)
display(eq)
dsolve(eq, y(x))
この例題をあえて積分因子法の公式に従ってやってみます。
公式にあてはめると，$\displaystyle P(x) = \frac{1}{x}, \ Q(x) = \frac{\sin x}{x}$ ですから積分因子 $g(x)$ は…
P = 1/x
Q = sin(x)/x
g = exp(integrate(P, x))
g
積分因子 $g(x)$ が上記のように求まったので，解は…
○練習：積分因子から解を求める
微分方程式
$$\frac{dy}{dx} + \frac{y}{x} = \frac{\sin x}{x}$$
の積分因子 $g(x)$ が上の上のように求められたとき，これを使って解を求めよ。
ヒント：
解は
$$y = \frac{1}{g(x)} \left\{\int^x g(x’) Q(x’) \, dx’ + C \right\}$$
と書くことができるから…
2階常微分方程式の例：単振動
定数係数2階線形常微分方程式の例として，力学の問題でよく出る，単振動の運動方程式として知られる式
$$ \frac{d^2 x}{dt^2} + \omega_0^2 \,x = 0, \quad (\omega_0 > 0)$$
を解いてみます。$x(t)$ は時間 $t$ の未知関数。
# まず x を関数として宣言
x = Function('x')
# 定数 omega0 > 0 を宣言
var('omega0', positive = True)
# 微分方程式の定義
eq = Eq(diff(x(t), t, 2) + omega0**2 * x(t), 0)
# 微分方程式の表示
display(eq)
dsolve(eq, x(t))
初期条件の設定例
上のセルで，dsolve() で解いた結果の $x(t)$ にあらわれる $C_1, C_2$ は積分定数です。（2階微分なので積分定数は2つ。）
初期条件を課すことで積分定数の値を決めることができます。同じ問題（$t$ での微分を $\dot{}$（ドット）で書いて）
$$ \ddot{x} + \omega_0^2 \,x = 0$$
を，$t = 0$ で $\displaystyle x(0) = x_0, \dot{x}(0) = v_0$ という初期条件で解きます。
初期条件は dsolve() のオプションとして ics={} のように辞書型のデータ形式で与えるのでした。ここでは，
$x(0) = x_0 \Rightarrow$ x(0):x0
$\dot{x}(0) = v_0 \Rightarrow$ x(t).diff(t).subs(t,0):v0
としています。（1階微分に対する初期条件の書き方がすぐには思いつきませんね。）
# （2文字以上の変数の場合は）変数として宣言
var('x0 v0')
dsolve(eq, x(t),
ics = {x(0):x0,
x(t).diff(t).subs(t,0):v0}
)
… というわけで，この初期条件によって $\displaystyle C_2 = x_0, \ C_1 = \frac{v_0}{\omega} $ と求まりました。
最も簡単な定数係数2階微分方程式
あらためて var('x') として独立変数を $x$，y = Function('y') として未知関数を $y(x)$ と表記して，最も簡単な形の
2階微分方程式
$$y^{\prime\prime} + K\, y = 0$$
を解いてみます。
$K$ をvar('K') で単に変数として宣言した場合は…
y = Function('y')
var('x')
var('K')
sahen = diff(y(x), x, 2) + K*y(x)
dsolve(sahen, y(x))
出力結果を見るに，平方根の中が正であることを前提にすれば $K < 0$ と仮定して解いているようです。そこで，以下のように $K$ の値については場合分けをして解いてみます。
$K > 0$ の場合
# K > 0 とする
var('K', positive = True)
sahen = diff(y(x), x, 2) + K*y(x)
dsolve(sahen, y(x))
$K < 0$ の場合
# K < 0 とする
var('K', negative = True)
sahen = diff(y(x), x, 2) + K*y(x)
dsolve(sahen, y(x))
$K = 0$ の場合
var('K')
sahen = diff(y(x), x, 2) + K*y(x)
dsolve(sahen.subs(K,0), y(x))
○練習：1つの解で3つの場合を表す
$\displaystyle y^{\prime\prime} + K\, y = 0$ の $K > 0$ の解は $\displaystyle C_2 = A, \ C_1 = \frac{B}{\sqrt{K}\ }$ と書き直すと
$$y = A \cos \left(\sqrt{K}x \right) + \frac{B}{\sqrt{K}\ } \sin \left(\sqrt{K}x \right) $$
のように書ける。この解を使って，$K < 0$ の場合と $K = 0$ の場合の解を直接求めてみよ。
# 念のために
var('A B K x')
y = A*cos(sqrt(K)*x) + B/sqrt(K)*sin(sqrt(K)*x)
y
この解をそのまま使って $K < 0$ の場合を表示されるには，$K$ の絶対値 $|K| = $ abs(K) を使って以下のように書けることを使う。
$$ \sqrt{K} = \sqrt{- |K|} = i \sqrt{|K|}\quad \mbox{for} \ \ K < 0$$
また，$K=0$ の場合は極限 $\displaystyle \lim_{K \rightarrow 0} y$ をとってみればよいであろう。
ということで，$K > 0$ の場合の解を1つ覚えておけば，$K < 0$ の場合も $K = 0$ の場合もそれから導くことができて別途覚えておく必要はないことが確認できます。
# 以下のように .subs() で一時的代入
sqrt(K).subs(K, -abs(K))
$K < 0$ の場合の $y$ は…
$K = 0$ の場合の $y$ は…
定数係数2階線形常微分方程式
同次方程式
定数係数2階線形微分方程式（同次方程式）は以下のように書ける。
$$y^{\prime\prime} + 2 b\, y^{\prime} + c\, y = 0$$
$b, c$ は定数。
SymPy でこの微分方程式を解くには以下のように…
y = Function('y')
var('b c x')
sahen = diff(y(x), x, 2) + 2*b*diff(y(x), x) + c*y(x)
display(sahen)
dsolve(sahen, y(x))
上の結果を見るに，$b^2 – c > 0$ を暗黙に仮定したときの解になっている。それ以外の場合の解はどのように書けるのであろうか…
ところで，この微分方程式は，以下のように両辺に $e^{b x}$ をかけると
\begin{eqnarray}
0 &=& e^{b x} \cdot \left\{ y^{\prime\prime} + 2 b\, y^{\prime} + c\, y\right\} \\
&=& \left( e^{b x} y \right)^{\prime\prime} + \left(c – b^2\right) \left( e^{b x} y \right) \\
&\equiv& Y^{\prime\prime} + K\, Y
\end{eqnarray}
と書き換えることができる。ここで $Y \equiv e^{b x} y, \ K \equiv c – b^2$ とした。
つまり
$$Y^{\prime\prime} + K\, Y = 0$$
の形の微分方程式を解けばよいことになり，これはすでに $K = c -b^2 > 0$ の場合には一般解が以下のように書けることがわかっているのであった。
\begin{eqnarray}
Y &=& A \cos \left(\sqrt{K}x \right) + \frac{B}{\sqrt{K}\ } \sin \left(\sqrt{K}x \right) \\
\therefore\ \ y = e^{-bx} Y &=& e^{-b x} \left\{A \cos \left(\sqrt{K}x \right) + \frac{B}{\sqrt{K}\ } \sin \left(\sqrt{K}x \right) \right\}
\end{eqnarray}
ためしにこれを関数 y1(K) として定義してみる。
def y1(K):
return exp(-b*x)*(A*cos(sqrt(K)*x)+B/sqrt(K)*sin(sqrt(K)*x))
display(y1(K))
○練習：同次方程式の一般解
同次方程式
$$ y^{\prime\prime} + 2 b\, y^{\prime} + c\, y = 0$$
の一般解を，
- $c – b^2 > 0$,
- $c – b^2 < 0$,
- $c – b^2 = 0$
それぞれの場合に y1(K) を使って表せ。
- $c – b^2 > 0$ の場合
y1(c-b**2)
- $c – b^2 < 0$ の場合
- $c – b^2 = 0$ の場合
非同次方程式
非同次項 $R(x)$ が一般にゼロでないときには
$$y^{\prime\prime} + 2 b\, y^{\prime} + c\, y = R(x)$$
の形になり，これを定数係数2階線形「非同次」方程式と呼ぶ。
非同次方程式の特殊解を人力で求めるには，あそこでまとめているように結構手間がかかるのだが，SymPy では同次方程式も非同次方程式も，微分方程式を解く際には dsolve() を使うだけであり，人間の手間としては同じ。
非同次方程式の例：強制振動
固有振動数 $\omega_0$ の調和振動子に，振動数 $\omega$ の周期的な外力 $F(\omega t)$ が働く場合，この運動は強制振動と呼ばれ，運動方程式は以下のような形になる。
$$\ddot{x} + \omega_0^2 \,x = \frac{F(\omega t)}{m}$$
簡単のために，$\displaystyle \frac{F(\omega t)}{m} = \sin \omega t$ とすると，
強制振動の運動方程式は
$$\ddot{x} + \omega_0^2 \,x = \sin \omega t$$
これは非同次項 $\sin \omega t$ をもつ非同次方程式である。
この微分方程式は SymPy では dsolve() を使って以下のように解ける。
x = Function('x')
var('omega0 omega', positive = True)
eq = Eq(diff(x(t), t, 2)+omega0**2 * x(t), sin(omega*t))
eq
dsolve(eq, x(t))
非同次方程式の例：共鳴
上の出力結果から，$\omega = \omega_0$ の場合は分母がゼロになるため，別途計算する必要があることがわかる。以下のようにして，
eq.subs(omega, omega0) で方程式 eq の中の omega を omega0 に入れ替えて…
display(eq.subs(omega, omega0))
この微分方程式を dsolve() で解きます。
dsolve(eq.subs(omega, omega0), x(t))
積分定数 $C_1, C_2$ がつくのは同次方程式の基本解。それ以外に時間 $t$ に比例して単調増加する振幅を持つ項が現れる。
調和振動子の固有振動数と同じタイミングの周期的外力を加えた強制振動に現れるこの現象は共鳴と呼ばれる。
○練習：極限としての共鳴
強制振動
$$\ddot{x} + \omega_0^2 \,x = \sin \omega t$$
の解は
$$x{\left(t \right)} = C_{1} \sin{\left(\omega_{0} t \right)} + C_{2} \cos{\left(\omega_{0} t \right)} + \frac{\sin{\left(\omega t \right)}}{ \omega_{0}^{2}- \omega^{2} }$$
であるので，$\omega \rightarrow \omega_0$ の場合は分母がゼロになって大変なことになると思いきや，
$$\ddot{x} + \omega_0^2 \,x = \sin \omega_0 t$$
の解は，何事もなかったかのように
$$x{\left(t \right)} = C_{2} \sin{\left(\omega_{0} t \right)} + \left(C_{1} -\frac{t}{2 \omega_{0}}\right) \cos{\left(\omega_{0} t \right)}$$
となる。これを $\omega \rightarrow \omega_0$ の極限として理解したいときには，以下のようにしたらどうか，という話。
まず，$\omega \neq \omega_0$ の場合の解を，積分定数 $C_1$ をちょっと変えて以下のようにする。
\begin{eqnarray}
x{\left(t \right)} &=& C_{1} \sin{\left(\omega_{0} t \right)} + C_{2} \cos{\left(\omega_{0} t \right)} + \frac{\sin{\left(\omega t \right)}}{ \omega_{0}^{2} -\omega^{2} } \\
&=& \left(C^{\prime}_1 -\frac{1}{\omega_0^2 -\omega^2}\right) \sin{\left(\omega_{0} t \right)}
+ C_{2} \cos{\left(\omega_{0} t \right)}
+ \frac{\sin{\left(\omega t \right)}}{ \omega_{0}^{2}- \omega^{2} }\\
&=& C^{\prime}_1 \sin{\left(\omega_{0} t \right)} + C_{2} \cos{\left(\omega_{0} t \right)}
+ \frac{\sin{\left(\omega t \right)} -\sin{\left(\omega_{0} t \right)}}{ \omega_{0}^{2} -\omega^{2} }
\end{eqnarray}
あらためて積分定数を $C_1^{\prime} \rightarrow C_1$ を置き直して書き直すと，$\omega \neq \omega_0$ の場合の強制振動の解は，
\begin{eqnarray}
x{\left(t \right)} &=&
C_1 \sin{\left(\omega_{0} t \right)} + C_{2} \cos{\left(\omega_{0} t \right)}
+ \frac{\sin{\left(\omega t \right)} -\sin{\left(\omega_{0} t \right)}}{ \omega_{0}^{2} -\omega^{2} }
\end{eqnarray}
と書いてもよいことになる。（積分定数は任意であったから）
この解の極限 $\omega \rightarrow \omega_0$ の極限をとると，共鳴の解が得られることを示せ。
$$\lim_{\omega \rightarrow \omega_0} \frac{\sin{\left(\omega t \right)} -\sin{\left(\omega_{0} t \right)}}{ \omega_{0}^{2}- \omega^{2} } = \cdots $$
人力でこの極限を求めるためには，
$$\sin A -\sin B = 2 \cos \frac{A+B}{2} \cdot \sin\frac{A -B}{2}$$
を使って，
\begin{eqnarray}
\frac{\sin{\left(\omega t \right)} -\sin{\left(\omega_{0} t \right)}}{ \omega_{0}^{2} -\omega^{2} } &=&
\frac{2}{\omega_{0}^{2} -\omega^{2} }
\cos \frac{(\omega+\omega_0) t}{2} \cdot
\sin \frac{(\omega -\omega_0) t}{2}\\
&=& -\frac{1}{\omega+\omega_0} \cos \frac{(\omega+\omega_0) t}{2} \cdot
\frac{2}{\omega -\omega_0} \sin \frac{(\omega -\omega_0) t}{2}
\end{eqnarray}
この形にすれば，$\omega \rightarrow \omega_0$ の極限で
\begin{eqnarray}
-\frac{1}{2 \omega_0} \cos \omega_0 t \cdot
\lim_{\omega \rightarrow \omega_0} \frac{2}{\omega -\omega_0} \sin \frac{(\omega -\omega_0) t}{2} &=& -\frac{t}{2 \omega_0}\cos \omega_0 t
\end{eqnarray}
となることがわかりやすいであろう。
○練習：非同次方程式の一般解
非同次方程式
$$y^{\prime\prime} + 2 b\, y^{\prime} + y = R(x) \equiv \sin x, \quad (b > 0)$$
の一般解は，右辺の非同次項を $R(x)=0$ とした同次方程式の一般解（同次方程式の基本解の線形和）に，非同次方程式の特殊解を加えたものです。
まずは SymPy の dsolve() で解かせてみましょう。
y = Function('y')
var('b c x')
eq = Eq(diff(y(x), x, 2)
+ 2*b*diff(y(x), x)
+ y(x), sin(x))
display(eq)
sol = dsolve(eq, y(x))
sol
上の出力結果をみると，$b^2 – 1 > 0$ を暗黙のうちに仮定して解いているようです。
以下のように，この解 sol.rhs を返す関数 y2(K) を定義して…
def y2(K):
return sol.rhs.subs(b**2-1, K)
y2(K)
$b^2-1>0$ の場合には K に b**2-1 を入れて y2(b**2-1) としてやるとよさそうですが…
y2(b**2-1)
$b^2-1<0$ の場合には y2(-abs(1-b**2)) としてやっても，虚数単位 $i$ を含んだ指数関数のままですし…
y2(-abs(1-b**2))
$b^2-1=0$ の場合として y2(0).subs(b, 1) を表示させると，間違った解が出てくる始末です。
y2(0).subs(b, 1)
上記の解が間違っていることは，以下のように，微分方程式の段階で
eq.subs(b, 1) としてから解いた解を見れば明らかです。
dsolve(eq.subs(b, 1), y(x))
ここはやはり，以下のように一旦微分方程式を変形してやります。
\begin{eqnarray}
y^{\prime\prime} + 2 b\, y^{\prime} + y &=& \sin x \\
e^{b x} \cdot \left\{y^{\prime\prime} + 2 b\, y^{\prime} + y \right\} &=& e^{b x}\sin x \\
\left(e^{b x} y\right)^{\prime\prime} + \left(1 -b^2\right) \left(e^{b x} y\right)&=& e^{b x}\sin x \\
\therefore\ \ Y^{\prime\prime} + K\, Y &=& e^{b x}\sin x
\end{eqnarray}
ここで，$\displaystyle Y \equiv e^{b x} y, \ K \equiv 1 -b^2$ とした。
これをまず以下のように $Y$ について解く。
まずは $0<b<1$ すなわち $K = 1 -b^2 >0$ の場合：
Y = Function('Y')
var('K', positive=True)
eq = Eq(diff(Y(x), x, 2)+ K*Y(x), exp(b*x)*sin(x))
sol1 = dsolve(eq, Y(x))
$Y$ が求まったら，$K=1-b^2$ であるから .subs(K, 1-b**2) で $b$ を使った表示に戻し，簡単化 .simplify()…
sol1.rhs.subs(K,1-b**2).simplify()
最後に $y = e^{-bx} Y$ であったから exp(-b*x) をかけて…
Eq(y(x), (exp(-b*x)*sol1.rhs.subs(K,1-b**2).simplify()).expand())
$b>1$ すなわち $K = 1-b^2 < 0$ の場合には…
var('K', negative=True)
eq = Eq(diff(Y(x), x, 2)+ K*Y(x), exp(b*x)*sin(x))
sol2 = dsolve(eq, Y(x))
Eq(y(x), (exp(-b*x)*sol2.rhs.subs(K,1-b**2).simplify()).expand())
$b = 1$ すなわち $K = 1-b^2 = 0$ の場合には…
sol3 = dsolve(eq.subs(K,0).subs(b,1), Y(x))
Eq(y(x), (exp(-x)*sol3.rhs).expand())