06. 時間対称スキーム, Runge-Kutta法
Photo by James Cheung on Unsplash
単純な常微分方程式を Euler 法「以外」で計算しようとすると、一般にはどんな技術が必要になるだろうか？:
Euler 法は簡単で扱いやすいんだけれども、パラメータを少し変えると不安定になってしまっていた． しかも、詳しく示してはいないけれども、一般に近似精度もあまり良くはないんだ．
そこで、もう少し良い近似解法が考えられないか、検討してみよう．
対象の微分方程式
対象として今回も前回扱ったものと同じ、人口 $u = u(t)$ ($t$: 時間) に対する常微分方程式
\begin{equation} \frac{du}{dt} = \gamma\, \left( \sin(u) \right)^{1.2}, \,\, \gamma > 0, t \in (0,\infty) \end{equation}
を考えよう． ただし、$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$ などで近似するという感じだ．
今回の対象としている微分方程式にこうしたシンプルな対称な近似を行ってみると、例えば
\begin{equation} \frac{u_{n+1} - u_n}{\Delta t} = \gamma \left( \sin( \frac{ u_{n+1} + u_n}{2} )\right)^{1.2} \end{equation}
というような近似式が得られる．
この近似式を採用したとしよう． すると数値計算のプロセスでは $u_n$ は既知で $u_{n+1}$ は未知量なので、この未知量をどう求めたら良いか、を考える必要がある．
そこでより分かりやすいように、 未知量を $u_{n+1} = u$, 既知量を $u_n = U$ と表記して上の近似式を整理しよう． すると、関数 $f(u, U)$ を
\begin{equation} f(u, U) := u - U - \Delta t \gamma\, \left( \sin( \frac{ u + U }{2} )\right)^{1.2} \end{equation}
と定義すれば、上の近似式の数値計算での「使い方」は、
与えられた $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$ に関する微分を求めよう． これは簡単で、
\begin{equation} \frac{df}{du} = 1 - 0.6 \,\Delta t \,\gamma \,\sin(\frac{u+U}{2})^{0.2} \, \cos(\frac{u+U}{2}) \end{equation}
となる．
これで Newton 法を適用できる． まず，準備段階は以下のようになるだろう．
1# 定数の設定
2γ = 1.0
3u0 = 0.1
4Δt = 0.1
5
6# 関数 f(u,U)
7f(u, U) = u - U - Δt * γ * (sin((u+U)/2))^(1.2)
8
9# 微分 df/du
10df(u, U) = 1 -
11 0.6 * Δt * γ * (sin((u+U)/2))^(0.2) * cos((u+U)/2)
そして，Newton 法で古い時間ステップの関数値 $U$ から新しい時間ステップの関数値 $u$ を求める関数を次のように定義する，という感じになるね．
1# Newton 法で、U から u を求める
2function newton(old_u)
3 u = old_u # とりあえず、初期値を古い値の U と一致させておいて…
4
5 for i in 1:100
6 if abs( f(u,old_u) ) < 1.0e-6 # 6桁あってれば
7 break # ループ(この場合は for )から脱出する．
8 end
9 u = u - f(u,old_u)/df(u,old_u) # Newton 法で u を改善する．
10 end
11
12 return u # f が十分にゼロに近くなったら、近似値 u を返す
13end
これで準備が整ったので、実際に数値計算をして時間発展をさせてみよう．
1u = u0 # 「今」の値は、もちろん最初は初期値そのもの．
2u_sq = [ u ] # 初期値を配列に入れておいて…
3
4for n in 1:100
5 u = newton(u) # Newton 法で新しい時間ステップの値が求まって出てくるはず．
6 push!(u_sq, u) # その値を配列に追加していく．
7end
計算がスムーズに終わっているようなので、早速、グラフにしてみてみよう．
1using Plots
2
3t_sq = [ n*Δt for n in 0:100 ]
4plot(t_sq, u_sq, marker = :circ,
5 xaxis = "time t", yaxis = "u", leg = false)
どうやら無事に計算がうまくいっているようだ．まずはめでたしめでたし．
Euler 法は不安定だったけどこの「対称なスキーム」はどうかな？
前回やってみたように、パラメータによっては Euler 法はふっとんでしまい、不安定だった． この数値スキームはどうだろう？ 試してみよう．
Euler スキームが動かなくなったパラメータ
1Δt = 1.2
を設定して、動かして、プロットしてみる．
1u = u0 # 「今」の値は、もちろん最初は初期値そのもの．
2u_sq = [ u ] # 初期値を配列に入れておいて…
3
4for n in 1:100
5 u = newton(u) # Newton 法で新しい時間ステップの値が求まって出てくるはず．
6 push!(u_sq, u) # その値を配列に追加していく．
7end
8
9t_sq = [ n*Δt for n in 0:100 ]
10plot(t_sq, u_sq, marker = :circ,
11 xaxis = "time t", yaxis = "u", leg = false)
おお、これはきちんと動作するぞ． どうやら、ちょっとばかり手間だけれども、「対称なスキーム」は Euler 法より安心して使えそうだ．
Newton法の代わりにライブラリに頼ってみる！
復習を兼ねて、以前学習した、NLsolve パッケージを使ってみよう． まず、このパッケージを毎回行うルーチン的な手続きは以下の通りだったので、そうしよう．
1using NLsolve
2
3function nls(func, params...; ini = [0.0])
4 if typeof(ini) <: Number
5 r = nlsolve((vout,vin)->vout[1]=func(vin[1],params...), [ini])
6 v = r.zero[1]
7 else
8 r = nlsolve((vout,vin)->vout .= func(vin,params...), ini)
9 v = r.zero
10 end
11 return v, r.f_converged
12end
あとは前やったように、以下のようにすれば良いはずだ．
1Δt = 0.1 # 時間刻み幅を元に戻しておいて、
2
3u_old = u0 # 「今」の値は、もちろん最初は初期値そのもの．
4u_sq = [ u_old ] # 初期値を配列に入れておいて…
5
6for n in 1:100
7 u = nls(f, u_old, ini = u_old)[1] # 関数 nls で新しい値をいきなり求める．
8 push!(u_sq, u) # その値を配列に追加していく．
9 u_old = u # 古い値を保持する変数に新しい値を入れて、次の時間へ
10end
ふむ、無事に動いたようだ．グラフにしてみてみよう．
1t_sq = [ n*Δt for n in 0:100 ]
2plot(t_sq, u_sq, marker = :circ,
3 xaxis = "time t", yaxis = "u", leg = false)
グラフを見ると、上の例とまったく同じようだ．内容的にも問題ないと見て良いだろう．
念のために、ライブラリを使った場合でもパラメータを変えてみよう
先と同じように、$\Delta t = 1.2$ として挙動を見てみよう．
1Δt = 1.2 # 時間刻み幅を大きくしてみる
2
3u_old = u0 # 「今」の値は、もちろん最初は初期値そのもの．
4u_sq = [ u_old ] # 初期値を配列に入れておいて…
5
6for n in 1:100
7 u = nls(f, u_old, ini = u_old)[1] # 関数 nls で新しい値をいきなり求める．
8 println(n, ": ", u, " ", sin((u+u_old)/2))
9 push!(u_sq, u) # その値を配列に追加していく．
10 u_old = u # 古い値を保持する変数に新しい値を入れて、次の時間へ
11end
…略…
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(x::Float64)
@ Base.Math .\math.jl:41
[2] ^
@ .\math.jl:1206 [inlined]
[3] f(u::Float64, U::Float64)
@ Main .\In[23]:7
…略…
Euler 法よりはずいぶん先まで計算できるという意味でかなりマシだけれども、この場合も途中で Euler法とまったく同じ理由で計算が停止してしまうことがわかる． ライブラリが万能ではない，ということだね．
このように、簡単な問題でも、手法次第で結果が得られるかどうかが容易に変わってしまう ので、現実には計算手法を複数知っておいてそれをうまく切り替えて使う必要がありそうだ．柔軟にやるしかないだろうね．
Runge-Kutta 法 (時間発展計算では定番)
さて，さらに新しい解法を紹介しよう．
ここで紹介する Runge-Kutta 法とは常微分方程式
\begin{equation} \frac{du}{dt} = r(u,t) \end{equation}
に対して $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$ として、
\begin{equation} \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. \end{equation}
と順番に $r_k$ を計算する． こうやって，$t$ から $t + Δt$ の間の $u$ の 時間変化 $r (= du/dt)$ の近似値を 4つ計算するのだ．
そしてこれらを用いて、
\begin{equation} u(t + \Delta t) \, { の近似値} := u + \Delta t \left( \frac{r_1 + 2r_2 + 2r_3 + r_4}{6} \right) \end{equation}
として、時間ステップを 1つ進める近似解法がこの方法だ．
早速、プログラムしてみよう． まず、この問題の場合は、常微分方程式の右辺 $r(u,t)$ は次のようになる．
この問題の場合右辺に $t$ が登場しないので，以下のプログラムでは右辺の $t$ に関連する計算を省略しているよ！
1r(u) = (sin(u))^(1.2)
そして、Runge-Kutta 法で 1ステップすすめる関数をプログラムで書くと次のような感じだな．
1function RK(u)
2 r1 = r(u)
3 r2 = r(u + Δt/2 * r1)
4 r3 = r(u + Δt/2 * r2)
5 r4 = r(u + Δt * r3)
6 return u + Δt * (r1 + 2*r2 + 2*r3 + r4)/6
7end
早速使ってみよう！ まずは安全そうなパラメータで．
1Δt = 0.1
2u0 = 0.1
3
4u = u0 # 最初の値はもちろん初期値
5u_sq = [ u ] # 初期値を配列に入れておいて…
6
7for n in 1:100
8 u = RK(u) # Runge-Kutta 法で新しい値をいきなり求める．
9 push!(u_sq, u) # その値を配列に追加していく．
10end
で、動いているようなのでプロットだ．
1t_sq = [ n*Δt for n in 0:100 ]
2plot(t_sq, u_sq, marker = :circ,
3 xaxis = "time t", yaxis = "u", leg = false)
ふむふむ、うまくいっているな．
では、$\Delta t$ を大きくしてみて、うまくいくかチェックしよう．
1Δt = 1.2 # この値を大きくした．
2u0 = 0.1
3
4u = u0 # 最初の値はもちろん初期値
5u_sq = [ u ] # 初期値を配列に入れておいて…
6
7for n in 1:100
8 u = RK(u) # Runge-Kutta 法で新しい値をいきなり求める．
9 push!(u_sq, u) # その値を配列に追加していく．
10end
今回も大丈夫なようなので、プロットしよう．
1t_sq = [ n*Δt for n in 0:100 ]
2plot(t_sq, u_sq, marker = :circ,
3 xaxis = "time t", yaxis = "u", leg = false)
どうやら、$\Delta t$ を大きくしてもこの問題では Runge-Kutta 法がうまく動作するようだ．
困ったらとりあえず Runge-Kutta 法で と覚えておくとよさそうな気がするね．
それぞれの方法の計算速度の比較をしてみよう
準備として、前回同様に初期値と計算回数を入力すると最後の解を出力する形で関数を作ろう．
まずは Euler 法だ． まず、前回の授業で定義した Euler 関数を再度以下のように定義して、
1function Euler(u)
2 return u + Δt * γ * (sin(u))^1.2
3end
やはり前回のように，以下のように関数をつくる．
1function Euler_all(init_u, N)
2 u = init_u # 最初の値
3 for n in 1:N
4 u = Euler(u) # Euler スキームで新しい値を求める．
5 end
6 return u
7end
次に、対称なスキームについてだ．非線形部分の解法にはライブラリを用いたものとしよう．
1function Symmetric_all(init_u, N)
2 u = init_u # 最初の値
3 for n in 1:N
4 u = nls(f, u, ini = u)[1] # 関数 nls で新しい値を求める．
5 end
6 return u
7end
それから Runge-Kutta 法についても以下のように同様につくろう．
1function RK_all(init_u, N)
2 u = init_u # 最初の値
3 for n in 1:N
4 u = RK(u) # Runge-Kutta 法で新しい値を求める．
5 end
6 return u
7end
念のために、どれもきちんと動くことを確認しよう．
1Δt = 0.1 # この値で試すことにしよう
2
3Euler_all(u0,100)
3.119420135870497
1Symmetric_all(u0, 100)
3.1184787108216545
1RK_all(u0, 100)
3.118421279457686
ほぼ同じ最終値1が出力されていることから見ても、無事に動作していることがわかる．
では、これを使って計算時間を比較しよう．
それには前回出てきた BenchmarkTools パッケージの @benchmark マクロを使おう．
BenchmarkTools を使って計算時間などを測定する！
さて早速，このマクロを使って計算速度を測定、比較してみよう． まずパッケージの利用宣言を以下のようにして、
1using BenchmarkTools
あとは測定したい関数を @benchmark マクロに渡すだけだ．
Euler法は前回計測したから，対称スキームと Runge-Kutta スキームの計算速度を測定しよう．
1@benchmark Symmetric_all(u0, 100)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min … max): 313.300 μs … 145.925 ms ┊ GC (min … max): 0.00% … 99.67%
Time (median): 327.500 μs ┊ GC (median): 0.00%
Time (mean ± σ): 374.796 μs ± 1.485 ms ┊ GC (mean ± σ): 10.34% ± 7.57%
▄▄▇▆█▇▅▄▂▂▂▁▁ ▁▁▁▂▁▁ ▂
███████████████▇██▇▆▆▆▆▅▄▅▁▄▃▄▅▄▄▃▅▁▅▇███████▇█▇█▇▆▇▆▅▅▅▅▅▅▁▅ █
313 μs Histogram: log(frequency) by time 499 μs <
Memory estimate: 340.64 KiB, allocs estimate: 9001.
1@benchmark RK_all(u0, 100)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min … max): 46.300 μs … 144.947 ms ┊ GC (min … max): 0.00% … 99.94%
Time (median): 60.700 μs ┊ GC (median): 0.00%
Time (mean ± σ): 73.300 μs ± 1.449 ms ┊ GC (mean ± σ): 19.86% ± 1.36%
▂ ▆▇▆▅▂ ▁ ▃▆▅▇██▇▆▅▄▃▃▂▂▁▁▁▁ ▁ ▁▁▁▁ ▃
█████████▇▇▆▇▇▄▅▅▃▁▁▃▁▃▆████████████████████████████████▇▆▆▇ █
46.3 μs Histogram: log(frequency) by time 76.8 μs <
Memory estimate: 35.95 KiB, allocs estimate: 2301.
すると、前回の結果とあわせて，この3つの解法では mean値 で比較するとおおよそ
| 解法 | mean値(μs) |
|---|---|
| Euler法 | 11 |
| 対称スキーム | 375 |
| Runge-Kutta スキーム | 73 |
となることがわかる．
計算式を眺めると Runge-Kutta スキームは Euler スキームの 4倍の計算時間がかかるはず． 実際に見てみると，Runge-Kutta スキームはそれ以上に時間がかかっているので(おそらく，GC に時間が取られている)，「理論による予測だけではなく，実際に測定してみることが重要だ」ということがよくわかる．
また, 非線形部分を毎回解かないといけない対称スキームは，容易に予想できるように，計算時間が大変大きいこともよく分かる．
今回は素直な結果になったが、プログラムが要する計算時間は直感による予想と異なることが多々ある．計算時間については、なるべくこのように「測定した値」をベースにして議論するほうが良い．
近似解の精度も見てみよう！
前回同様，time $t = 10$ での近似解値: $$u_{approx}(t=10) = 3.1184213196530965$$ を厳密解の代わりにして，近似解との差を誤差とみなしてそれぞれの解法を評価してみよう．
まず，様々な $\Delta t$ に対して計算して $t = 10$ のときの誤差がどうなるのか，記録を取ろう． 少し長いが，下記のようなプログラムになるだろう．
1# 全体条件
2T = 10.0 # この時間の…
3u_exact = 3.1184213196530965 # 厳密解もどき．
4
5ratio = 1.5 # 計算回数を毎回 1.5倍する <=> Δt を毎回 2/3 に．
6N_iter = 15 # この試行を何回やるか．
7
8# 初期化
9N = 100
10Δt = T/N
11u0 = 0.1
12
13# 最初に一回計算しておく．
14e_euler = Euler_all(u0,N) - u_exact
15e_sym = Symmetric_all(u0,N) - u_exact
16e_rk = RK_all(u0,N) - u_exact
17
18println("iter, N, Δt, e_euler, e_sym, e_rk") # 画面表示
19println("0, $N, $Δt, $e_euler, $e_sym, $e_rk")
20
21# 記録をとりはじめる．
22sq_dt = [ Δt ]
23sq_euler = [ e_euler ]
24sq_sym = [ e_sym ]
25sq_rk = [ e_rk ]
26
27# Δt を小さくしていったときの近似解をそれぞれ求める．
28for i in 1:N_iter
29 N = round(Int, ratio * N) # N を ratio 倍して，
30 Δt = T / N # Δt は 1/ratio 倍になるはず．
31 push!( sq_dt, Δt ) # 記録をとる
32
33 e_euler = Euler_all(u0,N) - u_exact # それぞれ計算
34 e_sym = Symmetric_all(u0,N) - u_exact
35 e_rk = RK_all(u0,N) - u_exact
36
37 push!(sq_euler, e_euler) # 記録をとる
38 push!(sq_sym, e_sym)
39 push!(sq_rk, e_rk)
40
41 println("$i, $N, $Δt, $e_euler, $e_sym, $e_rk") # 画面表示
42end
すると次のように次々と表示され，だんだん誤差が小さくなっていることや，それぞれの解法の違いがうっすらと見えてくる．
iter, N, Δt, e_euler, e_sym, e_rk
0, 100, 0.1, 0.000998816217400389, 5.739116855796311e-5, -4.0195410555554645e-8
1, 150, 0.06666666666666667, 0.0006564740322243523, 2.5518718715478172e-5, -7.922887679256974e-9
2, 225, 0.044444444444444446, 0.0004334550057212283, 1.1257736344472136e-5, -1.5630248206832675e-9
3, 338, 0.029585798816568046, 0.0002866709631739184, 4.869136503327098e-6, -3.0667823835983654e-10
4, 507, 0.01972386587771203, 0.00019028442774038368, 2.0658655817129556e-6, -6.054268197885904e-11
5, 760, 0.013157894736842105, 0.0001265708703082069, 8.775914257519446e-7, -1.198108279254484e-11
6, 1140, 0.008771929824561403, 8.421618377019158e-5, 3.9794948358817805e-7, -2.361666417982633e-12
7, 1710, 0.005847953216374269, 5.607101881643928e-5, 1.8996327044717987e-7, -4.60964599824365e-13
8, 2565, 0.003898635477582846, 3.734817716161132e-5, 1.1463056415195183e-7, -8.393286066166183e-14
9, 3848, 0.002598752598752599, 2.4881098828366532e-5, 5.093337351880223e-8, -1.021405182655144e-14
10, 5772, 0.0017325017325017325, 1.65809781838e-5, 2.2636969632117143e-8, 5.329070518200751e-15
11, 8658, 0.001155001155001155, 1.1051131361128341e-5, 1.0060840160974749e-8, 7.105427357601002e-15
12, 12987, 0.00077000077000077, 7.366152340981813e-6, 4.47147741056142e-9, 5.329070518200751e-15
13, 19480, 0.000513347022587269, 4.910330449359179e-6, 1.9874164536304306e-9, 6.217248937900877e-15
14, 29220, 0.00034223134839151266, 3.2733030281839604e-6, 8.832916620349351e-10, 1.0658141036401503e-14
15, 43830, 0.00022815423226100844, 2.1820906250802352e-6, 3.925770819535046e-10, 7.105427357601002e-15
もう少し状況を把握するために，グラフにしてみよう．
グラフの傾きを見やすくするために以下のように少し準備をして，
1a1 = 2.0e-2
2a2 = 2.0e-2
3a4 = 2.0e-3
4p1(t) = a1 * t
5p2(t) = a2 * t^2
6p4(t) = a4 * t^4
下記のようにグラフを描こう．
1default( marker = :auto, xaxis = ("Δt", :log10),
2 yaxis = ("error", :log10), legend = :outertopright )
3# 共通の設定は default 宣言で一回だけすれば良い．
4
5plot(sq_dt, abs.(sq_euler), label = "Euler")
6plot!(sq_dt, abs.(sq_sym), label = "Symmetric")
7plot!(sq_dt, abs.(sq_rk), label = "Runge-Kutta")
8
9default(marker = :none)
10
11plot!(sq_dt, p1.(sq_dt), label = "y = Cx")
12plot!(sq_dt, p2.(sq_dt), label = "y = Cx^2")
13plot!(sq_dt, p4.(sq_dt), label = "y = Cx^4")
このグラフを見て容易に以下のことがわかる．
- Euler 法よりも対称スキームのほうが誤差が小さく，さらに，Runge-Kutta法はもっと誤差が小さい．
- Euler 法では，$\Delta t$ と誤差 $\epsilon$ が比例関係にある．つまり，
$$\epsilon \propto \Delta t$$
これは $\Delta t$ を $1/2$ にすると誤差が $1/2$ になることを意味する.
つまり，計算量を倍にすると誤差が $1/2$ になるということだ．
ちなみに，こういう関係のある数値計算法を「1次精度の解法である」と呼ぶ． - 対称スキームでは $\Delta t$ と誤差 $\epsilon$ が二次関係にある．つまり，
$$\epsilon \propto \Delta t^2$$
これは $\Delta t$ を $1/2$ にすると誤差が $1/4$ になることを意味する.
つまり，計算量を倍にすると誤差が $1/4$ になるということなので，Euler 法よりも効率が良い．
こういう関係のある数値計算法を「2次精度の解法である」と呼ぶ． - 古典的 Runge-Kutta 法では $\Delta t$ と誤差 $\epsilon$ が 4次関係にある．つまり，
$$\epsilon \propto \Delta t^4$$
これは $\Delta t$ を $1/2$ にすると誤差が $1/16$ になることを意味する．
つまり，計算量を倍にすると誤差が $1/16$ になるということなので，Euler 法や対称スキームよりもさらに効率が良い．
こういう関係のある数値計算法を「4次精度の解法である」と呼ぶ．
これらからわかるのは，Runge-Kutta法などの結構複雑な解法は，誤差も小さいし，計算量を増やしたとき2の効率も良いので，本格的な計算をするのであれば積極的に利用していくべきだ，ということだね．
実際の計算時間も考慮に入れたいよね？
上のグラフは計算量を左右する $Δt$ と計算誤差の対応がはっきりわかる良いグラフだった． しかし，解法によって時間ステップを 1進めるための計算量は違っていたので，計算時間と誤差の関係はこのグラフからズレてくるはず．
そこで，実際の計算時間も測定してみようではないか．
というわけで，計算時間を測定するために @timed マクロを使用する．
使い方は簡単なのでヘルプを呼び出して読もう．
今回は，誤差と計算時間の両方をデータにしたいので，次のように関数を作って使うことにしよう．
1function error_and_time(scheme, N)
2 error, t = @timed (scheme(u0,N) - u_exact) # 値を求めて誤差の計算
3 return (error, t) # 誤差と計算時間を返す
4end
これで，例えば Eulerスキームを 100step 回したときの結果は下のようになるわけだ．
1# 間違いの無いように，念の為パラメータをあらためて設定しておく．
2T = 10.0
3N = 100
4Δt = T/N
5u_exact = 3.1184213196530965 # 厳密解もどき．
6
7error_and_time(Euler_all, 100) # こうすると，結果は下のように出るぞ．
(0.000998816217400389, 3.34e-5) # 誤差と計算時間
これを使って，誤差と計算時間の関係を探ろう． プログラムは以下のような感じになるだろう．
1# 全体条件
2T = 10.0 # この時間の…
3u_exact = 3.1184213196530965 # 厳密解もどき．
4
5ratio = 1.5 # 計算回数を毎回 1.5倍する <=> Δt を毎回 2/3 に．
6N_iter = 15 # この試行を何回やるか．
7
8# 初期化
9N = 100
10Δt = T/N
11u0 = 0.1
12
13# 最初に一回計算しておく．
14e_euler, t_euler = error_and_time(Euler_all, N)
15e_sym, t_sym = error_and_time(Symmetric_all, N)
16e_rk, t_rk = error_and_time(RK_all, N)
17
18println("iter, N, Δt") # 画面表示
19println("0, $N, $Δt")
20
21# 記録をとりはじめる．
22sq_dt = [ Δt ]
23sq_e_euler = [ e_euler ]
24sq_e_sym = [ e_sym ]
25sq_e_rk = [ e_rk ]
26sq_t_euler = [ t_euler ]
27sq_t_sym = [ t_sym ]
28sq_t_rk = [ t_rk ]
29
30# Δt を小さくしていったときの近似解をそれぞれ求める．
31for i in 1:N_iter
32 N = round(Int, ratio * N) # N を ratio 倍して，
33 Δt = T / N # Δt は 1/ratio 倍になるはず．
34 push!( sq_dt, Δt ) # 記録をとる
35
36 e_euler, t_euler = error_and_time(Euler_all, N)
37 e_sym, t_sym = error_and_time(Symmetric_all, N)
38 e_rk, t_rk = error_and_time(RK_all, N)
39
40 push!(sq_e_euler, e_euler) # 記録をとる
41 push!(sq_e_sym, e_sym)
42 push!(sq_e_rk, e_rk)
43
44 push!(sq_t_euler, t_euler)
45 push!(sq_t_sym, t_sym)
46 push!(sq_t_rk, t_rk)
47
48 println("$i, $N, $Δt") # 画面表示
49end
iter, N, Δt
0, 100, 0.1
1, 150, 0.06666666666666667
2, 225, 0.044444444444444446
3, 338, 0.029585798816568046
4, 507, 0.01972386587771203
5, 760, 0.013157894736842105
6, 1140, 0.008771929824561403
7, 1710, 0.005847953216374269
8, 2565, 0.003898635477582846
9, 3848, 0.002598752598752599
10, 5772, 0.0017325017325017325
11, 8658, 0.001155001155001155
12, 12987, 0.00077000077000077
13, 19480, 0.000513347022587269
14, 29220, 0.00034223134839151266
15, 43830, 0.00022815423226100844
そして，得られた結果をプロットする. 準備も含めて次のような感じにすれば良さそうだ．
1n1(t) = 1.0e-7 / t
2n2(t) = 1.0e-10 / (t^2)
3n4(t) = 2.0e-22 / (t^4)
4
5t_range_euler = minimum(sq_t_euler):0.0001:maximum(sq_t_euler)
6t_range_sym = minimum(sq_t_sym):0.0001:maximum(sq_t_sym)
7t_range_rk = minimum(sq_t_rk):0.0001:maximum(sq_t_rk)
8
9default( marker = :auto, xaxis = ("time (sec)", :log10),
10 yaxis = ("error", :log10), legend = :outertopright )
11
12plot(sq_t_euler, abs.(sq_e_euler), label = "Euler")
13plot!(sq_t_sym, abs.(sq_e_sym), label = "Symmetric")
14plot!(sq_t_rk, abs.(sq_e_rk), label = "Runge-Kutta")
15
16default(marker = :none)
17
18plot!(t_range_euler, n1.(t_range_euler), label = "y = C/x")
19plot!(t_range_sym, n2.(t_range_sym), label = "y = C/(x^2)")
20plot!(t_range_rk, n4.(t_range_rk), label = "y = C/(x^4)")
さて，このプロットを見るとわかるのは以下のような事柄だね．
- Runge-Kutta法が速くて誤差が小さく，明らかに他の解法より良い．
- 遅いと言っていた対称スキームだが，計算時間と誤差のバランスを勘案すると計算時間 $\geq 10^{-3}$ (← $\Delta t$ を細かくして誤差を小さくしようとしている状況)では Euler スキームよりマシと言って良さそうだ．
- Euler法, 対称スキーム, Runge-Kutta法は，それぞれ確かに 1次精度，2次精度, 4次精度だと思われる．
というわけで，Runge-Kutta法のような「高性能な数値解法」には大いに価値があることがよく分かる．
もう少し複雑な 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$ として、
\begin{equation} \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. \end{equation}
と上から順に計算して、$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$ の誤差を
\begin{equation} \epsilon_4 \cong \| U - \hat{U} \| \end{equation}
として、具体的に数値で見積もれることになる.
ベクトルの評価にノルムを使っているよ．
そして、そもそも $p$ 次近似値の局所誤差(1ステップだけで発生する誤差) $\epsilon_p$ は
\begin{equation} \epsilon_p \cong C\, (\Delta t)^{p+1} \end{equation}
となるように上の Fehlberg 公式は設計されている($C$ は時間と $u$ 等に依存する適当な係数)．
時間発展問題の近似式で「$p$ 次局所誤差」という場合は、通常はこの意味だ．
そして、この状況下で、$\Delta t$ を $\Delta t_{\rm opt}$ に修正して局所誤差を $C_R$ 程度以下に抑えたい、つまり、
\begin{equation} \epsilon_4 \lesssim C_R \end{equation}
と希望しているとしよう．
すると、現状と希望はそれぞれ以下のように要約できる．
現状: $\, \| U - \hat{U} \| \cong C\, (\Delta t)^5$,
希望: $\, C\, (\Delta t_{\rm opt} )^5 \lesssim C_R$
さらにこの二つの式から未知の定数 $C$ を消去すると、既知の値ばかりで右辺を構成できる．それが次の式だ．
理想的な時間ステップの推測値:
\begin{equation}
\Delta t_{\rm opt} \lesssim \left( \frac{C_R}{\| U - \hat{U} \|} \right)^{1 / 5} \, \Delta t
\end{equation}
そこで、この値の範囲に収まるように $\Delta t$ をどう変えるか、という方針というかアルゴリズム3を考えれば良い．
あとはこの基本方針にのっとって $\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) である． パラメータは webを参照してもらうとして、基本的な考え方は上の Fehlberg 法と同じである．
ちなみに，Dormand-Prince 法の5次スキームを実装すると次のような感じだ． なお，無駄な計算量がなるべく減るような書き方をしている4．
1function DP(u)
2 r1 = r(u)*Δt
3 r2 = r(u + (1/5)*r1 )*Δt
4 r3 = r(u + (3/40)*r1 + (9/40)*r2 )*Δt
5 r4 = r(u + (44/45)*r1 - (56/15)*r2 + (32/9) *r3 )*Δt
6 r5 = r(u + (19372/6561) *r1 - (25360/2187) *r2 + (64448/6561) *r3 - (212/729) *r4 )*Δt
7 r6 = r(u + (9017/3168) *r1 - (355/33) *r2 + (46732/5247) *r3 + (49/176) *r4 - (5103/18656) *r5 )*Δt
8
9 return u + (35/384)*r1 + (500/1113)*r3 + (125/192)*r4 - (2187/6784)*r5 + (11/84)*r6
10end
$\Delta t$ を調整しながら実際の問題の計算に使うなら、Fehlberg 法よりもこちらの Dormand-Prince 法の方が良いだろう．
レポート
下記要領でレポートを出してみよう．
- e-mail にて，
- 題名を 2024-numerical-analysis-report-06 として,
- 教官宛(アドレスは web の "TOP" を見よう)に，
- 自分の学籍番号と名前を必ず書き込んで，
- 内容はテキストでも良いし，pdf などの電子ファイルを添付しても良いので，
下記の内容を実行して，結果や解析，感想等をレポートとして提出しよう．
注意
近年はセキュリティ上の懸念から，実行形式のプログラムなどをメールに添付するとそのメールそのものの受信を受信側サーバが拒絶したりする．
そういうことを避けるため，レポートをメールで提出するときは添付ファイルにそういった懸念のあるファイルが無いようにしよう．
まあ要するに，レポートは pdf ファイルにして，それをメールに添付して送るのが良い ということだと思っておこう．
- 上のプログラムを少し手直しして，$\Delta t$ を一定の比率で小さくしていったときの近似解を計算して近似解を並べた数列を作り，その数列に Aitken加速5 を(何度か重ねがけして)適用して数列の収束先を予測する方法で厳密解もどきを自分で求めてみよう．
このときに用いる近似解法としては対称スキームか Runge-Kutta法を用いるとよいだろう．
ちなみに，この方法でも厳密解もどきとして十分な値が求まるだろう．おそらく，$u_{\rm exact} \cong 3.1184213196531$ ぐらいの値が求まるだろう． - 上で紹介した Runge-Kutta-Fehlberg 法か Dormand-Prince 法の誤差 vs $\Delta t$ のグラフを描いてみて，誤差が $\Delta t$ と何次の関係にあるかを確認しよう．
なお，この問では $\Delta t$ の自動調整機能は使わず，単なる普通の解法として扱おう．
ちなみに、パラメータの転記ミス等があると困るので、 英語 wikipedia の List of Runge-Kutta methods の該当する箇所を見てパラメータを確認しつつ、プログラムを進めよう． - 上の問題 2. で組んだプログラムが他の解法と比べてどれくらいの時間がかかるのか，授業で紹介した
BenchmarkToolsを使って測定しよう． - (チャレンジ問題) 上で紹介した Runge-Kutta-Fehlberg 法か Dormand-Prince 法の本来の使い方である $\Delta t$ の値の自動制御 を試みてみよう．
なお，この問題はなかなか難しいと思うので，無理に取り組まなくても良い．
-
最終結果が微妙に異なるのは誤差のせいだね．どれの誤差が大きいか，はまた別の議論が必要だね． ↩︎
-
同じ time t まで計算する状況下では，$\Delta t$ を小さくすることは全体の計算回数をそれに反比例して多くすることを意味する． ↩︎
-
あまり急に大きく値を修正すると良くないことが起きそうだし，多少は安全係数も掛けたいし，考えどころだね． ↩︎
-
含まれる多数の係数が分数なので，これをいったん先に実数に直しておけばもっと速いはず等々思うだろうが，コンパイラが賢くてそのあたりはよしなにやってくれる．むしろ，変数を先に計算して…なんてやるとかえって遅くなることも多い．気になる人は確かめてみると良い．賢いコンパイラよりうまくやるのはなかなか難しいぞ． ↩︎