06. 時間対称スキーム, Runge-Kutta法

Daisuke Furihata (降籏 大介)

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 法を適用できる． まず，準備段階は以下のようになるだろう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# 定数の設定
γ   = 1.0
u0   = 0.1
Δt  = 0.1

# 関数 f(u,U)
f(u, U) = u - U - Δt * γ * (sin((u+U)/2))^(1.2)

# 微分 df/du
df(u, U) = 1 - 0.6 * Δt * γ * (sin((u+U)/2))^(0.2) * cos((u+U)/2)

そして，Newton 法で古い時間ステップの関数値 $U$ から新しい時間ステップの関数値 $u$ を求める関数を次のように定義する，という感じになるね．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

# Newton 法で、U から u を求める
function newton(old_u)
    u = old_u # とりあえず、初期値を古い値の U と一致させておいて…
    
    for i in 1:100
      if abs( f(u,old_u) ) < 1.0e-6   # 6桁あってれば
        break   # ループ(この場合は for )から脱出する．
      end
      u = u - f(u,old_u)/df(u,old_u) # Newton 法で u を改善する．
    end
    
    return u # f が十分にゼロに近くなったら、近似値 u を返す
end

これで準備が整ったので、実際に数値計算をして時間発展をさせてみよう．

1
2
3
4
5
6
7

u = u0       # 「今」の値は、もちろん最初は初期値そのもの．
u_sq = [ u ] # 初期値を配列に入れておいて…

for n in 1:100
    u = newton(u) # Newton 法で新しい時間ステップの値が求まって出てくるはず．
    push!(u_sq, u) # その値を配列に追加していく．
end

計算がスムーズに終わっているようなので、早速、グラフにしてみてみよう．

1
2
3
4

using Plots

t_sq = [ n*Δt for n in 0:100 ]
plot(t_sq, u_sq, xaxis = "time t", yaxis = "u", leg = false)

svg

どうやら無事に計算がうまくいっているようだ．まずはめでたしめでたし．

Euler 法は不安定だったけどこの「対称なスキーム」はどうかな？

前回やってみたように、パラメータによっては Euler 法はふっとんでしまい、不安定だった． この数値スキームはどうだろう？ 試してみよう．

Euler スキームが動かなくなったパラメータ

1

Δt = 1.2

を設定して、動かして、プロットしてみる．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

u = u0       # 「今」の値は、もちろん最初は初期値そのもの．
u_sq = [ u ] # 初期値を配列に入れておいて…

for n in 1:100
    u = newton(u) # Newton 法で新しい時間ステップの値が求まって出てくるはず．
    push!(u_sq, u) # その値を配列に追加していく．
end

t_sq = [ n*Δt for n in 0:100 ]
plot(t_sq, u_sq, xaxis = "time t", yaxis = "u", leg = false)

svg

おお、これはきちんと動作するぞ． どうやら、ちょっとばかり手間だけれども、「対称なスキーム」は Euler 法より安心して使えそうだ．

Newton法の代わりにライブラリに頼ってみる！

復習を兼ねて、以前学習した、NLsolve パッケージを使ってみよう． まず、このパッケージを毎回行うルーチン的な手続きは以下の通りだったので、そうしよう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

using NLsolve

function nls(func, params...; ini = [0.0])
    if typeof(ini) <: Number
        r = nlsolve((vout,vin)->vout[1]=func(vin[1],params...), [ini])
        v = r.zero[1]
    else
        r = nlsolve((vout,vin)->vout .= func(vin,params...), ini)
        v = r.zero
    end
    return v, r.f_converged
end

あとは前やったように、以下のようにすれば良いはずだ．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

Δt = 0.1 # 時間刻み幅を元に戻しておいて、

u_old = u0  # 「今」の値は、もちろん最初は初期値そのもの．
u_sq = [ u_old ] # 初期値を配列に入れておいて…

for n in 1:100
    u = nls(f, u_old, ini = u_old)[1]  # 関数 nls で新しい値をいきなり求める．
    push!(u_sq, u) # その値を配列に追加していく．
    u_old = u # 古い値を保持する変数に新しい値を入れて、次の時間へ
end

ふむ、無事に動いたようだ．グラフにしてみてみよう．

1
2

t_sq = [ n*Δt for n in 0:100 ]
plot(t_sq, u_sq, xaxis = "time t", yaxis = "u", leg = false)

svg

グラフを見ると、上の例とまったく同じようだ．内容的にも問題ないと見て良いだろう．

念のために、ライブラリを使った場合でもパラメータを変えてみよう

先と同じように、$\Delta t = 1.2$ として挙動を見てみよう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

Δt = 1.2 # 時間刻み幅を大きくしてみる

u_old = u0  # 「今」の値は、もちろん最初は初期値そのもの．
u_sq = [ u_old ] # 初期値を配列に入れておいて…

for n in 1:100
    u = nls(f, u_old, ini = u_old)[1]  # 関数 nls で新しい値をいきなり求める．
    println(n, ": ", u, " ", sin((u+u_old)/2))
    push!(u_sq, u) # その値を配列に追加していく．
    u_old = u # 古い値を保持する変数に新しい値を入れて、次の時間へ
end

…略…
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. \]

と順番に $r_k$ を計算する． こうやって，$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$ が登場しないので省略しているよ！

1

r(u) = (sin(u))^(1.2)

そして、Runge-Kutta 法で 1ステップすすめる関数をプログラムで書くと次のような感じだな．

1
2
3
4
5
6
7

function RK(u)
    r1 = r(u)
    r2 = r(u + Δt/2 * r1)
    r3 = r(u + Δt/2 * r2)
    r4 = r(u + Δt * r3)
    return u + Δt * (r1 + 2*r2 + 2*r3 + r4)/6
end

早速使ってみよう！ まずは安全そうなパラメータで．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

Δt = 0.1
u0 = 0.1

u = u0  # 最初の値はもちろん初期値
u_sq = [ u ] # 初期値を配列に入れておいて…

for n in 1:100
    u = RK(u)  # Runge-Kutta 法で新しい値をいきなり求める．
    push!(u_sq, u) # その値を配列に追加していく．
end

で、動いているようなのでプロットだ．

1
2

t_sq = [ n*Δt for n in 0:100 ]
plot(t_sq, u_sq, xaxis = "time t", yaxis = "u", leg = false)

svg

ふむふむ、うまくいっているな．

では、$\Delta t$ を大きくしてみて、うまくいくかチェックしよう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

Δt = 1.2 # この値を大きくした．
u0 = 0.1

u = u0  # 最初の値はもちろん初期値
u_sq = [ u ] # 初期値を配列に入れておいて…

for n in 1:100
    u = RK(u)  # Runge-Kutta 法で新しい値をいきなり求める．
    push!(u_sq, u) # その値を配列に追加していく．
end

今回も大丈夫なようなので、プロットしよう．

1
2

t_sq = [ n*Δt for n in 0:100 ]
plot(t_sq, u_sq, xaxis = "time t", yaxis = "u", leg = false)

svg

どうやら、$\Delta t$ を大きくしてもこの問題では Runge-Kutta 法がうまく動作するようだ．

困ったらとりあえず Runge-Kutta 法で と覚えておくとよさそうな気がするね．

それぞれの方法の計算速度の比較をしてみよう

準備として、初期値と計算回数を入力すると最後の解を出力する形で関数を作ろう．

まずは Euler 法だ． まず、前回の授業で定義した Euler 関数を再度以下のように定義して、

1
2
3

function Euler(u)
    return u + Δt * γ * (sin(u))^1.2
end

以下のように全体をつくる．

1
2
3
4
5
6
7

function Euler_all(init_u, N)
  u = init_u  # 最初の値
  for n in 1:N
    u = Euler(u)  # Euler スキームで新しい値を求める．
  end
  return u
end

次に、対称なスキームについてだ．非線形部分の解法にはライブラリを用いたものとしよう．

1
2
3
4
5
6
7

function Symmetric_all(init_u, N)
  u = init_u  # 最初の値
  for n in 1:N
    u = nls(f, u, ini = u)[1]  # 関数 nls で新しい値を求める．
  end
  return u
end

それから Runge-Kutta 法についても以下のように同様につくろう．

1
2
3
4
5
6
7

function RK_all(init_u, N)
  u = init_u  # 最初の値
  for n in 1:N
    u = RK(u)  # Runge-Kutta 法で新しい値を求める．
  end
  return u
end

念のために、どれもきちんと動くことを確認しよう．

1
2
3

Δt = 0.1 # この値で試すことにしよう

Euler_all(u0,100)

3.119420135870497

1

Symmetric_all(u0, 100)

3.1184787108216545

1

RK_all(u0, 100)

3.118421279457686

ほぼ同じ最終値1が出力されていることから見ても、無事に動作していることがわかる．

では、これを使って計算時間を比較しよう． それには、BenchmarkTools というパッケージに入っている、@benchmark というマクロが便利だ．

おそらくみなさんの Julia 環境にはこのパッケージは入っていないだろうから，まずはインストールだ．

1
2

using Pkg
Pkg.add("BenchmarkTools")

インストールが済んだら，早速このマクロを使って計算速度を測定、比較してみよう． まずパッケージの利用宣言を以下のようにして、

1

using BenchmarkTools

あとは測定したい関数を @benchmark マクロに以下のように渡すだけだ．

1

@benchmark Euler_all(u0, 100)

すると，次のように計算時間に関する簡単な統計情報を教えてくれる．

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  11.400 μs …  1.437 ms  ┊ GC (min … max): 0.00% … 98.89%
 Time  (median):     13.700 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   13.349 μs ± 14.336 μs  ┊ GC (mean ± σ):  1.06% ±  0.99%

  █▇▆▅▂       ▂▂       ▆▇▅▃▄▂▁   ▁▅▆▅▄▄▃▂▂                    ▂
  ███████▆▅▇▇▇███▅▅▅▃▄████████▆▆▇█████████▇▇▆▅▅▆▇▇███▇▇▇▇▅▆▁▅ █
  11.4 μs      Histogram: log(frequency) by time      17.6 μs <

 Memory estimate: 7.83 KiB, allocs estimate: 501.

結果を見ると、10000回のサンプルで, 平均が 13 μsぐらいの計算時間だということや，その他諸々がわかる． 真ん中のグラフは計算時間の頻度グラフだ． 通常は 11.4μs 以上 17.6μsで済むということなども見て取れるし， Jupyter の画面だとこのグラフに「青い線」が入っていて，おそらくそれは median値を意味する箇所だろう．

Benchmark 結果の中の GC は garbage collection (ゴミ処理) の頭文字で，使わなくなったメモリの整理等の作業を意味している． 結果をみると GC が最大で計算時間の 98.89% を食ってしまうケースが有ると言っているので，GC が発生するとかなり時間を食うことがわかる． PC で使えるメモリが少ないほど GC が頻繁に発生するので，メモリが少ないとつらいことがわかるね．

同様に、対称スキームと Runge-Kutta スキームの計算速度も測定しよう．

1

@benchmark Symmetric_all(u0, 100)

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  356.000 μs …   6.075 ms  ┊ GC (min … max): 0.00% … 93.77%
 Time  (median):     375.600 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   423.466 μs ± 408.961 μs  ┊ GC (mean ± σ):  9.65% ±  9.09%

  █▃▁                                                           ▁
  ████▇▅▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃ █
  356 μs        Histogram: log(frequency) by time        4.2 ms <

 Memory estimate: 476.58 KiB, allocs estimate: 10601.

1

@benchmark RK_all(u0, 100)

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  49.800 μs …  1.158 ms  ┊ GC (min … max): 0.00% … 94.47%
 Time  (median):     50.800 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   53.551 μs ± 30.590 μs  ┊ GC (mean ± σ):  1.58% ±  2.68%

  ▄▆█▅▄▆▃▁▁                          ▂▂▁▂▃▁▂▂▁                ▂
  █████████▇▅▅▄▅▃▃▁▄▃▁▃▁▁▁▃▁▁▁▁▄▄▆▇▇▇█████████▇▆▅▅▅▄▅▄▄▄▄▃▄▁▅ █
  49.8 μs      Histogram: log(frequency) by time      71.5 μs <

 Memory estimate: 39.08 KiB, allocs estimate: 2501.

すると、対称スキームでは計算時間が 380μs ほどで、Runge-Kutta スキームでは 50μs ほどだということもわかる．

計算式を眺めると Runge-Kutta スキームは Euler スキームの 4倍の計算時間がかかるはずなので、Runge-Kutta スキームのこの計算時間測定値は想定通りと言ってよさそうだ．

また, 非線形部分を毎回解かないといけない対称スキームは，まあ予想通り，計算時間が大変大きいこともよく分かる．

今回は素直な結果になったが、プログラムが要する計算時間は直感による予想と異なることが多々ある．計算時間については、なるべくこのように「測定した値」をベースにして議論するほうが良い．

近似解の精度も見てみよう！

time $t = 10$ でのこの問題の近似解として，値 $$u_{approx}(t=10) = 3.1184213196530965$$ を厳密解の代わりにしよう2．

そして，近似解とこの厳密解(もどき)との違いをおおよそ誤差としてみなすことにして，それぞれの解法での近似解の誤差を評価してみよう．

具体的には，異なる値の $\Delta t$ に対する誤差をグラフにプロットして様子を見るのだ． これでいろいろわかるはずだ．

まず，様々な $\Delta t$ に対して計算して $t = 10$ のときの誤差がどうなるのか，記録を取ろう． 少し長いが，下記のようなプログラムになるだろう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

# 全体条件
T = 10.0 # この時間の…
u_exact = 3.1184213196530965 # 厳密解もどき．

ratio  = 1.5 # 計算回数を毎回 1.5倍する <=> Δt を毎回 2/3 に．
N_iter = 15  # この試行を何回やるか．

# 初期化
N = 100
Δt = T/N
u0 = 0.1

# 最初に一回計算しておく．
e_euler = Euler_all(u0,N) - u_exact
e_sym   = Symmetric_all(u0,N) - u_exact
e_rk    = RK_all(u0,N) - u_exact

println("iter, N, Δt, e_euler, e_sym, e_rk")  # 画面表示
println("0, $N, $Δt, $e_euler, $e_sym, $e_rk")

# 記録をとりはじめる．
sq_dt = [ Δt ] 
sq_euler = [ e_euler ]
sq_sym   = [ e_sym ]
sq_rk    = [ e_rk ]

# Δt を小さくしていったときの近似解をそれぞれ求める．
for i in 1:N_iter
  N = round(Int, ratio * N) # N を ratio 倍して，
  Δt = T / N                # Δt は 1/ratio 倍になるはず．
  push!( sq_dt, Δt )        # 記録をとる
  
  e_euler = Euler_all(u0,N) - u_exact     # それぞれ計算
  e_sym   = Symmetric_all(u0,N) - u_exact
  e_rk    = RK_all(u0,N) - u_exact

  push!(sq_euler, e_euler)        # 記録をとる 
  push!(sq_sym, e_sym)
  push!(sq_rk, e_rk)

  println("$i, $N, $Δt, $e_euler, $e_sym, $e_rk")  # 画面表示
end

すると次のように次々と表示され，だんだん誤差が小さくなっていることや，それぞれの解法の違いがうっすらと見えてくる．

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

もう少し状況を把握するために，グラフにしてみよう．

グラフの傾きを見やすくするために以下のように少し準備をして，

1
2
3
4
5
6

a1 = 2.0e-2
a2 = 2.0e-2
a4 = 2.0e-3
p1(t) = a1 * t
p2(t) = a2 * t^2 
p4(t) = a4 * t^4

下記のようにグラフを描こう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

default( marker = :auto, xaxis = ("Δt", :log10), 
    yaxis = ("error", :log10), legend = :outertopright )
# 共通の設定は default 宣言で一回だけすれば良い．

plot(sq_dt, abs.(sq_euler), label = "Euler")
plot!(sq_dt, abs.(sq_sym), label = "Symmetric")
plot!(sq_dt, abs.(sq_rk), label = "Runge-Kutta")

default(marker = :none)

plot!(sq_dt, p1.(sq_dt), label = "y = Cx")
plot!(sq_dt, p2.(sq_dt), label = "y = Cx^2")
plot!(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法などの結構複雑な解法は，誤差も小さいし，計算量を増やしたとき3の効率も良いので，本格的な計算をするのであれば積極的に利用していくべきだ，ということだね．

計算時間も考慮に入れたいよね？

上のグラフは計算量を左右する $Δt$ と計算誤差の対応がはっきりわかる良いグラフだった． しかし，解法によって時間ステップを 1進めるための計算量は違っていたので，計算時間と誤差の関係はこのグラフからズレてくるはず．

そこで，実際の計算時間も測定してみようではないか．

というわけで，まずは計算時間を測定するための下準備として，1つ小さなパッケージをインストールしよう．

1
2

using Pkg
Pkg.add("TickTock")

この TickTock パッケージを使用すると，

  • コマンド tick(): 時間測定を開始する．
  • コマンド tok() : 経過時間を測定し，「値として返して」終了．
  • コマンド tock(): 経過時間を測定し，「表示して」終了．
  • コマンド laptimer(): 経過時間を測定し，「表示して」測定を続行．

というようなコマンドが使える．

今回は，誤差と計算時間の両方をデータにしたいので，次のように関数を作って使うことにしよう．

1
2
3
4
5
6
7
8

using TockTock  # このパッケージの利用宣言を忘れずに．

function error_and_time(scheme, N)
  tick() # 時間測定開始
  error = scheme(u0,N) - u_exact # 値を求めて誤差の計算
  t = tok() # 時間測定(秒単位)
  return (error, t) # 誤差と計算時間を返す
end

これで，例えば Eulerスキームを 100step 回したときの結果は下のようになるわけだ．

1

error_and_time(Euler_all, 100) # こうすると，結果は下のように出るぞ．

(0.000998816217400389, 0.0001242) # 誤差と計算時間

これを使って，誤差と計算時間の関係を探ろう． プログラムは以下のような感じになるだろう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

# 全体条件
T = 10.0 # この時間の…
u_exact = 3.1184213196530965 # 厳密解もどき．

ratio  = 1.5 # 計算回数を毎回 1.5倍する <=> Δt を毎回 2/3 に．
N_iter = 15  # この試行を何回やるか．

# 初期化
N = 100
Δt = T/N
u0 = 0.1

# 最初に一回計算しておく．
e_euler, t_euler = error_and_time(Euler_all, N)
e_sym,   t_sym   = error_and_time(Symmetric_all, N)
e_rk,    t_rk    = error_and_time(RK_all, N)

println("iter, N, Δt")  # 画面表示
println("0, $N, $Δt")

# 記録をとりはじめる．
sq_dt = [ Δt ] 
sq_e_euler = [ e_euler ]
sq_e_sym   = [ e_sym ]
sq_e_rk    = [ e_rk ]
sq_t_euler = [ t_euler ]
sq_t_sym   = [ t_sym ]
sq_t_rk    = [ t_rk ]

# Δt を小さくしていったときの近似解をそれぞれ求める．
for i in 1:N_iter
  N = round(Int, ratio * N) # N を ratio 倍して，
  Δt = T / N                # Δt は 1/ratio 倍になるはず．
  push!( sq_dt, Δt )        # 記録をとる
  
  e_euler, t_euler = error_and_time(Euler_all, N)
  e_sym,   t_sym   = error_and_time(Symmetric_all, N)
  e_rk,    t_rk    = error_and_time(RK_all, N)
  
  push!(sq_e_euler, e_euler)  # 記録をとる 
  push!(sq_e_sym,   e_sym)
  push!(sq_e_rk,    e_rk)

  push!(sq_t_euler, t_euler)
  push!(sq_t_sym,   t_sym)
  push!(sq_t_rk,    t_rk)

  println("$i, $N, $Δt")  # 画面表示
end

そして，得られた結果をプロットすると，(準備も含めて)次のような感じにすれば良さそうだ．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

n1(t) = 1.0e-7 / t
n2(t) = 1.0e-10 / (t^2)
n4(t) = 2.0e-22 / (t^4)
n5(t) = 2.0e-26 / (t^5)

t_range_euler = minimum(sq_t_euler):0.0001:maximum(sq_t_euler)
t_range_sym = minimum(sq_t_sym):0.0001:maximum(sq_t_sym)
t_range_rk = minimum(sq_t_rk):0.0001:maximum(sq_t_rk)

default( marker = :auto, xaxis = ("time (sec)", :log10), 
    yaxis = ("error", :log10), legend = :outertopright )

plot(sq_t_euler, abs.(sq_e_euler), label = "Euler")
plot!(sq_t_sym, abs.(sq_e_sym), label = "Symmetric")
plot!(sq_t_rk, abs.(sq_e_rk), label = "Runge-Kutta")

default(marker = :none)

plot!(t_range_euler, n1.(t_range_euler), label = "y = C/x")
plot!(t_range_sym, n2.(t_range_sym), label = "y = C/(x^2)")
plot!(t_range_rk, n4.(t_range_rk), label = "y = C/(x^4)")
plot!(t_range_rk, n5.(t_range_rk), label = "y = C/(x^5)")

  ときどき変なプロット点があるのは，実際の計算時にディスクアクセスや他のアプリの割り込みなどで余計な時間が加わるばらつきがあるためだ．

さて，このプロットを見るとわかるのは以下のような事柄だね．

  • Runge-Kutta法が速くて誤差が小さく，明らかに他の解法より良い．

  • 遅いと言っていた対称スキームだが，計算時間と誤差のバランスを勘案すると Euler スキームよりマシと言って良さそうだ．

  • Euler法と対称スキームは，それぞれ確かに 1次精度，2次精度だ．

  • この場合，Runge-Kutta法は(理論予測より)速く誤差が小さくなっていて，5次精度ぐらいの性能が出ている．

というわけで，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$ として、

\[ \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$ をどう変えるか、という方針というかアルゴリズム4を考えれば良い．

あとはこの基本方針にのっとって $\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次スキームを実装すると次のような感じだ． なお，無駄な計算量がなるべく減るような書き方をしている5．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

function DP(u)
    r1 = r(u)*Δt
    r2 = r(u + (1/5)*r1 )*Δt
    r3 = r(u + (3/40)*r1 + (9/40)*r2 )*Δt
    r4 = r(u + (44/45)*r1 - (56/15)*r2 + (32/9) *r3 )*Δt
    r5 = r(u + (19372/6561) *r1 - (25360/2187) *r2 + (64448/6561) *r3 - (212/729) *r4 )*Δt
    r6 = r(u + (9017/3168) *r1 - (355/33) *r2 + (46732/5247) *r3 + (49/176) *r4 - (5103/18656) *r5 )*Δt

    return u + (35/384)*r1 + (500/1113)*r3 + (125/192)*r4 - (2187/6784)*r5 + (11/84)*r6
end

$\Delta t$ を調整しながら実際の問題の計算に使うなら、Fehlberg 法よりもこちらの Dormand-Prince 法の方が良いだろう．

レポート

下記要領でレポートを出してみよう．

  • e-mail にて，
  • 題名を 2021-numerical-analysis-report-06 として,
  • 教官宛(アドレスは web の "TOP" を見よう)に，
  • 自分の学籍番号と名前を必ず書き込んで，
  • 内容はテキストでも良いし，pdf などの電子ファイルを添付しても良いので，

下記の内容を実行して，結果や解析，感想等をレポートとして提出しよう．

  1. 上のプログラムを少し手直しして，$\Delta t$ を一定の比率で小さくしていったときの近似解を計算して近似解を並べた数列を作り，その数列に Aitken加速6 を(何度か重ねがけして)適用して数列の収束先を予測する方法で厳密解もどきを自分で求めてみよう．
    このときに用いる近似解法としては対称スキームか Runge-Kutta法を用いるとよいだろう．

    ちなみに，この方法でも厳密解もどきとして十分な値が求まるだろう．おそらく，$u_{\rm exact} \cong 3.1184213196531$ ぐらいの値が求まるだろう．

  2. 上で紹介した Runge-Kutta-Fehlberg 法か Dormand-Prince 法の誤差 vs $\Delta t$ のグラフを描いてみて，誤差が $\Delta t$ と何次の関係にあるかを確認しよう．

    なお，この問では $\Delta t$ の自動調整機能は使わず，単なる普通の解法として扱おう．

    ちなみに、パラメータの転記ミス等があると困るので、 英語 wikipedia の List of Runge-Kutta methods の該当する箇所を見てパラメータを確認しつつ、プログラムを進めよう．

  3. 上の問題 2. で組んだプログラムが他の解法と比べてどれくらいの時間がかかるのか，授業で紹介した BenchmarkTools を使って測定しよう．

  4. (チャレンジ問題) 上で紹介した Runge-Kutta-Fehlberg 法か Dormand-Prince 法の本来の使い方である $\Delta t$ の値の自動制御 を試みてみよう．
    なお，この問題はなかなか難しいと思うので，無理に取り組まなくても良い．

  1. 最終結果が微妙に異なるのは誤差のせいだね．どれの誤差が大きいか，はまた別の議論が必要だね． ↩︎

  2. これは古典的 Runge-Kutta 法と後述の 5次 Dormand-Prince 法 で $\Delta t$ をかなり小さくとることで得た近似値の収束の様子からおおよそこのあたりとして求めたものだ． ↩︎

  3. 同じ time t まで計算する状況下では，$\Delta t$ を小さくすることは全体の計算回数をそれに反比例して多くすることを意味する． ↩︎

  4. あまり急に大きく値を修正すると良くないことが起きそうだし，多少は安全係数も掛けたいし，考えどころだね． ↩︎

  5. 含まれる多数の係数が分数なので，これをいったん先に実数に直しておけばもっと速いはず等々思うだろうが，コンパイラが賢くてそのあたりはよしなにやってくれる．むしろ，変数を先に計算して…なんてやるとかえって遅くなることも多い．気になる人は確かめてみると良い．賢いコンパイラよりうまくやるのはなかなか難しいぞ． ↩︎

  6. 魔法の加速法で紹介している ↩︎