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( \frac{ u_{n+1} + u_n}{2} )\right)^{1.2} $$ となる.
これはつまり、数値計算のプロセス上では $u_n$ は既知で $u_{n+1}$ は未知量であるので、なにをしないといけないかというと、 関数 $f(u, U)$ ($u_{n+1} = u$, $u_n = U$ と表記)を $$ f(u, U) := u - U - \Delta t \gamma\, \left( \sin( \frac{ u + U }{2} )\right)^{1.2} $$ として与えられた $U$ に対して $f(u, U)$ がゼロになるような $u$ を求めるという「求解」を行なわないといけない、ということになる．
あとは Euler 法同様に、初期値は既知で $u_0 = u(0)$ とすることができるので、さらなる近似値を上のスキームを用いて $u_1, u_2, \cdots$ と順番に求めていけば好きなだけ先の近似値を計算できることがわかる．
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 法部分も含めて以下のようになるだろう．
# 定数の設定
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 を求める
function newton(old_u)
u = old_u # とりあえず、初期値を古い値の U と一致させておいて…
for i in 1:100
if abs( f(u,old_u) ) < 1.0e-6 # 6桁あってればいいね、という場合はこんな感じで判定する．
break # break というのは一番内側のループ(この場合は for )を脱出する… はず．
end
u = u - f(u,old_u)/df(u,old_u) # Newton 法で u を改善する．
end
return u # f が十分にゼロに近くなったら、近似値 u を返す
end
u_sq = [u0] # 初期値を配列に入れておいて…
u = u0 # 「今」の値は、もちろん最初は初期値そのもの．
for n in 1:100
u = newton(u) # Newton 法で新しい時間ステップの値が求まって出てくるはず．
push!(u_sq, u) # その値を配列に追加していく．
end
計算がスムーズに終わっているようなので、早速、グラフにしてみてみよう．
using Plots
gr()
plot(u_sq)
どうやら無事に計算がうまくいっているようだ．めでたしめでたし．
前回やってみたように、パラメータによっては Euler 法はふっとんでしまい、不安定だった． この数値スキームはどうだろう？ 試してみよう．
Δt = 1.2
u_sq = [u0] # 初期値を配列に入れておいて…
u = u0 # 「今」の値は、もちろん最初は初期値そのもの．
for n in 1:100
u = newton(u) # Newton 法で新しい時間ステップの値が求まって出てくるはず．
push!(u_sq, u) # その値を配列に追加していく．
end
plot(u_sq)
おお、これはきちんと動作するぞ． どうやら、ちょっとばかり手間だけれども、「対称なスキーム」は Euler 法より安心して使えそうだ．
復習を兼ねて、以前学習した、NLsolve パッケージを使ってみよう． まずは、使うとき毎回やるルーチン的な手続きは以下の通りだった．
using NLsolve
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
r = nlsolve((vin,vout)->vout[1]=func(vin[1],params...), [ini])
v = r.zero[1]
else
r = nlsolve((vin,vout)->vout .= func(vin,params...), ini)
v = r.zero
end
return v, r.f_converged
end
あとは前やったように、以下のようにすれば良いはずだ．
Δ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
ふむ、無事に動いたようだ．グラフにしてみてみよう．
plot(u_sq)
グラフを見ると、上の例とまったく同じようだ．内容的にも問題ないと見て良いだろう．
先と同じように、$\Delta t = 1.2$ として挙動を見てみよう．
Δ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
Euler 法よりは随分マシだけれども、この場合も途中で計算が停止してしまう． 今回は出力の $u$ は $\pi$ を超えていないのだけれども、$\sin$ が大変ゼロに近く、かなり微妙な状態になっているようだ．