Pythonで極める「微分方程式」

5月号日経ソフトウェア(3月24日(水)発売)
私の記事

Pythonで極める「微分方程式」
「厳密解」と「数値解」を完全理解

半年かけて制作した渾身の12ページ

プログラムファイル「dsol.py」
微分方程式の解法(解)には2つ
■公式による解法(厳密解)
□数値解法(数値解)
それぞれ、Pythonでは
■sympy.dsolve(微分方程式)
□scipy.integrate.odeint
を用いることで解けます。
厳密解だけ、数値解だけを得るのはそれぞれのコーディングですみます。
sympy.dsolve(微分方程式)は短く簡単
scipy.integrate.odeintを使いこなすには慣れが必要

そこで、数値解法に
オイラー法と4次ルンゲ・クッタ法を選び、解法を実際にコーディングすることにしました。

はたして、一つの微分方程式に対して
厳密解と2つの数値解、計3つの解曲線のグラフ同時に表示するプログラムをつくることにしました。
このアウトプットにより、数値解法の精度を見ることができます。

一つの微分方程式を、公式による解法と数値解法に渡して計算させるのに苦心しました。
それには、Pythonの変数の使い方がポイントになります
通常変数(数値計算)とシンボリック変数・シンボリック関数
数ヶ月のコーディングのおかげで、これらの変数の使い方に慣れることができます。
ちなみに、日本語のPython解説本でこれらを解説したものはほとんどありません。
Webドキュメントを頼りにしました。

驚異のsympy.dsolve()!
公式で解ける微分方程式のほとんどがこれで解けるような気がするほど強力です。
なかでもひっくり返るほど驚いたのが
sym.solvers.classify_ode()
これに微分方程式を渡してあげると、微分方程式の型を判別し返してきます
‘separable'
, ‘1st_exact'
, ‘Bernoulli'
, ‘1st_power_series'
, ‘lie_group'
, ‘separable_Integral'
, ‘1st_exact_Integral'
, 'Bernoulli_Integral'

大学で微分方程式を習ったことがある人ならば読めるはずです

分離形
線型
ベルヌーイ型
級数展開型
リー型

ここが面白いところなのですが、難しい微分方程式を解かせとsympy.dsolve()でもまともな結果を返してこない場合があります。
そのときに、sym.solvers.classify_ode()を用いて型を判別させ、その結果を
sympy. dsolve()の引数に「hint="Bernoulli"」として渡してあげることで解けます。
だったら自動でやってくれよとツッコミたくなるのですが、
ここはコーディングする人の経験からくる勘との合わせ技でsolutionを得ることに成功します。

モデル1:マルサス・人口モデル
モデル2:技術革新の普及モデル
モデル3:ベルタランフィの魚の成長モデル
を例にしてありますが、このdsol.pyを用いれば、
微分方程式+初期条件
を設定しさいすれば、
公式のよる解法と2つの数値解法により、厳密解表示および3つの解曲線のグラフの同時表示
ができるプログラムが完成しました。
グラフ表示にも細かいテクニックをたくさん投入しています。これにも時間がかかりました。

微分方程式を解いてみたい人
Pythonの通常変数(数値計算)とシンボリック変数・シンボリック関数の使い方を知りたい人
微分方程式の初学者
にはうってつけのプログラムです。