УПРАВЛЯЕМЫЙ ПЕРЕХОД В МОДЕЛИ ДИНАМИКИ БИОМАСС КОРНЕВЫХ РАСТЕНИЙ

5. Liou K.N., Ои S. С. The role of cloud microphysical processes in climate: An assessment from a one-dimensional perspective //J. Geophys. Res.: Atmospheres. 1989. 94, D6. 8599-8607.

6. Renno N. 0., Emanuel K.A., Stone P.H. Radiative-convective model with an explicit hydrologic cycle: 1. Formulation and sensitivity to model parameters //J. Geophys. Res.: Atmospheres. 1994. 99, D7. 14429-14441.

7. Renno N. 0., Stone P.H., Emanuel K.A. Radiative-convective model with an explicit hydrologic cycle: 2. Sensitivity to large changes in solar forcing //J. Geophys. Res.: Atmospheres. 1994. 99, D8. 17001-17020.

8. Kelly M.A., Randall D.A., Stephens G.L. A simple radiative-convective model with a hydrological cycle and interactive clouds // Quart. J. Roy. Meteorol. Soc. 1999. 125, N 555. 837-869.

9. Larson V.E. Stability properties of and scaling laws for a dry radiative-convective atmosphere // Quart. J. Roy. Meteorol. Soc. 2000. 126, N 562. 145-171.

10. Larson V.E. The effects of thermal radiation on dry convective instability // Dynamics Atmos. and Oceans. 2001. 34, N 1. 45-71.

11. Emanuel K., Wing A.A., Vincent E.M. Radiative-convective instability //J. Adv. Model. Earth Systems. 2014. 6, N 1. 75-90.

12. Fuchs Z., Raymond D.J. Large-scale modes in a rotating atmosphere with radiative-convective instability and WISH К // J. Atmos. Sci. 2005. 62, N 11. 4084-4094.

13. Beucler Т., Cronin Т., Emanuel K. A linear response framework for radiative-convective instability //J. Adv. Model. Earth Systems. 2018. 10, N 8. 1924-1951.

14. Нигматулин P.И. Уравнения гндро-н термодинамики атмосферы при малых силах инерции по сравнению с силой тяжести // Прикл. матем. и механ. 2018. 82, № 4. 472-484.

15. Монин А. С. Теоретические основы геофизической гидродинамики. М.: Гидрометеоиздат, 1988.

16. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. Т. VI. Гидродинамика. М.: Наука, 1986.

17. Holton J.R., Hakim G.J. An Introduction to Dynamic Meteorology. 5th ed. Waltham: Academic Press, 2012.

18. Марчук Г.И. Численное решение задач динамики атмосферы и океана. М.: Гидрометеоиздат, 1974.

Поступила в редакцию 21.02.2020

УДК 531.396

УПРАВЛЯЕМЫЙ ПЕРЕХОД В МОДЕЛИ ДИНАМИКИ БИОМАСС КОРНЕВЫХ РАСТЕНИЙ

В. В. Александров1, Т. Б. Александрова2, Л. Крусадо Лима3, X. А. Эскамилья Рейна4

В статье показана возможность решения проблемы перехода между периодическим и точечным аттракторами в бистабильной модели Розенцвейга-Макартура динамики корневых растений.

Ключевые слова: гемипаразитическое растение, модель Розенцвейга-Макартура, периодический аттрактор, множество достижимости.

The article shows the possibility of solving the problem of the transition between periodic

1 Александров Владимир Васильевич — доктор физ.-мат. наук, проф., зав. каф. прикладной механики и управления мех.-мат. ф-та МГУ; проф. физ.-мат. ф-та Автономного ун-та штата Пуэбла (Мексика), e-mail: vladimiralexandrov366Qhotmail.com.

2 Александрова Тамара Борисовна — канд. биол. наук, ст. науч. сотр. Ин-та математических исследований сложных систем МГУ, e-mail: tamara366Qyahoo.com.mx.

3 Крусадо Лима Лаура — асп. физ.-мат. ф-та Автономного ун-та штата Пуэбла (Мексика), e-mail: lcl.nfmaQgmail.com.

4Эскамилья Рейна Хуан Альберт,о — проф. физ.-мат. ф-та Автономного ун-та штата Пуэбла (Мексика), e-mail: j escami Qfcfm. buap. mx.

Aleksandrov Vladimir Vasilievich — Doctor of Fhysical and Mathematical Sciences, Frofessor, Lomonosov Moscow State University, Faculty of Mechanics and Mathematics, Head of the Chair of Applied Mechanics and Control; Professor, Autonomous University of Puebla, Mexico, Faculty of Physics and Mathematics.

Aleksandrova Tamara Borisovna — Candidate of Biological Sciences, Senior Researcher, Lomonosov Moscow State University, Institute of Complex Systems Mathematical Research.

Cruzado Lima Laura — Postgraduate, Autonomous University of Puebla, Mexico, Faculty of Physics and Mathematics.

Escamilla Reyna Juan Alberto — Professor, Autonomous University of Puebla, Mexico, Faculty of Physics and Mathematics.

and point attractors in the bistable Rosenzweig-MacArthur model with modifications for the dynamics of root hemiparasitic plants and their hosts.

Key words: hemiparasitic plant, Rosenzweig-MacArthur model, periodic attractor, reachable

set.

1. Введение. Корневые гемипаразитические (полупаразитические) растения паразитируют под землей, получая большую часть воды и питательных веществ от своих хозяев (растений — не паразитов) через гаусторий, и одновременно являются автотрофными [1]. Поскольку и гемипаразиты, и хозяева автотрофны, они конкурируют в борьбе за свет. В условиях высокой продуктивности усиление тени может снизить конкурентоспособность гемипаразитических растений по сравнению с их растениями-хозяевами, которые часто имеют преимущество в этой борьбе [2]. Подземные взаимодействия более сложны, поскольку все растения-хозяева сражаются и за почвенные ресурсы [3]. Среди корневых паразитарных растений выделяются растения Orobanchaceae, которые паразитируют на сельскохозяйственных культурах, овощах и масличных культурах, имеющих экономическое значение. Они являются самыми разрушительными вредителями сельскохозяйственных растений в мире [4].

2. Описание модифицированной модели Розенцвейга-Макартура. В модифицированной модели Розенцвейга-Макартура [5] построена взаимосвязь между этими двумя типами растений, которая расширяет хорошо известную модель хищник-жертва Розенцвейга-Макартура [6]. Модель, основанная на логистическом росте добычи, но с некоторыми изменениями, учитывающими конкуренцию в борьбе за солнечный свет и внутривидовой паразитизм среди гемипаразитов, выглядит следующим образом:

dyi (i_V1 + СУЛ аУ1У2

у 1 I

dt V K / y 1 + s

dy2 , . , . aed2yiy2 r=- = -{m + miy2)y2 +

2..

& ' (V2! + d2)(Уl + 5)'

где VI — биомасса хозяина (у1 > 0), у2 — биомасса гемипаразита (у2 > 0). Биомасса — это вес живого

2

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

Таблица1

Описание и единицы измерений переменных и параметров модели

Параметр Описание Единица измерения

У1 Биомасса растений-хозяев г/м2

У2 Биомасса растений-гемипаразитов г/м2

К Емкость среды, обусловленная сочетанием ограниченных ресурсов (свет, вода и питательные вещества) в окружающей среде г/м2

г Собственная скорость роста растения-хозяина 1/день

с Конкурентоспособность гемипаразитов за свет по отношению к своим хозяевам безразмерный

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

s Биомасса хозяина, необходимая для достижения средней максимальной скорости г/м2

то Снижение биомассы гемипаразитарных растений 1/день

ТО-1 Мера негативного взаимодействия между гемипаразитами м2/

е Эффективность использования воды и питательных веществ гемипаразитами для создания собственной биомассы безразмерный

d Биомасса хозяина, необходимая для достижения половинного значения функции д(у1) г/м2

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

a m m1

t = tr, ai =-, n = —, n\ = —.

r r r

Следовательно,

Лу, , 1 ЛУ2

У1 + су 2 к

= — (п + П1У2)У2 +

аЩЗ/2

У\ +8 ' ^1в^2у1у2

м ' (у2 + ¿2)(у1 + в)"

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

3. Качественный анализ. Критические точки модели (2) находятся на пересечении изоклин хозяев и гемипаразитов. Уравнение хозяев имеет две изоклины, одна из которых является тривиальной: у1 = 0, т.е. полное отсутствие хозяина, что биологически неинтересно. Вторая изоклина

%1,К) = У2 =

-у\ - зу\ + (2/1 + в) К

у1с + вс + а1К

пересекает ось хозяина у\ в точке у\ = К и ось гемипаразитов у2 в точке Л,(0,К) = У2 = Уравнение гемипаразитов также имеет две изоклины. Первая тривиальна: у2 = 0. Вторая имеет вид

р(У1) = У2 =

у1а1с12в

— п ) /п1.

Существуют два пересечения с осью у1, которые невозможно выразить аналитически, поэтому обозначим их А и В, предполагая, что А < В (рис. 1). Обращаем внимание на то, что р не зависит К

Всегда существуют два тривиальных положения равновесия: Ео = (0, 0) и Ек = (К, 0). Равновесие вымирания Ео описывает ситуацию без растений-хозяев и гемипаразитов, в то время как Ек описывает максимальную общую биомассу хозяина, которую может поддерживать окружающая среда. Если К принадлежит (А, В), то точка равновесия Ек является седловой точкой. Если точка К не принадлежит (А, В), то Ек — устойчивый узел. Таблица2

Наиболее интересен случай сосуществования, когда К > В. Параметры модели (2) Тогда есть две возможности: изоклина Н пересекает изоклину р в двух точках; Н нжодится выше р и они не пересекаются. В первом случае Е1 — неустойчивый фокус с расположенным вокруг него предельным циклом у0. При начальных условиях, принадлежащих области притяжения предельного цикла, траектории направляются к последнему или, если начальные условия находятся вне этой

Ек

притяжения предельного цикла через Вс, а область притяжения устойчивого узла Ек через Век (рис. 1).

Таким образом, можно сделать вывод, что система (4) имеет два аттрактора, разделенные сепаратрисой седловой точки Е2. В [7] продемонстрирована возможность перехода из области притяжения одного аттрактора в область притяжения другого аттрактора под действием небольшого постояннодействующего возмущения или управления. Для решения задачи перехода с параметрами е = 1 а = 0,6, й = 7 п = 0,1, п1 = 0,006, с = 0,1, Ь = 5 К = 15 необходимо перейти от области притяжения предельного цикла сосуществования к области притяжения устой-Ек

в правой части второго уравнения системы (2) добавляется управление которое представляет собой ручную прополку биомассы гемипаразитов в течение нескольких дней. Доказано, что прополка снижает производство новых семян, поэтому их популяция может сократиться в долгосрочной перспективе [8]. После нормализации времени при Ь2 = 1/г получаем

Параметр Единица измерения

У1 г/м2

У2 г/м2

К г/м2

с безразмерный

а\ безразмерный

в г/м2

п безразмерный

П1 м2/г

е безразмерный

(1 г/м2

£ = »11-

У1 + СУ2

к

^Г = -(П + П1У2)У2 +

атУ2 . 1 У1 + в

а1ей2у1у2

+ Ь2 ф),

, О € V = (О € КС| V- < ф) < 0}.

Система (5) может быть переписана в векторной форме:

y = f

где b = (0,b)T, € C 1(X), X С R2. Система (6) имеет периодическое решение y0(t + T) = y0(t) с периодом T, когда v(t) = 0 тождественно.

Чтобы решить эту проблему перехода, необходимо синтезировать множество достижимости Dtl для системы в отклонениях, построенной в окрестности периодического аттрактора нелинейной системы (6). Согласно [9] в окрестности предельного цикла y0(t) рассматривается система в вариациях по координате x = y — y0. Поэтому система (6) принимает вид

x = A

где b = (0, b2)T и a(t) = dip°^°{t)) при A(t + T) = A(t).

Система (2) имеет безразмерный период Т = 37.88, однако при г = 0.126 получается приблизительный период 10 месяцев. Периодическое решение y0 = (y0(t),y0(t)) нелинейной системы дифференциальных уравнений (2) может быть получено только численно в виде массива точек y0(ti) = (y0(ti), y0(¿¿)), i = 1, • • •, 500. Далее вычисляется якобиан системы (2) и в него подставляются численные массивы значений y0(ti). Матрица A(t) аппроксимируется кубическими сплайнами [10]. Находим нормированную фундаментальную матрицу решений X (T), интегрируя систему X(t) = A(t)X(t) с начальными условиями X(0) = I. Матрица монодромии X(T) имеет множители Флоке pi = 1 и р2 œ 0.0323. Согласно теореме Демидовича и Андронова-Витта [11], периодическое решение y0(t) системы (2) является асимптотически орбитально устойчивым.

Алгоритм, предложенный в [9], позволяет решить задачу построения множества достижимости линейной системы с периодическими коэффициентами при наличии возмущения. Применяя этот метод к системе (7), осуществляем переход XN(t) = X(t)M от системы координат (xi,x2) к новой системе координат (xN, x^ ). Матри ца M должна быть невырожденной для того, чтобы выполнялось условие M-1X(T)M = diag (1, P2). Для M

рами матрицы монодромии X(T), а именно M1 œ (—0.999, 0.049)T и M2 œ (0.421, — 0.907)T, соответствующими множителям Флоке pi и p^. Вычисления производятся в системе координат (xN ,xN), угол в между осями которой составляет в œ 70°. Чтобы вернуться к координатам нелинейной си-

„^i (—0.978 —0.454 \ стемы (yi, y2), во-первых, осуществляется переход X = X1 M , где M i œ nr0n n77n Ь и>

0.0529 —1.0779

во-вторых, осуществляется переход y = xN + y0.

С другой стороны, поскольку система (7) при v(t) = 0 имеет первый множитель Флоке |pi| = 1 и второй | p21 < 1, специальная матрица решений XN также может быть представлена в виде [11]

XN(t) = Ф(£) diag (1,

где Ф является непрерывно дифференцируемой T-периодической (n x п)-матрицей. Матрица Ф(t) может быть получена как В соответствии с [12] вектор-функция

ti

Xj1

0

где ei = (1,0)T, e2 = (0,1)т и G(ti,s) = Ф(íl)diag(l, e^1"^2^1-8^ Ф_1(з), 0 ^ s ^ ii, является решением системы (7) при нулевых начальных условиях. Из (8) следует, что для нахождения максимального отклонения xN, j = 1, 2, необходимо и достаточно, чтобы управление v(t) принимало экстремальные значения 0 или v-, совпадающие со знаком функции ejG(ti,t)b [13].

После получения максимальных отклонений для каждой координаты отдельно используется итерационный метод условного градиента поиска максимального отклонения по обеим координатам [13]. Выбирается любой начальный вектор cj, затем шаг за шагом ищется максимум функционала J = x2 + x2 до момента, пока не будет выполняться условие max [xi(ti, v0)xi(ti, v) +

v€V

x2(ti, v0)x2(ti, v)] = x2(ti,v0) + x2(ti,v0). Варьируя начальный вектор cj, возможно получить точки границы множества достижимости Dtl.

Период колебания T системы (2) делим на определенное количество точек и для каждой из них отстраиваем множество достижимости Dtl линейной системы (7). Данная последовательность множеств достижимости для линейной системы представляет собой аппроксимацию множества достижимости Dtl для нелинейной системы (5).

4. Численное решение задачи перехода. Интегрируем систему (2) с ti = T, v- = -0.05 и начальными условиями yo = (7.258, 9.295), принадлежащими орбите предельного цикла.

Период колебания Т делится на 76 точек, для каждой из которых строим множество достижимости. То есть существует последовательность Di(t) для г = 1,76 и i € [0, Т]. Тогда возможны два случая.

1. Если 0 ^ i < 13, то последовательность множеств достижимости Di(t), г = 1, 25, не пересекает область притяжения Bek.

2. Если 13 ^ i ^ Т , то последовательность множеств Di(t), i = 26,76, пересекает область притяжения Bek.

Таким образом, можно сделать вывод, что для второго случая возможен переход из области притяжения периодического аттрактора Be (сосуществование) в область притяжения стабильного узла Bek (вымирание гемипаразитов).

Рис. 1. Фазовое пространство случая, когда существует бистабилыгость между двумя аттракторами — точечным и периодическим. Изоклины Н (3) и р (4) обозначены пунктиром, траектории нелинейной системы (2) сплошными линиями со стрелками

Далее выбирается точка й0, которая является максимальным отклонением по двум координатам, принадлежит границе множества достижимости дБ (¿) и удовлетворяет условию € сШг(£) П Век для всех г = 26,76. Каждая точка й0 достигается за время движения от точки, принадлежащей предельному циклу у0(£), ДО точки Это движение происходит под действием управления V0 (¿), которое демонстрирует периоды, когда требуется ручная прополка. В [14] подтверждено, что ручная прополка или химическая борьба применяется после появления всходов урожая и повторно с интервалом в 2 недели во избежание заражения ОгоЬапскасеае.

Таким же образом рассчитываем приблизительное время вымирания гемипаразитов , интегрируя систему без управления (2), выби-

Рис. 2. На фазовом пространстве черным цветом показано оптимальное множество достижимости Бт, мпожество Бм — серым. Оптимальная траектория от начальной точки ^ к ^ представлена звездочками, максимальная траектория от начальной точки ¿М к ^ — серой штриховой линией

-0.01

-0.02

-0.03

-0.04

-0.05

_L

J_

0 50 100

Рис. 3. Оптимальное управление v0

150 t, дни

m для достижения точки d0 € dDj(t) П BEk

рая начальными условиями точки i = 26,76, до момента достижения сферы W = U(Ek,tk) радиуса г к = 2. Поэтому общее время tf достижения сферы W вычисляется как сумма tf = tld + tlw для i = 26, 76. Среди всех времен tj находится минимальное t™ и максимальное t^.

На рис. 2 отображены множества достижимости для минимального и максимального времени. Минимальное время tm вычисляется на первом этапе как время движения от точки, что лежит

на орбите предельного цикла y0(t), до точки = (0.004, 9.403) с управлением и составляет tm = 19.5 ( приблизительно 5.2 месяца). Ручная прополка v^ применяется 4 раза с интервалом в 8 дней (рис. 3). На втором этапе время tm вычисляется от тачального условия dm системы (2) до достижения W и составляет 300.67 (приблизительно 9.5 месяца). Таким образом, минимальное время tm ~ 1.2 года. Максимальное время tf вычисляется на первом этапе как tM до достижения точки df = (7.52, 7.96) с управлением vf и составляет ~ 37.8 (приблизительно 10 месяцев) . На втором этапе tf 2.9 года и траектория перехода от пробки df до WW, Таким образом, t^f ~ 3.74 года.

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

СПИСОК ЛИТЕРАТУРЫ

1. Press М.С. Autotrophy and heterotrophy in root hemiparasites // Trends Ecol. Evol. 1989. 4. 258-263.

2. Press M.C., Phoenix G.K. Impacts of parasitic plants on natural communities // New Phytol. 2005. 166. 737-751.

3. Cameron D.D., Hwangbo J.K., Keith A.M., Geniez J.M., Kraushaar D., Rowntree J., Seel W.E. Interactions between the hemiparasitic angiosperm Rhinanthus minor and its hosts: from the cell to the Ecosystem // Folia Geobot. 2005. 40. 217-229.

4. Parker C. Observations on the current status of Orobanche and Striga problems worldwide // Pest Manage Sci. 2009. 65.453-459.

5. Fibich P., Leps J., Berec L. Modelling the population dynamics of root hemiparasitic plants along a productivity-gradient 11 Folia Geobot. 2010. 45. 425-442.

6. Rosenzweig M.L., MacArthur R.H. Graphical representation and stability conditions of predator-prey interactions // Amer. Natur. 1963. 97. 209-223.

7. Simiu E. Chaotic Transitions in Deterministic and Stochastic Dynamical Systems. Princeton: Princeton University Press, 2002.

8. FAO (Food and Agriculture Organization of the United Nations). Progress on farmer training in parasitic weed management. Rome, Italy, 2008.

9. Коноваленко И. С. О построении множества достижимости в окрестности периодического аттрактора // Вестн. Моск. ун-та. Матем. Механ. 2020. № 3. 67-71.

10. Reithmeier Е. Periodic Solution of Nonlinear Dynamical Systems: Numerical Computation, Stability, Bifurcation and Transition to Chaos. Germany, Springer-Verlag, 1991.

11. Farkas M. Periodic motions. USA, Springer-Verlag, 1994.

12. Демидович Б.П. Лекции по математической теории устойчивости. М.: Изд-во МГУ, 1998.

13. Александров В.В., Болтянский В.Г., Лемак С.С., Парусников Н.А., Тихомиров В.М. Оптимальное управление движением. М.: Физматлит, 2005.

14. Habimana S., Nduwumuremyi A., Chinama R. Management of orobanche in field crops //J. Soil Sci. and Plant Nutrition. 2014. 14, N 1. 43-62.

Поступила в редакцию 20.1Í.2020