2020.11 大学から学生諸氏への通知

Volterra モデルの保存量を調べてみよう

Volterra モデルが保存系らしいことはこれまでの線形安定性解析や数値計算から推測されるところだ．なにせ相図でぐるっと回ってもとと同じ状態に戻ってくるのだ．いろいろ保存していると思って間違いない．

そこで少し真面目に考えてみよう．こうした考察は問題の系の「大域的性質」を調べることにあたり、問題の理解にも数値解析にもたいへん役に立つ．

保存量の探索: 方針

まず、保存系ということは、系になにか保存「量」(= 時間変化しない量)がある可能性が高い．そこで、こうした量が存在するかどうか、手探りで探してみよう．

そうした量、つまり保存量を仮に $G(u(t), v(t))$ と書くとして、一番簡単な考え方は次のようなものだろう．まず、$G$ が保存量であるということは時間で微分するとゼロということで、さらにその微分を連鎖律で分解して、

\[ \displaystyle 0 = \frac{d}{dt} G = \frac{\partial G}{\partial u} \frac{du}{dt} + \frac{\partial G}{\partial v} \frac{dv}{dt} = \displaystyle \left( \begin{array}{c} \displaystyle \frac{\partial G}{\partial u} \cr\cr \displaystyle \frac{\partial G}{\partial v} \end{array} \right) \cdot \left( \begin{array}{c} \displaystyle \frac{du}{dt} \cr\cr \displaystyle \frac{dv}{dt} \end{array} \right) = \nabla G \cdot \frac{d}{dt} \boldsymbol{u} \]

となることは間違いないので、

$\nabla G \perp \displaystyle \frac{d}{dt} \boldsymbol{u}$ となるような $G$ を探せばそれは保存量だ

ということになる．ではこの方針で、$G$ を探してみよう．

方程式を整理しておく

今回の Volterra モデルを改めて書くと(漁業効果ありで)

\[ \left\{\begin{array}{rcl} \displaystyle \frac{du}{dt} & = & u \left( -C_1 v + (D_1 - E) \right), \cr \displaystyle \frac{dv}{dt} & = & v \left( C_2 u - (D_2 + E) \right) \end{array}\right. \]

なわけだが、この右辺を見ると、それぞれ $u,v$ で割るときれいになりそうだ．そこでそうしてみると、

\[ \left\{\begin{array}{rcl} \displaystyle \frac{1}{u}\frac{du}{dt} & = & -C_1 v + (D_1 - E), \cr \displaystyle \frac{1}{v}\frac{dv}{dt} & = & C_2 u - (D_2 + E) \end{array}\right. \]

となり、右辺は扱いやすいかっこうになった．しかも、よく考えると、 \[ \displaystyle \frac{1}{u}\frac{du}{dt} = \frac{d}{dt} \log(u) \]

なので、これを用いて左辺を書き直すと

\[ \left\{\begin{array}{rcl} \displaystyle \frac{d}{dt} \log(u) & = & -C_1 v + (D_1 - E), \cr \displaystyle \frac{d}{dt} \log(v) & = & C_2 u - (D_2 + E) \end{array}\right. \]

となる．そこでこれをじっと眺めてみて，そうだなあ，

\[ \left\{\begin{array}{rcl} \displaystyle U & := \log(u), \cr \displaystyle V & := \log(v) \end{array}\right. \]

と変数変換してみると筋が良さそうだ．そこで、この変数変換を採用して、Volterra モデルを書き直すと、

\[ \left\{\begin{array}{rcl} \displaystyle \frac{dU}{dt} & = -C_1 e^V + (D_1 - E), \cr \displaystyle \frac{dV}{dt} & = C_2 e^U - (D_2 + E) \end{array}\right. \]

となる．だいぶ扱いやすくなった気がするので、この式をベースに、いろいろ考えよう．

$\nabla G \perp \displaystyle \frac{d}{dt} \boldsymbol{U}$ となるような $G$ を探す: 今回はラッキー

変数変換して扱いやすくなった式をベースに考えることにしたので、上の方針は

\[ \left( \begin{array}{c} \displaystyle \frac{\partial G}{\partial U} \cr \displaystyle \frac{\partial G}{\partial V} \end{array} \right) \cdot \left( \begin{array}{c} \displaystyle \frac{dU}{dt} \cr \displaystyle \frac{dV}{dt} \end{array} \right) = 0 \label{eq:perp} \]

であるような $G(U,V)$ を探そうという方針に翻訳できる．

本来この探索はここから難しいわけで、いろんなやり方を試してみないといけないし，それらも簡単ではないが，まあ今回はなんとかなるのだ．

さて，こういうときに最初に試すべき作業は、以下のようなものだろう．

まず、「仮に」 \[ \left( \begin{array}{c} \displaystyle \frac{dU}{dt} \cr \displaystyle \frac{dV}{dt} \end{array} \right) = \left( \begin{array}{c} \displaystyle - \frac{\partial G}{\partial V} \cr \displaystyle \frac{\partial G}{\partial U} \end{array} \right) \label{eq:hypo} \] という等式が成り立つと仮定してみよう．この右辺はあきらかに $ \displaystyle \left( \begin{array}{c} \displaystyle \frac{\partial G}{\partial U} \cr \displaystyle \frac{\partial G}{\partial V} \end{array} \right) $ と直交する(内積を計算してみるとすぐわかる)ので，仮定 $(\ref{eq:hypo})$ が成り立つならば狙っていた直交性 $(\ref{eq:perp})$ は満たされるということになる．

だから，この仮定の等式 $(\ref{eq:hypo})$ が成り立つような $G$ を見つけることができれば それは保存量だ、ということになる．そこで「式 $(\ref{eq:hypo})$ を満たす $G$ を見つけようと努力すること」が最初に試すべきことなのだ．

では、早速だが，まずはこの仮定の等式を成り立たせる $G$ が見つけられるか考えよう．この式 $(\ref{eq:hypo})$ を使ってもとの Volterra 方程式を書き直してみると、

\[ \left\{\begin{array}{rcl} \displaystyle - \frac{\partial G}{\partial V} & = -C_1 e^V + (D_1 - E), \cr \displaystyle \frac{\partial G}{\partial U} & = C_2 e^U - (D_2 + E) \end{array}\right. \]

となる．そしてわれわれは，これを満たす $G(U,V)$ が見つけられるか，を考えるわけだ. そして明らかに、大変ラッキーにもこれは可能で、

\[ G(U,V) = C_2 e^U - (D_2 + E) U + C_1 e^V - (D_1 - E) V \label{eq:invariant-by-log} \]

とすればよいことがみてすぐわかる．つまりこの $G$ は仮定の等式 $(\ref{eq:hypo})$ を満たす．

ということは，保存量 $G$ が見つかった¹のである！！いや～，ラッキーだな．

さて，$U, V$ の表記のままだとわかりにくいので，この保存量 $G$ をあらためて元の変数 $u,v$ を使ったものに書き直そう．すると

\[ G(u,v) = C_2 u - (D_2 + E) \log(u) + C_1 v - (D_1 - E) \log(v) \label{eq:invariant} \]

となる.

さて、以降はこの量が数値計算でどういう挙動をみせるかを調べていこう．

まずは Euler 法ではどうなのか、再チェックだ

というわけで、Euler スキームを準備して、計算してみよう．

1
2
3
4
5


# 各種定数を前回同様に設定して…
C1 = C2 = D1 = D2 = 1.0
u0 = v0 = 0.7
Δt = 0.05
E = 0.5

1
2
3
4
5
6


# Euler スキームを定義する．
function volterra_euler(u,v,E)
    u_new = (1 + Δt * (  D1 -E - C1*v ) ) * u
    v_new = (1 + Δt * ( -D2 -E + C2*u ) ) * v
    return u_new, v_new 
end

1
2


# 念のために、1ステップだけでも動かしてみる．
volterra_euler(u0,v0,E)

(0.693, 0.6719999999999999)

また、関数 $G$ を定義して、時間発展計算中にこの値を計算できるようにしよう．念のため、漁業係数も変えられるようにしておこう．

1

G(u,v,E) = C2*u - (D2+E)*log(u) + C1*v - (D1-E)*log(v)

では、いつものように動かしてみよう．あとでグラフにしてみたいので、$G$ の値も取っておくとしよう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10


# 大丈夫そうなので、500ステップほど動かしてみよう
u,v = u0,v0
uv_sq = [u0 v0]
G_sq = [ G(u,v,E) ]

for i in 1:500
    u,v = volterra_euler(u,v,E)
    uv_sq = vcat(uv_sq, [ u v ])
    push!( G_sq, G(u,v,E) )
end

無事に動いているか、簡単に確認してみる．

1

uv_sq

501×2 Array{Float64,2}:
0.7       0.7
0.693     0.672
0.68704   0.644885
0.682063  0.618672
0.678016  0.59337
⋮
0.39047   0.487357
0.390716  0.46032
0.391492  0.434789
0.392768  0.41069
0.394522  0.387954

とりあえず結果に全て数字が入っているので、計算はできているようだ．

では早速、グラフで状況を見ていこう．まずいつものお約束．

1

using Plots

そしてプロットだ．

1
2
3


t_sq = Δt * [0:500]
plot(t_sq, uv_sq, 
  xaxis = "time t", yaxis = "u and v", leg = false)

svg

先にやったとおり、Euler 法だとだんだん値の振れが大きくなっているなあ．一応、相図もみておこう．

1
2
3
4
5
6


u_sq = uv_sq[:,1] # データの 1行目 = u の近似値列
v_sq = uv_sq[:,2] # データの 2行目 = v の近似値列

plot(u_sq, v_sq, aspect_ratio = 1, xaxis = "u", yaxis = "v")

plot!((u0,v0), marker = :circle, lab = "initial value")

svg

解の軌道がだんだん大きくなっていることが相図からもわかるね．いや～、これはどうなんだろう．

ちなみに、この初期値からの解の軌道を、保存量 $G$ が一定の軌道として描いてみよう．まず、一定であるべき $G$ の値は初期値を使って

1

iniG = G(u0, v0, E)

2.113349887877465

と求まるので、あとはグラフを適当に描いてみる．今回は次のようにしてみよう．

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


U = 0.5:0.01:3.0   # 解の値をみて適当に
V = 0.01:0.01:1.5

Z = [ G(u,v,E) for u in U, v in V ]
# U, V の範囲全部の G を計算してしまえ．

contour(U, V, Z', aspect_ratio = 1, 
        clims = (iniG-0.001, iniG+0.001)) 
# その代わり，G が iniG に近い値の場合だけグラフ出力する．
# clims で、出力値の描画を制限できるのだ．

plot!((u0,v0), marker = :circle, lab = "initial value")

svg

Euler 法の解の軌道とずいぶん違うね．やっぱり Euler法はあんまり良くないね…

この問題ではこの $G = $constant. という式が，相図上での解軌道の式そのものである ことがここで(ほぼ)確認されたことになる！！
先に述べたように，相図上の解軌道の厳密な式が求まることには大きな価値があるが実際に求めることは通常は困難なので，今回のこれは大変ラッキーだ．

さて、この Euler法の結果で、本来は保存されるはずの量 $G(u,v)$ がどう変化しているかをみてみよう．

1

plot(t_sq, G_sq, xaxis = "time t", yaxis = "G", leg = false)

svg

…やっぱり、保存量のはずなのに $G$ が時間がたつにつれ大々的に増加している．これはまずい． Euler 法があまり信用ならないなあ、ということがこのケースではよく分かる．

時間対称なスキームならどうなの？

時間対称なスキームは無数のバリエーションがあるので、どう作るかは問題なんだけれども、まあ、とりあえずすごく単純に作って試してみよう．

\[ \left\{\begin{array}{rcl} \displaystyle \frac{ u_{n+1} - u_n}{\Delta t} & = & \displaystyle - C_1 \left( \frac{ u_{n+1} v_{n+1} + u_n v_n }{2} \right) + (D_1 - E) \left( \frac{ u_{n+1} + u_n }{2} \right), \cr \displaystyle \frac{ v_{n+1} - v_n}{\Delta t} & = & \displaystyle C_2 \left( \frac{ u_{n+1} v_{n+1} + u_n v_n }{2} \right) - (D_2 + E)\left( \frac{v_{n+1} + v_n}{2} \right) \end{array}\right. \]

これは次のように書き直せるので、

\[ \left\{\begin{array}{rcl} 0 & = & \displaystyle \left( 1 + \frac{\Delta t}{2} \left( - D_1 + E + C_1 v_{n+1} \right) \right) u_{n+1} + \left( -1 + \frac{\Delta t}{2} \left( - D_1 + E + C_1 v_{n} \right) \right) u_{n} , \cr 0 & = & \displaystyle \left( 1 + \frac{\Delta t}{2} \left( D_2 + E - C_2 u_{n+1} \right) \right) v_{n+1} + \left( -1 + \frac{\Delta t}{2} \left( D_2 + E - C_2 u_{n} \right) \right) v_{n} \end{array}\right. \]

$u_n, v_n$ を既知として、この連立非線形方程式を解いて $u_{n+1}, v_{n+1}$ を求めれば時間ステップを一つ進めることが出来る．

というわけで、まずはゼロになるべきこの右辺を関数として定義しよう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17


function symm_zero(uv_new, uv_curr, E) 
# uv_new が新しい u,v で、uv_curr は古い u,v
    
    up, vp = uv_new
    u, v   = uv_curr
    half = Δt/2

    eq1 = ( 1 + half * ( - D1 + E + C1*vp ) ) * up +
          ( -1 + half * ( - D1 + E + C1*v ) ) * u
    # eq1 = ... の右端の "+" に注意．
    # この式がそこの改行で終わっていないことを julia に伝えている.

    eq2 = ( 1 + half * (   D2 + E - C2*up ) ) * vp +
          ( -1 + half * (   D2 + E - C2*u ) ) * v

    return [ eq1, eq2 ]
end

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

そして、nlsolve を使ってこのスキームを作る．以下のように，u, v の今の値を渡すと新しい値が返ってくるように書くだけだ．

1
2
3
4
5
6


function volterra_symm(uv_curr, E)
    
  uv_new = nls(symm_zero, uv_curr, E, ini = uv_curr)[1]
  return uv_new
    
end

まずは、動くかどうか簡単にチェック．nlsolve をセッション中で初めて動かすときはいつものように少し待たされる．

1
2


uv_0 = [u0, v0]
volterra_symm(uv_0, E)

2-element Array{Float64,1}:
 0.6935102416466368
 0.6724420591556107

うん，無事に動くようですねえ．では、もう少し動かしてみる…前に、引数が楽に書ける G の定義をもう一つ作っておくと楽そうだ．

1

G(uv,E) = G(uv[1], uv[2], E)

この定義は、Julia の multi-dispatch 機能、つまり、「同じ関数名でも引数の型が違うときは違う定義を与えることができる」ことをうまく使っている

では、500ステップほど動かしてみよう

1
2
3
4
5
6
7
8
9


uv = copy(uv_0)
uv_sq = copy(uv_0') # 横向きにしておく
G_sq = [ G(uv, E)]

for i in 1:500
    uv = volterra_symm(uv, E)
    uv_sq = vcat( uv_sq, uv' )
    push!(G_sq, G(uv, E))
end

結果をプロットしてみる

1
2


plot(t_sq, uv_sq, 
     xaxis = "time t", yaxis = "u and v", leg = false) 

svg

…おお、振幅が一定だ．Euler 法よりは信用できるっぽい？

相図も描いておこう．

1
2
3
4
5
6


u_sq = uv_sq[:,1] # データの 1行目 = u の近似値列
v_sq = uv_sq[:,2] # データの 2行目 = v の近似値列

plot(u_sq, v_sq, aspect_ratio = 1, xaxis = "u", yaxis = "v")

plot!((u0,v0), marker = :circle, lab = "initial value")

svg

おお！ Euler 法と異なり、「元へ戻る」ことが明確にわかるぞ．

ちなみに、理論的に求めた $G$ = 一定の線と重ねて描画してみよう．

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


U = 0.5:0.01:3.0   # 解の値をみて適当に
V = 0.01:0.01:1.5

Z = [ G(u,v,E) for u in U, v in V ]

contour(U, V, Z', aspect_ratio = 1, 
  clims = (iniG-0.001, iniG+0.001)) 
# clims で、出力値の描画を制限できる

plot!((u0,v0), marker = :circle, lab = "initial value")

plot!(u_sq, v_sq, aspect_ratio = 1, xaxis = "u", yaxis = "v")

svg

黄色い線と青い線が重なって緑色に見えている．つまり，目で見る限り、もう区別がつかないということだな．この時間対称スキームによる近似解はだいぶ良いと言ってよさそうだな．

さて、この近似解が保存量 $G$ をどれくらい時間変化させてしまうかもチェックしよう．

1

plot(t_sq, G_sq)

svg

おお、Euler法よりは格段に良いな！おおまかに 3桁ぐらいは保たれている．長期的に見ればもとへ戻ってきている点も素晴らしい．

さて、格段に良いことは良いのだけれども、本来は $G$ は不変量であるのだから、「ブレ」をもうちょっと小さくできないか、もうちょっと欲張ってもよいだろう．なんとかならないものだろうか？

Runge-Kutta 法だとどうだろう？

古典的 Runge-Kutta 法を試してみよう．4次精度なので、だいぶ良いだろうことが期待できる．

まず、常微分方程式そのものの「右辺」を次のように定義して、

 1
 2
 3
 4
 5
 6
 7
 8
 9
10


function r(uv) 

    u,v = uv

    eq1 = u * ( -C1*v + (D1-E)) 
    eq2 = v * (  C2*u - (D2+E)) 

    return [ eq1, eq2 ]

end

前もそうしたが、Runge-Kutta 法を直接定義してしまおう．

1
2
3
4
5
6
7


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

早速動かしてみよう！

1
2
3
4
5
6
7
8
9


uv = copy(uv_0)
uv_sq = copy(uv_0')
G_sq = [ G(uv, E)]

for i in 1:500
    uv = RK(uv)
    uv_sq = vcat(uv_sq, uv')
    push!(G_sq, G(uv, E))
end

もう面倒なので、いきなりグラフを出してみる．

1
2


plot(t_sq, uv_sq, 
     xaxis = "time t", yaxis = "u and v", leg = false) 

svg

ふむ、結果は悪くなさそうだ．相図はもう面倒なので描かない(うまくいっていることは想像に難くない).

では、保存量である $G$ はどうかな．

1

plot(t_sq, G_sq, leg = false)

svg

少しずつ増えているが、ブレは小さそうだ．ただ、このグラフだと目盛りが変わらなくてわからないので、全体から初期値を引いて、そのグラフを描いてみよう．

1

plot(t_sq, G_sq .- G_sq[1], leg = false)

svg

こうしてみると、ブレはかなり小さくて、7～8桁程度は $G$ が保たれている．すごいね！

贅沢を言うならば，ズレは蓄積していくようなので長期計算には向かないかなあ，というところか．

構造保存解法の出番か

その理論については次回以降に解説するが、こうした問題に対して、保存量は保存するように、散逸量は散逸するように、という数値スキームを構成する方法がある．こうしたスキームを「構造保存数値解法」と呼んだりする²．

さて、どうやって導出するかは後述するが、今回の問題に対して、最もシンプルな構造保存スキームは以下のような物になる感じだ．

まず、以下のように変数変換しておく．

\[ \left\{\begin{array}{rcl} \displaystyle U & := & \log(u), \cr \displaystyle V & := & \log(v) \end{array}\right. \]

そして、これに対して、 $U_n \cong U(n\Delta t)$, $V_n \cong V(n\Delta t)$, として、

\[ \left\{\begin{array}{rcl} \displaystyle \frac{ U_{n+1} - U_n}{\Delta t} & = & \displaystyle - C_1 \left( \frac{ d (\exp) }{ d( V_{n+1}, V_n ) } \right) + D_1 - E, \cr \displaystyle \frac{ V_{n+1} - V_n}{\Delta t} & = & \displaystyle C_2 \left( \frac{ d (\exp) }{ d(U_{n+1}, U_n) } \right) - D_2 - E \end{array}\right. \label{eq:sp-scheme} \]

が、それである．ただし、$d/d(\cdot, \cdot)$ は差分商とよばれるもので、

\[ \displaystyle \frac{ d(f)}{d(a,b)} := \left\{\begin{array}{lcl} \displaystyle \frac{ f(a) - f(b) }{ a-b } & & {\rm when }\,\, a \neq b , \cr f'(a) & & {\rm when }\,\, a = b \end{array}\right. \]

で定義される．

そして、これは次のように書き直せるので、

\[ \left\{\begin{array}{rcl} 0 & = & \displaystyle U_{n+1} - U_n + \Delta t C_1 \left( \frac{ d(\exp) }{ d(V_{n+1}, V_n) } \right) + \Delta t ( - D_1 + E), \cr 0 & = & \displaystyle V_{n+1} - V_n - \Delta t C_2 \left( \frac{ d(\exp) }{ d(U_{n+1}, U_n) } \right) + \Delta t (D_2 + E) \end{array}\right. \]

$U_n, V_n$ を既知として、この連立非線形方程式を解いて $U_{n+1}, V_{n+1}$ を求めれば時間ステップを一つ進めることが出来るというスキームの出来上がりである³．

さて、まずは今回使う差分商を与えておこう．

1
2
3
4
5


function dedx(a,b)
  EPS = 1.0e-16
  return ifelse( abs(a-b) > EPS, (exp(a)-exp(b))/(a-b), exp(a) )
    # ifelse 関数を調べておこう
end

あとは、上の時間対称なスキームと同様の作業をしていけば良い．とちゅうで $\log$ が出てくることを忘れないようにしよう．

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


function sp_zero(uv_new, uv_curr, E) 
# uv_new は新しい u,v で、uv_curr は古い u,v

    up, vp = log.(uv_new)   
    # ドットをつけた計算は、各要素に作用するよ．
    u, v   = log.(uv_curr)
    
    eq1 = up - u + Δt*C1* dedx(vp,v) + Δt*(-D1 + E)
    eq2 = vp - v - Δt*C2* dedx(up,u) + Δt*( D2 + E)

    return [ eq1, eq2 ]
end

これを用いて1ステップ計算するための関数を定義する．

1
2
3
4
5
6


function volterra_sp(uv_curr, E)

  uv_new = nls(sp_zero, uv_curr, E, ini = uv_curr)[1]
  return uv_new

end

試しておこう．

1

volterra_sp(uv_0, E)

 2-element Array{Float64,1}:
  0.6935156878695156
  0.6724434211370223

動くようなので、これも 500ステップほど動かしてみよう．

1
2
3
4
5
6
7
8
9


uv = copy(uv_0)
uv_sq = copy(uv_0')
G_sq = [ G(uv, E)]

for i in 1:500
    uv = volterra_sp(uv,E)
    uv_sq = vcat(uv_sq, uv')
    push!( G_sq, G(uv, E) )
end

早速プロットしてみよう．

1
2


plot(t_sq, uv_sq, 
     xaxis = "time t", yaxis = "u and v", leg = false)

svg

ふむ、なかなか美しい．では、保存量 $G$ がどうなっているか見てみよう．

1

plot(t_sq, G_sq)

svg

GKS: Possible loss of precision in routine SET_WINDOW

おや、プロットと一緒になにか警告が出ている．これは、plot 命令にとって「値があまりにも変わらないので、図の上下幅の設定が怪しいかも」というような話のようだ．実際，この図も縦軸の目盛りの数字が全く変化していないしね．

$G$ の値を数字で直接見てみよう．

1

G_sq

501-element Array{Float64,1}:
2.113349887877465 
2.113349887877454 
2.1133498878774444

2.113349887876992 
2.113349887876999 
2.1133498878770056

う～む，上12桁ぐらいは変化していないな．これをグラフに描くのは確かに大変だ．

では、はっきりと見るために、$G$ の「初期値との差」をグラフで見てみよう．

1

plot(t_sq, G_sq .- G_sq[1], leg = false)

svg

図としては上と同じだが，縦軸の目盛りに意味ある数字が書かれたので，上下幅がはっきりしたぞ．

上下のズレの桁が大変小さい！

確かに $G$ が12桁程度保たれることがわかる．長期的に見たらズレは蓄積していくが、それも大変小さいことが見て取れる．

びっくりするぐらい、保存量が数値的に「良く」保存されている ことになる．

4次精度の Runge-Kutta 法に比べ、この解法はせいぜい 2次精度で精度は良くないはずなのに、保存量の保存性に関して言えば遥かに良いことになる．

さて、このからくりはどうなっているのだろう？考えておこう．

レポート

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

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

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

もとの常微分方程式に従うと保存量 $G$ が確かに保存されることを, 自力で確認しよう．

つまり，$u,v$ が元の Volterra モデルである常微分方程式の厳密解であるとき、 $ \displaystyle \frac{dG}{dt} $ の時間変化を自分で手で計算し、確かにこれがゼロになることを確認しよう, ということだね．
最後に紹介した構造保存解法 $(\ref{eq:sp-scheme})$ がなぜこんなによく $G$ を保つのか．自分なりに説明してみよう．

ヒントとしては，そうだなあ，厳密解 $u(t), v(t)$ に対しては $dG/dt = 0$ が成り立つわけだ．同じように考えて，構造保存解法 $(\ref{eq:sp-scheme})$ の解 $U_n, V_n$ に対して $\displaystyle \left\{G( U_{n+1}, V_{n+1}) - G( U_n, V_n) \right\} / \Delta t$ がどうなるかを手で計算してみてはどうかな．ただしこの $U_n, V_n$ は $\log(u), \log(v)$ に相当するのだから， $G$ の式は $(\ref{eq:invariant})$ ではなくて $(\ref{eq:invariant-by-log})$ を使うことを忘れないようにね．

普通、こんな簡単には見つからない．これは超々ラッキーなケース． ↩︎
英語だと "structure-preserving method" だね．少し古い名前だと "geometric integration" なんていうのもある． ↩︎
よく見ると，このスキームも「時間対称スキームの一種」である． ↩︎

数値計算法基礎 (2020)

08. 微分方程式の不変量