python - Lotka Volterra 和 Runge Kutta 不需要精确度

Question

人口随着时间的推移(每个高峰的高度应该相同

我编写了一个代码，使用 Runge-Kutta 4 阶来模拟小鼠和狐狸种群。但结果并不像想要的那样.. 每个峰应该几乎在相同的高度我不认为这是步长的问题..你有想法吗？

import matplotlib.pyplot as plt
import numpy as np

#function definition
def mice(f_0, m_0):
    km = 2.      #birthrate mice
    kmf = 0.02  #deathrate mice
    return km*m_0 - kmf*m_0*f_0

def fox(m_0, f_0):
    kf = 1.06   #deathrate foxes
    kfm = 0.01  #birthrate foxes
    return kfm*m_0*f_0 -kf*f_0

def Runge_kutta4( h, f, xn, yn):
    k1 = h*f(xn, yn)
    k2 = h*f(xn+h/2, yn + k1/2)
    k3 = h*f(xn+h/2, yn + k2/2)
    k4 = h*f(xn+h, yn + k3)
    return yn + k1/6 + k2/3 + k3/3 + k4/6

h = 0.01
f = 15.
m = 100.
f_list = [f]
m_list = [m]

for i in range(10000):
    fox_new = Runge_kutta4(h, fox, m, f)
    mice_new = Runge_kutta4(h, mice, f, m)
    f_list.append(fox_new)
    m_list.append(mice_new)
    f = fox_new
    m = mice_new

time = np.linspace(0,100,10001)
#Faceplot LV
fig = plt.figure(figsize=(10,10))
fig.suptitle("Runge Kutta 4")
plt.grid()
plt.xlabel('Mice', fontsize = 10)
plt.ylabel('Foxes', fontsize = 10)
plt.plot(m_list, f_list, '-')
plt.axis('equal')
plt.show()
fig.savefig("Faceplott_Runge_Kutta4.jpg", dpi=fig.dpi)

fig1 = plt.figure(figsize=(12,10))
fig1.suptitle("Runge Kutta 4")
plt.grid()
plt.xlabel('Time [d]', fontsize = 10)
plt.ylabel('Populationsize', fontsize = 10)
plt.plot(time, m_list , label='mice')
plt.plot(time, f_list , label='fox')
plt.legend()
plt.show()
fig1.savefig("Fox_Miceplot_Runge_Kutta4.jpg", dpi=fig.dpi)

Answer 1

在 Runge-Kutta 实现中，xn 是时间变量，yn 是标量状态变量。 f 是标量 ODE y'(x)=f(x,y(x)) 的标量 ODE 函数。但是，这不是您应用 RK4 过程的方式，您的 ODE 函数是自主的，不包含时间变量，而是包含两个耦合状态变量。实现时，结果应该是一个没有特定类型的复杂的一阶方法。

---

您需要将耦合系统作为一个耦合系统来求解，也就是说，必须以相同的增量同时计算两个变量的阶段。

kf1 = h*fox(mn, fn)
km1 = h*mice(fn, mn)
kf2 = h*fox(mn+0.5*km1, fn+0.5*kf1)
km2 = h*mice(fn+0.5*kf1, mn+0.5*km1)
kf3 = h*fox(mn+0.5*km2, fn+0.5*kf2)
km3 = h*mice(fn+0.5*kf2, mn+0.5*km2)
kf4 = h*fox(mn+km3, fn+kf3)
km4 = h*mice(fn+kf3, mn+km3)

等等

---

另见 Runge Kutta problems in JS对于 JavaScript 中的相同问题

另一种方法是对系统进行矢量化，使得Runge-Kutta过程可以保持不变，但在积分循环中必须构造和解包状态向量，

def VL(x,y): f, m = y; return np.array([fox(m,f), mice(f,m)])

y = np.array([f,m])
time = np.arange(x0,xf+0.1*h,h)

for x in time[1:]:
    y = Runge_kutta4(h, VL, x, y)
    f, m = y
    f_list.append(f)
    m_list.append(m)

其他一切保持不变。