ロトカヴォルテラ方程式 １捕食者ー２被食者系 自然災害シュミレーション
解決したいこと
t=200で大規模な自然災害が起きた想定をしたシュミレーションです。
（個体数が初期値の10⁻⁷以下になったら絶滅としています。）
3種すべて絶滅するはずなのですが、下のグラフのように被食者２（x2）がマイナスへ
捕食者１（ｙ１）がプラスへ発散してしまっています。
なぜこのような結果になってしまうのか分かりません。
解決方法を教えていただきたいです。
よろしくお願いいたします。
発生している問題・エラー
dydt = r2*(1- (a21/kx2*a22)* x - (y/kx2) - (a23/kx2*a22) * z) * y
C:\Users\s2180389\.ipython\code\runge1-2 Kkannsu.py:26: RuntimeWarning: overflow encountered in double_scalars
dzdt = (-r3 + a13*e* x+ a23*e* y) * z
C:\Users\s2180389\.ipython\code\runge1-2 Kkannsu.py:40: RuntimeWarning: overflow encountered in double_scalars
dydt = r2*(1 - (y/kx2) - (a23/kx2*a22) * z) * y
C:\Users\s2180389\.ipython\code\runge1-2 Kkannsu.py:41: RuntimeWarning: overflow encountered in double_scalars
dzdt = (-r3 + a23*e* y) * z
C:\Users\s2180389\.ipython\code\runge1-2 Kkannsu.py:55: RuntimeWarning: invalid value encountered in add
newXs = Xs + (k1 + 2 * k2 + 2 * k3 + k4) / 6
該当するソースコード
import numpy as np
import matplotlib.pyplot as plt
# パラメータ
e=2.6
a11=1
a12=1
a13=1
a21=1.5
a22=1
a23=0.2
a31=e*a13
a32=e*a23
r1=1
r2=1
r3=1
kx1=r1/a11
kx2=r2/a22
# 微分方程式の右辺
def func(t, Xs):
x, y, z = Xs
dxdt =r1* (1 - (x/kx1) - (a12/kx1*a11)* y -(a13/kx1*a11)* z) * x
dydt = r2*(1- (a21/kx2*a22)* x - (y/kx2) - (a23/kx2*a22) * z) * y
dzdt = (-r3 + a13*e* x+ a23*e* y) * z
if y<=0.0000003:
dxdt=r1* (1 - (x/kx1) -(a13/kx1*a11)* z) * x
dydt =0
dzdt = (-r3 + a13*e* x) * z
if z<=0.0000003:
dxdt =r1* (1 - (x/kx1) - (a12/kx1*a11)* y ) * x
dydt = r2*(1- (a21/kx2*a22)* x - (y/kx2) ) * y
dzdt = 0
if x<=0.0000003:
dxdt =0
dydt = r2*(1 - (y/kx2) - (a23/kx2*a22) * z) * y
dzdt = (-r3 + a23*e* y) * z
return np.array([dxdt,dydt,dzdt])
# 4次ルンゲクッタ
def rk4(dt, t, Xs, func):
k1 = dt * func(t, Xs)
k2 = dt * func(t + dt / 2, Xs + k1 / 2)
k3 = dt * func(t + dt / 2, Xs + k2 / 2)
k4 = dt * func(t + dt, Xs + k3)
new_t = t + dt
newXs = Xs + (k1 + 2 * k2 + 2 * k3 + k4) / 6
return new_t, newXs
# 計算の実行
dt = 0.05
t0 = 0
te = 1000
ts = np.arange(t0, te + dt, dt)
Xs = np.array([0.3, 0.3, 0.3])
xs, ys, zs = [], [], []
for t in ts
x, y, z = Xs
xs.append(x)
ys.append(y)
zs.append(z)
if t>=200:
ddt=10
dk=0.98
kx1=(r1/a11)-dk/(1+((t-200)**2/ddt**2))
kx2=(r2/a22)-dk/(1+((t-200)**2/ddt**2))
_, Xs = rk4(dt, t, Xs, func)
# 描画
fig, ax = plt.subplots()
ax.plot(ts, xs, label="x1")
ax.plot(ts, ys, label="x2")
ax.plot(ts, zs, label="y1")
ax.set_xlim(0,1000)
ax.set_ylim(0,1.0)
ax.set_xlabel("t")
ax.set_ylabel("x1, x2, y")
ax.legend()
plt.show()
災害の規模が小さい場合このような結果が得られます。
0
