Цитата(Дедъ @ Пятница, 13 Октября 2023, 13:50)
Содержать ВСЕ хорошие там не сказано. Главное - возлюбить пчелу
Здесь на многие вопросы дискуссии есть ответ:
Математическая биология и биоинформатика
2017. Т. 12. № 1. С. 224–236. doi: 10.17537/2017.12.224
================== МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ =================
УДК 51-76
Математическое моделирование динамики численности биологической популяции при
изменяющихся внешних условиях на примере
бурзянской бортевой пчелы (Apis mellifera L., 1758)
Ибрагимова Л.С. 1, Юмагулов М.Г. 2, Ишбирдин А.Р. 2, Ишмуратова М.М. 2
1Башкирский государственный аграрный университет, Уфа, Россия
2Башкирский государственный университет, Уфа, Россия
Аннотация. В статье рассмотрены некоторые вопросы использования неавтономных дискретных моделей для описания динамики численности биологической популяции. Обсуждаются вопросы определения коэффициентов модели в ситуации, когда эти параметры зависят от времени. В качестве приложения рассмотрена задача математического моделирования динамики численности семей природной популяции бурзянской бортевой пчелы, обитающей на территории Республики Башкортостан. Для описания численности пчелы в качестве базовой модели использовалось неавтономное уравнение Рикера. Предложена соответствующая дискретная математическая модель, позволяющая прогнозировать динамику численности бурзянской бортевой пчелы с учетом изменяющихся природных факторов. Разработаны соответствующие программы в средах MatLab и Mathcad. Приведены некоторые численные и графические результаты исследования предложенной модели.
Ключевые слова: математическая модель, биологическая популяция, модель Рикера, неавтономное уравнение, бортевая пчела.
1. ВВЕДЕНИЕ
В задачах математического моделирования динамики численности биологических популяций с неперекрывающимися поколениями, с четко выраженной периодичностью жизненных циклов, с небольшой численностью естественным является использование дискретных моделей, описываемых разностным уравнением (см., например, [1–5])
xn+1 = f x( n, ), n = 0,1,2,... (1)
В этом уравнении xn обозначает численность популяции в момент времени t = n, при этом единицей времени может быть, например, один год; μ – коэффициенты модели.
Уравнение (1) является автономным в том смысле, что функция f x( ,) и коэффициенты μ не зависят от времени n. Поэтому в нем численность популяции при всех t +n 1 полностью определяется только ее численностью в момент t = n. Такое предположение естественно, если считать, что внешние и внутренние условия для популяции постоянны во времени.
В действительности эти условия могут сильно изменяться и, как следствие, существенно сказываться на динамике популяции. В такой ситуации более естественным является использование неавтономного разностного уравнения
xn+1 = f x( n, ( )), n n = 0,1,2,..., (2)
в котором μ(n) – коэффициенты модели, изменяющиеся во времени.
В большинстве работ по популяционной динамике принимается предположение о постоянстве внешних и внутренних условий и, как следствие, используются автономные модели. Это часто приводит к расхождению моделирования и данных наблюдений, что связано, в первую очередь, с влиянием внешних изменяющихся факторов. В то же время задачи исследования динамики популяций при изменяющихся условиях изучены совершенно недостаточно, хотя и здесь получен ряд нетривиальных результатов (см., например, [6–9]). Это связано, в первую очередь, с возникающими при этом серьезными математическими трудностями.
При построении модели (2) одной из основных задач является задача определения характера зависимости коэффициентов μ от времени n. В настоящей работе изучается вопрос определения коэффициентов модели при изменяющихся внешних условиях на примере задачи моделирования динамики численности семей бурзянской бортевой пчелы, обитающей на территории заповедника «Шульган-Таш» (Республика Башкортостан).
2. ЗАДАЧА ОПРЕДЕЛЕНИЯ КОЭФФИЦИЕНТОВ МОДЕЛИ ПРИ ИЗМЕНЯЮЩИХСЯ УСЛОВИЯХ
Динамика автономной системы (1) определяется функцией f x( ,) . В математической экологии (см., например, [1–5, 10–12]) используются различные виды функций; большую популярность получили, в частности, модели, в которых функция f x( ,) имеет вид: f x( , =) x (модель Мальтуса), f x( , = −) x(1 x) (логистическая модель), f x( , =) xexp(r1− Kx ) (модель Рикера); здесь μ, r, K – некоторые коэффициенты; содержательный смысл этих коэффициентов раскрывается ниже на примере модели Рикера.
Использование этих моделей или их модификаций для описания динамики численности соответствующей биологической популяции требует определения соответствующих коэффициентов. Обычно эти коэффициенты определяются путем выделения в реальной системе характеристик, оказывающих основное влияние на динамику численности популяции; если таких характеристик много, то применяются те или иные варианты усреднения и т.п. операций для того, чтобы уменьшить количество базовых параметров модели. Для определения численных значений коэффициентов или диапазонов их изменения обычно используются методы параметрической идентификации (см., например, [13, 14]); здесь важное значение имеют экспериментальные данные.
Рассмотрим теперь неавтономную систему (2). Поведение решений этой системы определяется не только свойствами функции f x( ,) , но и характером зависимости коэффициента μ от n. Чем более сложной является эта зависимость, тем более сложно ведут себя решения системы (2).
Важный частный случай представляет ситуация, когда коэффициент μ(n) является периодической функцией времени n. Модели с периодически изменяющимся параметром часто применяются для решения биологически содержательных задач (см., например, [15–18]) для изучения влияния периодических воздействий на динамику численности популяций. Отметим также, что такие модели могут быть сведены к автономным системам.
В случае, когда коэффициент μ не является периодической функцией, задача исследования системы (2) усложняется. Здесь следует выделить две принципиально различные ситуации. Первая связана с предположением, что характер зависимости коэффициента μ от n известен для всех n. В этой ситуации модель (2) позволяет провести качественное и приближенное исследование динамики численности популяции и, в частности, прогнозировать изменение этой численности.
Вторая ситуация связана с предположением, что значение коэффициента μ известно лишь для данного момента n времени (или до данного момента n времени). В этой ситуации модель (2) позволяет определить численность популяции лишь в следующий момент времени n + 1. Тем не менее, и в этой ситуации, модель (2) полезна, так как она позволяет рассматривать различные сценарии изменения численности популяции в зависимости от возможных вариантов изменения внешних и внутренних условий.
Как и в случае автономных моделей, использование неавтономной модели вида (2) для описания динамики численности реальной биологической популяции требует определения соответствующих коэффициентов и (дополнительно) характера зависимости этих коэффициентов от времени. Здесь также коэффициенты определяются путем выделения в реальной системе характеристик, оказывающих основное влияние на динамику численности популяции. Так как в неавтономной системе коэффициенты меняются во времени, то для определения таких характеристик можно использовать методы регрессионного и корреляционного анализа (см., например, [19]), метод наименьших квадратов и т.д. Для определения численных значений коэффициентов или диапазонов их изменения здесь также эффективны методы параметрической идентификации, использование которых требует наличия достаточно большой базы экспериментальных данных.
Рассмотрим задачу определения коэффициентов μ модели вида (2) c выбранной функцией f x( ,) в следующей постановке. Пусть коэффициент μ является скалярной или векторной функцией времени n. Требуется определить значение функции μ(n) в момент n0 так, чтобы это значение учитывало бы реальные данные о внешних и внутренних условиях до момента n0 . Так как эти условия с течением времени изменяются (возможно, с некоторой цикличностью), то функция μ(n) будет непостоянной, т.е. получим неавтономное уравнение вида (2).
Решение указанной задачи при выбранной функции f x( ,) (и, соответственно, модели вида (2)) предлагается провести по следующей схеме.
1. На первом этапе среди многообразия реальных данных об изменяющихся внешних и внутренних условиях выделяются те данные 1, 2,..,m , которые оказывают наиболее существенное влияние на динамику системы.
2. Значение функции μ = μ(n) в момент времени n = n0 определяется как линейная комбинация (n0) =1 1(n0)+ 2 2(n0)+ +... m m(n0). При этом коэффициенты 1, 2,..,m (не зависящие от времени n) могут быть найдены, например, методом наименьших квадратов путем решения аналогичной задачи для предыдущих моментов времени n = −n0 k0 , n = − +n0 k0 1, n = n0 −1 при некотором k0 1.
Естественно, полученную модель следует затем протестировать, например, методами регрессионного и корреляционного анализа. Модель будет зависеть и от выбора функции f x( ,) в (2); поэтому естественно проанализировать ее свойства с учетом различных вариантов выбора функции f x( ,) . Указанный подход применяется ниже в задаче моделирования динамики численности бурзянской бортевой пчелы.
3. МОДЕЛИРОВАНИЕ ДИНАМИКИ ЧИСЛЕННОСТИ БУРЗЯНСКОЙ БОРТЕВОЙ ПЧЕЛЫ
В качестве приложения рассмотрим задачу математического моделирования динамики численности бурзянской бортевой пчелы, обитающей на территории заповедника «Шульган-Таш» (см. [20–22]). Наблюдения за численностью и состоянием популяции этой пчелы проводятся регулярно с 1960 года, со времени учреждения Прибельского филиала Башкирского заповедника (в данное время – заповедник «Шульган-Таш»). Учеты численности семей пчел проводятся во время весенних и осенних ревизий. Для характеристики погодно-климатических условий используются данные наблюдений метеостанции Башкирского заповедника.
На рисунке 1 представлен график изменения осенней численности бортевых семей за период 1960–2012 гг. (данные «Летописи природы» заповедника «Шульган-Таш»).
Рис. 1. Динамика численности бурзянской пчелы за период 1960–2012 гг. По горизонтали указаны годы наблюдений, а по вертикали – численность бортевых пчелиных семей.
Этот график наглядно иллюстрирует тот факт, что динамика численности бурзянской бортевой пчелы достаточно сложна, она характеризуется как периодами стремительного возрастания численности и, не менее стремительного, ее убывания, так и периодами стабильной численности. При этом можно говорить и о циклических явлениях в изменении численности; особо проявляется цикл с периодом примерно 11– 12 лет (см. [20, 21]). Естественно, говоря о циклах, мы имеем в виду лишь весьма приблизительную повторяемость.
Возникает естественный вопрос о том, какие причины и факторы являются основными с точки зрения влияния на изменение численности бурзянской бортевой пчелы. Эффективным аппаратом в изучении таких вопросов, как с качественной точки зрения, так и с точки зрения прогнозирования динамики численности биологических популяций, является разработка и анализ соответствующей математической модели.
3.1. Выбор математической модели
Бурзянскую бортевую пчелу можно рассматривать как популяцию с неперекрывающимися поколениями, с четко выраженной периодичностью жизненных циклов и с небольшой численностью (изменяющейся в пределах 50–190 семей). На численность пчелы сильно влияют внешние факторы, в первую очередь, погодные условия, которые нельзя назвать постоянными из года в год. Поэтому для описания динамики ее численности естественно использовать неавтономное разностное уравнение вида (2).
В настоящей работе в качестве базовой выбрана модель Рикера [12]
xn+1 =xn expr1− xKn . (3)
В этой модели положительный коэффициент r характеризует скорость роста популяции в отсутствии лимитирования, а параметр K показывает емкость экологической ниши данной популяции. Коэффициент K обычно связывают с максимальной численностью популяции xmax формулой
K = rxmax exp(1−r). (4)
Так как коэффициенты K и r связаны равенством (3), то фактически модель Рикера зависит только от одного коэффициента, например, от r. Приведенные в [20, 21] данные о динамике численности популяции бурзянской пчелы (см. также рис. 1) говорят о том, что в нашей задаче можно принять xmax =190 .
Замечание 1. Как было отмечено выше, основной целью настоящей статьи является изложение подхода, позволяющего определить коэффициенты дискретной модели (2) при изменяющихся внешних условиях по имеющейся информации о реальных данных о внешних и внутренних условиях до момента n0 . В этой связи выбор модели
Рикера в качестве базовой носит, в первую очередь, иллюстративный характер. Вместе с тем, выбор этой модели определялся также предоставляемой ею богатым многообразием вариантов качественного поведения дискретной системы.
Уравнение (3) является автономным. Оно имеет две точки равновесия: x = 0 и x = K. Первая из них неустойчива при всех значениях r 0, а вторая – устойчива при 0 r 2 и неустойчива при r 2. Уравнение (3) при r 2 может иметь циклы различных периодов (в частности, периодов 2, 4, 8, 16 и т.д.), некоторые из которых могут быть устойчивыми, а остальные неустойчивы. Это уравнение при r 14.77 может иметь решения самой различной сложности, в поведении которых не наблюдается никакой закономерности (см., например, [1, 12, 15]).
Уравнение (3) не учитывает тот факт, что в реальности внешние и внутренние условия изменяются во времени (часто непредсказуемым образом) и могут сильно изменить динамику системы. Для описания динамики численности пчел с учетом изменяющихся внешних и внутренних условий предлагается от автономной модели (3) перейти к неавтономной модели вида
xn+1 = xn exp(rn 1− xKn ) , n =0,1, 2,... , (5)
в которой коэффициент rn зависит от времени n. Здесь значение K в соответствии с (4) определяется равенством:
K rx= max exp(1−r), (6)
в котором число r определяется экспериментальным путем, обеспечивающим наилучшее среднеквадратичное приближение модельных значений к реальным данным (более детально об этом говорится ниже). В настоящей работе будет изучаться модель вида (5) в предположении, что на динамику численности пчел влияют только погодные и природные условия.
Для определения коэффициента rn в модели (5) воспользуемся приведенными результатами наблюдений в период с 1960 по 2012 гг. за рядом климатических и природных параметров, которые в большей или меньшей степени могут влиять на динамику численности пчел [20–22]. Здесь мы ограничимся исследованием влияния таких факторов как среднемесячные температуры воздуха, месячные суммы осадков и среднемесячные значения чисел Вольфа (показатель солнечной активности), а также среднегодовые температуры.
Как было отмечено выше, можно говорить о некоторой цикличности в изменении численности пчел, в частности, о наличии определенной цикличности с периодом примерно 11–12 лет. Отметим, что среди выбранных факторов числа Вольфа также имеют цикличность в 11–12 лет (закон Швабе-Вольфа [23]). Для определенности в качестве периода возьмем период T =12.
Замечание 2. Цикличность в динамике численности бурзянской бортевой пчелы с периодом примерно 11–12 лет отмечалась также в работах [20, 21]. Конечно, с целью получения более адекватной модели следует поварьировать длину периода, выбрав затем наиболее подходящий период. Однако, в настоящей статье мы для простоты ограничиваемся рассмотрением только периода 12 лет; именно для этого периода иллюстрируется подход определения параметров неавтономной модели (5).
Рассматриваемый интервал времени 1960–2012 гг. разобьем на циклы по 12 лет: 1960–1971 гг., 1972–1983 гг. и т.д. Примем следующие обозначения.
1) Пусть индекс n =1,2,...,12 обозначает номер года из выбранного цикла, индекс j =1,2,...,12 обозначает номер месяца выбранного года.
2) Пусть tn( )j , pn( )j , wn( )j и stn – это значения среднемесячной температуры, месячной суммы осадков, среднемесячное значение числа Вольфа в j -месяце n -го года и среднегодовая температура n -го года соответственно.
3) Положим tmin ( )j = min tn( )j , pmin ( )j = min pn ( )j , wmin ( )j = min wn ( )j ,
1 n 12 1 n 12 1 n 12
tmax ( )j = 1max ( ) n 12tn j , pmax ( )j = max pn ( )j , wmax( )j =1max wn( )j , stmin = 1min n 12 stn,
1 n 12 n 12
stmax = maxstn , t0( )j = t1( )j +t2( )j + +... t12( )j и st0 = st1 + + +st2 ... st12 . Естественно
1 n 12 12 12
считать, что tmin ( )j tmax ( )j , pmin ( )j pmax ( )j , wmin ( )j wmax ( )j и stmin stmax . 4) Определим безразмерные величины:
Tn( )j = 2 tn( )j −t0( )j , P jn( ) = pn( )j − pmin ( )j , tmax ( )j −tmin ( )j pmax ( )j − pmin ( )j
Wn( )j = wn( )j − wmin ( )j , STn = 2 stn − st0 . (7) wmax ( )j − wmin ( )j stmax − stmin
Тогда 0 P jn( ) 1 и 0 Wn( )j 1, а величины Tn( )j удовлетворяют оценкам: a Tn( )j b, где числа a и b таковы, что −2 a b 2 и b a− = 2 (аналогичным оценкам удовлетворяют величины STn ). Отметим также, что в случае нормального или равномерного распределения значений температур имеем a =−1 и b =1, т.е. −1Tn( )j 1. Таким образом, для года с номером n получим 37 величин вида (7). Все эти величины известны.
Естественным является предположение о том, что чем больше отклонения какойлибо из величин (7) от нуля, тем больше соответствующий фактор (температура, осадки или числа Вольфа) влияет на жизнедеятельность пчел и, соответственно, изменяет скорость роста популяции, т.е. коэффициент r в модели Рикера (3). Влияние значений величин (7) на коэффициент r в году с номером n определим линейной функцией вида:
12 12 12
rn = +r n( )j Tn( )j +n( )j P jn( )+n( )j Wn( )j + n STn , (8)
j=1 j=1 j=1
где суммирование ведется по месяцам, n( )j , n( )j , n( )j , n и r – некоторые коэффициенты, о вычислении которых будет указано ниже. Если эти коэффициенты известны, то неавтономная модель (5) позволяет спрогнозировать численность пчел в году с номером n + 1.
Замечание 3. На динамику численности пчел существенное влияние могут оказать не только природные факторы, но и ряд других факторов – таких как, например, болезни пчел, конкуренция, человеческий фактор и т.д. Изучение влияния таких факторов может составить предмет полноценного исследования. В настоящей статье мы ограничиваемся рассмотрением частной задачи – исследованием воздействия внешних изменяющихся природных факторов на динамику численности биологической популяции.
3.2. Расчет параметров модели
Таким образом, возникает задача определения 37 коэффициентов n( )j , n( )j , n( )j и n , соответствующих величинам (7). Определить их можно, используя имеющуюся информацию за предыдущие годы наблюдений. Однако, учет влияния всех 37 факторов на динамику численности популяции пчел приведет к весьма сложной задаче (в том числе с вычислительной точки зрения). В связи с этим отметим, что анализ корреляционной зависимости между величинами (7) и численностью пчелиных семей показывает, что только 6 из 37 величин (7) имеют более или менее выраженную корреляцию. Такими величинами являются среднемесячные температуры марта, апреля и июня, осадки февраля, числа Вольфа за май, а также среднегодовые температуры. Соответствующие коэффициенты корреляции также невелики, но для других факторов они практически нулевые. Например, коэффициент корреляции за все годы наблюдений между среднемесячной температурой апреля и численностью пчелиных семей равен 0.272, а для среднемесячной температуры июня он равен 0.016. Таким образом, анализируя все данные, перейдем к более простой модели, в которой будем учитывать только указанные шесть факторов.
Для простоты шесть соответствующих величин из (7) обозначим как
(1)n , (2)n , (3)n ,...,(6)n , где (1)n ,(2)n ,(3)n – соответствуют среднемесячным значениям температуры марта, апреля и июня, (4)n – среднемесячному значению чисел Вольфа за май, (5)n – сумме осадков февраля, (6)n – среднегодовой температуре. Обозначим далее через 1, 2, 3, 4, 5, 6 – соответствующие коэффициенты формулы (8), так что эта формула примет вид
rn = +r + + + + +1 (1)n 2 (2)n 3 (3)n 4 (4)n 5 (5)n 6 (6)n , (9)
где индекс n соответствует году наблюдения. Подставляя это равенство в (6), получим неавтономную модель, используя которую можно прогнозировать численность пчел с учетом изменяющихся погодных и природных условий.
Определение неизвестных коэффициентов r, 1, 2, 3, 4, 5, 6 из формулы (9) проводилось по следующей схеме. Задается значение r из интервала r (0;2) и из формулы (6) определяется значение K (в предположении, что xmax =190 ). Затем методом наименьших квадратов определяются значения 1, 2, 3, 4, 5, 6 (об этом ниже). Вычисления показали, что наилучшее среднеквадратичное приближение модельных значений (по формуле (5)) к эмпирическим данным достигается при r =1.6, в этом случае K =167. Полученное значение K является устойчивой точкой равновесия автономной модели Рикера (3) при r =r и K = K. Числа r =1.6 и K =167 можно трактовать как значения коэффициентов автономной модели (3), описывающей динамику численности бурзянской бортевой пчелы в оптимальных природных и климатических условиях, не изменяющихся во времени. В реальности же внешние и внутренние условия изменяются во времени и могут сильно изменить динамику системы.
По модели (3) определим функцию f x r( , ) = xexp(r1− Kx ) . Пусть известны значение численности xn пчелиных семей и коэффициент rn . Тогда значения xn+1 = f x r( n n, ) дадут нам (в соответствии с (5)) модельные значения численности в год
с номером n + 1.
Приведем теперь схему вычисления коэффициента rn , определенного формулой (9). Для этого вычислим значения коэффициентов 1, 2, 3, 4, 5, 6 на основе метода наименьших квадратов для указанных значений r и K . Пусть x x1, 2,..,x12 – это численность пчелиных семей в годы с номерами n =1,2,...,12. Численность xn определяется на начало года n; значения xn нам известны. Метод наименьших квадратов в рассматриваемой постановке состоит в подборе коэффициентов 1, 2, 3, 4, 5, 6 таким образом, чтобы минимизировать сумму:
( (f x r1, )1 − x2)2 +( (f x r2, 2)− x3)2 + +... ( (f x11,r11)− x12)2 → min, (10)
в которой rn определяется формулой (9).
К сожалению, задача (10) приводит к нелинейным уравнениям для неизвестных коэффициентов 1, 2, 3, 4, 5, 6 . Для того, чтобы получить линейную задачу, прологарифмируем функцию f x r( , ) :
g x r( , ) = ln f x r( , ) = ln x + r1− Kx ,
и положим
y1 = ln x1, y2 = ln x2, …, y12 = ln x12 .
Тогда получим равносильную (10) задачу:
( (g x r1, )1 − y2)2 +( (g x r2, 2)− y3)2 + +... ( (g x11,r11)− y12)2 → min ,
которая приводит к системе линейных алгебраических уравнений (11)
Ax =b
относительно неизвестных 1, 2, 3, 4, 5, 6 ; здесь: (12)
1 1
A=(aij ), x =...2 , b =−...2 ,
6 6
2 2 2 (13)
( )i ( )j x ( )i ( )j x ( )i ( )j x
aij = 1 1 1− K1 + 2 2 1− K2 + +... 11 11 1− K11 , (14)
=j y1 − y2 +r1− Kx1 1− Kx1 +1( )j y2 − y3 + r1− xK2 1− xK2 +( )2j ...
(15)
+ y11 − y12 +r1− xK11 1− xK11 11( )j .
Решая систему (12), найдем значения 1, 2, 3, 4, 5, 6 и, следовательно, получим искомую неавтономную модель (5).
3.3. Алгоритм использования модели
Приведем алгоритм использования предложенной модели (5) для описания динамики численности пчелиных семей в следующей постановке.
Пусть индекс n соответствует последовательным годам какого-либо выбранного промежутка времени. Пусть известна следующая информация:
– Численность x x1, 2,..,x12 пчелиных семей в годы с номерами n=1,2,...,12.
Значения xn определяются по состоянию на начало года с номером n.
– Климатические и природные параметры tn(3) , tn(4), tn(6), wn(5), pn(2) , stn в годы с номерами n =1,2,...,12. Здесь числа tn(3) , tn(4), tn(6) – это среднемесячные значения температуры марта, апреля и июня, wn(5) – среднемесячное значение чисел Вольфа за май, pn(2) – сумма осадков февраля, stn – среднегодовая температура. По этой информации требуется определить численность xn пчелиных семей в год с номером n = 13.
Предлагается следующий алгоритм решения этой задачи.
1. Для n =1,2,...,11 определим числа =(1)n Tn(3) , =(2)n Tn(4) , =(3)n Tn(6),
=(4)n Wn(5) , =(5)n Pn(2) и =(6)n STn, где Tn(3) , Tn(4), Tn(6), Wn(5), Pn(2) и
STn вычисляются в соответствии с формулами (7).
2. Используя числа (1)n , (2)n ,…,(6)n (n =1,2,...,11) и x x1, 2,..,x12 по формулам (11),
(14) и (15) определим элементы матрицы A и вектора b в соответствии с формулами
(13). При этом будем полагать: r =1.6 и K =167.
3. Найдем решение системы (12) и тем самым определим коэффициенты
1, 2, 3, 4, 5, 6 .
4. По формуле (9) определим число:
r12 = +r + + + + +1 12(1) 2 12(2) 3 12(3) 4 12(4) 5 12(5) 6 12(6) .
5. По формуле (5) найдем искомое значение x13 = x12exp(r121− xK12 ) .
Предложенная схема реализована в средах MatLab и Mathcad.
3.4. Численные результаты
Приведем для иллюстрации некоторые из полученных результатов. Данные, взятые в качестве исходных, приведены в таблице 1. Полученные результаты представлены на рисунке 2.
Таблица 1. Данные о численности пчелиных семей и погодных условиях за период 1961– 1971 гг. Здесь X – число пчелиных семей, tn(3), tn(4), tn(6) – температура марта, апреля и июня соответственно, wn(5) – числа Вольфа за май, pn(2) – осадки февраля, stn – среднегодовая температура
Год 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971
n 1 2 3 4 5 6 7 8 9 10 11
X 75 99 113 127 98 49 75 98 88 89 93
tn(3) –4.6 –3 –12 –11.1 –5.7 –5.8 –6.7 –7.2 –13.1 –7.1 –11.1
tn(4) 1.5 3.8 0.8 –1.4 –0.2 2.5 5.5 0.5 1.8 3 –0.3
tn(6) 15 13.2 14.5 13.4 14 12.8 14.7 12.4 13 12.2 13.3
wn(5) 51 43.7 43 9.5 24.1 45.3 96.5 127.2 120 127.5 57.5
pn(2) 7.5 11.5 18.1 10.3 33.6 75.3 4.4 30.5 51.2 32.7 29.8
stn 1.22 1.51 0.51 –0.33 0.77 0.1 1.06 –0.83 –2.72 –0.52 0.46
Определим коэффициенты 1, 2, 3, 4, 5, 6 . После вычислений по указанной выше схеме получим:
=1 0.8647, =−2 0.4906, =3 0.6033, =−4 2.4222, =−5 2.47 , =−6 0.0048.
Если теперь провести расчет по модели (5), то получим динамику численности пчел за эти годы, приведенные в таблице 2.
Таблица 2. Динамика численности пчелиных семей за период 1961–1971 гг. Здесь XX – результаты численного расчета по модели (5), а X – реальные данные
Год 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971
n 1 2 3 4 5 6 7 8 9 10 11
XX 78 88 111 127 110 50 79 87 87 102 88
X 75 99 113 127 98 49 75 98 88 89 93
Рис. 2. Динамика численности пчелиных семей за период 1960–1971 гг. Красным цветом представлен график изменения численности, полученный по предложенной модели, черным цветом – график, характеризующий реальные данные.
Проведя аналогичные вычисления по предложенной модели для всех годов наблюдений, получим численные результаты, графическое изображение которых приведено на рисунке 3.
Рис. 3. Динамика численности пчелиных семей за период 1960–2012 гг. Красным цветом представлен график изменения численности, полученный по предложенной модели, черным цветом – график, характеризующий реальные данные.
Для изучения вопроса о пригодности предложенной модели (5) в задаче аппроксимации данных были проведены соответствующие расчеты, показавшие, в частности, следующее. Среднее арифметическое отклонений эмпирических значений от модельных (деленное на 100, т.е. на среднее число пчелиных семей) равно – 0.01596. Средняя ошибка аппроксимации составляет 11.48 %. Наличие в последовательности отклонений сериальной корреляции проверялось с помощью критерия Дарбина– Уотсона [19]. Критические значения критерия Дарбина–Уотсона для 5 % уровня значимости (при выборке n=50 ) dL =1.291 и dU =1.822 . Для нашего ряда D=1.982; следовательно, сериальная корреляция в последовательности остатков отсутствует. Таким образом, можно говорить о том, что предложенная модель позволяет получить удовлетворительную аппроксимацию данных.
4. ЗАКЛЮЧЕНИЕ
В статье рассмотрены некоторые вопросы использования неавтономных разностных моделей вида (2) для описания динамики численности биологической популяции. Обсуждались вопросы определения коэффициентов модели в ситуации, когда эти коэффициенты зависят от времени. Предложен подход, позволяющий определять эти коэффициенты так, чтобы их значения учитывали бы реальные данные о внешних и внутренних условиях.
В качестве приложения рассмотрена задача математического моделирования динамики численности бурзянской бортевой пчелы, обитающей на территории Республики Башкортостан. В качестве базовой модели использовалось неавтономное уравнение Рикера. Предложена соответствующая дискретная математическая модель, позволяющая прогнозировать динамику численности бурзянской бортевой пчелы с учетом изменяющихся природных факторов. Разработаны соответствующие программы в средах MatLab и Mathcad. Приведен ряд численных и графических результатов, полученных с использованием предложенной модели. Эти результаты убедительно демонстрируют тот факт, что на динамику численности бурзянской бортевой пчелы существенное влияние оказывает совокупность природных факторов. Здесь особо значимы сумма осадков февраля (в частности, увеличение осадков отрицательно влияет на численность пчел) и значения температур марта, апреля и июня.
Предложенная модель может использоваться не только для прогнозирования численности пчел в конкретном году. Она может также использоваться для исследования вопроса о возможных сценариях изменения этой численности в зависимости от вариантов изменения внешних, в первую очередь, климатических условий, в частности, для изучения влияния глобального потепления климата. СПИСОК ЛИТЕРАТУРЫ
1. Братусь А.С., Новожилов А.С., Платонов А.П. Динамические системы и модели биологии. М.: Физматлит, 2010. 400 с.
2. Ризниченко Г.Ю. Лекции по математическим моделям в биологии. Часть 1. Ижевск: НИЦ «Регулярная и хаотическая динамика», 2002. 232 с.
3. Свирежев Ю.М., Логофет Д.О. Устойчивость биологических сообществ. М.:
Наука, 1978. 352 с.
4. Плюснина Т.Ю., Фурсова П.В., Тёрлова Л.Д., Ризниченко Г.Ю. Математические модели в биологии. М.–Ижевск: НИЦ «Регулярная и хаотическая динамика», 2014. 136 с.
5. Недорезов Л.В. Анализ динамики численности сосновой пяденицы с помощью дискретных математических моделей. Математическая биология и
биоинформатика. 2010. Т. 5. № 2. С. 114–123. doi: 10.17537/2010.5.114
6. Сидорин А.П. Поведение популяции с несколькими стационарными состояниями в случайной среде. В: Математические модели в экологии и генетике. М.: Наука. 1981. С. 71–74.
7. Неверова Г.П., Фрисман Е.Я., Математическое моделирование динамики локальных однородных популяций с учетом эффектов запаздывания. Математическая биология и биоинформатика. 2015. Т. 10. № 2. С. 309–324.
doi:
10.17537/2015.10.309
Фрисман Е.Я., Неверова Г.П., Кулаков М.П., Жигальский О.А. Смена
8. динамических режимов в популяциях видов с коротким жизненным циклом: результаты аналитического и численного исследования. Математическая биология и биоинформатика. 2014. Т. 9. № 2. С. 414–429. doi: 10.17537/2014.9.414
9. Li J. Discrete-time models with mosquitoes carrying genetically-modified bacteria. Mathematical Biosciences. 2012. V. 240. № 1. Р. 35–44.
10. Базыкин А.Д. Нелинейная динамика взаимодействующих популяций. МоскваИжевск: Институт компьютерных исследований, 2003. 368 с.
11. Юмагулов М.Г. Введение в теорию динамических систем. СПб.: Лань, 2015.
272 с.
12. Ricker W.E. Stock and recruitment. J. Fish. Res. Board of Canada. 1954. V. 11. № 5. P. 559–623. 13. Дейч A.M. Методы идентификации динамических объектов. М.: Энергия, 1979.
240 с.
14. Льюнг Л. Идентификация систем. Теория для пользователя. М.: Наука, 1991.
432 с.
15. Шлюфман К.В., Неверова Г.П., Фрисман Е.Я. 2-циклы уравнения Рикера с периодически изменяющимся мальтузианским параметром: устойчивость и мультистабильность. Нелинейная динамика. 2016. Т. 12. № 4. С. 553–565. doi: 10.20537/nd1604001
16. Ашихмина Е.В., Израильский Ю.Г., Фрисман Е.Я. Динамическое поведение модели Рикера при циклическом изменении одного из параметров. Вестник ДВО РАН. 2004. № 5. С. 19–28.
17. Sacker R.J., von Bremen H.F. A conjecture on the stability of the periodic solutions of Ricker’s equation with periodic parameters. Appl. Math. Comput. 2010. V. 217. № 3. P. 1213–1219.
18. Шлюфман К.В., Фишман Б.Е., Фрисман Е.Я. Особенности динамических режимов одномерной модели Рикера. Известия высших учебных заведений. Прикладная нелинейная динамика. 2012. Т. 20. № 2. С. 12–28.
19. Дрейпер Н., Смит Г. Прикладной регрессионный анализ. М.: Финансы и статистика. 1986. Т. 1. 366 с.
20. Косарев М.Н., Юмагужин Ф.Г., Нугуманов Р.Г. О динамике численности семей пчел башкирской популяции, заселяемости бортей и колодных ульев в государственном природном заповеднике «Шульган-Таш». В: Использование биологически активных продуктов пчеловодства в животноводстве и ветеринарной медицине. Сб. научных трудов. Москва-Уфа: ВГНКИ. БГАУ, 1999. С. 118–121.
21. Шарипов А.Я., Ишбирдин А.Р. Популяционная динамика бортевой пчелы (Apis mellifera mellifera L.) в заповеднике «Шульган-Таш» за полвека наблюдений. В:
Проблемы популяционной экологии: Материалы всероссийского семинара «Гомеостатические механизмы биологических систем». Под ред. Розенберга Г.С. Тольятти: Кассандра, 2015. С. 336–341.
22. Шарипов А.Я. Влияние изменений климата на состояние бурзянских бортевых пчел. Вестник Оренбургского государственного университета. 2010. № 12. С. 78–81.
23. Витинский Ю.И., Копецкий М., Куклин Г.В. Статистика пятнообразовательной деятельности Солнца. М.: Наука, 1986. 296 c.
Рукопись поступила в редакцию 16.04.2017, переработанный вариант поступил 20.06.2017. Дата опубликования 30.06.2017.
Однажды пробовал с 3-х сторон зайти, и копал даже, но в итоге разбросал от бессилия поленницу и ушёл.
За 5 лет использования пока ни одной потери(тьфу-тьфу-тьфу)
