CC BY
17
УДК 551.46.02
И.Е. Тимченко, Е.М. Игумнова
Процессы реакции - диффузии в адаптивной модели морской экосистемы
Рассмотрено применение метода адаптивного баланса влияний для построения моделей морских экосистем с учетом диффузии моделируемых параметров. Обсуждаются особенности адаптивных экосистем, связанные с наличием в их уравнениях отрицательных обратных связей второго порядка (логистических функций), учитывающих ресурсные ограничения развития процессов в экосистемах. Приведены результаты сравнения двух вариантов адаптивных моделей: с оператором диффузии в структуре логистической функции и вне ее структуры. На примерах уравнения Колмогорова - Фишера и одномерной модели морской экосистемы показано, что включение оператора диффузии в структуру логистической функции обеспечивает приспособление процессов реакции к влиянию диффузии.
Ключевые слова: уравнения реакции - диффузии, адаптивная модель, морская экосистема.
Введение. Адаптивные модели морских экосистем основаны на принципе динамического баланса моделируемых процессов, который обусловлен стремлением экосистемы к состоянию равновесия. Баланс осуществляется путем непрерывной подстройки процессов друг к другу, а также к тем изменениям, которые происходят в окружающей среде [1 - 4]. Внутрисистемное взаимодействие процессов в экосистемах подчинено законам выживания живых организмов, популяционной динамики, химической кинетики и другим законам. Все взаимодействия этого типа обычно относят к реакциям, которые существуют в данном объеме морской среды при определенных ресурсных ограничениях. Влияние окружающей среды проявляется в виде адвекции и диффузии концентраций моделируемых параметров в этот объем, которые не зависят от внутрисистемных процессов, но существенно изменяют условия прохождения реакций. Одновременный учет внутрисистемных и внешних влияний на процессы приводит к построению уравнений реакции - диффузии. Классическими примерами подобных уравнений могут служить уравнения реакции - диффузии Колмогорова - Фишера [5, 6], химического морфогенеза Тьюринга [7], теории диссипативных систем Пригожина [8] и многие другие.
Рост концентрации любого компонента экосистемы всегда ограничен имеющимися ресурсами развития реакции. Поэтому для моделирования эффекта насыщения концентрации в уравнениях моделей экосистем обычно используется наиболее простая логистическая функция [9]. Эта функция связывает отрицательной обратной связью скорость изменения концентрации и саму концентрацию: рост концентрации автоматически замедляется по мере приближения к некоторому максимуму, обеспеченному ресурсами. Что касается внутрисистемных и внешних влияний на концентрацию, то обычно они
© И.Е. Тимченко, Е.М. Игумнова, 2013
ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 2
51
присутствуют в правых частях уравнений экосистемы отдельно от логистической функции в виде функций источников и стоков, которые не охвачены отрицательной обратной связью, регулирующей скорость реакции.
Иной подход применяется в методе адаптивного баланса влияний (АВС-методе [10]). В нем используется управление всеми влияниями на динамику моделируемого процесса - как внутрисистемными, так и внешними. Сама логистическая функция появляется в уравнениях АВС-метода как частный случай общей базовой функции влияния, реализующей системный принцип адаптивного баланса влияний. Внутрисистемные взаимодействия в уравнениях АВС-моделей включаются в аргументы базовых функций влияния. Что касается адвекции и диффузии, то они могут быть включены в структуру базовой функции наряду с другими внутрисистемными влияниями или представлены отдельными слагаемыми в правых частях уравнений. В работе [11] АВС-методом была построена интегральная модель морской экосистемы, в которой процессы диффузии учитывались отдельно от базовой функции влияния. Вариант включения оператора диффузии в аргумент базовой функции до настоящего времени не рассматривался.
В данной статье сравниваются между собой оба возможных варианта построения адаптивных моделей: сначала на примере одномерного уравнения Колмогорова - Фишера (КФ), которое при включении оператора диффузии в аргумент базовой функции влияния переходит в адаптивное уравнение Колмогорова - Фишера (АКФ), а затем на примере интегральной АВС-модели морской экосистемы.
Учет влияющих факторов при моделировании процессов адаптации в экосистемах. Рассмотрим множество взаимосвязанных параметров морской среды, представляющих в некотором объеме концентрации тех химических веществ и биологических объектов, для изучения которых необходимо разработать соответствующую модель экосистемы. Функции, представляющие изменения параметров во времени, будем называть процессами реакции - диффузии, развивающимися в морской экосистеме. Для каждого процесса ui существует некоторая ресурсная емкость морской среды Ci (Current capacity [9]), которая определяет общие ресурсные возможности развития процесса. Удобно отнести Ci к серединам интервалов изменчивости процессов. Тогда внутрисистемные и внешние влияния, изменяющие ресурсные возможности развития процесса ui , проявляются в том, что они отклоняют значения процесса от Ci в пределах его интервала изменчивости (0, max ui). Если влияния постоянны во времени, экосистема находится в стационарном состоянии равновесия, при котором значения концентраций отклонены на определенные величины от центров соответствующих интервалов Ci . Так как
max ui = 2 Ci, максимально возможные отклонения составляют ± Ci.
В моделях математической биологии [9], популяционной динамики [12], химической кинетики диссипативных систем [8] ресурсные емкости среды Ci традиционно учитываются с помощью логистической функции роста с
52
ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 2
насыщением, предложенной Верхалстом [13]. Уравнение динамики Верхал-ста с логистической функцией в правой части имеет вид
ёи „ и ч
-Г = Ггиг (1 ~)> (1)
т Ci
где гI - удельная скорость роста концентрации (отношение скорости изменения концентрации к концентрации).
Стационарными решениями уравнения (1) служат значения и* = 0 и и*2 = Ci, при которых правая часть обращается в нуль. Устойчивым решением этого уравнения является значение ui 2 = Ц, в окрестности которого небольшим приращениям ui соответствуют противоположные по знаку приращения скорости изменения щ . При отрицательной обратной связи, обусловленной структурой логистической функции, решение стремится в стационарную точку ui2 = ц , т.е. в середину интервала изменчивости и.
При объединении в систему нескольких уравнений типа (1) в их правых частях появляются внешние и внутрисистемные влияния на процесс ui со
стороны других процессов и^ и внешних факторов. Поэтому имеет значение способ учета этих влияний. Обозначим функцию, представляющую все возможные влияния, Ai (и) . Обычно эта функция аддитивно добавляется к логистической функции в правой части уравнения (1) [5 - 9]:
^ = ги (1-Ц-)+А и). (2)
т Ci
Стационарными решениями уравнения (2) тогда становятся значения
которые являются физически обоснованными. когда функции влияния A (Uj) определены на полусегментах:
- ^ < A (Uj) <¥ . (4)
4
Отметим. что в отличие от уравнения Верхалста решения уравнения (2) могут принимать значения. выходящие за пределы интервала изменчивости (0. 2С). обеспечиваемого ресурсной емкостью С. так как неконтролируемое внешнее влияние A (Uj) может превысить скорость реакции ru. и тогда концентрация
U будет испытывать колебания или неограниченно возрастать.
Другой вариант учета как внутрисистемных. так и внешних влияний получается. если включить все влияния в структуру логистической функции.
ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 2 53
Подобный способ реализован в АВС-методе [10]. Идея метода заключается в том, что равновесие в модели любой динамической системы должно быть обеспечено балансом положительных и отрицательных влияний, формируемых обратными связями. Для сохранения баланса все происходящие в системе процессы непрерывно адаптируются к сумме приложенных к ним внутренних и внешних влияний.
Исходя из принципа адаптации, правая часть уравнения для любого процесса может быть выражена в виде [10]
^ = ru{1 - 2F (+)[аг (u - А, (и,))]}, (5)
dt
где F(aiüi) - базовая функция влияния, в качестве которой может быть использована любая монотонно растущая функция. Положим, что базовая функция равна своему аргументу и в качестве параметра at выберем 0,5 C-1. Тогда вместо уравнения (5) имеем
dU~ = rUi {1 - 1[и1 - Аг (и,)]}. (6)
dt Ci
Из уравнения АВС-метода (6) непосредственно следует, что устойчивое стационарное решение этого уравнения имеет вид
и*2 = C + А, (иj). (7)
Так как процессы в экосистеме не должны принимать отрицательные значения, необходимо выполнение условий
-Ct < А, (и,) < C, 0 < и*2 < 2C,. (8)
При этих условиях включение всех влияний в аргумент базовой функции F(aini) приводит к симметричным отклонениям значений процессов от
состояния невозмущенного равновесия и12 = Ci, совпадающего с серединами интервалов их изменчивости. Знаки отклонений определяются функциями влияния А (и,) .
Для проверки полученных выводов были проведены вычислительные эксперименты с уравнениями (2) и (6) при Ct = 0,5; rt = 1. Их результаты приведены на рис. 1.
В отсутствие влияющей функции в уравнениях (2) и (6) они переходят в уравнение АВС-метода (2) (или в уравнение Верхалста (1)). Так как ресурсная емкость приходится на середину интервала возможных значений процесса, обусловленных влияниями, при А = 0 наблюдается быстрая сходимость решений уравнения (2) в точку и{ = 0,5. Процесс сходимости иллюстрируют
стрелки на рис. 1, а. На рис. 1, б приведен пример отрицательного влияния А = - 0,2. Штриховой линией показано решение уравнения (6), которое соот-
54 ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 2
ветствует этому отклонению: точка равновесия перешла из положения 0,5 в положение 0,41, что указывает на сокращение С, вследствие отрицательного влияния на ресурсные возможности развития процесса иг . В то же время
уравнение (2) при влиянии А = - 0,2 не имеет решения, так как кривая графика 2 не пересекается с линией 3. Этот пример подтверждает полученное выше ограничение (4): А в уравнении (2) не должно быть меньше - 0,125 при гг = 1, Сг = 0,5.
а б
в г
Р и с. 1. Графики левых и правых частей уравнений (2) и (6), представленных в конечных разностях, при различных значениях влияющих функций А (а - г): 1 - уравнение (6); 2 - уравнение (2); 3 - линия и = и; 4, 5 - касательные в точке решения
Сходимость к устойчивому решению обеспечивается углом наклона касательных (см., например, штриховые линии 4 и 5 на рис. 1, г) к графикам 1 и 2 в точках их пересечений с линией 3. Когда этот угол превышает по модулю значение п/4, возможно появление периодических решений или вычислительной неустойчивости. Этот случай иллюстрирует рис. 1, в, на котором штриховыми стрелками отмечено формирование периодических решений уравнения (2) при А = 0,5. Как показали эксперименты, значение А = 0,5 является предельным по величине для устойчивости решения уравнения (6), тогда как для уравнения (2) проявления неустойчивости наступают при более низких значениях А. Такой вывод вытекает из анализа графиков 2 (рис. 1, б, в, г), так как они начинаются со значений оси ординат, равных А, и имеют большую крутизну в точках пересечения с линией 3, чем кривые 1. ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 2 55
Проведенный анализ показывает, что распространение отрицательной обратной связи на все влияния (внутрисистемные и внешние) в уравнении (6) усиливает управляющую роль базовой функции, так как она непосредственно влияет на текущую (эффективную) величину ресурсной емкости развития процесса. Диапазон влияний, не нарушающих условия устойчивости решений, в этом случае оказывается симметричным относительно первоначальной ресурсной емкости Сi. Этим подтверждается целесообразность включения А (и^) в аргументы базовых функций, как это имеет место в АВС-методе моделирования морских экосистем [6].
Адаптивный аналог уравнения реакции - диффузии Колмогорова -Фишера. Рассмотрим использование полученных результатов для случая, когда концентрация некоторого элемента экосистемы распределена вдоль отрезка прямой (0, Ь), лежащего на оси х, а в качестве внешнего влияния фигурирует диффузия концентрации, обусловленная условиями окружающей среды. Без потери общности мы можем перейти к безразмерным независимым переменным х и а также функции и путем соответствующей нормировки на масштабные множители. Если заменить функцию А (и;-) в уравнении (2) диффузионным членом, то это уравнение становится широко известным уравнением реакции - диффузии Колмогорова - Фишера [5, 6]
du „ u . ^ Э 2u
---- (1--) + D —
dt C Эх
= ru (1 + D^-y. (9)
где D - коэффициент диффузии. В большом количестве работ. посвященных этому уравнению. рассматриваются условия существования решений при различных постановках задач для функции u(х. t) . заданной на отрезке (0. L) или на всей оси х. В частности. для задачи Дирихле с начальными и краевыми условиями вида
u(0. х) - f (х). (10)
u(t.0) - 0. u(t. L) - 0 (11)
решения. которые должны обращать в нуль правую часть уравнения (9) и удовлетворять этим краевым условиям. находятся из уравнения
Э 2u r .„ 1 .
ЭХ2-- Du(1 - Cu). (l2)
Для выяснения свойств решений уравнения (9) методом фазовой плоскости [14] вводится система обыкновенных дифференциальных уравнений
du , .
- у - g(u. v).
Эх
dv r 1 (13)
— ---u(1--u) - h(u. v) .
Эх D V C 1 V ' '
56 ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 2
При краевых условиях (11) эта система имеет две точки равновесия: Ых (0,0) и М2 (С, 0) , соответствующие стационарным решениям уравнения (9). Для анализа устойчивости решений в окрестностях точек равновесия рассматривается якобиан J(и,V):
3(и, V) :
(а*; ( 0 1 ^
ди' дv
дк ; дк — 2г -и - Г; 0
V ди' дv ) V СБ В )
(14)
Линеаризация 3(и, V) в окрестности точки М1(0,0) приводит к нахождению двух чисто мнимых собственных чисел
4,2 = ±1\1Б •
(15)
Из этого следует, что точка М1(0,0) представляет собой центр на фазовой плоскости [14]. Подобная точка свидетельствует о наличии устойчивых периодических решений уравнения КФ.
Аналогичный анализ в окрестности второй точки М2(С,0) дает пару действительных собственных чисел
42 =±л-4,2 л| в
(16)
и эта точка равновесия представляет собой седло. Ей соответствует неустойчивое решение апериодического характера.
Анализ уравнения КФ позволяет сделать выводы о характере решений уравнения в зависимости от значений параметров г, В, С . В частности, было показано, что для обеспечения устойчивости решений длина Ь отрезка оси х, на котором рассматривается процесс реакции - диффузии, должна быть больше некоторой критической величины Ь*, определяемой параметрами уравнения КФ [15]:
ь > ь*=4в
(17)
Если это условие будет нарушено, концентрация и(х, I) с учетом поглощающих граничных условий (11) будет стремиться к нулю. Физически это означает, что геометрический размер той части среды, в которой развиваются процессы реакции, ограничивает возможности роста этих процессов. Параметр Ь* , зависящий от отношения интенсивности диффузии в к параметру скорости реакции г, лимитирует рост концентрации так же, как и ресурсная емкость среды С.
0233-7584. Мор. гидрофиз. журн., 2013, № 2
57
Перейдем теперь к рассмотрению адаптивного аналога уравнения Колмогорова - Фишера (АКФ). Для этого используем уравнение (6), в котором диффузионное слагаемое непосредственно включено в базовую логистическую функцию Ai и) с коэффициентом В соответствующей размерности.
Как уже было отмечено выше, смысл перехода к уравнению АКФ состоит в том, чтобы моделируемая концентрация и(х, I) автоматически подстраивалась под влияние диффузии за счет отрицательной обратной связи, обеспечиваемой базовой логистической функцией. Уравнение АКФ принимает вид
^ - ru[1 -l(u -Dd2u)]. dt L CK Эх
(18)
Из условия стационарности следует, что правая часть уравнения АКФ обращается в нуль при тех же значениях и* = 0 и и* = С , что и в уравнении КФ (9). Система уравнений (13) в данном случае переходит в систему
^ - у - g(u.у). Эх
dv 1 C
— - — u--- h(u. v).
Эх D D
(19)
Якобиан системы уравнений (19)
J(u. v) -
0;
D-1
(20)
оказывается не зависящим от и, что не дает возможности исследовать поведение решений уравнения АКФ в окрестностях точек равновесия. Однако условие стационарности АКФ представляет собой линейное уравнение
Э 2u 1
C
^ 2---u + — - 0.
Эх2 D D
решение которого. удовлетворяющее условиям (11). имеет вид
u (х) - С,еГах + C2е-^х + Ъ
(21)
C1 -
а -1
i2eFaL -e-VaL - 1Л
e^VaL - e^L
. C2 -
а-1
e
а -1
■jaL - ^
-1
->/aL s-sfaL
e
1
C
(22)
а - —. Ъ - — . D D
Поэтому некоторые свойства стационарных решений уравнения АКФ могут быть установлены по формуле (22) путем вычислительных экспериментов.
Рассмотрим зависимость между коэффициентом диффузии D и длиной интервала L для уравнения АКФ (18). которая вытекает из условий стационарности решений. На рис. 2 приведен рассчитанный по этому уравнению
58 ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 2
\
Ъ
Ъ
график кривой стационарного решения и (В, Ь) для варьируемых значений коэффициентов диффузии В и длин интервала Ь. По оси х заданы значения параметра N с которым коэффициенты диффузии и длины интервалов были связаны линейными зависимостями: В = 6,5N и Ь = 0,29Ж На этом же рисунке показано отношение Ь/л/В как функция параметра Ы, которое (при г = 1) входит в формулу (17), выведенную для уравнения КФ. Это отношение не должно быть таким, при котором значения концентрации и(В, Ь) становятся отрицательными, т.е. теряют физический смысл. Как следует из рисунка, критическим значением для уравнения АКФ служит величина отношения
Ь/л/В = 0,14, когда стационарные решения для концентраций и(х,I) обращаются в нуль и далее переходят в область отрицательных значений. Поэтому эмпирическая формула, связывающая коэффициент диффузии с длиной интервала, для уравнения АКФ имеет вид
Ь* @ 0,14 В. (23)
Р и с. 2. Стационарные решения адаптивного уравнения Колмогорова - Фишера и(Ь, В) и индекс реакции Ь/\1~0 при варьировании параметров Ь и В
Включение диффузии в базовую логистическую функцию в уравнении АКФ уменьшает по сравнению с уравнением КФ предельную длину интервала Ь* оси х, для которого процессы реакции затухают и перестают существовать в силу поглощающих граничных условий (11).
При переходе к конечно-разностным представлениям уравнений КФ и АКФ возникают дополнительные вопросы устойчивости решений этих уравнений. Известно, что пространственная дискретизация способствует появлению разнообразных периодических решений, характер которых зависит от соотношения интенсивностей реакции и диффузии [9, 15]. Поэтому для оценки устойчивости решений уравнений КФ и АКФ при конкретных значениях
0233-7584. Мор. гидрофиз. журн., 2013, № 2
59
параметров уравнений целесообразно прибегнуть к вычислительным экспериментам.
В проведенных экспериментах в качестве пространственной области Ь был выбран отрезок оси х, содержавший 10 узлов сеточной области, в которых были заданы следующие начальные условия: х1 = 1; х2 = 3; х3 = 5; х4 = 7; х5 = 9; х6 = 9; х7 = 7; х8 = 5; х9 = 3; х10 = 1. На отрезке Ь рассматривались задачи Дирихле для уравнений КФ и АКФ с нулевыми граничными условиями х0 = хЬ = 0 . При г =1, С = 5, Аt = 1, Ах = 1
были использованы следующие конечно-разностные аппроксимации по схеме Эйлера:
для уравнения КФ
DАt
(Ах)2
и для уравнения АКФ
DАt
(Ах?
u\ = 2uk- [1 - 0, l(uk-)] + —2 К- - 2uk- + ) (24)
u\ = 2uk_i{1 - 0,1[uk_i) - —— (uk-_1 - 2uk_i + uk+1)]}. (25)
Известно, что для устойчивости решения уравнения КФ (24) требуется выполнение следующего условия [9, 15]:
^ £ (26)
(Дх)2 2
в котором Dt и Дх обозначают шаги дискретизации задачи по времени и по пространству. Структура уравнения АКФ позволяет сформулировать аналогичное условие. Процесс итераций по формуле (25) остановится, когда правая часть формулы станет сколь угодно близкой к левой. Для этого необходимо, чтобы выражение в фигурных скобках стремилось к значению 0,5 и, следовательно, выражение в квадратных скобках стремилось к 5. Так как
0 £ u'k-1 £ 2С = 10, а (uk-- - 2u'k-1 + uk- ) £ 20 , диффузионный член в квадратных скобках не должен превышать 5 по модулю. Поэтому из условия устойчивости численного решения уравнения АКФ вытекает оценка
DDt 1
-7 £ - . (27)
(Дх)2 4
При решении уравнений (24) и (25) было поставлено дополнительное условие, чтобы в ходе итераций значения процессов не выходили за пределы интервала (0, 10), обусловленного заданной ресурсной емкостью С = 5:
u(х) = IF{u(х) < 0;0,01;IF[u(х) > 10;9,99; u(х)]}. (28)
В соответствии с оценкой (26) для обеспечения устойчивых решений уравнения КФ при Дх = 1 коэффициент диффузии не должен был превышать значено ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 2
ние В = 0,5. Для уравнения АКФ эта величина из формулы (27) не должна была превышать значение В = 0,25. Результаты экспериментов показаны на рис. 3.
108 -6 -
4 -2 i 0
7 6
5 4 3 210
U(Xi,t)
i—1—i—1—i—1—i—1—i—1—i—1—i—1—i—1—i—1—i
а
u(x)
АКФ
D = 0,25 N = 300
б
1—1—I—1—I—'—I—1—I—1—Г
x
"I—1—I
u(x)
1—1—I—1—I—1—I
D = 0,5 N = 300
"Т—1—I—1—I—1—I—1—Г
"I—1—I—1—Г
1210-8 -6 -4 -2 -0
D = 1,0 N = 300
123456789 10 123456789 10 д е
x
в
г
Р и с. 3. Значения концентрации и(х) в узлах сетки: а - сходимость к равновесию в различных узлах (цифры - номера узлов); б - сглаживание начального распределения после первых 4 итераций (цифры - номера итераций); в - е - распределения в узлах после 300 итераций при различных коэффициентах диффузии В (КФ - решения уравнения Колмогорова - Фишера, АКФ - решения адаптивного уравнения Колмогорова - Фишера)
Как следует из рис. 3, а, в отсутствие диффузии решения уравнений КФ и АКФ быстро сходятся к устойчивому стационарному значению u* = С = 5. Быстрая сходимость наблюдается во всех точках области, независимо от начальных условий. Поэтому уже после первых 5 итераций концентрации достигают равновесных значений (см. рис. 3, б). При малых коэффициентах ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 2 61
диффузии (меньше 0,25) решения уравнений КФ и АКФ практически не отличаются (см. рис. 3, в). С увеличением коэффициента возникают волны (см. рис. 3, г), амплитуда которых растет существенно быстрее для уравнения КФ, чем для АКФ. Уже при D = 0,75 решения уравнения КФ выходят за пределы удвоенной ресурсной емкости 2С = 10 (см. рис. 3, д) и в дальнейшем остаются ограниченными условием (28) (см. рис. 3, е). Следует отметить, что при больших коэффициентах диффузии распределения концентраций, полученные из уравнения АКФ после 300 итераций, не испытывают значительных колебаний (см. рис. 3, д, е). В этом случае диффузия моделируемой субстанции привела к резкому понижению ее концентрации ввиду поглощения ее на границах области. Результаты эксперимента свидетельствуют о том, что учет диффузии в структуре базовой функции влияния уравнения АКФ как фактора, управляющего ростом концентрации, приводит к повышению устойчивости решений, что выглядит физически вполне обоснованным.
Процессы реакции - диффузии в одномерной адаптивной модели морской экосистемы. Рассмотрим множество взаимосвязанных концентраций параметров экосистемы, распределенных вдоль некоторого горизонтального отрезка прямой линии в морской среде. В каждой точке отрезка происходят процессы взаимодействия химических и биологических объектов между собой, относящиеся к категории реакций, на которые влияют динамические процессы диффузии, развивающиеся в морской среде вдоль направления отрезка. В целях упрощения мы не рассматриваем влияние адвекции.
Одномерная адаптивная модель морской экосистемы может быть основана на использовании системы взаимосвязанных уравнений АВС-метода, взятых в форме КФ или АКФ. Таким образом, речь снова идет о двух вариантах построения адаптивной модели с операторами диффузии: стоящими отдельно от базовых функций влияния в уравнениях экосистемы или включенными в аргументы этих функций. Адаптивную модель экосистемы, у которой операторы диффузии не включены в базовые функции влияний, будем обозначать символом ABC+D, а в случае включения - ABCD.
Биологические объекты экосистемы будем характеризовать тремя интегральными концентрациями: фитопланктона PP, зоопланктона ZP и биоресурса BR, понимая под биоресурсом концентрацию всех гидробионтов, расположенных на более высоких, чем зоопланктон, уровнях пищевой цепи. Дополнительное усложнение задачи моделирования морской экосистемы связано с ресурсными ограничениями процессов реакции, которые в отличие от общих ресурсных емкостей морской среды Ci должны быть функциями времени. Так, например, ресурсами развития зоопланктона ZP служат: фитопланктон PP, кислород OX и биогенные элементы NT. В каждый момент времени рост концентрации зоопланктона будет зависеть только от одного из этих видов ресурсов, концентрация которого имеет минимальное значение. Динамика среды и, в частности, диффузия непрерывно влияют на соотношение ресурсов развития реакций. Поэтому в модели экосистемы должны быть предусмотрены агенты управления AG, следящие за лимитирующими значениями ресурсов. Так как целью исследований является анализ влияния диффузии на процессы реакции в одномерной модели, мы не будем рассматри-
62 ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 2
вать другие внешние влияния. Примем, что ресурсом развития фитопланктона является только концентрация углекислого газа СБ, а для замыкания окислительно-восстановительного цикла введем в структуру модели концентрацию детрита БТ. Концептуальная модель экосистемы с соответствующими агентами управления в ее структуре показана на рис. 4.
Р и с. 4. Концептуальная модель внутрисистемных связей в морской экосистеме с учетом внешнего влияния диффузии
Система уравнений адаптивной модели экосистемы ABC+D с операторами диффузии, стоящими отдельно от базовых логистических функций, может быть представлена таким образом:
^ = г1Щ(x){1->г(x) -4(и})]} + Dr ЩЩ . (29) dt ц dx
Система уравнений адаптивной модели экосистемы ABCD с операторами диффузии, включенными в аргументы базовых функций, имеет общий вид
du.(x) . . f1 l г . . . . . ^ д2и.(x)..-.
= ru(x){1 - —[u..(x)- A(uj)-D—. (30)
dt C. dx
С учетом концептуальной модели экосистемы, изображенной на рис. 4, в уравнениях (29) и (30) используем следующие обозначения и формулы:
ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 2 63
Щ (х) = [РР( х); 7Р( х); БЯ( х); ОХ (х); ИТ (х); БТ (х); СБ( х)],
РР / со
СБ,
= а2Р/ВКБЯ - ^С2Р[РР, ОХ, ИТ],
Л07Р [ РР, ОХ, ИТ ] = ПМр (г) = РР(г); -а7Р / РРРР(г); 0] + + ЩМгр (г) = ИТ (г); -а2р / тЫТ (г); 0] + ЩМгр (г) = ОХ (г); -а2р / охОХ (г); 0];
М2Р = а^ шт[ РР (г); ОХ (г); ИТ (г)],
ЛвК = ЛОвК [7Р, ОХ, ИТ ],
ЛОш [7Р ОХ ИТ] = Т¥[М^ (г) = 7.Р(1 )• -а„„ ^7Р(г); 0] +
+ ТБ[М7Р (г) = ОХ (г); ^ / охОХ (г); 0]; МвК = агв ш1п[ 7Р (г); ОХ (г); ИТ (г)],
Лит = ЛСдг [ОХ, БТ ],
ЛОш [ОХ, БТ ] = ПМ. (г) = ОХ (г); -а^т / охОХ (г); 0] + + ЩМм (г) = БТ (г); -аит / ШБТ (г); 0];
МИТ = аг§ шт[ОХ (г); БТ (г)],
В численных экспериментах с моделями (29) и (30) решения этих уравнений находились при тех же значениях параметров: г = 1, С = С = 5, Аг = 1, Ах = 1, а также начальных и краевых условиях, которые были использованы выше в экспериментах с уравнениями КФ и АКФ. Для упрощения задачи концентрации всех элементов экосистемы были представлены в безразмерной форме путем нормировки на ресурсные емкости С и приведения их значений к интервалу изменчивости (0, 10). В каждой из моделей коэффициенты диффузии считались одинаковыми для всех переменных: Д = Б. Кроме того, было поставлено условие (28), удерживающее решения в пределах интервала (0, 10), обусловленного заданной ресурсной емкостью С = 5. Коэффициенты взаимных влияний (реакций) а, входящие в соотношения (31), выбирались в диапазоне значений (0,1; 0,5) с таким расчетом, чтобы обеспечить компромисс между чувствительностью модели к изменениям концентраций пара-
(31)
64
1ББИ 0233-7584. Мор. гидрофиз. журн., 2013, № 2
метров экосистемы и устойчивостью вычислительной схемы Эйлера. Расчеты проводились на 370 безразмерных шагов по времени. Результаты расчетов представлены на рис. 5 - 7.
Как следует из рис. 5, а - з, начальные пространственные распределения концентраций фитопланктона, зоопланктона и биоресурса, изображенные штриховыми линиями, трансформировались после 300 итераций в распределения, согласованные между собой и с нулевыми краевыми условиями. При малых коэффициентах диффузии (см. рис. 5, а и б) распределения, полученные по моделям АВС+D и ABCD, практически совпадают. С увеличением влияния диффузии различия в решениях становятся заметными (см. рис. 5, в - е). Предельными значениями коэффициентов диффузии для модели АВС+D стали величины порядка 0,6, для модели ABCD - порядка 0,3 (рис. 5, ж и з). При дальнейшем увеличении коэффициентов диффузии величины концентраций выходят за пределы допустимого интервала изменчивости, который был ограничен условиями (28). На всех рисунках заметно влияние нелинейных агентов управления на концентрации зоопланктона и биоресурса, которое привело к более высокой изменчивости распределений ZP и BR по сравнению с концентрацией фитопланктона PP.
Сходимость итерационных процессов для моделей АВС+D и ABCD при одинаковых по величине коэффициентах диффузии после 300 итераций показана на рис. 6. Из сравнения сценариев, приведенных на рис. 6, а и б, видно, что при малых коэффициентах (D = 0,2) наблюдается различие в установившихся значениях концентраций PP, ZP и BR. При D = 0,3 модель ABCD демонстрирует установившиеся распределения концентраций ZP и BR, в то время как модель АВС+D воспроизводит колебания концентраций этих параметров.
Динамика гидрохимических процессов представлена на рис. 7. Цель эксперимента заключалась в том, чтобы сравнить роль агентов управления процессами AG и роль диффузии в формировании временных сценариев и пространственных распределений концентраций кислорода OX, биогенных элементов NT, детрита DT и углекислого газа CD. На рис. 7, а - г показана изменчивость их концентраций в 5-м узле расчетной сетки, в котором начальные значения концентраций составляли 9 безразмерных единиц. Влияние агентов управления при отсутствии диффузии заметно из сравнения сценариев, приведенных на рис. 7, а и б. При включении агентов увеличивается количество причинно-следственных связей между процессами в экосистеме (см. рис. 4) и происходит ресурсное лимитирование их взаимодействий. Поэтому стационарные состояния концентраций на рис. 7, б существенно иные, чем на рис. 7, а.
При включении диффузии с коэффициентом D = 0,2 адаптация гидрохимических процессов происходит более быстро в модели АВС+D (см. рис. 7, в и г). Однако в этой модели возникли периодические колебания концентрации углекислого газа CD, которые установились и продолжались весь период вычислений (см. рис. 7, д). Модель ABCD продемонстрировала более медленный пере6ход к стационарному состоянию, но результирующие значения концентраций параметров установились постоянными по величине (см. рис. 7, е).
ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 2 65
10-1 8
6 Н 4 2 1 0
108 6 -4 -2 0 108 -6 -4
2 Н 0
-"-1-"-1-"-1-"-1-1-г
а
ЛБО+Э Э = 0,3
п—"—I—"—I—"—I—"—I—'—I—I—I—I—I—'—I—'—I
в
т—I—I—I—I—I—I—I—I—I—I—I—I—I—I—>—I—I—I
б
ЛБОЭ Э = 0,15
н—'—I—'—I—'—I—'—I—|—I—|—I—|—I—1—I—|—I
ЛБО+Э Э = 0,4
\
ЛБОЭ Э = 0,2
г
N
12108 -6 -4 Н 2 0
ЛБО+Э Э = 0,6
—I—I—I—I—I—I—I—I—I—I—I—I—I—I—'—I—I—I
е
ЛБОЭ Э = 0,3
1 23456789 10 1 23456789 10
X
д
з
Р и с. 5. Пространственные распределения концентраций биологических элементов экосистемы при различных коэффициентах диффузии на 300-м шаге итераций (а - з): ЛВСБ - оператор диффузии в аргументе базовой функции (уравнение (30)), АВС+Б - вне ее (уравнение (29)); штриховыми линиями показаны начальные распределения всех параметров экосистемы
66
0233-7584. Мор. гидрофиз. журн., 2013, № 2
7 -6 -
5 -43 2
Н 0
10 -
8 -
6
4 Н 2 0
12 -| 10 8
6 Н 4
2 Н 0
ZP
в-в-в-в-в-в-в-в-в-в-в
PP
ZP
В-в-в-в-в-в-в-в-в-в-в
PP
BR
ABC+D D = 0,2
N + 360
I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I
а
ABC+D D = 0,3
ZP
PP
BR
N + 360
ABCD D = 0,2
BR
N + 360
I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I
б
BR
PP
ABCD D = 0,3
ZP
-а-в-в-
N + 360
-а-в-н
BR N + 360
123456789 10 11 123456789 10 11
I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I
г
ZP
□-в-в-в-в-в-в-в-в-в-□
ABCD D = 0,4
PP
д
е
Р и с. 6. Распределения концентраций в 5-м узле сетки при различных коэффициентах диффузии после 360 итераций (а - з): АВС+Б - диффузия вне базовых функций, АВСБ - диффузия включена в их аргументы
Рис. 7, ж и з дают представление о финальных распределениях концентраций гидрохимических элементов на 370-м шаге итераций. К этому времени они уже не зависели от начальных распределений, показанных штриховой линией, и оказались весьма разными для двух сравниваемых моделей. В модели АВС+Б распределение концентраций углекислого газа СБ имеет характер волны с пространственным периодом 2Дх (см. рис. 7, ж) и временным периодом 2Д1 (см. рис. 7, д). Распределения других элементов под влиянием диффузии оказались сглаженными и слегка убывают по мере приближения к поглощающим границам области (рис. 7, ж).
0233-7584. Мор. гидрофиз. журн., 2013, № 2
67
108 ■ 6 ■ 4 ■ 2 i 0
12 108 -6 -4 2 0
i 1 1 1 1 i 1 1 1 1 i 1 1 1 1 i 1 1 1 1 i
а
ABC+D, D = 0.2, AG*0
i 1 1 1 1 i 1 1 ' 1 i 1 1 1 1 i 1 1 1 1 i
б
ABCD, D = 0.2, AG*0
8 6 4 2
0
12108 -6 -4 -2
0
т-1-1-1-1-1-1-r
д
ABC+D, D = 0,2 AG*0 N=370
0-—<
-|-1-1-1-1-Г
X »
"1-1
ABCD, D = 0,2 AG*0 N= 370
360+N 1-1-1-1
X
-|-1
-1-1-1-1-1-1-
1 23456789 10 1 23456789 10
ж 3
Р и с. 7. Динамика и распределения концентраций химических элементов в 5-м узле сетки: а -без учета диффузии и агентов управления, б - без учета диффузии с агентами управления, в, г - с учетом диффузии и агентов управления на первых 50 шагах вычислений, д, е - с учетом диффузии и агентов управления на последних 10 шагах вычислений, а также финальные распределения концентраций химических элементов (ж, з) во всех узлах сетки с учетом диффузии и агентов управления на 370-м шаге вычислений (штриховая линия - начальные распределения концентраций)
68 ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 2
в
г
е
В модели ЛВСБ картина пространственных распределений оказалась более сложной, что подтверждает более высокую чувствительность этой модели к значениям коэффициента диффузии. Наибольшую изменчивость имеет распределение детрита БТ (см. рис. 7, з). Концентрации остальных элементов распределены более равномерно, причем, поглощающая роль граничных условий выражена весьма слабо.
Проведенные эксперименты дают основания для следующих выводов.
1. Метод адаптивного баланса влияний (АВС-метод), используемый для построения интегральных моделей морских экосистем, может быть распространен на пространственно-временные модели экосистем, учитывающие динамику гидрофизических процессов путем включения операторов переноса и диффузии в уравнения интегральной модели экосистемы в качестве дополнительных влияющих функций.
2. Базовые функции влияний, применяемые в модульных уравнениях АВС-метода, естественным образом учитывают ресурсные ограничения процессов развития в экосистемах, устанавливая отрицательные обратные связи между динамикой процесса и ресурсной емкостью окружающей среды.
3. Существует два способа построения интегральных моделей адаптивных пространственно-временных экосистем: по типу уравнения Колмогорова - Фишера без включения диффузии (и переноса) в аргументы базовых функций влияний и по типу АВС-метода, или адаптивного уравнения Колмогорова - Фишера, когда отрицательная обратная связь управляет подстройкой процессов реакции под внешние процессы диффузии (и переноса в общем случае).
4. Сравнение этих двух способов на примере упрощенной модели адаптивной морской экосистемы показало, что включение оператора диффузии в аргументы базовых функций влияний в уравнениях модели повышает ее чувствительность к изменениям ресурсной емкости развития процессов, которые обусловлены наличием диффузии.
Заключение. Вычислительные эксперименты, проведенные с уравнениями Колмогорова - Фишера и пространственными одномерными моделями морской экосистемы, подтвердили тот факт, что в адаптивных вариантах этих моделей происходит активная подстройка процессов реакции к процессам диффузии. Подобные результаты могут быть полезны для обоснования структуры адаптивных моделей морских экосистем, в которых описание химико-биологических процессов реакции и физических процессов переноса и диффузии объединяются в одну общую систему уравнений.
Традиционно в моделях морских экосистем процессы реакции и диффузии аддитивно входят в правые части уравнений динамики концентраций веществ в качестве функций источников и стоков [16]. Баланс процессов в модели экосистемы достигается за счет связей между переменными, входящими в различные уравнения модели. Ограничение функций, представляющих решения уравнений, происходит также за счет установления поглощающих условий для концентраций веществ и живых объектов на границах области.
Адаптивные модели экосистем отличаются от традиционных моделей тем, что они описывают отклонения параметров экосистемы от некоторого стацио-
ВБК 0233-7584. Мор. гидрофиз. журн., 2013, № 2 69
нарного состояния равновесия, в котором экосистема должна находиться, когда все внешние влияния равны нулю. В адаптивных моделях реализуется принцип управляемого баланса процессов: отрицательные обратные связи второго порядка присутствуют в каждом из уравнений динамики концентраций благодаря базовым функциям влияний, которые заставляют концентрации подстраиваться под сумму влияний, включенных в их аргументы. Базовая функция контролирует текущую ресурсную емкость в конкретном объеме морской среды и ограничивает пределы роста параметров экосистемы.
Поскольку концентрации живых объектов экосистем ограничены наличием ресурсов, необходимых для их существования, в уравнениях экосистем должны присутствовать агенты управления, следящие за минимальными количествами ресурсов. Перенос и диффузия объектов способны существенно изменять условия ресурсного лимитирования в данном объеме морской среды. Так, например, вселение вследствие диффузии (и переноса) дополнительного количества живых объектов уменьшает ресурсные возможности существования резидентов в этом объеме. Когда оператор диффузии стоит в уравнении отдельно от базовой функции влияния, он не влияет на эффективную ресурсную емкость среды обитания организмов. Включение его в аргумент базовой функции позволяет учесть влияние динамики морской среды на ресурсные возможности роста концентраций биологических элементов экосистемы.
СПИСОК ЛИТЕРАТУРЫ
1. Тимченко И.Е. Системные методы в гидрофизике океана. - Киев.: Наукова думка, 1988. - 180 с.
2. Тимченко И.И., Игумнова Е.М. Ассимиляция данных наблюдений и адаптивный прогноз природных процессов // Морской гидрофизический журнал. - 2009. - № 6. -С. 47 - 70.
3. Еремеев В.Н., Игумнова Е.М., Тимченко И.Е. Моделирование эколого-экономических систем. - Севастополь: ЭКОСИ-Гидрофизика, 2004. - 320 с.
4. Иванов В.А., Игумнова Е.М., Латун В. С., Тимченко И.Е. Модели управления ресурсами прибрежной зоны моря. - Севастополь: ЭКОСИ-Гидрофизика, 2007. - 258 с.
5. Kolmogoroff A., Petrovsky I., Piscounoff N. Etude de l'equation de la diffusion avec croissance de la quantité de matière et son application a un problème biologique // Moscow University. Bull. Math. - 1937. - № 1. - Р. 1 - 25.
6. Fisher R.A. The wave of advance of advantageous genes // Ann. Eugenics. - 1937. - 7. -P. 353 - 369.
7. Turing A.M. The chemical basis of morphogenesis // Phil. Trans. Roy. Soc. Lond. B. - 1952. - 237. - P. 37 - 72.
8. Николис Г., Пригожин И. Самоорганизация в неравновесных системах. - М.: Мир, 1979. - 512 с.
9. Murray J.D. Mathematical Biology. II: Spatial Models and Biomedical Applications. - Springer, 2008. - 736 p.
10. Тимченко И.Е., Игумнова Е.М., Тимченко И.И. Системный менеджмент и АВС-тех-нологии устойчивого развития. - Севастополь: ЭКОСИ-Гидрофизика, 2000. - 225 с.
70
ISSN 0233-7584. Мор. гидрофиз. журн, 2013, № 2
11. Романовский Е.В., Тимченко И.Е. Адаптивные процессы в модели морской экосистемы, основанной на уравнениях реакции - диффузии // Системы контроля окружающей среды. - Севастополь: ЭКОСИ-Гидрофизика, 2010. - С. 243 - 246.
12. Sharov A.A. Life-system approach: a system paradigm in population ecology // Oikos. -1992. - 63 - Р. 485 - 494.
13. Verhulst P.-F. Recherche mathémathiques sur le loi d'accroissement de la population // Nouveau Mémoires de l'Académie Royale des Sciences et Belles Lettres de Bruxelles. - 1845. -№ 18. - P. 3 - 38.
14. Немыцкий В.В., Степанов В.В. Качественная теория дифференциальных уравнений. -M.: Гостехиздат, 1950. - 450 с.
15. Cantrell R.S., Cosner C. Spatial Ecology via Reaction - Diffusion Equations // Ser. Math. Comput. Biol. - Chichester, UK: John Wiley and Sons, 2003. - 421 p.
16. Perez-Munuzuri V., Huhn F. The role of mesoscale eddies time and length scales on phytop-lankton production // Nonlin. Process. Geophys. - 2010. - № 17. - P. 177 - 186.
Морской гидрофизический институт НАН Украины, Материал поступил
Севастополь в редакцию 28.10.11
После доработки 06.12.11
АНОТАЦ1Я Розглянуто застосування методу адаптивного балансу впливiв для побудови моделей морських екосистем з урахуванням дифузп параметрiв, яш моделюються. Обговорюють-ся особливост адаптивних екосистем, пов'язаш з наявнютю в Тх рiвняннях негативних зворот-них зв'язюв другого порядку (лопстичних функцш), яю враховують ресурсш обмеження роз-витку процеав в екосистемах. Наведено результати порiвняння двох варiантiв адаптивних моделей: з оператором дифузп в структурi лопстичноТ функцп i поза И структурою. На прикладах рiвняння Колмогорова - Фшера та одновишрно'Т моделi морськоТ екосистеми показано, що включення оператора дифузп до структури лопстичноТ функцп забезпечуе пристосування процешв реакцп до впливу дифузп.
Ключовi слова: рiвняння реакцп - дифузп, адаптивна модель, морська екосистема.
ABSTRACT Application of the method of adaptive balance of causes for constructing marine ecosystems models taking into account diffusion of the simulated parameters is considered. The features of adaptive ecosystems are discussed in view of presence (in their equations) of the second order negative feedbacks (logistic functions) taking into account the resource constrains of the processes' development. The results of comparison of two versions of adaptive models - the diffusion operator is within the logistic function structure and outside it - are represented. The examples of the Kolmogorov-Fisher equation and the one-dimensional model of the marine ecosystem show that inclusion of the diffusion operator in the logistic function structure provides adaptation of the reaction processes to the diffusion influence.
Keywords: equations of reaction-diffusion, adaptive model, marine ecosystem.
ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 2
71
