Динамика популяции с учетом вылова¶
В этом примере мы покажем, как легко создать и решить дифференциальное уравнение, которое позволяет моделировать динамику популяции рыб с учетом вылова.
Описание модели¶
Мы будем моделировать популяцию рыб в водоеме, подверженную следующим процессам:
- влияние рождаемости и смертности описано квадратичной функцией $ax + bx^2$
- задание по вылову ежегодно уменьшает популяцию на $h$ тонн
- исходный размер популяции равен $x_0$.
Таким образом мы создаем модель следующей задачи:
$$\frac{dx}{dt} = (a + bx) x - h$$
$$x(0) = x_0$$
Общий вид блок-схемы, при помощи которой можно моделировать эту задачу, следующий:
Все параметры легко доступны, и их можно задавать как из скрипта (через переменые), так и на самом холсте. Решить эту систему нам поможет численное интегрирование, но стоит помнить:
есть комбинации входных условий, которые делают решение нестабильным (популяция уходит в $\pm \infty$);
при нестабильности решатель с постоянным шагом может не сообщить об ошибке, поэтому в настройках модели выбран решатель с переменным шагом;
максимальный размер шага симуляции установлен в 0.01, чтобы на интервале 0..1 мы всё-таки получили плавный график, даже если по условиям задачи можно делать шаги побольше.
Запуск и анализ результатов¶
Запустим модель при помощи средств программного управления:
# Загрузим модель, если она еще не открыта на холсте
if "fish_population_diff_eq_model" ∉ getfield.(engee.get_all_models(), :name)
engee.load( "$(@__DIR__)/fish_population_diff_eq_model.engee");
end
Выведем график зависимости динамики популяции от изменения ее начального размера:
# Запустим модель
model_data = engee.run( "fish_population_diff_eq_model" );
# Прочитаем нужные параметры для вывода графика
param_vector = engee.get_param("fish_population_diff_eq_model/Начальное количество рыбы (т)", "Value");
plot( model_data["x"].time, hcat( model_data["x"].value... )',
label=reshape(param_vector, 1, :), lw=2,
title="Зависимости изменения популяции от начального количества (т)", titlefont=font(10),
xlabel="Время, г", ylabel="Популяция, т" )
Теперь выставим одно значение x0
и построим график зависимости от изменения годового плана вылова:
param_1 = engee.get_param("fish_population_diff_eq_model/Годовой план вылова (т)", "Value" );
# Поменяем параметры модели
engee.set_param!( "fish_population_diff_eq_model/Начальное количество рыбы (т)", "Value"=>"5.0" );
engee.set_param!( "fish_population_diff_eq_model/Годовой план вылова (т)", "Value"=>"0.5:0.2:2" );
# # Запустим модель
model_data = engee.run( "fish_population_diff_eq_model" );
# # Получим параметры модели для графика
param_2 = engee.get_param("fish_population_diff_eq_model/Годовой план вылова (т)", "Value" );
plot( model_data["x"].time, hcat( model_data["x"].value... )',
label=reshape(param_2, 1, :), lw=2,
title="Зависимости изменения популяции от задания на вылов (т)", titlefont=font(10),
xlabel="Время, г", ylabel="Популяция, т" )
И вернем модель в исходное состояние.
engee.set_param!("fish_population_diff_eq_model/Годовой план вылова (т)", "Value"=>param_1 );
engee.set_param!("fish_population_diff_eq_model/Начальное количество рыбы (т)", "Value"=>string(param_vector) );
Заключение¶
Решение дифференциальных уравнений при помощи блок-схемы и интегратора упрощает передачу знаний между представителями разных специальностей, в чем мы убедились на примере решения дифференциального уравнения роста популяции.