Читать книгу Математические модели в естественнонаучном образовании. Том I - Денис Владимирович Соломатин - Страница 7
Глава 1. Динамическое моделирование разностными уравнениями
1.4. Вариации на тему логистической модели
ОглавлениеПредставляя дискретную логистическую модель в предыдущих разделах, старались делать модель максимально простой, чтобы сосредоточиться на разработке основных идей. Теперь, когда концепции равновесия и стабильности, а также техника построения паутинных диаграмм были разработаны, можно уделить больше внимания созданию более реалистичной модели.
Рассматривая график функции от на рисунке 1.9, для модели , одной из очевидных, но реально невозможных особенностей динамического поведения моделируемой численности, является тот факт, что парабола опускается ниже горизонтальной оси, когда отклоняемся достаточно далеко вправо. Это означает, что большие популяции становятся отрицательными на следующем временном этапе. Хотя можно интерпретировать отрицательную популяцию как вымершую, либо как долг, кредитное плечо, в экономических приложениях, но это может быть не то поведение, которое на самом деле произойдет и которое хотели бы, чтобы модель описала.
Рисунок 1.9. Модель с нереалистичными начиная с некоторого .
Возможно, более реалистичная модель допускала бы сколь угодно большие , от которых значения дают очень маленькие, но все же положительные, значения . Таким образом, популяция, значительно превышающая свою пропускную способность, может немедленно упасть до очень низких уровней, но, по крайней мере, часть популяции выживет. Графически должен зависеть от так, как показано на рисунке 1.10.
Рисунок 1.10. Новая модель с .
Функция с таким графиком имеет вид . Экспонента в этой формуле обеспечивает экспоненциальное убывание, когда движемся по графику горизонтально отдаляясь от начала координат, в то время как коэффициент вызывает начальный подъем на графике вблизи начала координат.
Модель иногда называют дискретной логистической моделью или моделью Рикера. Такая модель роста популяции, названная в честь её первооткрывателя Билла Рикера, была предложена в далёком 1954 году. Легко вычислить точки равновесия модели, ими являются и . Можно дополнительно проанализировать эту модель, нарисовав паутинную диаграмму и вычислив стабильность равновесий, как делалось неоднократно в предыдущих разделах.
Можно возразить против подхода к моделированию в формате «кролик из шляпы»; без объяснений, откуда взялось уравнение модели Рикера. Но ниже будет дано одно пояснение, важно понимать, что действительно важно, так это то, какие качественные изменения демонстрирует функция на графике, насколько реалистично такое поведение. Если странная формула дает нужный график, то этого уже достаточно для оправдания ей использования.
Для более полного обоснования адекватности модели Рикера вернемся к графику функции изменения численности населения на душу населения как функции от , что в свою очередь стимулировало развитие логистической модели. Единственная причина выбора формулы
заключалась в моделировании нисходящей тенденции, показанной на рисунке 1.1.
Как улучшить такую модель? Во-первых, изменение численности населения на душу населения не может быть меньше −1, потому что это будет означать более одной смерти на душу населения, но «Расстреливать два раза уставы не велят». Это означает, что график должен больше походить на рисунок 1.11.
Рисунок 1.11. Темпы роста на душу населения для новой модели.
Поскольку график выглядит как экспоненциально убывающая кривая, перемещенная вниз на одну единицу, это приводит к следующей формуле: , при некоторых положительных значениях и . Чтобы получить классическую формулу из модели Рикера, выполним замену переменных. Пусть и , тогда с новыми параметрами и модель принимает вид . Теперь элементарными преобразованиями можно прийти к формуле Рикера: . В этой формуле как и прежде следует интерпретировать как пропускную способность или грузоподъёмность логистической модели, потому что если , то ; а если , то . Конечная внутренняя скорость роста, однако, равна , а не просто , хотя для достаточно малого эти величины примерно одинаковы.
Конечно, кривая не обязана быть экспоненциально убывающей. Чтобы точнее смоделировать динамику популяции, нужно собрать данные о том, как численность популяции в момент времени зависит от численности популяции в момент времени . Тогда можно будет построить точки и найти формальное выражение функции, график которой через них проходит. Поскольку модель Рикера имеет два параметра, и , то изменяя каждый из них можно сделать так, чтобы теоретическая кривая достаточно хорошо покрывалась эмпирическими данными.
Другая часто используемая модель имеет вид . Физическое значение чисел , и в этой модели неочевидно, просто уравнения с тремя параметрами позволяют иметь больше свободы в выборе формы кривой и лучше соответствовать эмпирическим данным.
Представленные на рисунке 1.12 графики демонстрируют функциональную зависимость модели при двух различных вариантах значений параметров. Эти два графика описывают совершенно разную динамику населения. График слева, который асимптотически стремится к горизонтальной оси, представляет собой чистую конкуренцию за ресурсы между людьми, где каждый человек просто получает меньше ресурсов, если популяция очень велика. Таким образом, в рамках данной модели члены популяции страдают от наличия большой популяции вокруг. Следовательно, большое значение , вероятно, приведет к гораздо меньшему значению для и чем больше , тем меньше будет .
Рисунок 1.12. Две модели с разными значениями параметров.
График справа, который асимптотически стремится к прямой параллельной горизонтальной оси, представляет собой своеобразное соревнование, по условиям которого если численность популяция превышает свою пропускную способность, то некоторые особи получают все ресурсы, а другие не получают ничего. Поэтому любое большое значение может привести примерно к тому же значению . Конечно, многие популяции демонстрируют поведение, сочетающее аспекты этих двух крайних типов конкуренции, а поэтому описываются графиками промежуточного варианта.
Задачи для самостоятельного решения:
1.4.1. Для дискретной популяционной модели относительный темп роста определяется как .
а. Заполните пропущенные места: начиная с некоторого значения , если относительный темп роста окажется больше 1, то популяция будет _________ в течение следующего временного интервала, тогда как если он будет меньше 1, то популяция _________.
б. Какой смысл имеет относительный темп роста равный нулю? А отрицательный?
в. Приведите выражения для вычисления относительного темпа роста для геометрических и логистических моделей населения, а также других моделей пройденного раздела.
г. Постройте график каждой из относительных скоростей роста, которые выведите в части (в) как функции от . Возможно, придется задать несколько конкретных значений параметров, чтобы нарисовать графики.
1.4.2. На графиках (б), (в) и (г) из задачи 1.2.9 раздела 1.2 видно, что , когда достаточно мало. Объясните влияние этой особенности на динамику популяции. Почему это может оказаться важным с прикладной точки зрения? Обнаруженный эффект иногда называют эффектом Алле.
1.4.3. Постройте простую модель, показывающую эффект Аллее в следующих вариантах.
а. Объясните, что для некоторых параметров , средняя скорость роста , когда или , и , когда . Изобразите возможный график зависимости от .
б. Объясните, почему имеет нужные характеристики.
в. Исследуйте полученную модель, используя программы onepop.m из задачи 1.2.4, cobweb.m и cobweb2.m в MATLAB для некоторых вариантов значений и .
% cobweb.m
%
% Паутинная диаграмма для моделирования одной популяции разностным уравнением.
%
% У пользователя запрашивается уравнение, определяющее модель. Затем по
% щелчку на начальной численности популяции на графике будет отображаться
% «паутина» будущих численностей популяции.
%
p=0; % инициализация популяции
%
disp (' Введите формулу, определяющую модель популяции, используя "p" для')
disp ('обозначения численности: (Например: next_p = p+1.8*p*(1-p/10) )')
next_p=input ('next_p = ','s');
if isempty(next_p) next_p='p+1.8*p*(1-p/10)'; end;
eval( [next_p ';']); % проверяемая формула
%
disp (' ');
disp (' Введите верхний и нижний пределы значения численности в момент времени t,')
disp ('чтобы задать границы изображения на графике:')
limits=input('(По умолчанию [pmin pmax]=[0 20]) [pmin pmax]= ');
if isempty(limits) limits=[0 20]; end;
pinc=(limits(2)-limits(1))/50;
x=limits(1):pinc:limits(2);
%
p=limits(1); y=eval (next_p);
for i=x(2):pinc:limits(2); % цикл для создания вектора значений P
p=i;
p=eval (next_p);
y=[y p];
end;
%
figure % настройка графика
plot(x,y,x,x)
axis([limits(1),limits(2),limits(1),limits(2)]);
xlabel ('P_t');
ylabel ('P_{t+1}');
title (['следующее\_p=',next_p]);
%
continueb=1; % логическое значение продолжения цикла
while continueb % цикл пока кнопку не нажали
[p,q,button]=ginput(1); % получить начальную численность
if button==1
%
plot (x,y,x,x); % построение графика
axis([limits(1),limits(2),limits(1),limits(2)]);
hold on
xlabel ('P_t');
ylabel ('P_{t+1}');
title (['следующее\_p=',next_p]);
%
for i=1:50; % цикл построения секций паутины
w=p;
p=eval (next_p);
plot([w,w],[w,p],'k','EraseMode','none'); % рисуем вертикальный фрагмент
pause(.1);
if p<0; break; end; % фильтрация отрицательных значений P
plot([w,p],[p,p],'k','EraseMode','none'); % рисуем горизонтальный фрагмент
pause(.1);
end;
hold off;
else continueb=0; % конец цикла
end
end
%
% cobweb2.m
%
% Паутинная диаграмма для моделирования одной популяции разностным уравнением.
%
% У пользователя запрашивается уравнение, определяющее модель. Затем по
% щелчку на начальной численности популяции на графике будет отображаться
% «паутина» будущих численностей популяции. Старые линии постепенно стираются
% с течением времени.
%
m=[];
s=16; % количество линий для рисования
p=0; % инициализируем начальное значение
% % численности популяции
disp (' ')
disp (' Введите формулу, определяющую модель популяции, используя "p" для')
disp ('обозначения численности: (По умолчанию: next_p = p+2.5*p*(1-p/10) ) ')
next_p=input ('next_p = ','s');
if isempty(next_p) next_p='p+2.5*p*(1-p/10)';
end;
p=eval (next_p); % проверяем корректна ли формула
%
disp (' ')
disp ('Введите верхний и нижний пределы P в момент времени t, чтобы задать')
disp ('границы изображения графика:')
plimits=input ('(По умолчанию [pmin pmax]=[0 20]) [pmin pmax]= ');
if isempty(plimits) plimits=[0 20]; end;
%
% Формируем данные для построения модели
pinc=(plimits(2)-plimits(1))/20; % устанавливает интервал между
% % соседними значениями
h=[plimits(1):pinc:plimits(2)];
for k=1:21; % цикл создания вектора значений
p=h(k);
p=eval (next_p);
m=[m p];
end;
% начало построения нового графика с изображением функции модели и диагональной линии
figure;
hold on;
axis([plimits plimits]);
curve=plot(h,m,'Color','b');
diag=plot(h,h,'Color','g');
xlabel ('P_t');
ylabel ('P_{t+1}');
title (['следующее\_p=',next_p] );
% создаём вектор фрагментов для ступенек
stephan=ones(1,2*s);
button=1;
% получаем начальное значение численности популяции от пользователя
disp(' ')
disp(' Щелкните левой кнопкой на начальном значении или правой, чтобы выйти.')
[p,x,button]=ginput(1);
%
while(button==1)
%
x=p;
for i=1:s; % цикл для начала создания паутины
p=eval (next_p);
stephan(2*i-1)=plot([x;x],[x;p],'k','EraseMode','background');
pause(.1);
stephan(2*i)=plot([x;p],[p;p],'k','EraseMode','background');
pause(.1);
x=p;
end
%
for i=1:64; % цикл удаления первого элемента
p=eval(next_p); % вычисляем следующий член
delete(stephan(1)) % удаляем вертикальную линию
stephan(1:2*s-1)=stephan(2:2*s); % и указатель на неё
for k=1:2*s-1
set(stephan(k),'EraseMode','background');% перерисовываем линии
end;
set(curve,'Color','b'); % перерисовываем кривые
set(diag,'Color','g');
stephan(2*s)=plot([x;x],[x;p],'k','EraseMode','background');% добавляем линию
pause(.1);
delete(stephan(1)) % стираем горизонтальную линию
stephan(1:2*s-1)=stephan(2:2*s); % и указатель на неё
for k=1:2*s-1
set(stephan(k),'EraseMode','background');% перерисовываем линии
end;
set(curve,'Color','b'); % перерисовываем кривые
set(diag,'Color','g');
stephan(2*s)=plot([x;p],[p;p],'k','EraseMode','background');% добавляем линию
x=p; % сохраняем новую популяцию
pause(.1);
end
% получаем начальную популяцию от пользователя
disp(' ')
disp('Щелкните левой кнопкой на начальной численности или правой, чтобы выйти.')
[p,x,button]=ginput(1);
if (button==1) delete(stephan); end;
%
end
Является ли обнаруженная динамика популяции интуитивно ожидаемой?
г. Какие особенности этого уравнения кажутся нереалистичными? Как можно улучшить модель?
Проектные работы:
1. Исследуйте модель Рикера 1954 года более детально.
Рекомендации
Используйте калькулятор или компьютер для построения графика функции от для разных значений и . Сравните результат с соответствующим графиком логистической модели. Какие обнаруживаются качественные сходства и различия между графиками?
Найдите все точки равновесия модели.
Используйте программу onepop.m в MATLAB из задачи 1.2.4 для исследования динамического поведения этой модели при и различных . Обнаруживается ли стабильное равновесие? А 2-циклы? 4-циклы? Хаотичное поведение?
Используйте программу longterm.m в MATLAB из проектной работы 1.3.1 для создания диаграммы бифуркации этой модели по мере изменения .
2. Повторение из предыдущего проекта для модели , которая часто используется для моделирования популяций в живой природе. Для различных параметров можно сначала зафиксировать , и варьировать положительные значения . Затем зафиксируйте , и варьируйте и так далее.
3. Интересная модель популяции елового почкового червя была предложена Людвигом с соавторами в 1978 году. Исследуйте её. Авторы модели использовали дифференциальное уравнение и предполагали логистический рост популяции почкового червя, но вводили дополнительный параметр для учета влияния хищных птиц на моделируемую численность. Формализовалось явление «хищничества» функцией , где обозначало количество почковых червей, а параметры и могли быть выбраны для изменения графика в соответствии с экспериментальными данными.
Рекомендации
Изобразите график функции и подумайте, чем можно объяснить наблюдаемое количество почковых червей, потребляемых хищными птицами при разных размерах популяции почковых червей. В частности, увеличивается ли численность популяции и стабилизируется ли, как должно быть, согласно интуитивному представлению? Как значения α и β влияют на форму графика?
Изучите полную модель с помощью MATLAB для различных значений параметров, но выбирайте достаточно малым, чтобы избежать циклов или хаоса в логистической части модели. Найдите значения параметров, которые выглядят реалистичными.
Что можно сказать об устойчивых состояниях данной модели и типе их стабильности?
