06. 時間対称スキーム, Runge-Kuttaスキーム
Photo by James Cheung on Unsplash
単純な常微分方程式を Euler 法「以外」で計算しようとすると、一般にはどんな技術が必要になるだろうか?:
Euler 法は簡単で扱いやすいんだけれども、パラメータを少し変えると不安定になってしまっていた. しかも、詳しく示してはいないけれども、一般に「近似精度」もあまり良くはないんだ.
そこで、もう少し「良い」近似解法が考えられないか、検討してみよう.
対象の微分方程式
対象として今回も前回扱ったものと同じ、人口 $u = u(t)$ ($t$: 時間) に対する常微分方程式
\[ \frac{du}{dt} = \gamma\, \left( \sin(u) \right)^{1.2}, \,\, \gamma > 0, t \in (0,\infty) \]
を考えよう. ただし、$u$ の初期値 $u(0)$ については $0 < u(0) < 1$ としておく.
まずは、「対称な」スキームを扱ってみよう
時間刻み幅 $\Delta t$ を一定にして、微分を含め、全体を「時間方向に対称に」近似する, いいかえれば、Taylor 展開の中心点を $t = (n + 1 / 2)\Delta t$ にとることがもっとも素直になるような近似をする、ということで近似計算式(スキームと呼ぶ)を作ることを考えよう.
以降、$u_n \cong u(n\Delta t)$ と想定した近似値 $u_n$ を用いて話を進めよう.
こうした時間方向に対称なスキームの具体形は唯一ではなく無数に考えられるが、ここではなるべくシンプルに近似を作ってみよう. つまり、時間微分項 $du / dt$ を時間 $t = (n+ 1 / 2)\Delta t$ で中心対称な差分 $(u_{n+1} - u_n) / \Delta t$ で近似し、 他の部分に出てくる $u(t)$ は時間方向の平均 $(u_{n+1} + u_n) / 2$ などで近似するという感じだ.
今回の対象としている微分方程式にこうしたシンプルな対称な近似を行ってみると、例えば
\[ \frac{u_{n+1} - u_n}{\Delta t} = \gamma \left( \sin( \frac{ u_{n+1} + u_n}{2} )\right)^{1.2} \]
というような近似式が得られる.
この近似式を採用したとしよう. すると数値計算のプロセスでは $u_n$ は既知で $u_{n+1}$ は未知量なので、この未知量をどう求めたら良いか、を考える必要がある.
そこでより分かりやすいように、 未知量を $u_{n+1} = u$, 既知量を $u_n = U$ と表記して上の近似式を整理しよう. すると、関数 $f(u, U)$ を
\[ f(u, U) := u - U - \Delta t \gamma\, \left( \sin( \frac{ u + U }{2} )\right)^{1.2} \]
と定義すれば、上の近似式の数値計算での「使い方」は、
与えられた $U$ に対して $f(u, U)$ がゼロになるような $u$ を求める
ことで時間ステップを 1つ進める、ということが明確になる. これなら,以前の授業 03. Newton法! や 04. ゼロ点求解法その他 の知識が使えるね.
あとは Euler 法同様に、初期値は既知で $u_0 = u(0)$ とすることができるので、さらなる近似値を上のスキームを用いて $u_1, u_2, \cdots$ と順番に求めていけば好きなだけ先の近似値を計算できるはずだ.
Newton 法を用いてみよう!
以前 Newton 法について学習したが、復習とプログラムへの慣れを兼ねて、もう一度この問題を対象として丁寧にやってみよう.
まず、$f(u,U)$ の $u$ に関する微分を求めよう. これは簡単で、
\[ \frac{df}{du} = 1 - 0.6 \,\Delta t \,\gamma \,\sin(\frac{u+U}{2})^{0.2} \, \cos(\frac{u+U}{2}) \]
となる.
これで Newton 法を適用できる. まず,準備段階は以下のようになるだろう.
|
|
そして,Newton 法で古い時間ステップの関数値 $U$ から新しい時間ステップの関数値 $u$ を求める関数を次のように定義する,という感じになるね.
|
|
これで準備が整ったので、実際に数値計算をして時間発展をさせてみよう.
|
|
計算がスムーズに終わっているようなので、早速、グラフにしてみてみよう.
|
|
どうやら無事に計算がうまくいっているようだ.まずはめでたしめでたし.
Euler 法は不安定だったけどこの「対称なスキーム」はどうかな?
前回やってみたように、パラメータによっては Euler 法はふっとんでしまい、不安定だった. この数値スキームはどうだろう? 試してみよう.
Euler スキームが動かなくなったパラメータ
|
|
を設定して、動かして、プロットしてみる.
|
|
おお、これはきちんと動作するぞ. どうやら、ちょっとばかり手間だけれども、「対称なスキーム」は Euler 法より安心して使えそうだ.
非線形方程式を解く部分を、ライブラリに頼ってみよう!
復習を兼ねて、以前学習した、NLsolve パッケージを使ってみよう. まず、このパッケージを毎回行うルーチン的な手続きは以下の通りだったので、そうしよう.
|
|
あとは前やったように、以下のようにすれば良いはずだ.
|
|
ふむ、無事に動いたようだ.グラフにしてみてみよう.
|
|
グラフを見ると、上の例とまったく同じようだ.内容的にも問題ないと見て良いだろう.
念のために、ライブラリを使った場合でもパラメータを変えてみよう
先と同じように、$\Delta t = 1.2$ として挙動を見てみよう.
|
|
…略…
39: 3.14158103394284 1.2396717852940844e-5
40: 3.1415823736216253 1.0949807560126888e-5
41: 3.1415835327860404 9.700385960235457e-6
DomainError with -3.910542430870352e-7:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
[1] throw_exp_domainerror(::Float64) at .\math.jl:37
[2] ^ at .\math.jl:888 [inlined]
[3] f(::Float64, ::Float64) at .\In[1]:7
…略…
Euler 法よりはずいぶん先まで計算できるという意味でかなりマシだけれども、この場合も途中で Euler法とまったく同じ理由で計算が停止してしまうことがわかる. ライブラリが万能ではない,ということだね.
このように、簡単な問題でも、手法次第で結果が得られるかどうかが容易に変わってしまう ので、現実には計算手法を複数知っておいてそれをうまく切り替えて使う必要がありそうだ.柔軟にやるしかないだろうね.
Runge-Kutta 法を試してみよう
さて,さらに新しい解法を紹介しよう.
ここで紹介する Runge-Kutta 法とは常微分方程式
\[ \frac{du}{dt} = r(u,t) \]
に対して $u(t)$ から $u(t+Δt)$ の近似値を計算する手法の一つで、 $t$ から $t + Δt$ の間の $u$ の 時間変化 $r$ の近似値をいくつか計算して、それを使って $u(t+Δt)$ の近似値を作るもの、だな.
その中でも有名なのは、以下の「古典的 4次 Runge-Kutta 法」と呼ばれるもので、なかなか「頑丈 = 安定して計算できる」ことと、代入計算だけで済むので「計算が速い」ことが知られている. ただし、時間方向の対称性とかは失われてしまうのでそのあたりはバーターだね.
古典的 4次 Runge-Kutta 法
古典的 4次 Runge-Kutta 法というものでは,まず, 時間 $t$ での $u(t)$ の値を $u$ として、
\[ \left\{\begin{array}{rcl} r_1 & = & r(u, t), \cr r_2 & = & r(u + (\Delta t / 2) \, r_1, \, t+\Delta t / 2), \cr r_3 & = & r(u + (\Delta t / 2) \, r_2, \, t+\Delta t / 2), \cr r_4 & = & r(u + \Delta t \, r_3, \, t+\Delta t) \end{array} \right. \]
と順番に計算することで、$t$ から $t + Δt$ の間の $u$ の 時間変化 $r$ の近似値を 4つ計算する.
そしてこれらを用いて、
\[ u(t + \Delta t) \, { の近似値} := u + \Delta t \left( \frac{r_1 + 2r_2 + 2r_3 + r_4}{6} \right) \]
として、時間ステップを 1つ進める近似解法がこの方法だ.
早速、プログラムしてみよう. まず、この問題の場合は、常微分方程式の右辺 $r(u,t)$ は次のようになる.
この問題の場合,右辺に $t$ が登場しないので省略しているよ!
|
|
そして、Runge-Kutta 法で 1ステップすすめる関数をプログラムで書くと次のような感じだな.
|
|
早速使ってみよう! まずは安全そうなパラメータで.
|
|
で、動いているようなのでプロットだ.
|
|
ふむふむ、うまくいっているな.
では、$\Delta t$ を大きくしてみて、うまくいくかチェックしよう.
|
|
今回も大丈夫なようなので、プロットしよう.
|
|
どうやら、$\Delta t$ を大きくしてもこの問題では Runge-Kutta 法がうまく動作するようだ.
困ったらとりあえず Runge-Kutta 法で と覚えておくとよさそうな気がするね.
それぞれの方法の計算速度の比較をしてみよう
準備として、初期値と計算回数を入力すると最後の解を出力する形で関数を作ろう.
まずは Euler 法だ. まず、前回の授業で定義した Euler 関数を再度以下のように定義して、
|
|
以下のように全体をつくる.
|
|
次に、対称なスキームについてだ.非線形部分の解法にはライブラリを用いたものとしよう.
|
|
それから Runge-Kutta 法についても以下のように同様につくろう.
|
|
念のために、どれもきちんと動くことを確認しよう.
|
|
3.119420135870497
|
|
3.1184787108216545
|
|
3.118421279457686
ほぼ同じ最終値1が出力されていることから見ても、無事に動作していることがわかる.
では、これを使って計算時間を比較しよう.
それには、BenchmarkTools というパッケージに入っている、@benchmark
というマクロが便利だ.
おそらくみなさんの Julia 環境にはこのパッケージは入っていないだろうから,まずはインストールだ.
|
|
インストールが済んだら,早速このマクロを使って計算速度を測定、比較してみよう. まずパッケージの利用宣言を以下のようにして、
|
|
あとは測定したい関数を @benchmark
マクロに以下のように渡すだけだ.
|
|
BenchmarkTools.Trial:
memory estimate: 7.83 KiB
allocs estimate: 501
--------------
minimum time: 15.499 μs (0.00% GC)
median time: 16.599 μs (0.00% GC)
mean time: 17.222 μs (0.69% GC)
maximum time: 1.204 ms (98.38% GC)
--------------
samples: 10000
evals/sample: 1
結果を見ると、10000回測定していて統計的な結果までわかる. 大雑把に言えば Euler 法で 100ステップ計算させるとおおよそ 17μs 程度の計算時間がかかることがわかるね.
Benchmark 結果の中の GC は garbage collection (ゴミ処理) の頭文字で,使わなくなったメモリの整理等の作業を意味している.これが発生するとかなり時間を食うことがわかるね.garbage collection は PC で使えるメモリが少ないほど頻繁に発生するので,メモリが少ないとつらいことがわかるね.
同様に、対称スキームと Runge-Kutta スキームの計算速度も測定しよう.
|
|
BenchmarkTools.Trial:
memory estimate: 485.95 KiB
allocs estimate: 11001
--------------
minimum time: 451.700 μs (0.00% GC)
median time: 471.500 μs (0.00% GC)
mean time: 527.247 μs (9.83% GC)
maximum time: 6.514 ms (91.65% GC)
--------------
samples: 9475
evals/sample: 1
|
|
BenchmarkTools.Trial:
memory estimate: 39.08 KiB
allocs estimate: 2501
--------------
minimum time: 69.599 μs (0.00% GC)
median time: 70.101 μs (0.00% GC)
mean time: 73.282 μs (1.29% GC)
maximum time: 1.291 ms (94.17% GC)
--------------
samples: 10000
evals/sample: 1
すると、対称スキームでは計算時間が 400-500μs ほどで、Runge-Kutta スキームでは 70μs ほどだということもわかる.
計算式を眺めると Runge-Kutta スキームは Euler スキームの 4倍の計算時間がかかるはず. なので、Runge-Kutta スキームのこの計算時間測定値は想定通り!!と言っていいかな.
また, 非線形部分を毎回解かないといけない対称スキームは想像できたように計算時間が大変大きいこともよく分かる.
今回は素直な結果になったが、プログラムが要する計算時間は直感による予想と異なることが多々ある.計算時間については、なるべくこのように「測定した値」をベースにして議論するほうが良い.
もう少し複雑な Runge-Kutta 法も紹介
詳しく知りたければ、英語だけれども気軽に読める wikipedia の List of Runge-Kutta methods を見ておくと良いかな.
4(5)次 Runge-Kutta 法 (Runge-Kutta-Fehlberg 法)
これは、計算途中で $\Delta t$ をうまく調整していくことを想定したタイプの Runge-Kutta 法の一つで、有名なものである. 具体的には、まず、時間 $t$ での $u(t)$ の値を $u$ として、
\[ \left\{ \begin{array}{rcl} r_1 & = & \displaystyle r(u, t), \cr r_2 & = & \displaystyle r\left(u + \frac{1}{4} \Delta t \, r_1 , t+ \frac{1}{4} \Delta t \right), \cr r_3 & = & \displaystyle r\left( u + \Delta t \left( \frac{3}{32} \, r_1 + \frac{9}{32}\, r_2 \right) , t + \frac{3}{8} \Delta t \right), \cr r_4 & = & \displaystyle r\left( u + \Delta t \left( \frac{1932}{2197} \, r_1 - \frac{7200}{2197}\, r_2\ + \frac{7296}{2197}\, r_3 \right) , t + \frac{12}{13} \Delta t \right), \cr r_5 & = & \displaystyle r\left( u + \Delta t \left( \frac{439}{216} \, r_1 - 8\, r_2 + \frac{3680}{513}\, r_3 - \frac{845}{4104}\, r_4 \right) , t + \Delta t \right), \cr r_6 & = & \displaystyle r\left( u + \Delta t \left( -\frac{8}{27} \, r_1 + 2\, r_2 -\frac{3544}{2565}\, r_3 + \frac{1859}{4104}\, r_4 - \frac{11}{40}\, r_5 \right) , t + \frac{1}{2} \Delta t \right) \end{array} \right. \]
と上から順に計算して、$t$ から $t + Δt$ の間の $u$ の 時間変化 $r$ の近似値を 6つ計算する(かなり複雑に見えると思うが).
そしてこれらのうち、4つの値 ($r_1, r_3, r_4, r_5$)を用いて、
$u(t + \Delta t)$ の 4次近似値 $U$ $\displaystyle = u + \Delta t \, \left( \frac{25}{216}\,r_1 + \frac{1408}{2565}\, r_3 + \frac{2197}{4104}\, r_4 - \frac{1}{5}\, r_5 \right) $
を、そして、同様に 5つの値 ($r_1, r_3, r_4, r_5, r_6$)を用いて、
$u(t + \Delta t)$ の 5次近似値 $\hat{U}$ $\displaystyle = u + \Delta t \, \left( \frac{16}{135}\,r_1 + \frac{6656}{12825}\, r_3 + \frac{28561}{56430}\, r_4 - \frac{9}{50}\, r_5 + \frac{2}{55}\,r_6 \right) $
を求める.
こうすると、 低次近似値からみたら高次近似値は随分「良い」はずなのでこれを真値の代わりに使うというアイディアが考えられる. このアイディアに基づくと、低次の方の 4次近似値 $U$ の誤差を
\[ \epsilon_4 \cong \| U - \hat{U} \| \]
として、具体的に数値で見積もれることになる.
ベクトルの評価にノルムを使っているよ.
そして、そもそも $p$ 次近似値の局所誤差(1ステップだけで発生する誤差) $\epsilon_p$ は
\[ \epsilon_p \cong C\, (\Delta t)^{p+1} \]
となるように上の Fehlberg 公式は設計されている($C$ は時間と $u$ 等に依存する適当な係数).
時間発展問題の近似式で「$p$ 次局所誤差」という場合は、通常はこの意味だ.
そして、この状況下で、$\Delta t$ を $\Delta t_{\rm opt}$ に修正して局所誤差を $C_R$ 程度以下に抑えたい、つまり、
\[ \epsilon_4 \lesssim C_R \]
と希望しているとしよう.
すると、現状と希望はそれぞれ以下のように要約できる.
現状: $\, \| U - \hat{U} \| \cong C\, (\Delta t)^5$,
希望: $\, C\, (\Delta t_{\rm opt} )^5 \lesssim C_R$
さらにこの二つの式から未知の定数 $C$ を消去すると、既知の値ばかりで右辺を構成できる.それが次の式だ.
理想的な時間ステップの推測値:
\[
\Delta t_{\rm opt} \lesssim \left( \frac{C_R}{\| U - \hat{U} \|} \right)^{1 / 5} \, \Delta t
\]
そこで、この値の範囲に収まるように $\Delta t$ をどう変えるか、という方針というかアルゴリズム2を考えれば良い.
あとはこの基本方針にのっとって $\Delta t$ を調整しつつ、計算を進めていけば良い.
ちなみに、次のステップの近似解として 4次近似値 $U$ と 5次近似値 $\hat{U}$ のどちらを採用するかであるが、「高次の方が良いんだろうからそちらを採用する」としたいのが通常だろう.実際、他の多くの解法はそう想定して設計されている.
しかし、この Fehlberg 公式は 4次近似が良くなるように片寄って設計されているため、原理的には、この Fehlberg 公式では次のステップの近似解として 4次近似値 $U$ を採用したほうが良いだろう.
5(4)次 Runge-Kutta 法 (Dormand-Prince 法)
Fehlberg 公式とは逆に、次のステップの近似式としてより高次な近似値である 5次近似値 $\hat{U}$ を採用した方が良いように素直に設計されているのが Dormand-Prince 法(wikipedia) である. パラメータ等はさすがに省略するが、基本的な考え方は上の Fehlberg 法と同じである.
$\Delta t$ を調整しながら実際の問題の計算に使うなら、Fehlberg 法よりもこちらの Dormand-Prince 法の方が良いだろう.
レポート
下記要領でレポートを出してみよう.
- e-mail にて,
- 題名を 2020-numerical-analysis-report-06 として,
- 教官宛(アドレスは web の "TOP" を見よう)に,
- 自分の学籍番号と名前を必ず書き込んで,
- 内容はテキストでも良いし,pdf などの電子ファイルを添付しても良いので,
下記の内容を実行して,結果や解析,感想等をレポートとして提出しよう.
- 上で紹介した Runge-Kutta-Fehlberg 法か Dormand-Prince 法を用いてみよう.
ちょっと面倒だけれども、上に紹介した Runge-Kutta-Fehlberg 法か Dormand-Prince 法を実際にプログラムして、今回の問題に適用してみよう. ステップ毎に $\Delta t$ を調整するようにプログラムを組めれば理想的だが、難しいようならば $\Delta t$ は固定値のままでもまあ良い.
ただし、パラメータの転記ミス等があると困るので、 英語 wikipedia の List of Runge-Kutta methods の該当する箇所を見てパラメータを確認しつつ、プログラムを進めよう. - 上の問題 1. で組んだプログラムが他の解法と比べてどれくらいの時間がかかるのか,授業で紹介した
BenchmarkTools
を使って測定しよう.
そしてその結果が,Runge-Kutta-Fehlberg 法等の計算式から予測されるものと整合するか評価してみよう.