Juliaのパッケージを使って常微分方程式を解きます。
はじめに
RSSリーダーに登録しているThe Julia Language Blogから新着記事が届いていた。
記事タイトルにあるComposabilityとは組み合わせられること、という感じだろうか。 冒頭のパラグラフで言われているように同一のツール（パッケージと呼ぶべきか）が種々の問題に適用可能になっているのだと理解した。 ODEとは常微分方程式 (ordinary differential equation) の意であり、しかし略語で示されると一瞬なんのことかわからなくなる。 Neuralと形容されていて、また本文中にもneural networksとあり、そうした分野では常微分方程式を取り扱うということなのだろう。ニューラルネットはほとんどわかっていない。
記事の書き出しにあるSciML Common Interfaceからリンクを辿ってみると、このパッケージは本当に色々な微分方程式、あるいは問題設定を取り扱うことができて便利そう。 最近の解説動画としてSciMLのエコシステムを説明しているものを見つけた。 ここに出てくる例題、climate modelは自分にとって多少理解しやすいものであったから、少し探してみたらarXivで文字になった内容を見つけることができた（ので、いま読んでいる途中）。
チュートリアルをやってみる
前述のSciMLにはチュートリアルが存在するということで、SciMLTutorialsパッケージを追加してみた。その後、簡単な気持ちでusing SciMLTutorialsからのSciMLTutorials.open_notebooks()を実行したら、Precompilingになんか5時間くらいかかって疲れた……事情はよくわからないまま、しかし例題と解説の書かれたJupyter notebooksが手に入った。
ノートブックのうち "01-ode_introduction.ipynb" というもの、内容としてはAn Intro to DifferentialEquations.jlに相当するところは、DifferentialEquationsパッケージを追加することで無事に動作させることができた：
add DifferentialEquationsでパッケージを追加したさいにはProgress途中で失敗したらしい挙動が見られて、27/30の表示。NGになった依存関係は以下：SundialsBoundaryValueDiffEqDiffEqFinancial
- 日を改めて
up DifferentialEquationsを実行してみたら特にエラーは出ることなく "67 dependencies successfully precompiled" で67/67の表示。なぜか依存数が増えて、エラーも解消されているのはよくわからない…
An Intro to DifferentialEquations.jlの冒頭で取り扱っているのは単純な常微分方程式である：
数値解析に慣れた身でこうした方程式を解こうとするとすると反射的に離散化してしまうけれど、そうではなく元の式のまま扱うことができるのはとても便利。 Juliaでこういったパッケージが提供されていることは意識の片隅に置いておきたい。
問題の検算をしようと思って方程式を離散化しはじめたら、つい凝ってしまった：
function main(n::Int)
u = Array{Float64,1}(undef,n)
t = Array{Float64,1}(undef,n)
u .= 0.0
t .= 0.0
u0 = 1.0
dt = 1.0 / n
c = 0.98
u[1] = u0
t[1] = 0.0
for i = 2:n
t[i] = dt * i
u[i] = c * u[i-1] * dt + u[i-1]
end
return n, u[n], abs(u[n] - exp(c * 1.0))
end
println(main(10))
println(main(100))
println(main(1000))
これくらいのsnippetならばJuliaでさくっと書けるようになったのでうれしい。
この常微分方程式の解はネイピア数を底とする指数関数 になることに気づいたので、戻り値として解析解との誤差の絶対値を返してみる。
実行したときの様子：
$ julia main.jl (10, 2.3196425204343, 0.3448137214951168) (100, 2.6260395923548363, 0.038416649574580575) (1000, 2.6605705529590007, 0.00388568897041619)
分割数nを10倍に増やすと、誤差はざっくり1/10になってるので良い感じ。
精度を上げるのにもう少しやりようはあるけれども、そこにリソースを割くよりかは素直にDifferentialEquationsパッケージを使うほうが早そうだ。
パッケージの使い方はもう少しおさらいしていきたい。
Follow-up
前回記事で触れたクエット流れの支配方程式にはナビエ・ストークス方程式由来の偏微分記号が出てくるけど、実質的には2階常微分方程式になっている。 なので積分すればたやすく解を得ることができて、速度uは壁からの距離に応じた直線分布となる。 記事ではフォローできてなかったけど、この偏微分方程式は分類としては放物型になる。 Computational Fluid Dynamics: The Basics With Applications (McGraw-Hill Series in Aeronautical and Aer)の第3章では、偏微分方程式の場合分けとその性質が解説されているので読み進めているところ。 英語なのでしっかり読もうとすると骨が折れるけど、この著者の書く英語はとても読みやすいので助かる。
おわりに
The Julia Language Blogの新着記事をきっかけに微分方程式を取り扱うパッケージSciMLに入門した。
冒頭で触れたRSSリーダーにはReederを使っていて、これはRebuildのエピソード287、Taro Minowaさんがゲストの回で教わったもの。 数多くのRSSリーダーがサービスを終了していく中で（おそらく収益化しにくいのだろう）、これだけ高品質なアプリを作り続けてくれる製作者には素直に頭が下がる。 気に入った読み物があればRSSリーダーにそっと登録して読み続ける、その書き手の静かなファンになる感覚がよい。 ポッドキャストもまた購読して配信を待つという意味で、広義のRSSフィードといえるのかもしれない。 登録してファンになる感覚はそう違うものではなく、ポッドキャストが好きな理由にはそういう側面もあるのかなと思う。■
