КОМПЬЮТЕРНЫЙ ПРАКТИКУМ ПО КУРСУ¶

Системный анализ и основы компьютерного моделирования экосистем¶

Cоставил: к.э.н. Варюхин А.М.¶

В пособии приведены методические указания и порядок выполнения работ, обеспечивающих компьютерную поддержку авторского курса А.М. Варюхина "Системный анализ и основы моделирования экосистем". В каждой работе сформулированы цели ее выполнения, кратко изложены теоретические основы, на базе которых она построена, приведены контрольные вопросы и рекомендуемая литература.Пособие будет полезно студентам и аспирантам узучающим методы компьютерного моделирования динамических экологических, социально-экомических и технических систем.

Саратов 2025¶

МЕТОДИЧЕСКИЕ УКАЗАНИЯ И ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТ¶

"Исследование моделей роста популяций"¶

Основная цель работ: закрепить теоретические сведения по математической теории динамики популяций и получить, практические навыки по построению и применению математических моделей роста популяций при различных предположениях об их функционировании и с использованием программных средств Julia , Engee и Markdown.

Введение

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

Методические указания при использовании программных средств.

После ознакомления с теорией предлагается изучить несколько типов математических моделей. Логика работы с каждой из моделей следующая, пользователь строит модель, задает соответствующие параметры модели и получает в графическом представлении различные типы поведения системы. Наборы параметров, обеспечивающие характерные режимы динамики численности приведены в заданиях, приложенных к каждой из моделей. После выполнения очередного задания пользователь отвечает на контрольные вопросы. Модельные параметры, которые должны задаваться пользователем, следующие: начальная численность популяции, скорость роста популяции, емкость среды, коэффициент смертности, коэффициент внутривидовой конкуренции и.т.д..Работы основывается на последовательном выполнении заданий, в каждом из которых изучается влияние одного или группы параметров на режим динамики численности популяции.

В качестве программного инструментария предлогается использовать язык Julia и отечественную среду моделирования Engee. Язык программирования Julia разрабатывается с 2009 года в Массачусетском технологическом институте (MIT). Распространяется бесплатно по лицензии MIT. Российская среда моделирования Engee в которую интегрирован язык Julia является прекрасной альтернативой системе MATLAB и ,прежде всего, это связано с тем в ней реализована технология визуального моделирования, которая не уступает ни по простоте использования ни по функциональности хорошо известным Simulink Matlab и Wolfram SystemModeler.

Язык программирования Julia довольно прост. Его наиболее сильная сторона- возможность неограниченного расширения с помощью пакетов. В Julia реализованы практически все актуальные средства универсальных статистических вычислений, методы оптимизации и решения дифференциальных уравнений, а также множество алгоритмов машиного обучения, включая построение и обучение нейронных сетей. Кроме того, он является кросс-платформенным, компилируемым языком, что позволяет создавать программы, быстродействие которых сопоставимо с быстродействием программ, написанных на C, Fortran. Ещё одна особенность языка — возможность создания качественной графики типографского уровня, которая может быть экспортирована в распространённые графические форматы и использована для презентаций или публикаций. Использование Julia в среде Engee открывает дополнительные возможности, связанные с технологией визуального моделирования. Все это обусловило выбор языка Julia как базового программного инструментария для изучения математических моделей развития и взаимодействия популяций.

Для выполнения работ небходимо установить на свой компьютер актуальную версию языка Julia и Microsoft VS Code .Для реализации технологии визуального моделирования необходимо зарегистрироваться в Engee https://engee.com/account/login Изучить основы программной среды Microsoft VS Code, Engee , языка программирования Julia и языка разметки теста Markdown, который должен быть использован при создании отчетов о проделанной работе.

Лабораторная работа 1 "Модель экспоненциального роста".¶

В основе этой модели, предложенной Мальтусом в 1798 г., лежит предположение, что прирост численности вида за время t пропорционален этой численности и интервалу времени, за который произошел прирост:

$$\Delta X = r X \Delta t $$

Здесь r- константа собственной скорости роста популяции. Совершив предельный переход, получим линейное дифференциальное уравнение:

$$\frac{{dX}}{{dt}} = r \cdot X $$

Примером применения модели Мальтуса может служить описание развития однородной популяции в условиях неограниченных ресурсов питания (рост клеточной культуры до начала истощения питательной среды)

Порядок выполнения работы¶

В данном задании исследуют влияние начальной численности и параметра r на динамику роста популяции. Цель задания - изучить влияние начальной численности и скорости роста на характер кривой динамики численности и освоить общие принципы работы с программными средствами Julia и Markdown

ЭКСПОНЕНЦИАЛЬНАЯ МОДЕЛЬ¶

скрипт 1

using DifferentialEquations, Plots, PlotThemes
f(x, r, t) = 0.2 * x
x0 = 4
tspan = (0.0, 20.0)
prob = ODEProblem(f, x0, tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
theme(:dao)
plot(sol, linewidth=3, title="График решения по модели Мальтуса",
    xaxis="Time (t)", yaxis="x(t)",label="x(t)-численность популяции")

1.Скопируйте к себе на компьютер, выше привенный скрипт.

Запуститу его в VSC . Задайте r – скоростЬ роста популяции (r < 0,5) и начальное значение численности $x_0$ $(x_0 < 5)$. Задайте интервал моделирования 20. Запустите скрипт на исполнение. Посмотрите на результаты. Поэкспериментируйте с моделью, меняйте параметры, и посмотрите, как меняются графики. Предложенный выше скрипт, позволяют анализировать отклик модели на задаваемые вручную параметры $r$ и $x_0$.

Необходимо разработать скрипт, позволяющий автоматически задавать параметры и получать оклики модели. Вариант скрипта позволяющего задавать набор из нескольких значений $r$ и $x_0$ приведен ниже.

скрипт 2¶

using DifferentialEquations, Plots, PlotThemes

f(x, r, t) = r * x

r=0.5
x0=200
tspan = (0.0, 1.5)

prob = ODEProblem(f, x0, tspan, r)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
theme(:dao)
plot(sol, linewidth=3, label="x0= 200.0 r= 0.5")


#План 'эксперимента по x0
x0exp = range(100,1200,length=3)
#План эксперимента по r
rexp = (0.1:0.2:0.5)

#План 'эксперимента по x0
for i in (1:3)
    x0 = x0exp[i]
    #План эксперимента по r 
    for j in (1:3)
        r = rexp[j]
        x0str=string(x0)
        rstr=string(r)
        labstr = "x0= "*x0str * " r= "*rstr
        println(labstr)
        prob = ODEProblem(f, x0, tspan, r)
        sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
        plot!(sol, linewidth=3, label=labstr)
    end
end
xlabel!("Time (t)")
ylabel!("x(t)")
title!("Семейство графиков при изучение модели Мальтуса")

x0= 100.0 r= 0.1
x0= 100.0 r= 0.3
x0= 100.0 r= 0.5
x0= 650.0 r= 0.1
x0= 650.0 r= 0.3
x0= 650.0 r= 0.5
x0= 1200.0 r= 0.1
x0= 1200.0 r= 0.3
x0= 1200.0 r= 0.5

2.Скопируйте к себе на компьютер, данный скрипт.

Запустите его в $VSC$ . Используйте этот скрипт для исследования поведения модели Мальтуса. Поэкспериментируйте с программой, наблюдая отклики модели при различных планах изменения параметров $r$ и $x_0$.

3.Создайте $maltus.ipynb$ документ, который содержит:

Наименование доумента и фамилию студента;
Постановку задачи (историческая справка, математическая модель Мальтуса и т. д.);
Всроенные скрипты и результаты их работы.
Интерпретацию полученных результатов, ответы на контрольные вопросы и выводы.

Контрольные вопросы:

чем отличаются полученные кинетические кривые?
как получить эти кривые в логарифмическом масштабе?

4.Сформируйте на основе созданного Вами $maltus.ipynb$ файла $HTML$ файл

Лабораторная работа.2 " Непрерывная модель логистического роста".¶

Модель логистического роста была предложена Ферхюльстом в 1838 г. для описания развития популяции в условиях ограниченных ресурсов питания. В основу модели положено уравнение:

$dX/dt = r\cdot X -b\cdot X^{2}$ , которое приводится к виду: $dX/dt = r\cdot X\cdot(1-X/K)$

Член - $b\cdot X^{2}$, пропорциональный количеству встреч между особями учитывает "самоотравление" популяции, объяснимое многими причинами (конкуренция за ресурсы питания, выделение в среду вредного метаболита и др.). Коэффициент $b$ называется коэффициентом внутривидовой конкуренции. Величина $К = г/b$ соответствует устойчивому стационарному состоянию с максимально возможной в данных условиях численностью популяции и называется "емкостью среды". Параметр $r$ называется скоростью роста и характеризует способность популяции к увеличению численности. Решением дифференциального уравнения является функция:

$$X(t)=X_0/(X_0+(K-X_0)\cdot exp(-r\cdot t))$$

$X_0$ - начальная численность популяции. Характер логистической кривой зависит oт величин параметров $r$ и $К$ и начальной численности $X_0$

ЛОГИСТИЧЕСКАЯ МОДЕЛЬ¶

скрипт 3

using DifferentialEquations, Plots, PlotThemes
f!(u, p, t) = r * u - u^2 * b
#-----------ПАРАМЕТРЫ МОДЕЛИ--------------------
r=0.1
k=1000.0
b=r/k
tspan = (0.0, 80.0)
u0 = 200

prob = ODEProblem(f!, u0, tspan)
sol = solve(prob)
theme(:dao)
plot(sol, linewidth=3, label="x0= 200.0 r= 0.1 k=1000")
xlabel!("Time (t)")
ylabel!("x(t)")
title!("Логистическая модель Ферхюльста")

Порядок выполнения работы¶

Цель задания - изучить влияние параметров: скорости роста, начальной численности, емкости среды на режим динамики численности популяции.

1.Скопируйте к себе на компьютер, выше привенный скрипт.Задайте значения $r$ – скорости роста популяции ( $0 < r < 0,5$), начальное значение численности $X_0$ ($10 < X_ < 100$) и емкость среды $K$ ($100 < K < 10000$). Задайте интервал моделирования в диапазоне (20-200). Запустите модель на исполнение в $VSC$. Поэкспериментируйте с моделью, оцените полученные результаты.

2.Самостоятельно разработайте свой скрипт, который позволяет получать два семейства графиков.Семейство графиков зависимости $Xср(r)$ при разных значениях $K$ и семейство графиков $X(t)$ при разных значениях $r$ и заданном значении $K$. Поэкспериментируйте со скриптом, наблюдая отклики модели при различных планах изменения параметров.

скрипт 4

using DifferentialEquations, Plots, PlotThemes
using StatsBase
f!(u, p, t) = r * u - u^2 * b
#-ПАРАМЕТРЫ МОДЕЛИ
#План эксперимента по r
rexp = (0.1:0.1:1.2)
n_rexp = length(r)
bexp = rexp / k

k = 100.0
r = 0.1
b = r / k
tspan = (0.0, 20.0)
u0 = 50.0
#построение графиков  x(t) при разных r и k=100,x0=50
prob = ODEProblem(f!, u0, tspan)
sol = solve(prob)
theme(:dao)
plot(sol, linewidth=3, label="x0= 50.0 r= 0.1 k=100.0")

#План эксперимента по r 
for j in (2:12)
    r = rexp[j]
    b = bexp[j]
    rstr = string(r)
    labstr = "x0= 50.0"* " r= " * rstr*" k=100.0"
    #println(labstr)
    prob = ODEProblem(f!, u0, tspan)
    sol = solve(prob)
    plot!(sol, linewidth=3, label=labstr)
end

xlabel!("Time (t)")
ylabel!("x(t)")
title!("Семейство графиков по модели Ферхюльста")

#построение графиков xсрt(r,k)

kexp = (100.0:50.0:300.0)
D = zeros(5, 12)
u0 = 50.0
for i in (1:5)
    k = kexp[i]
    for j in (1:12)
        r = rexp[j]
        b = rexp[j]/kexp[i]
        prob = ODEProblem(f!, u0, tspan)
        sol = solve(prob)
        D[i, j] = mean(sol[:])
    end
end

plot(rexp[2:12], D[1, 2:12], label="k=100", legend=:outertopright, wsize=(800, 500), linewidth=3)
plot!(rexp[2:12], D[2, 2:12], label="k=150", linewidth=3)
plot!(rexp[2:12], D[3, 2:12], label="k=200", linewidth=3)
plot!(rexp[2:12], D[4, 2:12], label="k=250", linewidth=3)
plot!(rexp[2:12], D[5, 2:12], label="k=300", linewidth=3)
ylabel!("xср(t)")
xlabel!("r")
title!("Зависимость средней численности популяции от r и k")

3.Создайте $logistNp.ipynb$ документ, который содержит:

Наименование доумента и фамилию студента;
Постановку задачи (историческая справка, математическая модель Ферхюльста и т. д.);
Всроенные скрипты и результаты их работы.
Интерпретацию полученных результатов, ответы на контрольные вопросы и выводы.

Контрольные вопросы:

Какие характеристики логистической кривой определяются константой скорости, емкостью среды и начальной численностью?
Какие начальные условия приводят к решению с точкой перегиба?
Как меняется форма кривой, если $X_0< 0.5\cdot K$; $0.5\cdot K< X_0 < K$; $X_0>K$; $X_0 = K$

4.Сформируйте на основе созданного Вами $logistNp.ipynb$ файла $HTML$ файл

Лабораторная работа 3 "Дискретная модель логистического роста".¶

Дискретные модели применяются для описания развития популяций, численность которых в момент времени t зависит численности в k предшествующих моментов времени:

$$X(t)= F(X(t-1),X(t-2),...,X(t-k))$$

В простейшем случае численность каждого следующего поколения зависит лишь от численности предыдущего поколения, и говорят, что поколения в популяции не перекрываются. Это справедливо для многих видов насекомых, а также для некоторых синхронных культур микроорганизмов. В качестве примера дискретной модели, рассмотрим разностный аналог логистического уравнения (см. непрерывную логистическую модель):

$$dX/dt = r\cdot X\cdot(1-X/K)$$

Заменив $dX/dt$ на $\Delta X/\Delta t$ получим $\Delta X =X(t+1)-X(t)$ и $\Delta t =1$

$$X(t+1)=X(t)\cdot (1+r\cdot (1-X(t)/K))$$

Учитывая биологические соображения, преобразуем данное уравнение к виду:

$$X(t+1)=X(t)\cdot exp(r\cdot (1-X(t)/K))$$

Это уравнение можно считать разностным аналогом логистического уравнения. При различных соотношениях параметров $r$ и $К$, пользуясь этой моделью, можно получать различные режимы динамики численности популяции:

$0 < r < 1$ - монотонное приближение численности к стационарной;
$1 < r < 2$ - затухающие колебания;
$2 < г < 3$ - стационарные колебания;
$r > 3,1$ - нерегулярное поведение ( хаос ).

Цель задания - получение характеристик колебательных режимов динамики численности популяции и демонстрация возможностей и области применимости дискретного моделирования по сравнению с непрерывным. При изучении этой модели варьируется параметр скорости роста при заданных параметрах начальной численности и емкости среды.

1.Разработайте скрипт исследования дискретной модели логистического роста с построением соответствующих графиков. Скрипт для этой модели может иметь следующий вид:

скрипт 5

using Plots, PlotThemes
using StatsBase
# дискретная  логистическая модель

global n = 50
global k = 200.0 # емкость среды
global ts = (1:n)# интервал моделирования
global x0 = 50.0 # начальная численность популяции
global n_ts = length(ts)

function Y(n_r, ri, x, D_r)
    for i in (1:n_r1)
        r = ri[i]
        x[i] = x0
        D_r[i, 1] = x[i]
        for j in (1:n_ts-1)
            x[j+1] = x[j] * exp(r * (1 - x[j] / k))
            D_r[i, j+1] = x[j+1]
        end
    end
    [D_r[1, :] D_r[2, :] D_r[3, :] D_r[4, :] D_r[5, :]]
end

x1 = zeros(n_ts)
ri1 = (0.2:0.2:1.0)# первый диапозон изменения параметра  0.0<r<1.0
n_r1 = length(ri1)# число значений параметра r
D_r1 = zeros(n_r1, n_ts)
y1 = Y(n_r1, ri1, x1, D_r1)

x2 = zeros(n_ts)
ri2 = (1.2:0.2:2.0)# второй диапозон изменения параметра  0.0<r<1.0
n_r2 = length(ri2)# число значений параметра r
D_r2 = zeros(n_r2, n_ts)
y2 = Y(n_r2, ri2, x2, D_r2)

x3 = zeros(n_ts)
ri3 = (2.2:0.2:3.0)#третий диапозон изменения параметра  0.0<r<1.0
n_r3 = length(ri3)#число значений параметра r
D_r3 = zeros(n_r3, n_ts)
y3 = Y(n_r3, ri3, x3, D_r3)

x4 = zeros(n_ts)
ri4 = (4.0:1.0:8.0)#четвертый диапозон изменения параметра  0.0<r<1.0
n_r4 = length(ri4)#число значений параметра r
D_r4 = zeros(n_r4, n_ts)
y4 = Y(n_r4, ri4, x4, D_r4)

theme(:dao)
p1 = plot(ts, y1, title="0 <r< 1,k=200,x0=50", ylabel="x(t)", xlabel="t", lw=2)
p2 = plot(ts, y2, title="1 <r< 2,k=200,x0=50", ylabel="x(t)", xlabel="t", lw=3)
p3 = plot(ts, y3, title="2 <r< 3,k=200,x0=50", ylabel="x(t)", xlabel="t", lw=2)
p4 = plot(ts, y4, title="3 <r< 8,k=200,x0=50", ylabel="x(t)", xlabel="t", lw=2)
plot(p1, p2, p3, p4, layout=(2, 2), wsize=(1000, 700), legend=false)

Поэкспериментируйте со скриптом, наблюдая отклики модели при различных планах изменения параметров.

2.Создайте $logistD.ipynb$ документ, который содержит:

Наименование доумента и фамилию студента;
Постановку задачи (историческая справка, дискретная логистическая модель и т. д.);
Всроенные скрипт и результаты его работы.
Интерпретацию полученных результатов, ответы на контрольные вопросы и выводы.

Контрольные вопросы:

При каких значениях $r$ дискретная и непрерывная логистическая модели дают одинаковые решения?
При каких значениях $r$ получаются устойчивые решения и в чем их отличие друг от друга?
При каких значениях $r$ модель может описывать хаотические вспышки численности насекомых?

4.Сформируйте на основе созданного Вами $logistD.ipynb$ файла $HTML$ файл

Лабораторная работа 4 Исследование динамики двух видов популяций¶

Основная цель работы: исследование взаимоотношений «хищник-жертва» на фазовой плоскости и во времени методом математического моделирования.

Теоретическое введение. Биологические системы вступают во взаимодействие друг с другом на всех уровнях, будь то взаимодействие биомакромолекул в процессе биохимических реакций, или взаимодействие видов в популяциях.В популяционной динамике принято классифицировать взаимодействия по их результатам. Наиболее распространенными и хорошо изученными являются взаимодействия конкуренции (когда численность каждого из видов в присутствии другого растет с меньшей скоростью), симбиоза (когда виды способствуют росту друг друга) и типа хищник-жертва или паразит-хозяин (когда численность вида-жертвы в присутствии вида-хищника растет медленнее, а вида-хищника - быстрее). В природе также встречаются взаимодействия, когда один из видов чувствует присутствие второго, а другой - нет (аменсализм ** и **комменсализм), или виды нейтральны.

Первое глубокое математическое исследование закономерностей динамики взаимодействующих популяций дано в книге В Вольтерра "Математическая теория борьбы за существование" (1931). Крупнейший итальянский математик Вито Вольтерра - основатель математической биологии предложил описывать взаимодействие видов подобно тому, как это делается в статистической физике и химической кинетике, в виде мультипликативных членов в уравнениях (произведений численностей взаимодействующих видов). Тогда в общем виде с учетом самоограничения численности по логистическому закону система дифференциальных уравнений, описывающая взаимодействие двух видов, может быть записана в форме:

$$dy_1/dt = a_1\cdot y_1 + b_{12}\cdot y_1\cdot y_2 - c_1\cdot y_1^{2}$$ $$dy_2/dt = a_2\cdot y_2 + b_{21}\cdot y_1\cdot y_2 - c_2\cdot y_2^{2}$$

Здесь параметры $a_i$ - константы собственной скорости роста видов, $c_i$ - константы самоограничения численности (внутривидовой конкуренции),$b_{ij}$ - константы взаимодействия видов, $(i,j=1,2)$. Соответствие знаков этих последних коэффициентов различным типам взаимодействий приведено в таблице:

Таблица типов взаимодействий популяций

Исследование свойств моделей подобного типа приводит к некоторым важным выводам относительно исхода взаимодействия видов. Уравнения конкуренции $( b_{12}>0 , b_{21} <0)$ предсказывают выживание одного из двух видов, в случае если собственная скорость роста другого вида меньше не-которой критической величины. Оба вида могут сосуществовать, если произведение коэффициентов межпопуляционного взаимодействия меньше произведе-ния коэффициентов внутри популяционного взаимодействия: $b_{12}b_{21} < c_2c_1$ .

Вольтерр предположил по аналогии со статистической физикой, что интенсивность взаимодействия пропорциональна вероятности встречи (вероятности столкновения молекул), то есть произведению концентраций. Это и некоторые другие предположения позволили построить математическую теорию взаимодействия популяций одного трофического уровня (конкуренция) или разных трофических уровней (хищник-жертва).Системы, изученные Вольтерра, состоят из нескольких биологических видов и запаса пищи, который используют некоторые из рассматриваемых видов. О компонентах системы формулируются следующие допущения:

Пища либо имеется в неограниченном количестве, либо ее поступление с течением времени жестко регламентировано.
Особи каждого вида отмирают так, что в единицу времени погибает по-стоянная доля существующих особей.
Хищные виды поедают жертвы, причем в единицу времени количество съеденных жертв всегда пропорционально вероятности встречи особей этих двух видов, т.е. произведению количества хищников на количество жертв.
Если имеются пища в неограниченном количестве и несколько видов, ко-торые способны ее потреблять, то доля пищи, потребляемая каждым видом в единицу времени, пропорциональна количеству особей этого вида, взятого с некоторым коэффициентом, зависящим от вида (модели межвидовой конкуренции).
Если вид питается пищей, имеющейся в неограниченном количестве, прирост численности вида за единицу времени пропорционален численности вида.
Если вид питается пищей, имеющейся в ограниченном количестве, то его размножение регулируется скоростью потребления пищи, т.е. за единицу вре-мени прирост пропорционален количеству съеденной пищи.

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

Классическая модель Лотки и Вольтера.¶

Базовой моделью незатухающих колебаний служит классическое уравнение Вольтерра, описывающее взаимодействие видов типа хищник-жертва.

Рассмотрим модель взаимодействия хищников и их добычи, когда между особями одного вида нет соперничества. Пусть $y_1$ и $y_2$ число жертв и хищников соответственно. Предположим, что относительный прирост жертв $y_1^{'}/y_1$ = $\alpha-\beta x_2$, $\alpha>0$, $\beta>0$, где $\alpha$ — скорость размножения жертв в отсутствие хищников, $\beta y_2$ - потери от хищников. Развитие популяции хищников зависит от количества пищи (жертв), при отсутствии пищи ( $y_1=0$ ) относительная скорость изменения популяции хищников $y_2^{'}/x_2 =-\gamma$, $\gamma >0$ , наличие пищи компенсирует убывание, и при $y_1>0$ имеем: $$y_2^{'}/y_2 = (-\gamma + \delta y_1), \delta >0$$

Систему Вольтера - Лотка можно представить в следующем виде:

$$dy_1/dt =\alpha y_1-\beta y_2y_1$$ $$dy_2/dt =-\gamma y_2 + \delta y_1y_1$$

$\alpha,\beta , \gamma , \delta > 0$

Скрипт реализующий решение этой системы дифференциальных уравнений может иметь следующий вид:

скрипт 6

using DifferentialEquations, Plots, PlotThemes
using StatsBase
#Классическая модель Лотки и Вольтера.
function volterra!(du, u, p, t)
    alfa, beta, gamma,delta = p
    du[1] = alfa * u[1] - beta*u[2]*u[1]
    du[2] = -gamma * u[2] + delta*u[1]*u[2]
end

#параметры модели
alfa = 4.0
beta = 2.5
gamma = 2.0
delta = 1.0
u0=[3.0,1.0]
tspan=(0.0,40.0)
p = (alfa, beta, gamma,delta)

prob = ODEProblem(volterra!, u0, tspan, p)
sol = solve(prob,abstol=1e-8, reltol=1e-8)

theme(:dao)
#plotly()
y=[sol[1,:],sol[2,:]]
p1 = plot(sol.t, y, title="Численность жертв y1 и хищников y1", ylabel="y1(t),y2(t)", xlabel="t",lw=2)
p2 = plot(y[1], y[2], title="ФАЗОВЫЙ ПОРТРЕТ y1,y2", xlabel="y1", ylabel="y2", lw=3)
plot(p1, p2, layout=(2, 1), wsize=(850, 800), legend=false)

Видно, что процесс имеет колебательный характер. При заданном начальном соотношении числа особей обоих видов 3 : 1, обе популяции сначала растут. Когда число хищников достигает величины $\beta$, популяция жертв не успевает восстанавливаться и число жертв начинает убывать. Уменьшение количества пищи через некоторое время начинает сказываться на популяции хищников и когда число жертв достигает величины $y_1=\gamma/\delta$ (в этой точке $y_2^{'}=0$), число хищников тоже начинает сокращаться вместе с сокращением числа жертв. Сокращение популяций происходит до тех пор, пока число хищников не достигнет величины $y_2= \alpha/\beta$ (в этой точке $y_1^{'}=0$).С этого момента начинает расти популяция жертв, через некоторое время пищи становится достаточно, чтобы обеспечить прирост хищников, обе популяции растут, и процесс повторяется снова и снова. На графике четко виден периодический характер процесса. Периодичность процесса явственно видна на фазовой плоскости — фазовая кривая ($y_1(t), y_2(t)$) — замкнутая линия. Самая левая точка, этой кривой, - это точка, в которой число жертв достигает наименьшего значения. Самая правая точка, - точка пика популяции жертв. Между этими точками количество хищников сначала убывает, до нижней точки фазовой кривой, где достигает наименьшего значения, а затем растет до верхней точки фазовой кривой. Если в начальный момент система находилась в стационарной точке, то решения $y_1(t), y_2(t)$ не будут изменяться во времени, останутся постоянными. Всякое же другое начальное состояние приводит к периодическому колебанию решений. Неэллиптичность формы траектории, охватывающей центр, отражает негармонический характер колебаний.

Уравнения Вольтерра-Лотка с логистической поправкой.¶

Рассмотрим модель конкурирующих видов с "логистической поправкой":

$$dy_1/dt =\alpha y_1-\beta y_2y_1 - \mu y_1^{2}$$ $$dy_2/dt =-\gamma y_2 + \delta y_1y_1 -\mu y_2^{2}$$

Внесем в скрипт необходимые изменения и получим следующие результаты ($\mu =0.1$)

скрипт 7

using DifferentialEquations, Plots, PlotThemes
using StatsBase
#Модель Вольтера с логистической поправкой.
function volterra2!(du, u, p, t)
    alfa, beta, gamma, delta, mu = p
    du[1] = alfa * u[1] - beta * u[2] * u[1]- mu*u[1]^2
    du[2] = -gamma * u[2] + delta * u[1] * u[2]- mu*u[2]^2
end

#параметры модели
alfa = 4.0
beta = 2.5
gamma = 2.0
delta = 1.0
mu= 0.1
u0 = [3.0, 1.0]
tspan = (0.0, 20.0)
p = (alfa, beta, gamma, delta, mu)

prob = ODEProblem(volterra2!, u0, tspan, p)
sol = solve(prob, abstol=1e-8, reltol=1e-8)
theme(:dao)
#plotly()
y = [sol[1, :], sol[2, :]]
p1 = plot(sol.t, y, title="Численность жертв y1 и хищников y2", ylabel="y1(t),y2(t)", xlabel="t", lw=3)
p2 = plot(y[1], y[2], title="ФАЗОВЫЙ ПОРТРЕТ y1,y2", xlabel="y1", ylabel="y2", lw=3)
plot(p1, p2, layout=(2, 1), wsize=(800, 800), legend=false)

В этом случае поведение решений в окрестности стационарной точки меняется в зависимости от величины и знака параметра $\mu$ . Рассмотрим фазовый портрет системы Вольтерра—Лотка для $\mu =0.1, \alpha =4, \beta =2.5, \gamma =2, \delta =1$ и графики ее решения с начальным условием $y_1(0)=3, y_2(0)=1$.

Видно, что в этом случае стационарная точка превращается в устойчивый фокус, а решения — в затухающие колебания. При любом начальном условии состояние системы через некоторое время становится близким к стационарному и стремится к нему при $t\longrightarrow \infty$ . В случае при отрицательном значении параметра $\mu$, стационарная точка является неустойчивым фокусом и амплитуда колебаний численности видов растет. В этом случае как бы близко ни было начальное состояние к стационарному, с течением времени состояние системы будет сильно отличаться от стационарного. Рассмотренные модели могут описывать поведение конкурирующих фирм, рост народонаселения, численность воюющих армий, изменение эколо-гической обстановки, развитие науки и пр.

Порядок выполнения задания.¶

Цель задания: Разработка (использование) скриптов для решения соотвествующих систем дифференциальных уравнений и построение на этой основе графиков динамики переменных $y_1$ (жертв) и $y_2$ (хищников) во времени и фазовом пространстве.

Задайте в первом скрипте интервал моделирования в диапазоне $(0-20)$, значения вектора коэффициентов правых частей классической системы $(0< \alpha, \beta, \gamma, \delta, <10)$,начальные значения численности $y_1 ,y_2$ $(0 < y_1 < 5, 0 < y_2< 5)$. Запустите скрипт на исполнение. Просмотрите результаты. Поэкспериментируйте с моделью оцените полученные результаты.
Задайте во втором скрипте дополнительно параметр $\mu$ в диапазоне $(0.1-1.0)$. Запустите скрипт на исполнение. Просмотрите результаты. Поэкспериментируйте с моделью оцените полученные результаты.
Создайте $volterra.ipynb$ документ, который содержит:

Наименование доумента и фамилию студента;
Постановку задачи (историческая справка, математическая модель Вольтера, основные принципы взаимодействия популяций и т. д.);
Всроенные скрипты и результаты их работы.
Интерпретацию полученных результатов, ответы на контрольные вопросы и выводы.

Контрольные вопросы:

При каких значениях параметра кривые системы (2) и (3) подобны?
Как значения параметра $\mu$ для системы с поправкой влияют на форму кривых?
Как начальные значения $y_1$ и $y_2$ влияют на форму кривых?
Дайте определение константам $\alpha,\beta , \gamma , \delta$, как их значения влияют на формы кривых классической системы и системы с поправкой.

Сформируйте на основе созданного Вами $volterra.ipynb$ файла $HTML$ файл

Литература¶

Гиляров А.М. Популяционная экология, М.: Из-во МГУ, 1990.
Демидович Б.П. Лекции по математической теории устойчивости. М.: Наука, 1967
Зарядов И.С. Введение в статистический пакет R. Учебно-методическое пособие, РУДН, 2010.
Кабаков Р. R в действии. Анализ и визуализация данных в программе R/пер.с англ. М.: ДМК Пресс, 2016.-588 с.: ил.
Мастицкий С. Э.Визуализация данных с помощью ggplot2.М.: ДМК Пресс, 2017.
Пискунов Н.С. Дифференциальное и интегральное исчисления. Для втузов. Учеб-ник. М.: Наука, 1985.
Полуэктов PA., Пых ЮА., Швытов ИА. Динамические модели экологических систем. Л.: Гидрометеоиздат. 1980.
Ризниченко Г.Ю., Рубин А.Ь. Математические модели биологических продукционных процессов: Учебное пособие. М.: Из-во МГУ, 1993.
Смит Дж. М, Модели в экологии, М.: Мир, 1976.
Федоров В.Д., Гильманов Г. Г. Экология. М.: Из-во МГУ. 1980.
Yihui Xie, J.J. Allaire, Garrett Grolemund (2019).R Markdown The Definitive Guide. Published July 31, 2018 by Chapman and Hall/CRC 304 Pages.ссылка на книгу

Приложения¶

Элементы математические теории устойчивости систем¶

Динамическое описание объектов в виде систем обыкновенных дифференциальных уравнений часто используются для исследования реальных процессов в различных областях науки и, в частности, биологии, экологии, химии, экономики. Решения системы обыкновенных дифференциальных уравнений могут быть устойчивыми или неустойчивыми. Введем понятие устойчивого и неустойчивого решения системы дифференциальных уравнений:

$$\bar X^{'} = f(t,\bar X)\ldots(1)$$

с начальными условиями $\bar X(t_0)=\bar X_0$ где $\bar X=(X_1,X_2,\ldots,X_n)$

$n$ - мерный вектор; $t \in I = [t_0, + \infty)$ - независимая переменная, по которой производится дифференцирование;

$\bar X^{'}=(dX_1/dt,dX_2/dt,\ldots,dX_n/dt)$ - $n$ -мерный вектор производных по $t$

$\bar f(t,\bar X) = (f_1(t,\bar X),f_2(t,\bar X),\ldots,f_n(t,\bar X))$ $n$-мерная вектор - функция.

Если начальные данные $(t_0 ,\bar X_0)$ изменяются, то изменяется и решение. Тот факт, что решение зависит от начальных данных, обозначается следующим образом:$\bar X(t) =\bar X(t;t_0,\bar X_0)$ . Естественно, что в качестве математической модели пригодна лишь та задача Коши, которая устойчива к малым изменениям начальных данных. Определим понятие устойчивости, асимптотической устойчивости и неустойчивости в смысле Ляпунова. Для этого отклонение решения $\bar X(t) =\bar X(t;t_0,\bar X_0)$, вызванное отклонением $\bar \Delta X_0$ начального значения $\bar X_0$, будем записывать следующим образом: $$\mid \bar X(t;t_0,\bar X_0 +\bar \Delta X_0)-\bar X(t)\mid = \mid \bar X(t;t_0,\bar X_0 +\bar \Delta X_0)-\bar X(t;t_0,\bar X_0)\mid $$

Определение 1. Решение $\bar X(t) =\bar X(t;t_0,\bar X_0)$ системы (1) называется устойчивым по Ляпунову в положительном направлении (или устойчивым), если оно непрерывно по $\bar X_0$ на интервале $I = [t_0, + \infty)$, т.е. $\forall \epsilon >0$ $\exists \delta >0$ такое, что если для $\lor \bar \Delta X_0$ $\mid \bar \Delta X_0 \mid \leq \delta$ $\Rightarrow$ $\mid \bar X(t;t_0,\bar X_0 +\bar \Delta X_0)-\bar X(t)\mid \leq \delta$ $\ \forall t \geq t_0$. Если, кроме того, если отклонение решения $\bar X(t)$ стремится к нулю при $t → + ∞$ для достаточно малых $\bar \Delta X_0$ , т.е. $\exists \Delta > 0$, что $\forall \bar \Delta X_0$ $\ \mid \bar \Delta X_0\mid \leq \Delta \Rightarrow \mid \bar X(t;t_0,\bar X_0 +\bar \Delta X_0)-\bar X(t)\mid \longrightarrow 0 ,t → + ∞ $. То решение $\bar X(t)$ системы (1) называется асимптотически устойчивым в положительном направлении (или асимптотически устойчивым). Аналогично определяются различные типы устойчивости решения в отрицательном направлении.

Определение 2.Решение $\bar X(t)=\bar X(t;t_0,\bar X_0)$ системы (1) называется неустойчивым по Ляпунову в положительном направлении (или неустойчивым), если оно не является устойчивым в положительном направлении. Аналогично определяется неустойчивость в отрицательном направлении.

Геометрическая интерпретация

Геометрически устойчивость по Ляпунову решения $\bar X(t)$ можно интерпретировать следующим образом (рисунок 1): все решения ) $\bar X(t;t_0,\bar X_0 +\bar \Delta X_0)$, близкие в начальный момент $t0$ к решению $\bar X(t)$ (т.е. начинающиеся в пределах $δ$ - трубки ), не выходят за пределы $ε$ - трубки при всех значениях $t ≥ t_0$ .
Асимптотическая устойчивость есть устойчивость с дополнительным условием. Геометрически это означает, что любое решение $\bar X_1(t)$ , начинающееся в момент $t_0$ в $Δ$ - трубке, с течением времени неограниченно приближается к решению $\bar X(t)$. Трубка радиуса $Δ$ называется областью притяжения решения $\bar X(t)$.
Неустойчивость по Ляпунову геометрически означает, что среди решений, близких в начальный момент $t_0$ к решению $\bar X(t)$ найдется хотя бы одно, которое в некоторый момент $t_1$ (свой для каждого такого решения) выйдет за пределы $ε$ - трубки (рисунок 1)

Исследование устойчивости произвольного решения $\bar X(t)$ системы (1) всегда можно свести к исследованию устойчивости нулевого решения некоторой преобразованной системы. В дальнейшем будем предполагать, что система (1) имеет нулевое решение, т.е. $\bar f(t,0)= 0 \ ∀ t ≥ t_0$ , и ограничимся исследованием устойчивости нулевого решения. Переформулируем определения различных типов устойчивости для нулевого решения $\bar X(t) ≡ 0$ системы (1).

Определение 3. Нулевое решение $\bar X(t) ≡ 0$ системы (1) называется устойчи-вым по Ляпунову в положительном направлении (или устойчивым), если $∀ ε > 0 \ ∃ δ = δ ( ε ) > 0$ такое, что $∀ \bar X_0$

$\mid \bar \Delta X_0 \mid ≤ δ \Rightarrow \mid \bar X(t;t_0,\bar X_0)\mid ≤ ε$ для $∀ t ≥ t_0$

Если кроме того, $∃ Δ > 0$ , что для $∀ X_0$ $\mid \bar \Delta X_0 \mid ≤ δ \Rightarrow \mid \bar X(t;t_0,\bar X_0)\mid → 0,t → + ∞$, то решение $\bar X(t) ≡ 0$ системы (1) называется асимптотически устойчивым в положительном направлении ( или асимптотически устойчивым ).

Определение 4. Нулевое решение $\bar X(t) ≡ 0$ системы (1) называется неустойчивым по Ляпунову в положительном направлении (или неустойчиво), если оно не явля-ется устойчивым в положительном направлении.

Устойчивость решения стационарной системы линейных дифференциальных уравнений с постоянными коэффициентами.¶

Система обыкновенных дифференциальных уравнений называется автономной (или стационарной, или консервативной, или динамической), если независимая переменная не входит явно в систему уравнений. Нормальную автономную систему 2- го порядка можно записать в следующем виде:

$$dX/dt =P(X,Y)\ldots \;\;\; (2)$$

Эффективным методом исследования устойчивости динамической системы является графическое представление характеристик траекторий в фазовом пространстве, так называемый фазовый портрет системы, а также изучение зависимости фазовых пе-ременных от времени и зависимости структуры фазового портрета от параметров.

Такой подход допускает наглядное представление поведения переменных на фа-зовой плоскости $(X,Y)$. Совокупность траекторий точек $М(X,Y)$ плоскости $(X,Y)$, изображающих (представляющих) значения $X и Y$ в последовательные моменты времени $t$ соответствуют состояниям системы в процессе изменения согласно уравнениям (2).

Множество фазовых траекторий, построенных из различных начальных условий (задаваемых точкой с координатами $(Xо,Yо)$, образует так называемый фазовый портрет системы и позволяет анализировать характер изменений в системе без знания аналитических выражений решений системы уравнений (2). Особый интерес представляют стационарные состояния (точки покоя $X_{ст}Y_{ст}$.) системы, в которых производные переменных по времени равны нулю.

$$dX/dt \mid X_{ст}Y_{ст} =0 ;dY/dt \mid X_{ст}Y_{ст} =0$$

Точки покоя системы (2) могут быть устойчивыми и неустойчивыми по Ляпунову. Как известно, исследование устойчивости любого, а значит, и постоянного решения можно свести к исследованию устойчивости нулевого решения. Поэтому далее будем считать, что система (2) имеет нулевое решение $X(t)=Y(t)=0$ , т. е.$P(0,0=Q(0,0)=0$ , и точка покоя совпадает с началом координат фазового пространства. Таким образом, устойчивость нулевого решения системы (2) означает устойчивость начала координат фазового пространства системы (2), и наоборот. Рассмотрим типы стационарных состоянии и их устойчивости для нормальной системы двух линейных дифференциальных уравнений с постоянными коэффициентами:

$$dX/dt= P_1X+P_2Y \ldots \;\;\; (3)$$

$$dY/dt= P_3X+P_4Y$$

Точка $(0,0)$ является точкой покоя системы (3). Исследуем расположение траектории системы (3) в окрестности этой точки. Из теории известно, что решение можно представить в следующем виде:

$$X=\alpha_1 e^{k_t},Y=\alpha_2 e^{k_t}$$

Для определения $k$ получаем характеристическое уравнение:

$$\begin{vmatrix} P_1-k& P_2\\ P_3& P_4-k \end{vmatrix}=0\ldots \;\;\;\;(4)$$

Рассмотрим возможные случаи. $I.$ Корни характеристического уравнения действительные и различные.

$k_1 < 0, k_2 < 0.$ Точка покоя асимптотически устойчива (устойчивый узел).
$k_1 > 0, k_2 > 0.$ Точка покоя неустойчива (неустойчивый узел).
$k_1 > 0, k_2 < 0.$ Точка покоя неустойчива (седло).
$k_1 = 0, k_2 > 0.$ Точка покоя неустойчива.
$k_1 = 0, k_2 < 0.$ Точка покоя устойчива, но не асимптотически

$II.$ Корни характеристического уравнения комплексные: $k_1 = p + q_i, k_2 = p - q_i$

$p < 0 , q ≠ 0.$ Точка покоя асимптотически устойчива (устойчивый фокус).
$p > 0 , q ≠ 0.$ Точка покоя неустойчива (неустойчивый фокус).
$p = 0, q ≠ 0.$ Точка покоя устойчива (центр). Асимптотической устойчивости нет.

$III.$ Корни кратные: $k1 = k2$

$k_1 = k_2 < 0.$ Точка покоя асимптотически устойчива (устойчивый узел).
$k_1 = k_2 > 0.$ Точка покоя неустойчива (неустойчивый узел).
$k_1 = k_2 = 0.$ Точка покоя неустойчива. Возможен исключительный случай, когда все точки плоскости являются устойчивыми точками покоя. Для системы (3) двух линейных уравнений с постоянными действительными коэффициентами характеристическое уравнение (4) приводится к виду $k_2 + a_1 k + a_2 = 0.$
Если $a_1 > 0 , a_2 > 0$, то нулевое решение системы (3) асимптотически устойчиво.
Если $а_1 > 0 , a_2 = 0$, или $a_1 = 0 , a_2 > 0$ , то нулевое решение устойчиво, но не асимптотически.
Во всех остальных случаях нулевое решение неустойчиво; однако при $a_1 = a_2 = 0$ возможен исключительный случай, когда нулевое решение устойчиво, но не асимптотически. Разновидности особых точек в фазовом пространстве для системы (3) представлены на рисунке 2.

Аналитическая оценка устойчивости (вида особой точки) нулевого решения системы (3) для заданного набора значений коэффициентов правых частей (пример).¶

Создадим скрипт для численного решения системы (3) и выполним с его помощью расчеты и построение изображений динамики переменных $x$ и $y$ во времени и фазовом пространстве для набора праметров $p = (1.0, -2.0, 2.0, -1.0)$. Исследование решений системы (3) при заданных параметрах в окрестностях точки покоя (0,0) показывает, что фазовый портрет представляет собой замкнутый контур (центр) и, следовательно, система вблизи нулевой точки будет находиться в устойчивом колебательном режиме. Таким образом точка покоя устойчива, а асимптотической устойчивости нет

using DifferentialEquations, Plots, PlotThemes
using StatsBase
#Система дифференциальных уравнений 3.
function sysdif3!(du, u, p, t)
    P1, P2, P3, P4 = p
    du[1] = P1 * u[1] + P2 * u[2]
    du[2] = P3 * u[1] + P4 * u[2]
end

tspan = (0.0, 15.0)
P1 = 1.0
P2 = -2.0
P3 = 2.0
P4 = -1.0
p = (P1, P2, P3, P4)
u0 = [0.0015 0.0022]

theme(:dao)
#plotly()
prob = ODEProblem(sysdif3!, u0, tspan, p)
sol = solve(prob, abstol=1e-8, reltol=1e-8)

y = [sol[1, :], sol[2, :]]
p1 = plot(sol.t, y, title="p=(1.0  -2.0  2.0 -1.0)", ylabel=("x(t),  y(t)"), xlabel=("t"), lw=2)
p2 = plot(sol, vars=(1, 2), title="Фазовый портрет x,y", lw=2)
plot(p1, p2, layout=(1, 2), wsize=(950, 450), legend=false)

Покажем, что для заданной системы дифференциальных уравнений корни характерестического уровнения мнимые

using Symbolics
using Latexify
@variables P¹ P² P³ P⁴ k
b = [P¹-k P²
    P³ P⁴-k]
latexify(b)

(P¹ - k) * (P⁴ - k) - P³ * P²

Symbolics.expand((P¹ - k) * (P⁴ - k) - P³ * P²)

#значения коэффициентов
@variables P¹ P² P³ P⁴
ex = [:(P¹ ~ 1.0)
    :(P² ~ -2.0)
    :(P³ ~ 2.0)
    :(P⁴ ~ -1.0)]
eqs = parse_expr_to_symbolic.(ex, (Main,))

#решение характеристического уравнения: получены два мнимых корня
using Nemo, Groebner
@variables P¹ P² P³ P⁴ k
symbolic_solve(-(-2.0) * 2.0 + 1.0 * (-1.0) - 1.0 * k - (-1.0) * k + k^2)

2-element Vector{SymbolicUtils.BasicSymbolic{Complex{Real}}}:
 (0 + 1im)*√(3)
 (0 - 1im)*√(3)

ПРИЛОЖЕНИЕ 2¶

Эффект Олли и его учет в моделях роста популяций¶

Упрощенный учет минимальной численности популяции, поддерживающей ее существование.¶

В теоретической биологии достаточно долго обсуждается необходимость наличия минимальной численности особей L для поддержания существования популяции. Гипотеза популяционного L-минимума может быть упрощенно отражена в следующей модифицированной модели Ферхюльста:

$$\frac{{dx}}{{dt}} = r \cdot x \cdot \left( {1 - \left( {\frac{x}{k}} \right)^m } \right) \cdot \left( {x - L} \right)^n$$

$m = 1,\frac{1}{2},\frac{3}{4},\frac{5}{6}...{\rm n = 1}{\rm ,}\frac{{\rm 1}}{{\rm 3}},....$

Рассмотрим два случая:

$m = 1,\frac{1}{2},\frac{3}{4},\frac{5}{6}$ ; $n=1$ и $L=40 < x0=50$
$m = 1,\frac{1}{2},\frac{3}{4},\frac{5}{6}$ ; $n=1$ и $L=100 > x0=50$

Скрипт для исследования дифференциального уравнения представлен ниже. Как видно из графиков, в первом случае имеем:

$$L < x0{\rm \;\;\;\; x}\left( {\rm t} \right){\rm } \to {\rm k }{\rm , t } \to \infty$$

Во втором случаем имеет место следующее:

$$L > x0 \;\;\;x\left( t \right) \to 0{\rm , t} \to \infty$$

Параметр $m$ позволяет незначительно регулировать скорость роста(падения) численности популяции $x(t)$

using DifferentialEquations, Plots, PlotThemes
using StatsBase
f!(u, p, t) = r * u*(1-(u/k)^m)*(u-L) 
#-ПАРАМЕТРЫ МОДЕЛИ
#План эксперимента по L
Lexp = (10.0:0.5:100.0)
n_Lexp = length(Lexp)
#План эксперимента по m
mexp =[1 1/2 3/4 5/6] 
n_mexp = length(mexp)

k = 40.0
r = 0.01
tspan = (0.0, 20.0)
u0 = 50.0
D = zeros(n_mexp, n_Lexp)
m=0
L=40
for i in (1:n_mexp)
    m = mexp[i]
    for j in (1:n_Lexp)
        L = Lexp[j]
        prob = ODEProblem(f!, u0, tspan)
        sol = solve(prob, abstol=1e-8, reltol=1e-8)
        D[i, j] = mean(sol[:])
    end
end

prob = ODEProblem(f!, u0, tspan)
sol = solve(prob, abstol=1e-8, reltol=1e-8)
D

4×181 Matrix{Float64}:
 43.1087  43.1145  43.1205  43.1266  …  86.6539  86.735   87.1391  87.8717
 44.2339  44.2382  44.2449  44.2483     81.0401  81.3383  82.1973  82.7791
 43.5132  43.5197  43.526   43.5335     85.0135  85.3804  85.7347  86.0854
 43.3672  43.3738  43.3796  43.3867     85.5779  85.9436  86.3225  86.6867

theme(:dao)
plotly()
plot(Lexp[1:n_Lexp], D[1, 1:n_Lexp], label="m=1", linewidth=3)
plot!(Lexp[1:n_Lexp], D[2, 1:n_Lexp], label="m=1/2", linewidth=1)
plot!(Lexp[1:n_Lexp], D[3, 1:n_Lexp], label="m=3/4", linewidth=1)
ylabel!("xср(t)")
xlabel!("L")
title!("Зависимость средней численности популяции от L и m")
p1=plot!(Lexp[1:n_Lexp], D[4, 1:n_Lexp], label="m=5/6", linewidth=1)
p5 = plot(sol, label="x(t) L=100>x0  m=5/6", linewidth=2)
plot(p1, p5, layout=(1, 2), wsize=(900, 450), legend=:bottomleft)

plot(Lexp[1:n_Lexp], D[1, 1:n_Lexp], label="m=1",  linewidth=3)
plot!(Lexp[1:n_Lexp], D[2, 1:n_Lexp], label="m=1/2", linewidth=1)
plot!(Lexp[1:n_Lexp], D[3, 1:n_Lexp], label="m=3/4", linewidth=1)
ylabel!("xср(t)")
xlabel!("L")
title!("Зависимость средней численности популяции от L и m")
p2 = plot!(Lexp[1:n_Lexp], D[4, 1:n_Lexp], label="m=5/6", linewidth=1)
p6 = plot(sol, label="x(t) L=40<x0  m=5/6", linewidth=2)
plot(p2, p6, layout=(1, 2), wsize=(900, 450), legend=:bottomleft,)

Принцип Оли в модели Базыкина¶

Не вдаваясь в биологические подробности, приведем модель Базыкина, в которой реализуется принцип Оли:

$$\frac{{dx}}{{dt}} = \frac{{c \cdot x^2 }}{{M + x}} - r \cdot x - b \cdot x^2$$

Здесь первый член представляет собой коэффициенет рождаемости, а $r$ и $b=r/k$ параметры модели Ферхюльста. Исследуем правую часть дифференциального уравнения. В Wolfram Mathematica и методом подбора горафически найдем параметры дифференциального уравнения Базыкина (с=0.738, M=215, r=0.1 ). Зададим k= 1000 и b=r/k=0.0001. Расчитаем с точностью до трех знаков стационарные точки, где производная равняется нулю. Срипт для расчета стационарных точеприведен ниже.

#Поиск стационарных точек для модели Базыкина
xst=0.038
z=0.0
for x in (0.1:0.0002:6200.0)
    z=0.738*x^2/(215+x)-0.1*x -0.0001*x^2
    if (z>=-0.00002)&&(z <= 0.00002) 
        xst=x
        println("хстац=",round(xst, digits=3)) 
    end  
end

хстац=35.074
хстац=35.074
хстац=35.074
хстац=6129.926

using Roots
f(x) = 0.738 * x^2 / (215 + x) - 0.1 * x - 0.0001 * x^2
xct1=find_zero(f, (1.0, 40.0))
xct2=find_zero(f, (40.0, 6200.0))
println("хстац1=", round(xct1, digits=3))
println("хстац2=", round(xct2, digits=3))

Покажем стационарные точки графически. для того, чтобы увидеть все стационарные точки, построим графики с разным масштабом $x$

using DifferentialEquations, Plots, PlotThemes
x=(0.0:0.001:40)
z = @. 0.738 * x^2 / (215 + x) - 0.1 * x - 0.0001 * x^2
theme(:dao)
plotly()
plot(x, z, label="dx/dt=z(x)", linewidth=3)
plot!([0.0 35.074], [0.0 0.0],label=["xстац=0.0"  "хстац=35.074"], 
    seriescolor=:green,
    marker=:circle,
    markersize=6,
    markercolor=:green,
    markerstrokecolor=:green)
ylabel!("dx/dt")
xlabel!("x")

x = (1000.0:0.1:6200.0)
z = @. 0.738 * x^2 / (215 + x) - 0.1 * x - 0.0001 * x^2
theme(:dao)
plotly()
plot(x, z, label="dx/dt=z(x)", linewidth=3)
plot!([6129.926], [0.0], label="xстац=6129.926",
    seriescolor=:green,
    marker=:circle,
    markersize=6,
    markercolor=:green,
    markerstrokecolor=:green)
ylabel!("dx/dt")
xlabel!("x")

Таким образом, уравнение Базыкина имеет три стационарные точки: $x=0,x=k1=35.074$ и $x=k2=6129.926$. Итак, в случае, когда наблюдается эффект Олли, в популяции возможны два ненулевых равновесных значения численности $k1, k2$ при этом $k1<k2$, причем $k2$ - устойчивое равновесие, а $k1$ - неустойчивое. Это означает, что в популяции существует критический порог уровня численности, равный $k1$. Если начальное значение численности $x(0)$ окажется больше $k1$ , то с течением времени численность популяции приближается к равновесному значению $k2$. Если же $x(0)$ меньше $k1$, то численность популяции монотонно убывает до нуля и популяция вымирает.Следовательно, любое (в том числе антропогенное) снижение численности такой популяции ниже критического уровня чревато ее вырождением. Посмотрим, как это моделируется с помощью в julia

using DifferentialEquations, Plots, PlotThemes
using StatsBase
f!(u, p, t) = ((c*u^2)/(M+u)) - r * u - b*u ^2 
#-ПАРАМЕТРЫ МОДЕЛИ
#План эксперимента по L
c=0.738
M=215
k = 1000.0
r = 0.1
b=r/k
tspan = (0.0, 40.0)
u0 = 50.0
prob = ODEProblem(f!, u0, tspan)
sol = solve(prob, abstol=1e-8, reltol=1e-8);

p1=plot(sol.t, sol.u, title="x0 > k1",legend=false,lw=3)
ylabel!("x(t)")
xlabel!("t")

u0=35
tspan = (0.0, 120.0)
prob = ODEProblem(f!, u0, tspan)
sol = solve(prob, abstol=1e-8, reltol=1e-8);
p2 = plot(sol.t, sol.u, title="x0 < k1", lw=3)
xlabel!("t")
plot(p1, p2, layout=(1, 2), wsize=(900, 450), legend=false)

ПРИЛОЖЕНИЕ 3¶

Реализация классической модели Мальтуса с использованием инструментов визуального моделировани Engee

using RDatasets, Plots, PlotThemes
using CSV,DataFrames

#открытие модели
modelName = "Maltus01";
PID_model = modelName in [m.name for m in engee.get_all_models()] ? engee.open( modelName ) : engee.load( "/user/Vamscripts/$(modelName).engee");
engee.set_param!(modelName, "StopTime" => 6.0) # интервал моделирования

#Программа исследования классической модели Мальтуса 
r=0.5
#Начальная ЗНАЧЕНИЕ ДОХОДА - x0
x0=10.0
#Выполнение модели
engee.run(modelName) 

# ВЫДЫЧА РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ
#engee.get_results("Maltus01")
#modeling_variable = engee.get_results("Maltus01")
df_Xt1=CSV.read("Xt.csv",DataFrame)

#ГРАФИК
#plotly()
theme(:ggplot2)

plot(df_Xt1[!,"time"], df_Xt1[!,"1"], label="xt", legend=:outerbottomright, wsize=(800, 500), linewidth=3)
xlabel!("t")
ylabel!("Численность популяции хt")
title!("Классическая модель Мальтуса")

using RDatasets, Plots, PlotThemes
using CSV,DataFrames

#открытие модели
modelName = "Maltus01";
PID_model = modelName in [m.name for m in engee.get_all_models()] ? engee.open( modelName ) : 
engee.load( "/user/Vamscripts/$(modelName).engee");
engee.set_param!(modelName, "StopTime" => 6.0) # интервал моделирования

#Программа исследования классической модели Мальтуса 
#План експеримента по r
rexp=(0.2:0.1:0.6)
r=rexp[1]
#Начальная ЗНАЧЕНИЕ ДОХОДА - x0
x0=10.0
#Выполнение модели
engee.run(modelName) 
# ВЫДЫЧА РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ
df_Xt1=CSV.read("Xt.csv",DataFrame)
#ГРАФИК
#plotly()
theme(:ggplot2)
plot(df_Xt1[!,"time"], df_Xt1[!,"1"], label="x0 10.0 r=0.2", legend=:outerbottomright, wsize=(800, 500), linewidth=3)
#План эксперимента по r 
for j in (2:length(rexp))
    r = rexp[j]
    rstr = string(r)
    labstr = "x0= 50.0"* " r= " * rstr
    engee.run(modelName)
    df_Xt1=CSV.read("Xt.csv",DataFrame)
    plot!(df_Xt1[!,"time"], df_Xt1[!,"1"], label=labstr, linewidth=3)
end

xlabel!("t")
ylabel!("Численность популяции хt")
title!("Классическая модель Мальтуса")

ПРИЛОЖЕНИЕ 4¶

Реализация модели логистического роста Ферхюльста в Engee:

#Программа исследования непрерывной модели логистического роста по r
using RDatasets, Plots, PlotThemes
using CSV,DataFrames
#открытие модели
modelName = "Logrost";
PID_model = modelName in [m.name for m in engee.get_all_models()] ? engee.open( modelName ) : 
engee.load( "/user/Vamscripts/$(modelName).engee");
engee.set_param!(modelName, "StopTime" => 100.0) # интервал моделирования
#План експеримента по r
rexp=(0.0:0.1:0.6)
r=rexp[2]
#План 'эксперимента по K
k1 = 2500.0;
k2 = 5000.0;
kexp=(k1:500.0:k2)
k = kexp[1] 
x0 = 50.0;
emk = 1.0/k

#Выполнение модели
engee.run(modelName)
# ВЫДЫЧА РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ
df_Xt = CSV.read("Xt.csv", DataFrame)
theme(:ggplot2)
plot(df_Xt[!, "time"], df_Xt[!, "1"], label="x0 =50.0 r=0.1 k=2500.0", 
legend=:outerbottomright, wsize=(800, 500), linewidth=3)
#Реализация плана эксперимента по r 
for j in (3:length(rexp))
    r = rexp[j]
    rstr = string(r)
    labstr = "x0= 50.0" * " r= " * rstr*" "*"k=2500.0"
    engee.run(modelName)
    df_Xt = CSV.read("Xt.csv", DataFrame)
    plot!(df_Xt[!, "time"], df_Xt[!, "1"], label=labstr, linewidth=3)
end

xlabel!("t")
ylabel!("Численность популяции Xt")
title!("Модель логистического роста Ферхюльста")

#Программа исследования непрерывной модели логистического роста по k и r
using RDatasets, Plots, PlotThemes
using CSV, DataFrames

#открытие модели
modelName = "Logrost";
PID_model = modelName in [m.name for m in engee.get_all_models()] ? engee.open(modelName) :
            engee.load("/user/Vamscripts/$(modelName).engee");
engee.set_param!(modelName, "StopTime" => 100.0) # задание интервала моделирования

#План експеримента по r
rexp = (0.0:0.005:0.2)
#r = rexp[1]
#План 'эксперимента по K
k1 = 2500.0;
k2 = 5000.0;
kexp = (k1:500.0:k2)
k = kexp[1]
x0 = 50.0

D = zeros(length(kexp), length(rexp))

#Реализация плана эксперимента по k и r 
for i in (1:length(kexp))
    k = kexp[i]
    emk=1/k
    for j in (2:length(rexp))
        r = rexp[j]
        engee.run(modelName)
        df_Xt = CSV.read("Xt.csv", DataFrame)
        D[i, j] = df_Xt[end,2]
    end
end
D

6×41 Matrix{Float64}:
 0.0  81.3709  131.341  209.3    327.093  …  2500.0  2500.0  2500.0  2500.0
 0.0  81.5448  132.071  211.59   333.365     3000.0  3000.0  3000.0  3000.0
 0.0  81.6695  132.597  213.256  337.995     3500.0  3500.0  3500.0  3500.0
 0.0  81.7633  132.995  214.523  341.553     4000.0  4000.0  4000.0  4000.0
 0.0  81.8364  133.306  215.519  344.372     4500.0  4500.0  4500.0  4500.0
 0.0  81.8949  133.555  216.323  346.661  …  5000.0  5000.0  5000.0  5000.0

# ВЫДЫЧА РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ
theme(:ggplot2)
plot(rexp[2:length(rexp)], D[1, 2:length(rexp)], label="k=2500.0", legend=:outerbottomright, wsize=(800, 500), linewidth=3)

for  i   in  (2:length(kexp))
    k=kexp[i]
    kstr = string(k)
    labk ="k = "*kstr
    plot!(rexp[2:length(rexp)], D[i, 2:length(rexp)], label=labk, linewidth=3)
end

xlabel!("r")
ylabel!("Численность популяции Xtn")
title!("Модель логистического роста Ферхюльста (x0=50)")

ПРИЛОЖЕНИЕ 5¶

Реализация дискретной модели логистического роста Ферхюльста в Engee:

#=Программа исследования дискретой модели логистического роста по  r 
(с использванием функции пользователя и визуальной модели Engee)=#

using RDatasets, Plots, PlotThemes
using CSV, DataFrames

# дискретная  логистическая модель
#открытие модели
modelName = "LogDrost";
PID_model = modelName in [m.name for m in engee.get_all_models()] ? engee.open(modelName) :
            engee.load("/user/Vamscripts/$(modelName).engee");

#задание глобальных параметров модели
global k = 5000.0 # емкость среды
global tn = 30 #конец интервала
global h = 0.5 #шаг
global x0 = 100.0 # начальная численность популяции
global n_ts = tn / h + 1
 
#задание функции
function Y(ri,n_r,D_r)
    for i in (1:n_r)
    global r = ri[i]
    engee.run(modelName)
    df_Xdt = CSV.read("Xdt.csv", DataFrame)
    D_r[i, :] = (df_Xdt[!, 2])
    end
    y=[D_r[1, :] D_r[2, :] D_r[3, :] D_r[4, :] D_r[5, :]]
    return y
end

#задание параметров Y()
ri1 = (0.2:0.2:1.0)# первый диапозон изменения параметра  0.0<r<1.0
ri2 = (1.2:0.2:2.0)# второй диапозон изменения параметра  1.0<r< 2.0
ri3 = (2.2:0.2:3.0)# третий диапозон изменения параметра  2.0<r< 3.0
ri4 = (4.0:1.0:8.0)# четвертый диапозон изменения параметра  3.0<r< 8.0
n_r1 = length(ri1)# число значений параметра r1
D_r1 = zeros(n_r1, Int(n_ts))
n_r2 = length(ri2)# чиcло значений параметра r2 
D_r2 = zeros(n_r2, Int(n_ts))
n_r3 = length(ri3)# чиcло значений параметра r3 
D_r3 = zeros(n_r3, Int(n_ts))
n_r4 = length(ri4)# чиcло значений параметра r4 
D_r4 = zeros(n_r4, Int(n_ts))

#расчет координат графиков с помощью функции Y()
y1 = Y(ri1,n_r1,D_r1)
y2 = Y(ri2,n_r2, D_r2)
y3 = Y(ri3,n_r3, D_r3)
y4 = Y(ri4,n_r4, D_r4)

#вывод графиков
theme(:dao)
p1 = plot((1:n_ts), y1, title="0 <r< 1,k=200,x0=50", ylabel="x(t)", lw=2)
p2 = plot((1:n_ts), y2, title="1 <r< 2,k=200,x0=50", ylabel="x(t)", lw=1)
p3 = plot((1:n_ts), y3, title="2 <r< 3,k=200,x0=50", ylabel="x(t)", xlabel="t", lw=1)
p4 = plot((1:n_ts), y4, title="3 <r< 8,k=200,x0=50", ylabel="x(t)", xlabel="t", lw=1)
plot(p1, p2, p3, p4, layout=(2, 2), wsize=(1000, 700), legend=false)

ПРИЛОЖЕНИЕ 6¶

Реализация модели решения линейной системы обыкновенных дифференциальных уравнений в Engee:

#ИССЛЕДОВАНИЕ УСТОЙЧИВОСТИ НУЛЕВЫХ РЕШЕНИЙ СИСТЕМЫ ДИФУРАВНЕНИЙ
using RDatasets, Plots, PlotThemes
using CSV, DataFrames
#Задание коэффициентов правых частей
P = [1.0, -2.0, 2.0, -1.0]
#Начальные условия
U0 = [0.020, 0.020]

#открытие модели
modelName = "SysDifyr";
PID_model = modelName in [m.name for m in engee.get_all_models()] ? engee.open(modelName) :
            engee.load("/user/Vamscripts/$(modelName).engee");
engee.set_param!(modelName, "StopTime" => 10.0) # задание интервала моделирования

engee.run(modelName)
df_Xt = CSV.read("Xt.csv", DataFrame)
df_Yt = CSV.read("Yt.csv", DataFrame)
y = [df_Xt[!, "1"], df_Yt[!, "1"]]
#вывод графиков
theme(:dao)
p1 = plot(df_Xt[!, "time"], y, xlabel="t", ylabel="Xt,Yt", lw=3)
p2 = plot(df_Xt[!, "1"], df_Yt[!, "1"], xlabel="X", ylabel="Y", title="Фазовый портрет x,y", lw=3)
plot(p1, p2, layout=(2, 1), wsize=(400, 700), legend=false)

ПРИЛОЖЕНИЕ 7¶

Реализация модели Волтерра-Лотки в Engee:

using RDatasets, Plots, PlotThemes
using CSV, DataFrames
#Модель Вольтера с логистической поправкой.

#Задание коэффициентов правых частей уравнений и начальных условий
a = 4.0
b = 2.5
с = 2.0
d = 1.0
β = 0.1
x10 = 3.0
x20 = 1.0

v = [a, b, c, d, β, x10, x20]#входной вектор параметров модели

#открытие модели
modelName = "VoltLog";
PID_model = modelName in [m.name for m in engee.get_all_models()] ? engee.open(modelName) :
            engee.load("/user/Vamscripts/$(modelName).engee");
engee.set_param!(modelName, "StopTime" => 20.0) # задание интервала моделирования

#выполнение модели
engee.run(modelName)
df_X1t = CSV.read("X1t.csv", DataFrame)
df_X2t = CSV.read("X2t.csv", DataFrame)
y = [df_X1t[!, "1"], df_X2t[!, "1"]]
#вывод графиков
theme(:dao)
p1 = plot(df_X1t[!, "time"], y, xlabel="t", ylabel="X1t,X2t", lw=3)
p2 = plot(df_X1t[!, "1"], df_X2t[!, "1"], xlabel="X1", ylabel="X2", title="Фазовый портрет x1,x2", lw=3)
plot(p1, p2, layout=(2, 1), wsize=(400, 700), legend=false)

Классические динамические модели развития популяций

КОМПЬЮТЕРНЫЙ ПРАКТИКУМ ПО КУРСУ¶

Системный анализ и основы компьютерного моделирования экосистем¶

Cоставил: к.э.н. Варюхин А.М.¶

Саратов 2025¶

МЕТОДИЧЕСКИЕ УКАЗАНИЯ И ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТ¶

"Исследование моделей роста популяций"¶

Лабораторная работа 1 "Модель экспоненциального роста".¶

Порядок выполнения работы¶

ЭКСПОНЕНЦИАЛЬНАЯ МОДЕЛЬ¶

скрипт 2¶

Лабораторная работа.2 " Непрерывная модель логистического роста".¶

ЛОГИСТИЧЕСКАЯ МОДЕЛЬ¶

Порядок выполнения работы¶

Лабораторная работа 3 "Дискретная модель логистического роста".¶

Лабораторная работа 4 Исследование динамики двух видов популяций¶

Классическая модель Лотки и Вольтера.¶

Уравнения Вольтерра-Лотка с логистической поправкой.¶

Порядок выполнения задания.¶

Литература¶

Приложения¶

Элементы математические теории устойчивости систем¶

Устойчивость решения стационарной системы линейных дифференциальных уравнений с постоянными коэффициентами.¶

Аналитическая оценка устойчивости (вида особой точки) нулевого решения системы (3) для заданного набора значений коэффициентов правых частей (пример).¶

ПРИЛОЖЕНИЕ 2¶

Эффект Олли и его учет в моделях роста популяций¶

Упрощенный учет минимальной численности популяции, поддерживающей ее существование.¶

Принцип Оли в модели Базыкина¶

ПРИЛОЖЕНИЕ 3¶

ПРИЛОЖЕНИЕ 4¶

ПРИЛОЖЕНИЕ 5¶

ПРИЛОЖЕНИЕ 6¶

ПРИЛОЖЕНИЕ 7¶

Тип

Краткое описание

Связанные материалы

Теги

Категории

Языки

Форматы

Уровень

Статус

Опубликовано

Обновлено

Источник