Модель Вольтерры-Лотки

[Alm99]

Добрый день, уважаемые форумчане! Возник один вопрос, очень нужна консультация. Имеется всем известная модель Лотки-Вольтерры "хищник-жертва". Я написал скрипты на языке Python 3.0 для моделирования, и выбираю в качестве решателя системы ОДУ метод Рунге-Кутты 4 порядка. А вопрос заключается вот в чем: почему при решении методом Рунге-Кутты я на фазовой плоскости получаю такую картинку:

и можно ли получить такую картинку ( я просто не понгимаю, почему при использовании другого метода картина иная)может я что-то делаю не так:

Код программы выкладываю ниже:

сам решатель метод Рунге-Кутта 4 порядка

Используется синтаксис Python

def rk4(funceq, v, p, dt):
    """ Реализация метода Рунге-Кутта 4 порядка дял решения системы ДУ модели Вольтеррры - Лотки
    Функция получает массив переменных системы ОДУ funceq, а также временной шаг dt и переменные
    u,v."""

    vder1, pder1 = funceq(v, p)
    vder2, pder2 = funceq(v + vder1 * dt / 2, p + pder1 * dt / 2)
    vder3, pder3 = funceq(v + vder2 * dt / 2, p + pder2 * dt / 2)
    vder4, pder4 = funceq(v + vder3 * dt, p + pder3 * dt)
    vnew = v + (vder1 + 2 * vder2 + 2 * vder3 + vder4) * dt / 6
    pnew = p + (pder1 + 2 * pder2 + 2 * pder3 + pder4) * dt / 6
    return vnew, pnew
 

а это основной файл:

Используется синтаксис Python

import pylab as pp
import numpy as np
import RungeKutta4
import csv
import pandas as pd

# Параметры модели Волтерра-Лотки
a = 0.8
b = 0.2
c = 0.33
d = 0.08
# Шаг дискретизации по времени
dt = 0.05

# Расчет времени моделирования ( временной интервал для модели динамики популяций)
N = 500
t = np.arange(N) * dt

# Начальные значения популяций
v0 = 20
p0 = 7

# Уравнения модели
def funceq(v, p):
    vder = a * v - b * v * p
    pder = d * v * p - c * p
    return vder, pder

# Симуляция процесса (цикл повторений)
v = np.zeros(N)
p = np.zeros(N)
v[0] = v0
p[0] = p0
for k in range(N - 1):
    v[k + 1], p[k + 1] = RungeKutta4.rk4(funceq, v[k], p[k], dt)

# Вывод графиков
f1 = pp.figure()
pp.plot(t, v,'g--', label = 'Жертвы', linewidth = 1)
pp.plot(t, p,'b-.', label = 'Хищники', linewidth = 1)
pp.grid()
pp.legend(loc = 'best')
pp.xlabel('Время')

f2 = pp.figure()
pp.plot(v, p, 'r-->', linewidth = 0.3)
pp.grid()
pp.xlabel('Жертвы')
pp.ylabel('Хищники')
pp.title('Фазовый портрет')

pp.show()
 

Я изучил коды, написанные другими людьми, при использовании метода Рунге-Кутта 4 порядка у всех фазовый портрет имеет такой же вид как и на первой картинке. Прошу помощи, не могу разобраться, почему так?

[Vladimir-80]

А откуда вы взяли вторую картинку? Она какая-то странная.

[Alm99]

Vladimir-80
Вот здесь пример моделирвания:
https://scientific-python.readthedocs.io/en/latest/notebooks_rst/3_Ordinary_Differential_Equations/02_Examples/Lotka_Volterra_model.html

[Pphantom]

Вообще говоря, в модели Лотки-Вольтерры должен получаться цикл, и первая картинка этому вполне соответствует. Так что, по-видимому, проблема в "другом методе", из-за которого вторая картинка получается попросту неправильной.

P.S. Ну да - используется существенно менее точный метод, отсюда и результат.

[Alm99]

Pphantom
Вот да, колебательный процесс, который повтоярется, просто я думаю, решил задать вопрос, неужели так сильно влияет выбранный численный метод? Еще конечно есть вопрос, в данном случае как определяется точка равновесия (стационарное положение) - я знаю общую формулу, но не знаю как определить собсвтенные значения матрицы и коэффициентов системы ОДУ, как это организовать командой в python?

Используется синтаксис Python

lambda1, lambda2 = linalg.eigvals()

Как правильно передать эти значения, чтобы получить чисел $\lambda_{1,2}$

[Pphantom]

Alm99 в сообщении #1532329 писал(а):

Вот да, колебательный процесс, который повтоярется, просто я думаю, решил задать вопрос, неужели так сильно влияет выбранный численный метод?

Естественно. Иначе зачем вообще нужно было бы возиться с более сложными методами?..

Alm99 в сообщении #1532329 писал(а):

Еще конечно есть вопрос, в данном случае как определяется точка равновесия (стационарное положение) - я знаю общую формулу, но не знаю как определить собсвтенные значения матрицы и коэффициентов системы ОДУ, как это организовать командой в python?

А зачем? Точка равновесия - это ситуация, когда обе производные равны нулю. Запишите соответствующую систему и решите ее руками, это существенно проще.

[Alm99]

Pphantom
Руками то я ее решил) Ну ладно, по сути это и не требуется! Спасибо за оказанную помощь!

[dtn888]

Любопытная система, нелинейная и одновременно с очень простым решением для двух точек равновесия:
1. $x(t)=0,y(t)=0$

2. $x(t)=\frac{\delta}{\gamma},y(t)=\frac{\alpha}{\betta}$

Вот численное решение с помощью Mathematica:

Clear["Derivative"]; ClearAll["Global`*"];

pars = {\[Alpha] = 1, \[Beta] = 1, \[Gamma] = 1, \[Delta] = 1};

s[x0_?NumericQ] := 
  First@NDSolve[{x'[t] == \[Alpha] x[t] - \[Beta] x[t] y[t], 
     y'[t] == \[Gamma] x[t] y[t] - \[Delta] y[t], x[0] == x0, 
     y[0] == x0}, {x, y}, {t, 0, 30}, 
    Method -> {"TimeIntegration" -> {"ExplicitRungeKutta", 
        "DifferenceOrder" -> 4}}];

ParametricPlot[
 Evaluate[{x[t], y[t]} /. s[#] & /@ Range[0.1, 1, 0.05]], {t, 0, 30}, 
 PlotRange -> Full]

Интересно следующее:
1. Эта модель описывает реальные процессы в природе, или это просто "игрушка"?
2. Есть ли такие значения коэффициентов, при которых система находится в одной из точек равновесия?

[Alm99]

dtn888
Это стандартная модель Лотки - Вольтерра, описывающая процесс динамики популяции, хотя есть ее и более точная постановка...