logo

常微分方程式を時間を逆に解くトリック 📂数値解析

常微分方程式を時間を逆に解くトリック

メソッド

dydt=f(y) {\frac{ d y }{ d t }} = f(y) 上記のように与えられる常微分方程式を解くための数値的ソルバー KK は未来の時点 T>0T > 0 と初期値 y0y_{0} に対して K(f,T,y0)={y(t):t[0,T]} K \left( f , T , y_{0} \right) = \left\{ y(t) : t \in [0, T] \right\} のように動作すると仮定する。逆に、S>0S > 0 だけ過去の時点に対して {y(s):s[S,0]}\left\{ y(s) : s \in [-S, 0] \right\} を得たいならば、K(f,S,y0)K \left( -f , S , y_{0} \right) を計算すればよい。

説明

上記で紹介されたトリックは、ソルバーがどのような形であれ時間を過去に送る方法を扱っている。時間を過去に戻すということは、特に力学系の分析で不安定固定点を見つけるなどの作業で有用に使われる。

トリックの核心は、時間を逆向きに流れるようにする方法が非常に簡単に、ff の代わりに f-f を使うだけで済むという点だ。

モチーフ

f(y)=dydt=limh0y(t+h)y(t)h    y(t+h)y(t)+hf(y) \begin{align*} & f(y) = {\frac{ d y }{ d t }} = \lim_{h \to 0} {\frac{ y(t + h) - y(t) }{ h }} \\ \implies & y (t + h) \approx y(t) + h f(y) \end{align*} 例としてオイラー法を考えると、十分に小さい hh に対して上記のような近似で導かれる。しかし、微分係数 y˙\dot{y} の定義に従えば、その分子は必ずしも y(t+h)y(t)y(t + h) - y(t) である必要はなく、y(t)y(th)y(t) - y(t - h) でも全く問題ないはずであり、それなら次のような計算も可能となる。 f(y)=dydt=limh0y(t)y(th)h    y(th)y(t)hf(y)=y(t)+h(f(y)) \begin{align*} & f(y) = {\frac{ d y }{ d t }} = \lim_{h \to 0} {\frac{ y(t) - y(t - h) }{ h }} \\ \implies & y (t - h) \approx y(t) - h f(y) = y(t) + h \cdot \left( - f (y) \right) \end{align*} このような過去方向への解法はオイラー法を KK とした時に、実は K(f,h,y0)K \left( -f , h , y_{0} \right) を求めることと同じである。

変数置換

しかし、KK がオイラー法でない場合、上記のような直観だけで結果を導出することは難しい。この一般化を s=ts = -t のような変数置換によって説明しよう。ここで ss が正の方向 S>0S > 0 に向かって大きくなるということは、tt が負の方向に動くということ、つまり時間 tt が逆向きに流れることを意味する。そうすると、元の常微分方程式 y˙=f(y)\dot{y} = f(y)連鎖律に従って次のように書くことができる。 dydt=f(y)    dydsdsdt=f(y)    dyds=f(y)    y=f(y) \begin{align*} & {\frac{ d y }{ d t }} = f(y) \\ \implies & {\frac{ d y }{ d s }} {\frac{ d s }{ d t }} = f(y) \\ \implies & {\frac{ d y }{ d s }} = - f(y) \\ \implies & y ' = - f(y) \end{align*} これは形式的に g=fg = - f に対して y=g(y)y ' = g (y) を解くのと異ならず、結局 ff の符号だけを反転させて元のメソッド KK を適用することと同じである。