ただのメモ

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

Lotka Volterraの方程式で遊ぶ

ロトカ・ボルテラの方程式で遊んでみる。

まずは、普通の捕食・被捕食系。 dx/dt = ax-bxy, dy/dt = cxy-dy

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

#initial settings
N=30000;dt=0.01;T=N*dt
a=0.1;b=0.1;c=0.2;d=0.3


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


def main(p,q,c):
    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
        y[i+1]=y[i]+(yk1+2*yk2+2*yk3+yk4)/6

    plt.plot(x,y,color=c)

main(1,1,"blue")
main(2,1,"red")
main(1,0.5,"green")
main(0.5,0.5,"cyan")
main(1.5,1.1,"orange")
main(1.5,1.01,"purple")

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

plt.show()

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

次は、捕食者同士が共食いをするケース。 dx/dt = ax - bxy, dy/dt = cxy-dy^2

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

最終項を2乗から3乗に変えるとこうなる。

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

次は、最後の項が、 dxy/(1+x)である時。つまり、xが多い時は、普通に減り、xが少ない時は減りにくい場合。xが毒性のあるものだとこうなりそう。

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

次は、最終項が、 dy/(1+x)の時。つまり、xが多い時は、全然減らない(栄養過多だから死ににくい)が、xが少ない時は、ちゃんと死ぬ場合。

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

さっきの項を変えて dy/(1+x^2)にした場合。

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

3乗にするとこうなる。

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

最後は、最終項を dy/(1+e^s)にしてみた。

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

 

いろいろと弄ってみて面白い。

ここで湧き出た疑問は、果たして、この軌道の形が意味するのは、どういうことか、ということだ。

数理モデルにセンスが必要なら、そのセンスがある人とセンスが無い人の違いは何か。

そのセンスの違いを生み出すものは、実は数式と形の関係性を記述する上で非常に重要なものではないだろうか。

数理モデルは奥が深い。

バイバイ!