-1

我正在尝试使用 Lotka-Volterra 方程编写一个程序来进行捕食者与猎物的相互作用。使用 ODE 求解：

dx/dt = a*x - B*x*y
dy/dt = g*x*y - s*y

使用四阶龙格-库塔法

我需要绘制一个图表，显示 x 和 y 作为从 t = 0 到 t=30 的时间的函数。

a = alpha = 1
b = beta = 0.5
g = gamma = 0.5
s = sigma = 2
初始条件 x = y = 2

到目前为止，这是我的代码，但未在图表上显示任何内容。一些帮助会很好。

#!/usr/bin/env python

from __future__ import division, print_function
import matplotlib.pyplot as plt
import numpy as np

def rk4(f, r, t, h):
        """ Runge-Kutta 4 method """
        k1 = h*f(r, t)
        k2 = h*f(r+0.5*k1, t+0.5*h)
        k3 = h*f(r+0.5*k2, t+0.5*h)
        k4 = h*f(r+k3, t+h)
        return (k1 + 2*k2 + 2*k3 + k4)/6

def f(r, t):
        alpha = 1.0
        beta = 0.5
        gamma = 0.5
        sigma = 2.0
        x, y = r[2], r[2]
        fxd = x*(alpha - beta*y)
        fyd = -y*(gamma - sigma*x)
        return np.array([fxd, fyd], float)


tpoints = np.linspace(0, 30, 0.1)
xpoints = []
ypoints = []

r = np.array([2, 2], float)
for t in tpoints:
        xpoints += [r[2]]
        ypoints += [r[2]]
        r += rk4(f, r, t, h)

plt.plot(tpoints, xpoints)
plt.plot(tpoints, ypoints)
plt.xlabel("Time")
plt.ylabel("Population")
plt.title("Lotka-Volterra Model")
plt.savefig("Lotka_Volterra.png")
plt.show()
4

2 回答 2

4

tpoints运行脚本后对变量的简单检查表明它是空的：

In [7]: run test.py

In [8]: tpoints
Out[8]: array([], dtype=float64)

这是因为你使用np.linspace不正确。第三个参数是输出中所需的元素数量。您请求了一个长度为 0.1 的数组。

看看np.linspace's 的文档字符串。弄清楚如何调整代码不会有问题。

于 2015-04-21T02:50:55.770 回答
1

1）定义'h'变量。

2) 使用

tpoints = np.arange(30) #array([0, 1, 2, ..., 30])

不是

np.linspace()

并且不要忘记将时间步长设置为等于 h：

h=0.1
tpoints = np.arange(0, 30, h)

3）小心索引：

def f(r,t):
    ...
    x, y=r[0], r[1]
    ...

for t in tpoints:
    xpoints += [r[0]]
    ypoints += [r[1]]
    ...

并更好地使用 .append(x)：

for t in tpoints:
    xpoints.append(r[0])
    ypoints.append(r[1])
    ...

这是python 3.7的测试代码（我设置了h = 0.001以获得更多的presize）

import matplotlib.pyplot as plt
import numpy as np

def rk4(r, t, h):                    #edited; no need for input f
        """ Runge-Kutta 4 method """
        k1 = h*f(r, t)
        k2 = h*f(r+0.5*k1, t+0.5*h)
        k3 = h*f(r+0.5*k2, t+0.5*h)
        k4 = h*f(r+k3, t+h)
        return (k1 + 2*k2 + 2*k3 + k4)/6

def f(r, t):
        alpha = 1.0
        beta = 0.5
        gamma = 0.5
        sigma = 2.0
        x, y = r[0], r[1]
        fxd = x*(alpha - beta*y)
        fyd = -y*(gamma - sigma*x)
        return np.array([fxd, fyd], float)

h=0.001                               #edited
tpoints = np.arange(0, 30, h)         #edited
xpoints, ypoints  = [], []
r = np.array([2, 2], float)
for t in tpoints:
        xpoints.append(r[0])          #edited
        ypoints.append(r[1])          #edited
        r += rk4(r, t, h)             #edited; no need for input f
plt.plot(tpoints, xpoints)
plt.plot(tpoints, ypoints)
plt.xlabel("Time")
plt.ylabel("Population")
plt.title("Lotka-Volterra Model")
plt.savefig("Lotka_Volterra.png")
plt.show()

您还可以尝试绘制“循环”：

plt.xlabel("Prey")
plt.ylabel("Predator")
plt.plot(xpoints, ypoints)
plt.show()

https://i.stack.imgur.com/NB9lc.png

于 2019-05-28T01:07:09.850 回答