08. 微分方程式の不変量
Photo by Patrick Federi on Unsplash
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\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} \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. \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 \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. \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$ を探す: 今回はラッキー
変数変換して扱いやすくなった式をベースに考えることにしたので、 上の方針は
\begin{equation} \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} \end{equation}
であるような $G(U,V)$ を探そうという方針に翻訳できる．
本来この探索はここから難しいわけで、いろんなやり方を試してみないといけないし，それらも簡単ではないが，まあ今回はなんとかなるのだ．
さて，こういうときに最初に試すべき作業は、以下のようなものだろう．
まず、「仮に,とある歪対称な行列$S$を用いて」この式の右辺について
\begin{equation}
\left( \begin{array}{c} \displaystyle \frac{dU}{dt} \cr \displaystyle \frac{dV}{dt} \end{array} \right) =
S
\left( \begin{array}{c} \displaystyle \frac{\partial G}{\partial U} \cr \displaystyle \frac{\partial G}{\partial V} \end{array} \right)
\label{eq:hypo}
\end{equation}
という等式が成り立つと仮定してみよう．
ちなみに，$n\times n$行列 $S$ が歪対称行列(skew-symmetric matrix)もしくは交代行列(alternating matrix)であるとは，
任意の $n$次元ベクトル $\boldsymbol{x}$ に対して
\begin{equation}
\boldsymbol{x}^{\top} S \boldsymbol{x}
\left( = \boldsymbol{x} \cdot S \boldsymbol{x} \right)
= 0
\label{eq:def-skew}
\end{equation}
であることである．
これは，$S$ の $(i,j)$成分を $s_{ij}$ と書いたとき，
\begin{equation}
s_{ij} = - s_{ji}\,\, \mbox{ for } {}^{\forall}i,j
\end{equation}
であることと必要十分である．
さて，この仮定 $(\ref{eq:hypo})$ が成り立つならば式 $(\ref{eq:hypo})$ の右辺は歪対称行列の定義から $ \displaystyle \left( \begin{array}{c} \displaystyle \frac{\partial G}{\partial U} \cr \displaystyle \frac{\partial G}{\partial V} \end{array} \right) $ と直交する． なぜなら， $S$ と任意のベクトル $\boldsymbol{x}$ は式 $(\ref{eq:def-skew})$を 満たすので， ベクトル $ \displaystyle \left( \begin{array}{c} \displaystyle \frac{\partial G}{\partial U} \cr \displaystyle \frac{\partial G}{\partial V} \end{array} \right) $ に対しても
\begin{equation} \left( \begin{array}{c} \displaystyle \frac{\partial G}{\partial U} \cr \displaystyle \frac{\partial G}{\partial V} \end{array} \right)^{\top} S \left( \begin{array}{c} \displaystyle \frac{\partial G}{\partial U} \cr \displaystyle \frac{\partial G}{\partial V} \end{array} \right) = \left( \begin{array}{c} \displaystyle \frac{\partial G}{\partial U} \cr \displaystyle \frac{\partial G}{\partial V} \end{array} \right) \cdot S \left( \begin{array}{c} \displaystyle \frac{\partial G}{\partial U} \cr \displaystyle \frac{\partial G}{\partial V} \end{array} \right) = 0 \end{equation}
となり，これは
\begin{equation} \left( \begin{array}{c} \displaystyle \frac{\partial G}{\partial U} \cr \displaystyle \frac{\partial G}{\partial V} \end{array} \right) \perp S \left( \begin{array}{c} \displaystyle \frac{\partial G}{\partial U} \cr \displaystyle \frac{\partial G}{\partial V} \end{array} \right)\,\, \left( \Longleftrightarrow\,\, \nabla G \perp S \,\nabla G \right) \end{equation}
と同値だからである．
つまり，仮定 $(\ref{eq:hypo})$ が成り立つならば 狙っていた直交性 $(\ref{eq:perp})$ は自動的に満たされるということになる．
だから，この仮定の等式 $(\ref{eq:hypo})$ が成り立つような $S$, $G$ の組み合わせを見つけることができれば その $G$ は保存量だ、ということになる． そこで「式 $(\ref{eq:hypo})$ を満たす $S$, $G$ の組み合わせを見つけようと努力すること」が最初に試すべきことなのだ．
$S$, $G$ の両方が不明なのでこの探索は本来大変なのだが，実は $S$ については参考になる前例がある． それが以下の内容だ．
物理学などでよく見られる例として，系の状態を示す物理量を並べたベクトル $\boldsymbol{U}$ に対して
\begin{equation} \frac{d \boldsymbol{U}}{dt} = \left( \begin{array}{c} 0 & 1 \cr -1 & 0 \end{array} \right) \nabla_{\boldsymbol{U}} G \end{equation}
というモデル方程式が多々あり，この形で支配方程式が書かれる系を ハミルトン系 などと呼ぶ．
つまりこの場合は
\begin{equation}
S = \left( \begin{array}{c} 0 & 1 \cr -1 & 0 \end{array} \right)
\end{equation}
なのだ．今回の話ではこれを参考にしよう．
ちなみに，この方程式で支配される系の代表例はみんなもよく知る Newton運動方程式系だ．
つまり，上の $S$ を参考にして変形が可能ならば Volterra モデル方程式はハミルトン系であり，Newton運動方程式の親戚のようなものだと言える．
そして，ハミルトン系については解析力学という研究分野によっていろいろな性質が知られていて，量子力学をはじめ，現代物理学の基礎ともなっている．
気になる人は少し調べてみよう．
さて話を戻そう． 早速だが，今回は(あとの都合もあって符号だけ変えて)上の $S$ を参考にして \begin{equation} S = \left( \begin{array}{c} 0 & -1 \cr 1 & 0 \end{array} \right) \end{equation} とし，この $S$ を採用して話を進めよう．
そしてこの $S$ を代入した 仮定の等式 $(\ref{eq:hypo})$ を成り立たせる $G$ を見つけられるか考えよう． $S$ を代入した式 $(\ref{eq:hypo})$ を使ってもとの Volterra 方程式を書き直してみると、
\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. \label{eq:condhypo} \end{equation}
となる． そしてわれわれは，これを満たす $G(U,V)$ が見つけられるか，を考えるわけだ.
そしてこの式 $(\ref{eq:condhypo})$ を見ると，都合の良い特徴を持っていることに気づく． それは 第1式は($G$にとっての)独立変数 $V$ のみの関係式で, 第2式は独立変数 $U$ のみの関係式となっているので， $G$ の $U$ と $V$ への依存性を別々に考えて \begin{equation} G(U,V) = G_1(U) + G_2(V) \end{equation} とすれば(一種の変数分離だ)，
\begin{equation} \left\{\begin{array}{rcl} \displaystyle - \frac{\partial G_2}{\partial V} & = -C_1 e^V + (D_1 - E), \cr \displaystyle \frac{\partial G_1}{\partial U} & = C_2 e^U - (D_2 + E) \end{array}\right. \label{eq:condhypotwo} \end{equation}
という，異なる関数 $G_1(U)$, $G_2(V)$ それぞれに対する独立な関係式がある，という格好にできるのだ．
あとは この2つの式 $(\ref{eq:condhypotwo})$ それぞれを満たす $G_1(U)$, $G_2(V)$ を求めるだけだがこれは誰でもできる簡単な問題で， \begin{equation} \left\{\begin{array}{rcl} \displaystyle G_1(U) & = & C_2 e^U - (D_2 + E) U + \mbox{const.} , \cr \displaystyle G_2(V) & = & C_1 e^V - (D_1 - E) V + \mbox{const.} \end{array}\right. \end{equation} と解が求められる．
よって，大変ラッキーにも保存量 G が求まったのだ！1 具体的には以下のようになる．
\begin{equation} G(U,V) = C_2 e^U - (D_2 + E) U + C_1 e^V - (D_1 - E) V + \mbox{const.} \label{eq:invariant-by-log} \end{equation}
そしてもちろん，この議論の流れから明らかだが， この $G$ は仮定の等式 $(\ref{eq:hypo})$ を満たすのだ．
さて，$U, V$ の表記のままだとわかりにくいので，この保存量 $G$ をあらためて元の変数 $u,v$ を使ったものに書き直し2，以降はこれを用いよう． 具体的には
\begin{equation} G(u,v) = C_2 u - (D_2 + E) \log(u) + C_1 v - (D_1 - E) \log(v) \label{eq:invariant} \end{equation}
となる.
さて、以降はこの量が数値計算でどういう挙動をみせるかを調べていこう．
まずは Euler 法ではどうなのか、再チェックだ
というわけで、Euler スキームを準備して、計算してみよう．
1# 各種定数を前回同様に設定して…
2C1 = C2 = D1 = D2 = 1.0
3u0 = v0 = 0.7
4Δt = 0.05
5E = 0.5
1# Euler スキームを定義する．
2function volterra_euler(u,v,E)
3 u_new = (1 + Δt * ( D1 -E - C1*v ) ) * u
4 v_new = (1 + Δt * ( -D2 -E + C2*u ) ) * v
5 return u_new, v_new
6end
1# いつものように，念のために 1ステップ動かしてみる．
2volterra_euler(u0,v0,E)
(0.693, 0.6719999999999999)
また、関数 $G$ を定義して、時間発展計算中にこの値を計算できるようにしよう．念のため、漁業係数も変えられるようにしておこう．
1G(u,v,E) = C2*u - (D2+E)*log(u) + C1*v - (D1-E)*log(v)
では、いつものように動かしてみよう．あとでグラフにしてみたいので、$G$ の値も取っておくとしよう．
1# 大丈夫そうなので、500ステップほど動かしてみよう
2u,v = u0,v0
3uv_sq = [u0 v0]
4G_sq = [ G(u,v,E) ]
5
6for i in 1:500
7 u,v = volterra_euler(u,v,E)
8 uv_sq = vcat(uv_sq, [ u v ])
9 push!( G_sq, G(u,v,E) )
10end
無事に動いているか、簡単に確認してみる．
1uv_sq
501×2 Array{Float64}:
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
とりあえず結果に全て数字が入っているので、計算はできているようだ．
では早速、グラフで状況を見ていこう． まずいつものお約束．
1using Plots
そしてプロットだ．
1default( legend = :outertopright,
2 xlabel = "time t", ylabel = "number" )
3
4t_sq = Δt * [0:500]
5plot(t_sq, uv_sq, label = ["魚 🐟" "鮫 🦈"], fontfamily="" )
先にやったとおり、Euler 法だとだんだん値の振れが大きくなっているね．
さて、この Euler法による近似解で、本来は保存されるはずの量 $G(u,v)$ がどれくらい保たれているかをみてみよう．
1plot(t_sq, G_sq, yaxis = "G", leg = false)
やっぱり、保存量のはずなのに Euler法による近似解では $G$ が大きく増加している． これははっきり言って良くないねえ． Euler 法があまり信用ならないなあ、ということがこのことからもよく分かる．
さて次に，相図で保存量がどう見えるのか見ていこう． まず最初は Euler法の近似解を単独で見てみる．
1u_sq = uv_sq[:,1] # データの 1行目 = u の近似値列
2v_sq = uv_sq[:,2] # データの 2行目 = v の近似値列
3
4plot(u_sq, v_sq)
5plot!((u0,v0), marker = :circle, fontfamily="",
6 aspect_ratio = 1, legend = false,
7 xaxis = "魚 🐟", yaxis = "鮫 🦈" )
8# aspect_ratio は縦横の軸の比率．
9
10annotate!(u0,v0, ("初期値", :bottom, :left))
この相図からも解の軌道がだんだん大きくなっていることがわかるね． まあ，どちらへ移動しているかなどがわかりやすいという変なメリットもあるが．
さて，この相図上において，きちんとした計算方法で求めた近似解ならば，その軌道はぐるっと周を描いて戻ってくるだろうというのは前回得られた予想だった．
そして今われわれは保存量 $G$ を一つ求めることに成功した(2つ目以降があるのかどうかは知らないが)．
さらに，この予想と保存量 $G$ との間に直接の関係があるのでは？，というのが素直な疑問というものだ． そこでこの疑問をチェックしよう．
具体的には，解の軌道が保存量 $G$ が一定の$(u,v)$の集合と一致するのか，を調べてみるのが直接的で良い．
うまくいくならばそれを Euler法を始め，いろんな解法で求めた近似解の軌道と比べるなどして知見が得られそうではないか．
さて，ではまず、一定であるべき $G$ の値を求めておこう． これは簡単で，初期値 $(u_0, v_0)$ を使って直接求めれば良い． じっさいは
1iniG = G(u0, v0, E)
2.113349887877465
と求まる．
そこで次に，$G(u,v)$ が(ほぼ)この値を持つ点 $(u,v)$ をかき集めて相図上に描いてみる． こうした「条件を満たす点の集合」を描く方法はいろいろあるが，今回は強引に次のようにしてみよう．
1U = 0.5:0.01:3.0 # 解の値をみて適当に
2V = 0.01:0.01:1.5
3
4Z = [ G(u,v,E) for u in U, v in V ]
5# U, V の範囲全部の G を計算した結果を行列に！ 強引だなあ．
6
7contour(U, V, Z', aspect_ratio = 1, legend = false,
8 xaxis = "魚 🐟", yaxis = "鮫 🦈",
9 clims = (iniG-0.001, iniG+0.001))
10# そして，G が iniG に近い値の場合だけグラフ出力する．
11# clims で、出力値の描画を制限できるのだ．
12
13plot!((u0,v0), marker = :circle, fontfamily = "")
14annotate!(u0,v0, ("初期値", :bottom, :left) )
なにかそれっぽいね(あとでより良い解法での相図と比べるとよりはっきりする)． ということは，
Volterra 問題では $G = \mbox{constant.}$ が相図上での解軌道の式だと言えそうだ！
ということになり，保存量を求めることがすなわち相図での解軌道を求めることそのものになっていたのだ，ということになる． めでたいね3．
あと，念の為に，この図と Euler法による近似解の相図を重ねて見ておこう．
1contour(U, V, Z', aspect_ratio = 1, legend = false,
2 xaxis = "魚 🐟", yaxis = "鮫 🦈",
3 clims = (iniG-0.001, iniG+0.001))
4
5plot!((u0,v0), marker = :circle, fontfamily = "")
6annotate!(u0,v0, ("初期値", :bottom, :left) )
7
8plot!(u_sq, v_sq) # Euler法による近似解
黄色の線が(おそらく)正しい解軌道で，緑色の線が Euler法による近似解の軌道だ． 最初の 1/3周ぐらいはまあ一致しているが，それ以降は 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}$ を求めれば時間ステップを一つ進めることが出来る．
というわけで、まずは ゼロ になるべきこの右辺を関数として定義しよう．
1function symm_zero(uv_new, uv_curr, E)
2# uv_new が新しい u,v で、uv_curr は古い u,v
3
4 up, vp = uv_new
5 u, v = uv_curr
6 half = Δt/2
7
8 eq1 = ( 1 + half * ( - D1 + E + C1*vp ) ) * up +
9 ( -1 + half * ( - D1 + E + C1*v ) ) * u
10 # eq1 = ... の右端の "+" に注意．
11 # この式がそこの改行で終わっていないことを julia に伝えている.
12
13 eq2 = ( 1 + half * ( D2 + E - C2*up ) ) * vp +
14 ( -1 + half * ( D2 + E - C2*u ) ) * v
15
16 return [ eq1, eq2 ]
17end
NLsolve パッケージを使って計算するとしよう．いままでと同様に、まずは次のように手続きをして…
1using NLsolve
2
3function nls(func, params...; ini = [0.0])
4 if typeof(ini) <: Number
5 r = nlsolve((vout,vin)->vout[1]=func(vin[1],params...), [ini])
6 v = r.zero[1]
7 else
8 r = nlsolve((vout,vin)->vout .= func(vin,params...), ini)
9 v = r.zero
10 end
11 return v, r.f_converged
12end
そして、nlsolve を使ってこのスキームを作る． 以下のように，u, v の今の値を渡すと新しい値が返ってくるように書くだけだ．
1function volterra_symm(uv_curr, E)
2
3 uv_new = nls(symm_zero, uv_curr, E, ini = uv_curr)[1]
4 return uv_new
5
6end
まずは、動くかどうか簡単にチェック．nlsolve をセッション中で初めて動かすときはいつものように少し待たされるので少し待とう．
1uv_0 = [u0, v0]
2volterra_symm(uv_0, E)
2-element Array{Float64,1}:
0.6935102416466368
0.6724420591556107
うん，無事に動くようで安心安心． では、もう少し動かしてみる…前に、引数が楽に書ける G の定義をもう一つ作っておくと楽そうだ．
1G(uv,E) = G(uv[1], uv[2], E)
この定義は、Julia の multi-dispatch 機能、つまり、「同じ関数名でも引数の型が違うときは違う定義を与えることができる」ことをうまく使っている
では、500ステップほど動かしてみよう
1uv = copy(uv_0)
2uv_sq = copy(uv_0') # 横向きにしておく
3G_sq = [ G(uv, E)]
4
5for i in 1:500
6 uv = volterra_symm(uv, E)
7 uv_sq = vcat( uv_sq, uv' )
8 push!(G_sq, G(uv, E))
9end
結果をプロットしてみる
1plot(t_sq, uv_sq, label = ["魚 🐟" "鮫 🦈"], fontfamily="" )
おお、振幅がかなり一定にみえる．Euler 法よりは信用できるっぽいぞ．
早速，相図も描こう．
1u_sq = uv_sq[:,1] # データの 1行目 = u の近似値列
2v_sq = uv_sq[:,2] # データの 2行目 = v の近似値列
3
4plot(u_sq, v_sq)
5plot!((u0,v0), marker = :circle, fontfamily="",
6 aspect_ratio = 1, legend = false,
7 xaxis = "魚 🐟", yaxis = "鮫 🦈" )
8# aspect_ratio は縦横の軸の比率．
9
10annotate!(u0,v0, ("初期値", :bottom, :left))
おお！ Euler 法と異なり、「元へ戻る」ことが明確にわかるぞ．
ちなみに、保存量 $G$ = 一定 の線と重ねて描画してみよう．
1contour(U, V, Z', aspect_ratio = 1, legend = false,
2 xaxis = "魚 🐟", yaxis = "鮫 🦈",
3 clims = (iniG-0.001, iniG+0.001))
4
5plot!((u0,v0), marker = :circle, fontfamily = "")
6annotate!(u0,v0, ("初期値", :bottom, :left) )
7
8plot!(u_sq, v_sq) # 対称スキームによる近似解
黄色い線と青い線が重なって緑色に見えていて，目で見る限り区別がつかないね． ということは，
- この時間対称スキームによる近似解はだいぶ良いと言ってよさそう
- 保存量 $G = \mbox{const.}$ は相図の解軌道を表すとみてまず間違いない
と言えそうだ．
さて、この近似解が保存量 $G$ からどれくらいずれるか，その時間変化をより丁寧にチェックしよう． Euler法のとき同様，近似解 $(u,v)$ から計算される $G(u,v)$ をプロットしてみるのだ．
1plot(t_sq, G_sq, yaxis = "G", leg = false)
おお、Euler法よりは格段に良いな！ おおまかに 3桁ぐらいは保たれている．長期的に見ればもとへ戻ってきている点も素晴らしい．
さて、Euler法に比べて格段に良いことは良いのだけれども、本来は $G$ は不変量であるのだから、「ブレ」をもうちょっと小さくできないか、もうちょっと欲張ってもよいだろう． なんとかならないものだろうか？
Runge-Kutta 法だとどうだろう？
古典的 Runge-Kutta 法を試してみよう． 4次精度なのでかなり良いだろうと期待できる．
まず、常微分方程式そのものの「右辺」を次のように定義して、
1function r(uv)
2
3 u,v = uv
4
5 eq1 = u * ( -C1*v + (D1-E))
6 eq2 = v * ( C2*u - (D2+E))
7
8 return [ eq1, eq2 ]
9
10end
前もそうしたが、Runge-Kutta 法を直接定義してしまおう．
1function RK(uv)
2 r1 = r(uv)
3 r2 = r(uv + Δt/2 * r1)
4 r3 = r(uv + Δt/2 * r2)
5 r4 = r(uv + Δt * r3)
6 return uv + (Δt/6) * (r1 + 2*r2 + 2*r3 + r4)
7end
早速動かしてみよう！
1uv = copy(uv_0)
2uv_sq = copy(uv_0')
3G_sq = [ G(uv, E)]
4
5for i in 1:500
6 uv = RK(uv)
7 uv_sq = vcat(uv_sq, uv')
8 push!(G_sq, G(uv, E))
9end
そしてグラフを出してみる．
1plot(t_sq, uv_sq, leg = false)
ふむ、これも結果は良いことが一目で分かる． 相図はもう面倒なので描かない(うまくいっていることは想像に難くない).
では、肝心の，保存量である $G$ の時間変化はどうかな？
1plot(t_sq, G_sq, yaxis = "G", leg = false)
ほんの少しずつ増えているが、ブレは小さい，といったところか．
ただ、このグラフだと目盛りの数字が変わらなくてその程度がよくわからない． そこで全体から初期値を引いて、そのグラフを描いてみよう．
1plot(t_sq, G_sq .- G_sq[1], yaxis = "ΔG", leg = false)
こうしてみると、はっきりわかる． ブレはかなり小さくて、7～8桁程度は $G$ が保たれている．すごいね！
贅沢を言うならば，ズレは蓄積していくようなので長期計算には向かないかもしれない，という懸念があるぐらいかな．
構造保存解法の出番か
こうした問題に対して、保存量は保存するように、散逸量は散逸するように、という数値スキームを設計する方法がある．こうしたスキームを「構造保存数値解法」と呼んだりして4，数値解析の理論派では今かなりホットな分野だ．
さて、どうやって導出するかは後述するが、今回の問題に対して、最もシンプルな構造保存スキームは以下のような物になる．
まず、保存量の導出のところで扱ったように変数変換を考える．
\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. \label{eq:sp-scheme} \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 } & & {\rm when }\,\, a \neq b , \cr f'(a) & & {\rm when }\,\, a = b \end{array}\right. \end{equation}
で定義される5．
そして、NLsolve の利用を見越して，この計算式 $(\ref{eq:sp-scheme})$ を次のように書き直しておく．
\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. \label{eq:sp-scheme-zero} \end{equation}
これで，$U_n, V_n$ を既知として、この連立非線形方程式を解いて $U_{n+1}, V_{n+1}$ を求めれば時間ステップを一つ進めることが出来るという計算式の出来上がりだ6．
ではプログラムに移ろう．
まずは今回使う差分商を与えておこう．
1function dedx(a,b)
2 EPS = 1.0e-14 # この値はいいかげんに決めた．
3 return ifelse( abs(a-b) > EPS, (exp(a)-exp(b))/(a-b), exp(a) )
4 # ifelse 関数を調べておこう
5end
あとは、上の時間対称なスキームと同様の作業をしていけば良い(途中で $\log$ が出てくることを忘れないようにしよう)． まずはゼロになるべき式 $(\ref{eq:sp-scheme-zero})$ の値を返す関数を次のように定義する．
1function sp_zero(uv_new, uv_curr, E)
2# uv_new は新しい u,v で、uv_curr は古い u,v
3
4 up, vp = log.(uv_new)
5 u, v = log.(uv_curr)
6
7 eq1 = up - u + Δt*C1* dedx(vp,v) + Δt*(-D1 + E)
8 eq2 = vp - v - Δt*C2* dedx(up,u) + Δt*( D2 + E)
9
10 return [ eq1, eq2 ]
11end
そして，この値がほぼゼロになる $(u,v)$ を計算することで 1ステップ分の計算を行う関数を定義する．
1function volterra_sp(uv_curr, E)
2
3 uv_new = nls(sp_zero, uv_curr, E, ini = uv_curr)[1]
4 return uv_new
5
6end
念の為に 1ステップ試しておこう．
1volterra_sp(uv_0, E)
2-element Array{Float64,1}:
0.6935156878695156
0.6724434211370223
動くようなので、これも 500ステップほど動かしてみよう．
1uv = copy(uv_0)
2uv_sq = copy(uv_0')
3G_sq = [ G(uv, E)]
4
5for i in 1:500
6 uv = volterra_sp(uv,E)
7 uv_sq = vcat(uv_sq, uv')
8 push!( G_sq, G(uv, E) )
9end
早速プロットしてみよう．
1plot(t_sq, uv_sq, label = ["魚 🐟" "鮫 🦈"], fontfamily="" )
ふむ、なかなか美しい．では、保存量 $G$ がどうなっているか見てみよう．
1plot(t_sq, G_sq, yaxis = "G", leg = false)
これも Runge-Kutta法の近似解同様，だいぶ良さそうだ． そして値が目盛りの数字から読み取れない状況であることも同様だ．
このようなグラフ描画時に下記のような warning が表示されることがある．
GKS: Possible loss of precision in routine SET_WINDOW
これは、plot 命令にとって「値があまりにも変わらないので、図の上下幅の設定が怪しいかも」という話のようだ． 実際，この図も縦軸の目盛りの数字が全く変化していないしね． こういうときは，全体から適切な基準値を引くなどして，値の幅をはっきりさせてから描画するとよいだろう．
そこで，これも Runge-Kutta法のとき同様，はっきりと見るために $G$ の初期値との差をグラフで見てみよう．
1plot(t_sq, G_sq .- G_sq[1], yaxis = "ΔG", leg = false)
するとすぐわかるように，保存量 $G$ のズレが大変小さく，おおよそ 12桁程度は保たれることがわかる． 長期的に見たらズレは蓄積していくが、それも大変小さいことが見て取れる． つまり，
びっくりするぐらい、保存量が数値的に「良く」保存されている ことがわかる．
4次精度の Runge-Kutta 法に比べ、この解法はせいぜい 2次精度で精度は良くないはずなのに、保存量の保存性に関して言えば遥かに良いことになる．
さて、このからくりはどうなっているのだろう？ 考えておこう(→ レポート問題へ)．
レポート
下記要領でレポートを出してみよう．
- e-mail にて，
- 題名を 2024-numerical-analysis-report-08 として,
- 教官宛(アドレスは web の "TOP" を見よう)に，
- 自分の学籍番号と名前を必ず書き込んで，
- 内容はテキストでも良いし，pdf などの電子ファイルを添付しても良いので，
下記の内容を実行して，結果や解析，感想等をレポートとして提出しよう．
注意
近年はセキュリティ上の懸念から，実行形式のプログラムなどをメールに添付するとそのメールそのものの受信を受信側サーバが拒絶したりする．
そういうことを避けるため，レポートをメールで提出するときは添付ファイルにそういった懸念のあるファイルが無いようにしよう．
まあ要するに，レポートは pdf ファイルにして，それをメールに添付して送るのが良い ということだと思っておこう．
- もとの常微分方程式に従うと保存量 $G$ が確かに保存されることを, 自力で確認しよう．
つまり，$u,v$ が元の Volterra モデルである常微分方程式 の厳密解であるとき、 $ \displaystyle \frac{dG}{dt} $ の時間変化を自分で手で計算し、確かにこれがゼロになることを確認しよう, ということだね． - ここで出てきた 4つの解法(Euler法，時間対称なスキーム, Runge-Kutta法, 構造保存数値解法)が要する計算時間を
BenchmarkToolsパッケージか@timedマクロなどを用いて測定しておこう． - (チャレンジ問題) ここで出てきた 4つの解法(Euler法，時間対称なスキーム, Runge-Kutta法, 構造保存数値解法)の保存量の誤差と計算時間の関係
について，
授業資料 06. 時間対称スキーム, Runge-Kutta法
の最後に出てきたようなグラフを描いてみよう．
- (チャレンジ問題) 最後に紹介した構造保存解法 $(\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})$ を使うことを忘れないようにね．
-
普通、こんな簡単には見つからない．これは超々ラッキーなケース． ↩︎
-
あってもなくても良い $\mbox{const.}$ はゼロにして，実質的に無視してしまおう． ↩︎
-
先に述べたように相図上の解軌道の厳密な式が求まることには大きな価値があるが，実際に求めることは通常は困難なので，今回のこれは大変ラッキーだ． ↩︎
-
英語だと "structure-preserving method" だね．少し古い名前だと "geometric integration" なんていうのもある． ↩︎
-
これをより一般化させた，"離散勾配(discrete gradient)" という概念がある．かなり面白い概念で，拡張性も高いぞ． ↩︎
-
よく見ると，このスキームも「時間対称スキームの一種」である． ↩︎