Лотка Вольтерра с Рунге Куттой не нужной точности - PullRequest
Новая
       45

Лотка Вольтерра с Рунге Куттой не нужной точности

2 голосов
/ 06 октября 2019

Популяция с течением времени (должна быть одинаковой высоты на каждом пике

Я запрограммировал код для имитации популяции мышей и лис, используя 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)

1 Ответ

2 голосов
/ 06 октября 2019

В реализации Рунге-Кутты xn - это переменная времени, а yn - скалярная переменная состояния. f - скалярная функция ODE для скалярного ODE y'(x)=f(x,y(x)). Однако это не то, как вы применяете процедуру 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)

и т. д.

См. также Задачи Рунге-Кутты в JS для той же проблемы в JavaScript

enter image description here

Другой способ - это векторизация системы, чтобы процедура Рунге-Кутта могла оставаться прежней, но в цикле интегрирования вектор состояния должен быть построен и распакован. , 

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)

все остальное остается прежним.

...