Euler法「以外」の解法 for ODE
単純な常微分方程式を 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$ を一定にして、微分を含め、全体を「時間方向に対称に」近似することで近似計算式(スキームと呼ぶ)を作ることを考えよう. 以降、$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\left( \frac{ u_{n+1} + u_n}{2} \right) \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.0
|
|
f (generic function with 1 method)
|
|
df (generic function with 1 method)
|
|
newton (generic function with 1 method)
これで準備が整ったので、実際に数値計算をして時間発展をさせてみよう.
|
|
計算がスムーズに終わっているようなので、早速、グラフにしてみてみよう.
|
|
Plots.GRBackend()
|
|
どうやら無事に計算がうまくいっているようだ.めでたしめでたし.
Euler 法は不安定だったけどこの「対称なスキーム」はどうかな?
前回やってみたように、パラメータによっては Euler 法はふっとんでしまい、不安定だった. この数値スキームはどうだろう? 試してみよう.
|
|
1.2
|
|
|
|
おお、これはきちんと動作するぞ. どうやら、ちょっとばかり手間だけれども、「対称なスキーム」は Euler 法より安心して使えそうだ …が、しかし、次を見てみよう.
Newton 法の代わりにライブラリに頼ってみよう!
復習を兼ねて、以前学習した、NLsolve パッケージを使ってみよう. まずは、使うとき毎回やるルーチン的な手続きは以下の通りだった.
|
|
あとは前やったように、以下のようにすれば良いはずだ.
|
|
ふむ、無事に動いたようだ.グラフにしてみてみよう.
|
|
グラフを見ると、上の例とまったく同じようだ.内容的にも問題ないと見て良いだろう.
念のために、ライブラリを使った場合でもパラメータを変えてみよう
先と同じように、$\Delta t = 1.2$ として挙動を見てみよう.
|
|
1: 0.2445828206408548: 0.17144028204630682 2: 0.7241519610699485: 0.4656486127362301 3: 1.8709232707904873: 0.962896580041833 4: 2.7192167342635747: 0.7489808758544224 5: 2.985520115427732: 0.2852087610786438 6: 3.0725584588957093: 0.1123158743514151 7: 3.1069544379675222: 0.0518129943639024 8: 3.1225677124829754: 0.026828358987976905 9: 3.1304047960356782: 0.015105824781990605 10: 3.1346510039397333: 0.009064629461162748 11: 3.1370946477756876: 0.005719796543410794 12: 3.1385715278210715: 0.003759556934926038 13: 3.139501169967586: 0.0025563019113543084 14: 3.140106810914135: 0.0017886621951831016 15: 3.140513174918811: 0.00128266032161094 16: 3.1407928947061197: 0.000939618639065662 17: 3.140989805063468: 0.0007013036475126663 18: 3.141131196989348: 0.0005321525382685236 19: 3.1412345332852865: 0.0004097884410071393 20: 3.1413112628758575: 0.0003197555037724399 21: 3.1413690567791157: 0.00025249375962353975 22: 3.1414131562201306: 0.0002015470888057806 23: 3.141447206188787: 0.00016247238461971784 24: 3.1414737826971733: 0.00013215914642826653 25: 3.1414947332426224: 0.00010839561968291782 26: 3.1415114008269587: 8.958655488308557e-5 27: 3.1415247739178596: 7.456621731479871e-5 28: 3.1415355805694105: 6.247634611733667e-5 29: 3.1415443858855614: 5.267036228293089e-5 30: 3.1415516094845417: 4.46559047268447e-5 31: 3.1415575733329284: 3.8062181048994965e-5 32: 3.141562526596144: 3.260362525126794e-5 33: 3.141566663653864: 2.8058464785367872e-5 34: 3.141570137293549: 2.4253116084303685e-5 35: 3.1415730684652154: 2.105071040937595e-5 36: 3.1415755535587575: 1.8342577805956342e-5 37: 3.141577669879381: 1.6041870723227465e-5 38: 3.1415794798010395: 1.4078749582692603e-5 39: 3.14158103394284: 1.2396717852940844e-5 40: 3.1415823736216253: 1.0949807560126888e-5 41: 3.1415835327860404: 9.700385960235457e-6
DomainError: Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace: [1] nan_dom_err at ./math.jl:300 [inlined] [2] ^ at ./math.jl:699 [inlined] [3] f(::Float64, ::Float64) at ./In[2]:2 [4] #10 at ./In[13]:5 [inlined] [5] finite_difference_jacobian!(::Array{Float64,2}, ::##10#12{#f,Tuple{Float64}}, ::Array{Float64,1}, ::DiffEqDiffTools.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}}) at /home/jrun/.julia/v0.6/DiffEqDiffTools/src/jacobians.jl:138 [6] (::NLsolve.#fj!#1{##10#12{#f,Tuple{Float64}}})(::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}) at /home/jrun/.julia/v0.6/NLsolve/src/objectives/autodiff.jl:6 [7] value_jacobian!!(::NLSolversBase.OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1},Val{false}}, ::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}) at /home/jrun/.julia/v0.6/NLSolversBase/src/interface.jl:124 [8] trustregion(::NLSolversBase.OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1},Val{false}}, ::Array{Float64,1}, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::Bool, ::Float64, ::Bool, ::NLsolve.NewtonTrustRegionCache{Array{Float64,1}}) at /home/jrun/.julia/v0.6/NLsolve/src/solvers/trust_region.jl:118 [9] #nlsolve#38(::Symbol, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::Bool, ::LineSearches.Static{Float64}, ::Float64, ::Bool, ::Int64, ::Float64, ::NLsolve.#nlsolve, ::NLSolversBase.OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1},Val{false}}, ::Array{Float64,1}) at /home/jrun/.julia/v0.6/NLsolve/src/nlsolve/nlsolve.jl:26 [10] (::NLsolve.#kw##nlsolve)(::Array{Any,1}, ::NLsolve.#nlsolve, ::NLSolversBase.OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1},Val{false}}, ::Array{Float64,1}) at ./
:0 [11] #nlsolve#39(::Symbol, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::Bool, ::LineSearches.Static{Float64}, ::Float64, ::Bool, ::Int64, ::Float64, ::Symbol, ::Bool, ::NLsolve.#nlsolve, ::##10#12{#f,Tuple{Float64}}, ::Array{Float64,1}) at /home/jrun/.julia/v0.6/NLsolve/src/nlsolve/nlsolve.jl:59 [12] #nls#9(::Float64, ::Function, ::Function, ::Float64, ::Vararg{Float64,N} where N) at ./In[13]:5 [13] (::#kw##nls)(::Array{Any,1}, ::#nls, ::Function, ::Float64, ::Vararg{Float64,N} where N) at ./ :0 [14] macro expansion at ./In[16]:7 [inlined] [15] anonymous at ./ :?
Euler 法よりは随分マシだけれども、この場合も途中で計算が停止してしまう. 表示されているのはエラーの直前の数字なのでわかりにくいが、この場合も $\sin$ 関数の中身が次の計算ステップで $\pi$ を超えるために $\sin$ 関数の値がマイナスになり、1.2乗に失敗するという、同じ形の失敗のようだ.
結局この近似式では $\Delta t$ を大きくした場合、Newton 法だとうまくいったけれども、NLsolve パッケージだとうまくいかないという結果になった. このように、簡単な問題でも、手法次第で結果が得られるかどうかが容易に変わってしまう ので、実際の計算はある程度「予備計算で試行錯誤してみてから」柔軟に取り組む必要がありそうだ.
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 & = & \displaystyle r(u, t), \cr r_2 & = & \displaystyle r\left(u + \frac{1}{2} \,\Delta t\, r_1 ,\, t+ \frac{1}{2} \,\Delta t \right), \cr r_3 & = & \displaystyle r\left(u + \frac{1}{2} \,\Delta t\, r_2 ,\, t+ \frac{1}{2} \,\Delta t \right), \cr r_4 & = & \displaystyle r(u + \Delta t \, r_3 , \, t+\Delta t) \end{array} \right. $$
と計算して、$t$ から $t + Δt$ の間の $u$ の 時間変化 $r$ の近似値を 4つ計算する. そしてこれらを用いて、
$u(t + \Delta t)$ の近似値 $\displaystyle = u + \frac{1}{6} \,\Delta t \, (r_1 + 2r_2 + 2r_3 + r_4)$
として、時間ステップを 1つ進める近似解法がこの方法だ.
早速、プログラムしてみよう. まず、この問題の場合は、常微分方程式の右辺 $r(u,t)$ は次のようになる.
注: この問題の右辺には $t$ は登場しないので省略しているよ!
|
|
r (generic function with 1 method)
すると、Runge-Kutta 法で 1ステップすすめる関数をプログラムで書くと次のような感じだな.
|
|
RK (generic function with 1 method)
早速使ってみよう! まずは安全そうなパラメータで.
|
|
|
|
ふむふむ、うまくいっているな. では、$\Delta t$ を大きくしてみて、うまくいくかチェックしよう.
|
|
|
|
どうやら、$\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 の該当する箇所を見てパラメータを確認しつつ、プログラムを進めよう.