Lotka Volterra和Runge Kutta并不需要精确 [英] Lotka Volterra with Runge Kutta not desired precision

问题描述

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

我已经编写了一个代码，以使用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)

解决方案

在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)

等

---

有关JavaScript中的相同问题，另请参见 JS中的Runge Kutta问题

另一种方法是对系统进行矢量化处理，以使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)

 

其他所有内容保持不变.

Population over time (should be the same height at every peak

I've programmed a code to simulate a mice and fox population using Runge-Kutta 4th order. But the result is not as wanted to be.. each peak should nearly be at same height I don't think that it is a problem of step size.. Do you have an idea?

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)

解决方案

In the Runge-Kutta implementation, xn is the time variable and yn the scalar state variable. f is the scalar ODE function for the scalar ODE y'(x)=f(x,y(x)). However, this is not how you apply the RK4 procedure, your ODE functions are autonomous, contain no time variable but instead of it two coupled state variables. As implemented, the result should be a convoluted first order method of no specific type.

---

You need to solve the coupled system as a coupled system, that is, the stages for both variables have to be calculated simultaneously with the same increments.

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)

etc.

---

See also Runge Kutta problems in JS for the same problem in JavaScript

The other way is to vectorize the system so that the Runge-Kutta procedure can remain the same, but in the integration loop the state vector has to be constructed and unpacked,

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)

everything else remaining the same.

这篇关于Lotka Volterra和Runge Kutta并不需要精确的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！