時間方向に対称なスキーム, Runge-Kutta スキーム

単純な常微分方程式を 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つ進める、ということが明確になる．

あとは 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

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

# 関数 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_sq = [u0] # 初期値を配列に入れておいて…
u = u0       # 「今」の値は、もちろん最初は初期値そのもの．

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

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

1
2
3
4
5

using Plots
gr()

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

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

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

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

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

1

Δt = 1.2

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

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

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)

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

非線形方程式を解く部分を、ライブラリに頼ってみよう！

復習を兼ねて、以前学習した、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_sq = [u0] # 初期値を配列に入れておいて…
u_old = u0  # 「今」の値は、もちろん最初は初期値そのもの．

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)

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

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

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

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

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

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

1: DomainError: Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar. …略… 39: 3.14158103394284: 1.2396717852940844e-5 40: 3.1415823736216253: 1.0949807560126888e-5 41: 3.1415835327860404: 9.700385960235457e-6

Euler 法よりは随分マシだけれども、この場合も途中で計算が停止してしまう． 表示されているのはエラーの直前の数字なのでわかりにくいが、この場合も $\sin$ 関数の中身が次の計算ステップで $\pi$ を超えるために $\sin$ 関数の値がマイナスになり、1.2乗に失敗するという、同じ形の失敗のようだ．

このように、簡単な問題でも、手法次第で結果が得られるかどうかが容易に変わってしまう ので、計算には柔軟に取り組む必要がありそうだ．

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 法

時間 $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つ ($r_1, \cdots, 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_sq = [u0] # 初期値を配列に入れておいて…
u = u0  # 最初の値はもちろん初期値

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)

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

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

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

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

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)

どうやら、$\Delta t$ を大きくしてもこの問題では 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} \| \]

として、具体的に数値で見積もれることになる.

tips: ベクトルの評価にノルムを使っていることに慣れよう．$r$ がスカラー値とは限らないからね． $r$ がスカラー値の場合は絶対値でも良い．けど、ベクトル値や関数だったら、上の式のように例えばノルムにしないとね．そのあたりの配慮があれば、ノルムでなくても良いよ．

そして、そもそも $p$ 次近似値の局所誤差(1ステップだけで発生する誤差) $\epsilon_p$ は

\[ \epsilon_p \cong C\, (\Delta t)^{p+1} \]

となるように上の Fehlberg 公式は設計されている($C$ は時間と $u$ 等に依存する適当な係数)．

tips: 時間発展問題の近似式で「$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$ をどう変えるか、という(工夫とでもいうべき)アルゴリズムを考えれば良い(あまり急に大きく値を修正すると良くない、とされているし、多少は安全係数も掛けたい)．

あとはこの基本方針にのっとって $\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 法の方が良いだろう．

Report No.04: 上で紹介した Runge-Kutta-Fehlberg 法か Dormand-Prince 法を用いてみよう．

ちょっと面倒だけれども、上に紹介した Runge-Kutta-Fehlberg 法か Dormand-Prince 法を実際にプログラムして、上の問題に適用してみよう． ステップ毎に $\Delta t$ を調整するようにプログラムを組めれば理想的だが、難しいようならば $\Delta t$ は固定値のままでもまあ良い．

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