Pythonのodeintで微分方程式の数値解析をした話

というわけで

今の世の中いろいろな情報が検索で出てくるわけなのですが、とある講義の課題で微分方程式の数値解析をする必要に迫られた際にPython使おうとして、中々目当ての情報が出なかったんですね。いや、Pythonのodeintで確かに例を示しながら解いてるサイトはあるんですけど、(何故か)私がそれを読んで実装するまでに非常に時間がかかってしまい。

なので、いつか先の未来の私がまたodeintを使おうとしてハマらないよう、彼に向けての意味合いも含めて記録しておきます。

出来るだけ一般化して

出てくるサイトが皆ご丁寧に $x''+x=0$ とかの例を解くコードをそのまま載せてくださっているために、それの解釈でまず時間かかってしまったのが私の失敗で。というわけで、先ずは細かいコード云々よりももう少しユルく、こんな感じに書きたいよねっていう頭の中を再現してみます。

1.データをどうするか

先ずは一変数(例えば $x$ )の微分方程式を考えましょう。Python君に数値的に解かせてあげる際、解いた後のデータをどう整理して貰うかってなると、配列を使ってあげて
[
[ $x,~x',~x'',\cdots$ ] ${}_{t=0}$ ,
[ $x,~x',~x'',\cdots$ ] ${}_{t=1}$ ,
[ $x,~x',~x'',\cdots$ ] ${}_{t=2}$ ,
$\cdots$
]
とすれば便利そうですね(刻み幅は上では1ですが、別に例なので任意で)。

なので、こんな感じ([ $x,~x',~x'',\cdots$ ])にデータを扱ってやれば良さそうです。この一次元配列を、ここではdataという名前にしてあげましょう。
data[0]= $x$ で、もっと一般化すればdata[n]= $\dfrac{d^n}{dt^n}x$ です。

2.微分方程式の分割…の前に、理想的な"ODEINT君"

1.のような考え方を元に、オイラー法をfor文使って実装出来なくもないですが、ここでは先ずodeintを使ってサクッと解くことに主眼を置きましょう。odeintの方がfor文でやるよりアッサリ書けてしまうのですが、それでももう一段階だけお膳立てをしてあげなければなりません。それが、微分方程式の「分割」です。余談ですが、私が初めて数値解析をしてみようとした時、この辺りで「何故こうしたいのか」がよく分からなくって一旦odeintから尻尾巻いて逃げ出したんですね。そしてお手製のオイラー法をfor文で実装したりして。なので、お気持ちも出来る限り混ぜながら進めていきます。

一番人間に楽させてくれるodeint君の理想的偽物(ODEINT君)ってどんなやつでしょうか。ODEINT君に対して解きたい微分方程式 $f(x,~x',~x'',\cdots)=0$ と初期値を見たまま

#デタラメコード

t=np.arange(0, 10000000, 0.1) #時間(開始, 終了, 刻み幅)

data_0=[1, 2, 3] #data_0:t=0での[x, x', x'']

def func(data):
    #f(x, x', x'')=0をdataの要素で書く
    f(data[0], data[1], data[2])=0

result = ODEINT(func, data_0, t)
#ODEINT(微分方程式, 初期値, 時間)

(勿論上のコード(もどき)は嘘ばっかりで、funcの定義の部分もダメダメです。実際に書いて実行しようとしたら沢山エラー吐かれてしまいそう)

って感じに打ち込んであげると、resultが1.で(こうしたいなぁ)って思ってた二次元配列の形
[
[ $x,~x',~x'',\cdots$ ] ${}_{t=0}$ ,
[ $x,~x',~x'',\cdots$ ] ${}_{t=1}$ ,
[ $x,~x',~x'',\cdots$ ] ${}_{t=2}$ ,
$\cdots$
]
で得られる…なんて、非常に楽じゃないでしょうか。一部Pythonの構文さえ知っておけば、特に頭を使うこともなく書けてしまいます。

3.微分方程式の分割

では実際に、ODEINT君程では無いけども楽はさせてくれるodeint君を見ていきましょう。
といっても、上のデタラメコードとほぼ同じ書き方で書けてしまいます。直す必要があるのが、「def func(data):」の部分。微分方程式をPython君に理解して貰うわけですが、そもそもPythonで関数functionを定義する時って

def function({引数}):
    return {戻り値}

で大抵書かれます。ODEINT君を考えたのと同様に微分方程式をそのまま書こうにも、一体何が「戻り値」なのさ?ってなってしまいます。

odeint君に渡してあげる、引数と戻り値の存在する関数は、微分方程式を以下のように"分割"したものになります。
f:id:sGya_youtoo:20190108200902p:plain
…と、数式苦手な方(私もそう)が見るとこれだけで逃げたくなるような画像になってしまいましたが、そこまで慌てなくても大丈夫です。

大抵の場合(例外を私はまだ知らない)、 $y_n=\frac{d^n}{dt^n}x$ としてしまえば、勝手に分割出来ます。

こうして分割した微分方程式を、Python君に教えてあげましょう。関数名は、diffとしておきます。引数は、data=[ $x,~x',~x'',\cdots$ ]です。

def diff(data):
    return [φ_0, φ_1, φ_2, ...]

勿論、実際にコードを書く際は、 $\phi_n$ を文字通りそのまま戻り値に書くのではなく、 $\phi_n$ を $x,~x',~x'',\cdots$ の式で表したモノを書いて下さい。ちゃんと、 $x$ とかではなく、dataの要素として書くこともお忘れなく。

これで、odeint君へのお膳立て、微分方程式の分割が終了です。

先程までの $\phi_n$ への注意を踏まえて、お気持ちを優先して一般的にコードを書いてみると、こうなります。

t=np.arange(T_0, T_end, ⊿t)
    #時間(開始, 終了, 刻み幅)

data_0=[x_0, x'_0, x''_0, ...]
    #data_0:t=0での[x, x', x'', ...]

def diff(data):
    #分割した微分方程式
    #data=[x, x', x'', ...] : 引数
    #[φ_0, φ_1, φ_2, ...] : 戻り値
    return [φ_0, φ_1, φ_2, ...]

result = odeint(diff, data_0, t)
    #odeint(分割した微分方程式, 初期値, 時間)

実際に解いてみよう

$x''+(x^2-1)x'+x=0,~(x,x')_{t=0}=(1,1),~t:0～100,~\Delta t=0.01$
の微分方程式を数値的に解いてみましょう。

ちなみに、この微分方程式はVan der Pol方程式と呼ばれる微分方程式 $x''+\mu(x^2-1)x'+x=0$ の $\mu=1$ の場合です。

先ずは分割

というわけで、分割しましょう。
$\begin{align*} &y_0=x && \frac{d}{dt}y_0=y_1 \\ &y_1=x' && \frac{d}{dt}y_1=x''=-(y_0{}^2-1)y_1-y_0 \end{align*}$

いよいよコードを

さて、では実際にコードを書いてみましょう。Windows10 64bit, Python3.7(Anaconda), Jupyter Notebookを使用しています。

import

必要なのは、時間を設定してやるためのnumpyと、odeintを使うためのscipy.integrateなので、

import numpy as np
from scipy.integrate import odeint

と書いてやります。

時間

時間を以下のように指定してあげましょう。

t = np.arange(0, 100, 0.1)
     #時間(T_0, T_end, ⊿t)

初期条件

data=[ $x,~x'$ ]の初期条件を入れましょう。

data_0 = [1, 1]
     #初期値 [x, x'']_t=0

分割した微分方程式

今回は、分割した微分方程式の名前をdiffとしましょう。

def diff(data, dummy):
    return [data[1], -(data[0]**2-1)*data[1]-data[0]]

※diffの引数にdummyは要るのか?という話ですが、これをdiff(data)のようにすると

TypeError: diff() takes 1 positional argument but 2 were given

というエラーメッセージが出てくるので、一先ず何らかの引数をダミーで入れてあげないといけないようです(よくわかってない)。

最後に

odeint君に仕事してもらいましょう。

result = odeint(diff, data_0, t)
     #odeint(分割後の微分方程式, 初期値, 時間)

解いたとはいえ

折角Python君とodeint君で微分方程式を数値的に解いても、データがパソコンの中で眠っていては意味がないですね。出そうと思えばprint(result)なんかで

[[ 1.00000000e+00  1.00000000e+00]
 [ 1.00994948e+00  9.89850878e-01]
 [ 1.01979604e+00  9.79406878e-01]
 ...
 [ 2.00861912e+00  1.68830849e-03]
 [ 2.00853633e+00 -1.81464086e-02]
 [ 2.00825817e+00 -3.73867155e-02]]

みたいに出ますけど、そんな膨大な数値の羅列見ても大抵の人は眠くなるだけです。

というわけで、今度は解析結果を視覚的に捉えられるようにしていきましょう。Python君にグラフを描かせるのです…が、そろそろこの記事もかなりの分量になってきてしまいましたし、コレは次の機会に回すことにしましょう。書いたら[リンク]欄にリンク貼っときます。

最後に

以上で、Pythonのodeintを使って微分方程式を数値的に解く方法の紹介は終わりです。書いている人は、別にodeintの動作原理を詳しく学んだ人ではなくただ単に数あるネットの情報を拾いつつ何個か微分方程式を解いてみただけの人なので、厳密性は保証しかねる部分が多いです。

何かしらマズイ点を見つけた方は、コメントなり私のTwitterアカウント[@sGya_youtoo]までご連絡なり頂ければ幸いです。

参考

SciPyで微分方程式を解く
 odeintを使った常微分方程式の数値計算 - とつとつとしてろうとせず
 ファン・デル・ポール振動子 - Wikipedia

リンク

ちょっと待ってね。

思考のお道具箱