@lotkavposted at 2021-11-15

ロトカヴォルテラ方程式競争系３すくみの関係

Q&A

Python 数値計算シミュレーション数式

解決したいこと

ロトカヴォルテラ方程式競争系
３すくみの関係についてです。

上のようにプロットしたところ３種すべて個体数が０となってしまいます。

画像の左図のようにしたいのですが、方法がわかりません。
解決方法を教えていただきたいです。
よろしくお願いいたします。

該当するソースコード

import numpy as np
import matplotlib.pyplot as plt

a,b=1.6,0.5

x１=0.3
x2=0.3
x3=0.3

dt=0.1
n=10000

x=np.zeros(n)
x[0]=x1
y=np.zeros(n)
y[0]=x2
z=np.zeros(n)
z[0]=x3

for i in range(1,n):
    x[i]=(1-x[i]-a*y[i]-b*z[i])*x[i]*dt
    y[i]=(1-b*x[i]-y[i]-a*z[i])*y[i]*dt
    z[i]=(1-a*x[i]-b*y[i]-z[i])*z[i]*dt
    print(i,x[i],y[i],z[i])

t=np.arange(0,n*dt,dt)

plt.plot(t,x,label="x1")
plt.plot(t,y,label="x2")
plt.plot(t,z,label="x3")
plt.legend()
plt.show()

1Answer

@bizzpaper posted at 2021-11-15

微分方程式の差分解法は既知の値を使用して次の時間の値を計算します。

@loktav さんがfor文で計算している内容では、求めたい次の時間の値を使用して次の時間の値を計算しようとしています。

なぜ0になるかというと、上で述べたことに加えてx = np.zeros(n)と定義しているので、計算式から次の時間の値は0のままです。

また、for文内の計算式は合っていますでしょうか？
計算式が微分方程式の右辺のみで、差分で展開した左辺の微分項のうち既知の値を移項し忘れているように見えます。

以下に実装例とその結果を載せますので参考にしてください。

ご希望のようなグラフになるかは、パラメータ次第だと思います。

import numpy as np
import matplotlib.pyplot as plt

dt = 0.05
t0 = 0
te = 500

ts = np.arange(t0, te + dt, dt)

a = 2.9
b = 0.8
x = 0.3
y = 0.3
z = 0.3

xs, ys, zs = [], [], []

for t in ts:
	xs.append(x)
	ys.append(y)
	zs.append(z)

	x += (1. - xs[-1] - a * ys[-1] - b * zs[-1]) * xs[-1] * dt
	y += (1. - b * xs[-1] - ys[-1] - a * zs[-1]) * ys[-1] * dt
	z += (1. - a * xs[-1] - b * ys[-1] - zs[-1]) * zs[-1] * dt

	print(f"(x, y, z) = ({x:.3e}, {y:.3e}, {z:.3e}) ({t:.2e} sec)\r", end=" ")


fig, ax = plt.subplots()

ax.plot(ts, xs, label="x")
ax.plot(ts, ys, label="y")
ax.plot(ts, zs, label="z")

ax.set_xlabel("t")
ax.set_ylabel("x, y, z")
ax.legend()

plt.show()

0Like

Comments

@lotkav
Questioner
＠bizzpaperさん　ご回答ありがとうございます。

初歩的なところからご教授いただき助かりました。
０になる理由、計算式の改善点を理解することができました。

実装例を参考にさせていただきパラメータをいじって検討します。
ありがとうございました。
@bizzpaper
お力になれたようで何よりです！
引き続きがんばってください😎👍

Are you sure you want to delete the question?