2020/10/04 13:53 更新
巡回型Lotka-Volterra方程式の陽的Euler法による計算
83 いいね ブックマーク
目次

生物種の個体数の時間変化を記述する巡回型Lotka-Volterra方程式(式 1)を、陽的Euler法で離散化します(式 2)。

$$\tag{1} \frac{\mathrm{d}}{\mathrm{d} t}\left(\begin{array}{c}p(t) \\ q(t) \\ r(t)\end{array}\right)=\left(\begin{array}{c}p(t)(q(t)-r(t)) \\ q(t)(r(t)-p(t)) \\ r(t)(p(t)-q(t))\end{array}\right)$$$$\tag{2}\frac{1}{h}\left(\begin{array}{c}p_{n+1}-p_{n} \\ q_{n+1}-q_{n} \\ r_{n+1}-r_{n}\end{array}\right)=\left(\begin{array}{c}p_{n}\left(q_{n}-r_{n}\right) \\ q_{n}\left(r_{n}-p_{n}\right) \\ r_{n}\left(p_{n}-q_{n}\right)\end{array}\right)$$

離散化する前の巡回型Lotka-Volterra方程式では、$pqr=$一定であるため、この微分方程式の解を$pq$平面に射影すると、$pqr=$一定の曲面を$p+q+r=1$で切断して射影した軌跡と一致します。
一方、離散化された巡回型Lotka-Volterra方程式では、$p_nq_nr_n$は狭義単調減少であるため、解を$pq$平面に射影した軌跡は時間発展に伴い徐々に大きくなり、二等辺三角形に漸近します。
実際に、juliaで計算して、確かめてみました。

h = 0.5
num = 200
function explicit_euler(p,q,r)
    p_new = p+h*p*(q-r)
    q_new = q+h*q*(r-p)
    r_new = r+h*r*(p-q)
    return p_new, q_new, r_new
end
function LV(p,q,r)
    pqr = [p, q, r]
    for i in 1:num
        p,q,r = explicit_euler(p,q,r)
        pqr = hcat(pqr, [p,q,r])
    end    
    return pqr
end

実行結果は、こちらです。コンピュータの計算結果からも、二等辺三角形に漸近している様子が確認できました。
このように、単純に陽的Euler法を適用すると、保存するべき量が保存しないという、弊害が発生してしまう事があります。そのため数値計算には、さまざまな差分法がありますが、保存量が保存するかどうかというのは、一つの目安になります。