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

Volterra モデルが保存系らしいことはこれまでの線形安定性解析や数値計算から推測されるところだ． そこで少し真面目に考えてみよう． こうした考察は問題の系の「大域的性質」を調べることにあたり、数値解析にもたいへん役に立つ．

保存量の探索: 方針

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

そうした量、つまり保存量を仮に $G(u(t), v(t))$ と書くとして、一番簡単な考え方は次のようなものだろう． まず、$G$ が保存量であるということは、連鎖律で分解して、 \begin{equation} \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} \frac{\partial G}{\partial u} \cr \frac{\partial G}{\partial v} \end{array} \right) \cdot \left( \begin{array}{c} \frac{du}{dt} \cr \frac{dv}{dt} \end{array} \right) = \nabla G \cdot \frac{d}{dt} \boldsymbol{u} \end{equation} となることは間違いないので、

$\nabla G \perp \frac{d}{dt} \boldsymbol{u}$ となるような $G$ を探せば良い

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

方程式を整理しておく

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

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

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

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

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

なので、

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

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

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

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

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

変数変換して扱いやすくなった式をベースに考えることにしたので、 上の方針は $ \displaystyle \left( \begin{array}{c} \frac{\partial G}{\partial U} \\ \frac{\partial G}{\partial V} \end{array} \right) $ と $ \left( \begin{array}{c} \frac{dU}{dt} \\ \frac{dV}{dt} \end{array} \right) $ が直交するような $G(U,V)$ を探そうということに翻訳できる．

本来この探索はここから難しいわけで、いろんなやり方を試していかないといけない． そして、最初に試すべき作業は、例えば、「仮に」 \begin{equation} \left( \begin{array}{c} \frac{dU}{dt} \cr \frac{dV}{dt} \end{array} \right) = \left( \begin{array}{c} - \frac{\partial G}{\partial V} \cr \frac{\partial G}{\partial U} \end{array} \right) \end{equation} とするとこれはあきらかに $ \displaystyle \left( \begin{array}{c} \frac{\partial G}{\partial U} \\ \frac{\partial G}{\partial V} \end{array} \right) $ と直交するので、これが成り立つような $G$ があるか、というチェックだろう．

そこで、このチェックを行おう． これは書き直すと、

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

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

\begin{equation} G(U,V) = C_2 e^U - (D_2 + E) U + C_1 e^V - (D_1 - E) V \end{equation}

とすればよいことがわかる． つまり、保存量が見つかったのである(普通、こんな簡単には見つからない．これは超々ラッキーなケース)．

この保存量 $G$ をあらためて元の $u,v$ で書き直すと

\begin{equation} G(u,v) = C_2 u - (D_2 + E) \log(u) + C_1 v - (D_1 - E) \log(v) \end{equation}

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

自力でこの $G$ が確かに保存されることを確認しておこう(レポート問題)

自分の手で $ \displaystyle \frac{dG}{dt} $ を計算して、確かにこれがゼロになることを確認しておこう．

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

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

In [1]:
# 各種定数を前回同様に設定して…
C1 = C2 = D1 = D2 = 1.0
u0 = v0 = 0.7
Δt = 0.05
E = 0.5
Out[1]:
0.5
In [2]:
# Euler スキームを再び定義する．変えたくなるかもしれないので、漁業係数 E を引数にいれておく．
function volterra_euler(u,v,E)
    u_new = (1 + Δt*(D1-E) - Δt*C1*v) * u
    v_new = (1 - Δt*(D2+E) + Δt*C2*u) * v
    u,v = u_new, v_new
    return u, v
end
Out[2]:
volterra_euler (generic function with 1 method)
In [3]:
# 念のために、1ステップだけでも動かしてみる．
volterra_euler(u0,v0,E)
Out[3]:
(0.6929999999999998, 0.672)
In [4]:
# 大丈夫そうなので、500ステップほど動かしてみよう
u,v = u0,v0
uv_sq = [u0, v0]

for i in 1:500
    u,v = volterra_euler(u,v,E)
    uv_sq = hcat(uv_sq, [u,v])
end
In [5]:
# 早速、結果をプロットしてみよう．まずはその準備．
using Plots
gr()
Out[5]:
Plots.GRBackend()
In [6]:
plot(uv_sq') # プロットしてみると… やっぱり Euler 法は振幅がだんだん大きくなっている？
Out[6]:
<?xml version="1.0" encoding="utf-8"?> 100 200 300 400 500 0.5 1.0 1.5 2.0 2.5 3.0 3.5 y1 y2

なにはともあれ、保存量 $G$ の挙動をチェックしてみよう

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

In [7]:
G(u,v,E) = C2*u - (D2+E)*log(u) + C1*v - (D1-E)*log(v)
Out[7]:
G (generic function with 1 method)

Euler法の結果だと、この GG がどう変化するか、みてみよう． まず、次のようにすれば、GG の列を作ることが出来る．

In [8]:
G_euler_sq = [G(uv_sq[1,n], uv_sq[2,n],E) for n in 1:501]
Out[8]:
501-element Array{Float64,1}:
 2.11335
 2.11384
 2.11431
 2.11477
 2.11523
 2.11568
 2.11613
 2.11657
 2.11702
 2.11746
 2.1179 
 2.11835
 2.11879
 ⋮      
 2.64204
 2.64292
 2.64376
 2.64459
 2.64541
 2.64621
 2.64701
 2.64781
 2.64861
 2.64941
 2.65022
 2.65103

んん～？ 間違えたかな？ 値が変わっているのが明らかだな． 念のためプロットしてみよう．

In [9]:
plot(G_euler_sq)
Out[9]:
<?xml version="1.0" encoding="utf-8"?> 100 200 300 400 500 2.2 2.3 2.4 2.5 2.6 y1

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

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

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

\begin{equation} \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. \end{equation}

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

\begin{equation} \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. \end{equation}

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

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

In [10]:
function symm_zero(u_new, u_curr, E)
    up, vp, u, v = u_new[1], u_new[2], u_curr[1], u_curr[2]
    half = Δt/2
    r = similar(u_new)

    r[1] = ( 1 + half * ( - D1 + E + C1*vp ) ) * up + ( -1 + half * ( - D1 + E + C1*v ) ) * u
    r[2] = ( 1 + half * (   D2 + E - C2*up ) ) * vp + ( -1 + half * (   D2 + E - C2*u ) ) * v
    return r
end
Out[10]:
symm_zero (generic function with 1 method)

NLsolve パッケージを使って計算するとしよう．いままでと同様に、まずは次のように手続きをして…

In [11]:
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
Out[11]:
nls (generic function with 1 method)
In [14]:
# で、nlsolve を使ってこのスキームを作る．
function volterra_symm(u,v, E)
  u_curr = [u,v]
  u_new = nls(symm_zero, u_curr, E, ini = u_curr)[1]
  return u_new[1], u_new[2]
end
Out[14]:
volterra_symm (generic function with 1 method)
In [15]:
# まずは、動くかどうか簡単にチェック．
volterra_symm(u0,v0,E)
Out[15]:
(0.6935102416466368, 0.6724420591556107)
In [16]:
# 動くようなので、500ステップほど動かしてみよう
u,v = u0,v0
uv_symm_sq = [u0, v0]

for i in 1:500
    u,v = volterra_symm(u,v,E)
    uv_symm_sq = hcat(uv_symm_sq, [u,v])
end
In [17]:
# 結果をプロットしてみる
plot(uv_symm_sq')  # …おお、振幅が一定だ．Euler 法よりは信用できるっぽい？
Out[17]:
<?xml version="1.0" encoding="utf-8"?> 100 200 300 400 500 0.5 1.0 1.5 2.0 2.5 y1 y2

保存量 $G$ をチェックしよう．やり方はさっきと同じで、$G$ の変化データを作ってからプロットすれば良い．

In [18]:
# まず、G の時間変化のデータを作る
G_symm_sq = [G(uv_symm_sq[1,n], uv_symm_sq[2,n],E) for n in 1:501]
Out[18]:
501-element Array{Float64,1}:
 2.11335
 2.11336
 2.11336
 2.11337
 2.11337
 2.11337
 2.11338
 2.11338
 2.11338
 2.11339
 2.11339
 2.11339
 2.1134 
 ⋮      
 2.11341
 2.11341
 2.11341
 2.11341
 2.11341
 2.11341
 2.11342
 2.11342
 2.11342
 2.11342
 2.11342
 2.11343
In [19]:
# 保存量 G の時間変化を見てみると…
plot(G_symm_sq)
Out[19]:
<?xml version="1.0" encoding="utf-8"?> 100 200 300 400 500 2.1130 2.1131 2.1132 2.1133 2.1134 y1

う～む、Euler法よりは格段にマシだけれども． 保存されるはずの GG が微妙に振動していることが分かる． さて、どうしたものだろう？

構造保存解法の出番か

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

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

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

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

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

\begin{equation} \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. \end{equation}

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

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

で定義される．

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

\begin{equation} \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. \end{equation}

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

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

In [20]:
function dedx(a,b)
  EPS = 1.0e-16
  return ifelse( abs(a-b) > EPS, (exp(a)-exp(b))/(a-b), exp(a) )
    # ifelse 関数を調べておこう
end
Out[20]:
dedx (generic function with 1 method)

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

In [21]:
# ゼロになるべき関数を作って…
function sp_zero(u_new, u_curr, E)
    up, vp, u, v = map(log, (u_new[1], u_new[2], u_curr[1], u_curr[2])) 
    # map というのは、複数の対象に関数を一度に作用させる命令だ．便利なので覚えておこう．
    
    r = similar(u_new)
    r[1] = up - u + Δt*C1* dedx(vp,v) + Δt*(-D1 + E)
    r[2] = vp - v - Δt*C2* dedx(up,u) + Δt*( D2 + E)
    return r
end
Out[21]:
sp_zero (generic function with 1 method)
In [24]:
# スキームを関数として与える．
function volterra_sp(u, v, E)
  u_curr = [u,v]
  u_new = nls(sp_zero, u_curr, E, ini = u_curr)[1]
  return u_new[1], u_new[2]
end
Out[24]:
volterra_sp (generic function with 1 method)
In [25]:
# 試してみよう．
volterra_sp(u0, v0, E)
Out[25]:
(0.6935156878695156, 0.6724434211370223)
In [26]:
# 動くようなので、500ステップほど動かしてみよう
u,v = u0,v0
uv_sp_sq = [u0, v0]

for i in 1:500
    u,v = volterra_sp(u,v,E)
    uv_sp_sq = hcat(uv_sp_sq, [u,v])
end
In [27]:
# 結果をプロットしてみる
plot(uv_sp_sq')  # …おお、これも振幅が一定だ．
Out[27]:
<?xml version="1.0" encoding="utf-8"?> 100 200 300 400 500 0.5 1.0 1.5 2.0 2.5 y1 y2
In [28]:
# 早速、保存量がどうなっているか見てみよう．
G_sp_sq = [G(uv_sp_sq[1,n], uv_sp_sq[2,n], E) for n in 1:501]
Out[28]:
501-element Array{Float64,1}:
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 ⋮      
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
 2.11335
In [29]:
# 上の結果を見ると、ほぼ 2.11335 のままのようなので、この値をひいて、ズレをプロットしよう
plot(G_sp_sq - 2.11335)
Out[29]:
<?xml version="1.0" encoding="utf-8"?> 100 200 300 400 500 -0.000000112123 -0.000000112122 -0.000000112121 -0.000000112120 y1

ズレの桁が大変小さくて、びっくりするぐらい、保存量が数値的に保存されている ことがわかる．さて、このからくりはどうなっているのだろう？ 考えておこう．