ただのメモ

他の人に見せても良い方のメモ

差分がロトカボルテラ方程式な場合

ロトカ・ボルテラの方程式とは、以下のような微分方程式である。

 dx/dt = ax -bxy \\ dy/dt = cxy -dy

ここで、yが新たな変数zとxとの差、だとして、考える。

 dx/dt = ax - bx(z-x) \\ d(z-x)/dt = cx(z-x)-d(z-x)

少し式変形をすると、こうなる。

 dx/dt = ax - bx(z-x) \\ dz/dt = ax -bx(z-x)+cx(z-x)-d(z-x)

これが、差分がロトカ・ボルテラの方程式の形、というものだ。

もとが振動する解を持つので、こちらも振動して欲しいものだ。

確認する。

コードはこちら。

#%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import math
import random

#initial settings
N=20000;dt=0.01;T=N*dt
a=0.1;b=0.1;c=0.2;d=0.3;K0=5;K=20
tau=20

def f(s,t):
    return a*s-b*s*(t-s)

def g(s,t):
    return a*s-b*s*(t-s)+c*s*(t-s)-d*(t-s)

def main(p,q,col):
    x=[0]*N;y=[0]*N;z=[0]*N
    x[0]=p;y[0]=q
    for i in range(N-1):
        xk1=f(x[i],y[i])*dt
        yk1=g(x[i],y[i])*dt
        xk2=f(x[i]+xk1/2,y[i]+yk1/2)*dt
        yk2=g(x[i]+xk1/2,y[i]+yk1/2)*dt
        xk3=f(x[i]+xk2/2,y[i]+yk2/2)*dt
        yk3=g(x[i]+xk2/2,y[i]+yk2/2)*dt
        xk4=f(x[i]+xk3,y[i]+yk3)*dt
        yk4=g(x[i]+xk3,y[i]+yk3)*dt
        x[i+1]=x[i]+(xk1+2*xk2+2*xk3+xk4)/6#+0.05*(random.random()-0.5)
        y[i+1]=y[i]+(yk1+2*yk2+2*yk3+yk4)/6
    plt.plot(x,y,color=col)

#we add little noise!
main(3,4,"red")
main(1.5,2.5,"blue")
main(12,14,"green")
main(15,17,"lightgreen")
main(3,6,"purple")
main(8,9,"cyan")
main(2.6,3.1,"orange")
main(4,6,"pink")

plt.xlabel("x")
plt.ylabel("y")
plt.title("Example 403")

plt.show()

f:id:medical-science:20220403214031p:plain

今回は、差分がロトカ・ボルテラ方程式のような場合なので、普通のロトカ・ボルテラの方程式によって描かれる形とは違う。しかし、y座標からx座標との差分を取った値を新たなyの値として、描いてみると、おなじみの形となることだろう。

ちなみに、原点を通りそうな直線が、グラフ上に見えるが、それは y=xである。

そして、その後、xとyの重み付け和が一定の場所の後、xが0近くでyが急激に下がっていく。

これを繰り返す、という現象だ。

 

今回の試みを利用すれば、差分以外にも、和や、積や、商などや、非線形関数に咬ませた時の値にロトカ・ボルテラの方程式が見られる、という形に出来るかもしれない。

 

バイバイ!