微分方程式の不変量: Volterra eq.の場合

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\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\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 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\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$ を探す: 今回はラッキー

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

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

まあ、最初に試すべき作業は、以下のように考えると良さそうだ．まず、「仮に」 \[ \left( \begin{array}{c} \displaystyle \frac{dU}{dt} \cr\cr \displaystyle \frac{dV}{dt} \end{array} \right) = \left( \begin{array}{c} \displaystyle - \frac{\partial G}{\partial V} \cr\cr \displaystyle \frac{\partial G}{\partial U} \end{array} \right) \] という等式が成り立つと仮定すると、この状況下ではこの左辺( = 右辺)はあきらかに

$ \displaystyle \left( \begin{array}{c} \displaystyle \frac{\partial G}{\partial U} \cr\cr \displaystyle \frac{\partial G}{\partial V} \end{array} \right) $ と直交する．これは この等式が成り立つような $G$ が存在するか をチェックして、見つかるならばそれは保存量だ、ということを意味している．そこでこの「チェック」が最初に試すべきこと、だと考えてもいいだろう．

では、この等式を成り立たせる $G$ が存在するかのチェックを早速行ってみよう．この等式チェックを書き直して考えると、

\[ \left\{\begin{array}{rcl} \displaystyle - \frac{\partial G}{\partial V} & = -C_1 e^V + (D_1 - E), \cr\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 \]

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

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

\[ G(u,v) = C_2 u - (D_2 + E) \log(u) + C_1 v - (D_1 - E) \log(v) \]

となる.

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

まずは 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.674851 0.568983 0.672523 0.545508 0.670993 0.522938 0.670223 0.501262 0.670181 0.480465 0.670836 0.460531 0.672159 0.441438 0.674128 0.423166 ⋮ 0.406685 0.724975 0.40211 0.685343 0.398384 0.647722 0.395441 0.612045 0.393226 0.578243 0.391687 0.546243 0.390782 0.515973 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


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


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

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\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\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


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

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.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 = vcat( 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)

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


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

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

svg

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

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

1

plot(t_sq, G_sq)

svg

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

1

plot(t_sq, G_sq .- G_sq[1] )

svg

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

構造保存解法の出番か

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

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

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

\[ \left\{\begin{array}{rcl} \displaystyle U & := & \log(u), \cr\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\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. \]

が、それである．ただし、$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\cr f’(a) & & {\rm otherwise } \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\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


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

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

svg

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

1

plot(t_sq, G_sq)

svg

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.113349887877465 2.113349887877454 2.1133498878774444 2.113349887877436 2.113349887877429 2.1133498878774235 2.1133498878774195 2.1133498878774173 2.1133498878774164 2.113349887877417 2.113349887877419 2.113349887877422 2.113349887877427

2.113349887876912 2.113349887876921 2.11334988787693 2.1133498878769394 2.1133498878769488 2.113349887876958 2.113349887876967 2.113349887876976 2.113349887876984 2.113349887876992 2.113349887876999 2.1133498878770056

上12桁ぐらいは変化していないな．

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

1

plot(t_sq, G_sq .- G_sq[1] )

svg

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

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

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

Report No.06: 保存量について…

自力で保存量 $G$ が確かに(もとの常微分方程式に従うと)保存されることを確認しよう．
$u,v$ が元の Volterra モデルである常微分方程式の厳密解であるとき、 $ \displaystyle \frac{dG}{dt} $ の時間変化を自分で手で計算し、確かにこれがゼロになることを確認しよう．
最後の解法が、なぜこんなによく $G$ を保つのか．自分なりに説明してみよ．