微分方程式の不変量: Volterra モデルの場合

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} \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) = \nabla G \cdot \frac{d}{dt} \boldsymbol{u} \end{equation} となることは間違いないので、

$\nabla G \perp \displaystyle \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), \cr \displaystyle \frac{dv}{dt} & = v \left( C_2 u - (D_2 + E) \right) \end{array}\right. \label{eq:volterra} \end{equation}

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

\begin{equation} \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. \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), \cr \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), \cr \displaystyle \frac{dV}{dt} & = C_2 e^U - (D_2 + E) \end{array}\right. \end{equation}

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

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

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

本来この探索はここから難しいわけで、いろんなやり方を試していかないといけない．

まあ、最初に試すべき作業は、例えば、「仮に」 \begin{equation} \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:assumption} \end{equation} という等式が成り立つと仮定すると、この状況下ではこの左辺( = 右辺)はあきらかに $ \displaystyle \left( \begin{array}{c} \displaystyle \frac{\partial G}{\partial U} \cr \displaystyle \frac{\partial G}{\partial V} \end{array} \right) $ と直交することから、$(\ref{eq:assumption})$ が成り立つような $G$ が存在するか をチェックすることだろう．

注: この流れからわかるように、もしそうした $G$ が存在するならば、それは保存量である．

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

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

となる.

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

まずは 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) - Δt*C1*v) * u
    v_new = (1 - Δt*(D2+E) + Δt*C2*u) * v
    return u_new, v_new
end

volterra_euler (generic function with 1 method)

1
2

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

(0.6929999999999998, 0.672)

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

1

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

G (generic function with 1 method)

では、いつものように動かしてみよう．あとでグラフにしてみたいので、$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 = hcat(uv_sq, [u,v])
    push!( G_sq, G(u,v,E) )
end

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

1

uv_sq

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

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

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

1
2

using Plots
gr()

Plots.GRBackend()

そしてプロットだ．

1
2

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

先にやったとおり、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")

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

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

1

iniG = G(u0, v0, E)

2.113349887877465

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

1
2
3
4
5
6
7
8
9

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")

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

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

1

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

…やっぱり、保存量のはずなのに $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}$ を求めれば時間ステップを一つ進めることが出来る．

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

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
    eq2 = ( 1 + half * (   D2 + E - C2*up ) ) * vp + ( -1 + half * (   D2 + E - C2*u ) ) * v

    return [ eq1, eq2 ]
end

symm_zero (generic function with 1 method)

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

nls (generic function with 1 method)

で、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

volterra_symm (generic function with 1 method)

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

1
2

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

2-element Array{Float64,1}: 0.69351 0.672442

動くようだな．では、もう少し動かしてみる…前に、G の計算式をもう一つ作っておくと楽そうだ．

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

1

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

G (generic function with 2 methods)

では、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 = hcat(uv_sq, uv)
    push!(G_sq, G(uv, E))
end

結果をプロットしてみる

1

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

…おお、振幅が一定だ．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")

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

ちなみに、$G$ = 一定 の線と重ねて描画してみよう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

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")

目で見る限り、もう区別がつかないね(黄色い線と青い線が重なって緑色に見えている)．この近似解はだいぶ良いと言ってよさそうだ．

さて、保存量 $G$ をチェックしよう．

1

plot(t_sq, G_sq)

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

さて、格段に良いことは良いのだけれども、本来は $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

r (generic function with 1 method)

前もそうしたが、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

RK (generic function with 1 method)

早速動かしてみよう！

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 = hcat(uv_sq, uv)
    push!(G_sq, G(uv, E))
end

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

1

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

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

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

1

plot(t_sq, G_sq)

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

1

plot(t_sq, G_sq - 2.11335)

こうしてみると、ブレはかなり小さくて、8桁程度は $G$ が保たれている．しかし、ズレは蓄積していくようなので、長期計算には向かないかなあ．

構造保存解法の出番か

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

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

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

\begin{equation} \left\{\begin{array}{rcl} \displaystyle U & := & \log(u), \cr \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}$ を求めれば時間ステップを一つ進めることが出来るというスキームの出来上がりである．

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

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

dedx (generic function with 1 method)

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

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

sp_zero (generic function with 1 method)

これを用いて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

volterra_sp (generic function with 1 method)

試しておこう．

1

volterra_sp(uv_0, E)

2-element Array{Float64,1}: 0.693516 0.672443

動くようなので、これも 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 = hcat(uv_sq, uv)
    push!( G_sq, G(uv, E) )
end

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

1

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

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

1

plot(t_sq, G_sq)

GKS: Rectangle definition is invalid in routine SET_WINDOW GKS: Rectangle definition is invalid in routine SET_WINDOW GKS: Rectangle definition is invalid in routine SET_WINDOW GKS: Rectangle definition is invalid in routine SET_WINDOW

おや、プロットされないぞ． これは、plot 命令にとって「値があまりにも変わらないので、上下幅をどうしたらよいかわからない」ということによるようだ．

実際、$G$ の値を数字で見てみよう．

1

G_sq

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

少なくとも、表示される6桁程度では差が見えないということだな．

では、この値 2.11335 を引いて、プロットしてみよう．

1

plot(t_sq, G_sq - 2.11335 )

ズレをもう少しはっきりとみるために、もう少し丁寧にプロットし直してみよう．

とりあえず、$G$ の初期値との差を見てみようか．

1

plot(t_sq, G_sq - G_sq[1] )

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

つまり、びっくりするぐらい、保存量が数値的に保存されている ことになる． 4次精度の Runge-Kutta 法に比べ、この解法はせいぜい 2次精度で精度は良くないはずなのに、保存量の保存性に関して言えば遥かに良いことになる．

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

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

$u,v$ が元の Volterra モデルである常微分方程式 $(\ref{eq:volterra})$ の厳密解であるとき、 $ \displaystyle \frac{dG}{dt} $ の時間変化を自分で手で計算し、確かにこれがゼロになることを確認しよう．

Report No.08: 最後の解法が、なぜこんなによく $G$ を保つのか．自分なりに説明してみよ．