Численная реализация модели таксис - реакция - диффузия, описывающей динамику системы хищник - жертва Текст научной статьи по специальности «Математика»

УДК 519.6:574.3

ЧИСЛЕННАЯ РЕАЛИЗАЦИЯ МОДЕЛИ ТАКСИС - РЕАКЦИЯ - ДИФФУЗИЯ, ОПИСЫВАЮЩЕЙ ДИНАМИКУ СИСТЕМЫ ХИЩНИК - ЖЕРТВА

© 2010 г. А.Д. Загребнева1, Ю.В. Тютюнов1, Ф.А. Сурков1, А.И. Азовский

2

1Научно-исследовательский институт механики и прикладной математики Южного федерального университета, пр. Стачки, 200/1, г. Ростов н/Д, 344090, niimpm@ns.math.sfedu.ru

1Research Institute of Mechanics and Applied Mathematics of Southern Federal University, Stachki Ave, 200/1, Rostov-on-Don, 344090, niimpm@ns.math.sfedu.ru

2Московский государственный университет, Ленинские горы, 1, стр. 12, г. Москва, 119991, main@soil.msu.ru

2Moscow State University, Leninskie Gory, 1/12, Moscow, 119991, main@soil.msu.ru

Для одномерной (отрезок) и двумерной (прямоугольник) областей представлены схемы численного решения нелинейной системы дифференциальных уравнений в частных производных, описывающей вызванную трофотаксисом миграцию хищников. Переменными модели являются плотности популяции хищников и жертв, а также сытость хищников. Показано, что схемы сохраняют свойство консервативности исходной системы дифференциальных уравнений; с высокой точностью удовлетворяют условию устойчивости ненулевого стационарного режима, аналитически полученного для исходной непрерывной модели. Экспериментально продемонстрировано возникновение сложных пространственно-неоднородных структур.

Ключевые слова: разностная схема, трофотаксис, таксис-диффузия-реакция, хищник - жертва, пространственные структуры.

Difference schemes for numerical solution of partial differential equations, constructed to describe predator movements caused by tro-photaxis, are presented for one-dimensional (segment) and two-dimensional (rectangle) closed domains. Predator and prey population densities and predator satiety are the model variables. It is shown that schemes keep general integral of predator density; schemes satisfy to the nonzero equilibrium stability condition, obtained analytically for initial continuous system. The emergence of spatially heterogeneous patterns is experimentally demonstrated.

Keywords: difference scheme, trophotaxis, taxis-reaction-diffusion, predator-prey model, spatial patterns.

В [1, 2] для описания пространственно-временной динамики трофической системы веслоногие раки гар-пактициды - диатомовые водоросли предложена модель хищник - жертва, представляющая собой систему уравнений типа таксис - диффузия - реакция. После введения безразмерных величин соответствующая система уравнений, заданная в области О (в случае одномерного местообитания О = [0,1], двумерного -О = [0,1] х [0, Ь\), имеет следующий вид:

^ = Я(1 - Я)-я(( + п)ЛЯ,

дП = -V • ((S )(( + n)vs )+V • ((S )Vn),

(1)

dS

— = R - vS + 8S AS. dt S

Здесь R - плотность популяции жертв (диатомовых водорослей); (Щ = n/|q| - осредненное по пространству (области Q) значение плотности хищников (гар-пактицид) N; N = JQ Ndx - общая численность популяции хищников; n - отклонение плотности хищников N от осредненного значения (Щ, N = (ы) + n; S - средняя сытость, наполненность желудков хищников; v - коэффициент переваривания; 5r , ju(S) -коэффициенты диффузии плотности жертв и хищников соответственно; x(s) - коэффициент таксиса хищников, %(S) = - dj(S)/dS . В [1] показано, что для обеспечения положительного таксисного ответа, т.е. движения хищников по направлению к жертве, коэффициент диффузии j(S) должен удовлетворять следующим условиям: j(s) > 0 для любых S > 0 ; j(S) убывает с возрастанием сытости, т.е. dj(S)/dS < 0 .

В силу граничных условий, задающих отсутствие потоков через границу местообитания Q,

VRlx=aa=VHx=aa=VSx=aa= ^ (2)

N постоянно во времени. Таким образом, (Щ = const, а среднее по пространству значение отклонения n равно нулю, т.е. (n) = l/|Q| Jq ndx = 0 .

Предполагается, что в начальный момент заданы распределения:

Rlt=0 = *pR , nlt=0 =Vn, SL=0 =Vs .

Численная аппроксимация в случае одномерной области

Выпишем схему численного решения системы дифференциальных уравнений (1)-(3) для случая одномерной области. На отрезке О = [0,1] введем равномерную сетку {х, = 0% I = 0..м} с шагом к = 1/М; Я, (/) = Я(х,, !), щ (!) = п(х,, !), Б, (() = 5(х,, !) - значения плотностей популяции жертвы, хищников и сытости хищников в узлах сетки. Аппроксимация системы (1)-(3) во внутренних узлах сетки [3-7]:

^ = Я, (1 -< N - щ - Я, Ьдя^-^кЯ1^,

а! к

dni 1 ( Б,+1 — Б, Б,- — Б,- 1 —^ =--1 а, -- а^-^ | +

а! к I к к

+11 k - ^v-n-L h I г+L h h

^ = R,-vS, +SSS'-1 2Sj + S+1 , i = 1... M-1 _dt S h2

где коэффициенты a, и b, принимают значения:

схема I: центральная разностная [3-5],

a, = (x(Si )((W) + n, )+ x(S,-1 )((W) + n,-1 ))/2 ,

b, =(p(St )+M(St-1 ))/2;

схема II: противопотоковая разностная [6, 7], a [Ж-1 X(W) + пг-i ) S,-1 < S, ; ^ US XM + n, ) S,-i > S, ; \^(S,-1 \ nt-1 < nt ;

bi =

AS1)' n-1 ^ ni.

Граничные условия (2) приблизим центральными разностями. С учетом полученных разностных граничных условий аппроксимация системы (1) - (3) во внешних узлах сетки выглядит следующим образом:

^ = Я0 (1 -< ^ - П0 - Я0 )+ 2дКЯ^,

dt

dn0

—- =--1 aL

dt h I

h2

SL - S,

0 1+ 2 f bL nLZi h \ 1 h

dS0. = R0 - vS0 + 2^sSl S°-

it=0 rR'-"it=0 rn' ~lt=0 rs' (3)

Свойство консервативности N накладывает условия на начальное распределение n : \ayndx = 0.

Важно, чтобы при численном решении задачи (1)-(3) сохранялись свойства исходной системы, в частности, N = const.

В статье предложена аппроксимация нелинейных членов в уравнении для n, позволяющая сохранять во

времени величину N . Рассмотрены стандартные схемы решения [3-6]. Также использовался опыт работы [7], в которой исследуются волновые явления в системах таксис - реакция - диффузия.

dt

dRM dt

S h 2

- = RM (l - N - nM - RM ) + 2SrRm-1- Rm

dn

M 2 ( SM-1 SM

dt

+ 2 b

dS

M

dt

= RM -vSM + 28S

h ) h\

SM -1 - SM

h1

Утверждение 1. Схемы I и II сохраняют свойство

консервативности N исходной системы дифференциальных уравнений (1)-(3).

Доказательство. Для схем I и II N = const, если d (jQ Ndx) dt = 0. Проверим, что это условие выполняется. Аппроксимируем интеграл по формуле трапеций и учитываем, что Nt = (n} + nt. Тогда

0

h

n

n

M -1

h

h

3-1 J N (x, t )dx Ufjz

dt

d f MN + N -

dt

i =1

2

h I =

= h

M-1

= Z

i=1

M-1 dn, 1 f dn0 dnM Y

1 dt 2 i dt dt

f

_ "i+1

^ Si - ai ——>Si-1 I +

+ I bi

-I a.

+1 a.

ni+1 - ni

- bi

ni - ni-1

h

S1 - So | + fb n1 - no

h

sm sm-1 I _ J ь nM nM-1 I = o

h

Утверждение доказано.

Численная аппроксимация в случае двумерной области

Выпишем схему численного решения системы дифференциальных уравнений (1) - (3) для случая прямоугольной области О. На прямоугольнике О = [0,1] х [0, Ь] введем равномерную сетку { (,Уi) = К;,), г = 0...ых , у = 0...Му } с шагом кх = 1/Мх и ку = ЬМУ по осям X, Г. Яг у (г) =

= К(хг , Уг , ^), "г,у (() = , Уг , I), $,у (() = 4х, у, /) - значения плотности жертвы, хищников и сытости хищников в узлах сетки. Аппроксимация системы (1) - (3) во внутренних узлах сетки для г = 1. Мх —1,

у = 1.М., —1, имеет вид

^ у

'Ж

% /

dt

+ SR

= R

i,j (-<N - ni,j -Ri,j )+

f Ri+1, j- 2R/, j + Ri-1, j +R/, j+1 - 2R/, j +R j-1 ^

h

2

h

dt h*

2

у у

i

dni,J 1 f x Si+1,J - Si,J x Si,J - Si a*i ,-----a^j 'J

'i+1, j

Л+1, j

4 j+1

dS,

i, j

^+1

= R, , -vS, , +

К

ni+1,j - ni, j

к

Si, J+1 ~ "Si, j

hv

y

ni,/+1 - ni, J

hx

n -n

i, j

-a

Si-, j - Si', /-1

i, j

by.

i, j

ni, j - n

y У I

У У

dt i,j i,j

Si+1, j - 2Si', j + Si-1, j Si, j+1 - 2Si, j + Si, j-1

+ S

h

2

-+-

h

2

y

at,j = 0,5(^(Si, / ^ N + п,, /)+ Z(S^ j ( + n^, j )),

b*j = 0,5(4—, j )+ß—-1, j )),

a(j = 0,5(x(Si, j X( N + n,, /)+ M,/-1 X(N^> + nu)),

ЬУ/ = 0,5ß(—, j )+ß—,/-!)); схема IV: противопотоковая разностная [6, 7],

a* = |x(si-1,j tn) + i Si-1,j < Si,j,

a'J [ж(,/ + n.,j), Si-!,/ > Si,/;

f^(si-1,j ), ni-1,j < ni,j ,

[d—i,j ni-1,j > ni,j ; fx(Si,/-1 +/), Si,/-i < Si,/,

[х((,,/ + nh}), Si, j-1 > Sj;

id(si,J-1} ni,J-1 < ni,J,

biX/

4/ =

ß

j}

ni,/-1 > ni, / •

Аналогично случаю одномерной области Q, для аппроксимации системы (1) - (3) во внешних узлах сетки граничные условия (2) приблизим центральными разностями.

Утверждение 2. Схемы III и IV сохраняют свойство консервативности N исходной системы дифференциальных уравнений (1) - (3).

Доказательство. Аналогично доказательству утверждения 1.

Вычислительные эксперименты

С применением описанных выше аппроксимаций проведем ряд вычислительных экспериментов. Полученные системы (М +1) и (Mx + l) (My +1) обыкновенных дифференциальных уравнений проинтегрируем методом Рунге - Кутта 5-го порядка, улучшенным Ричардсоном [8], с контролем точности на шаге RKF45. Вычисления проведем при следующих значениях параметров и начальных данных: v = 0,1, SR = 0,005, ^(S) = e~3S,

S— = 1-10, e = 110-

(4)

где коэффициенты , afj, b*j, bfj принимают значения: схема III: центральная разностная [3-5],

одномерная область Q:

N = 0,82, M = 200, cpR = 1 - (N) - 0,01 cos(^x),

Vn = 0, ер— =(1 - N)/v; (5)

двумерная:

L = 4, N = 0,7975, Mx = 100, My = 300, pR = 1 N - 0,001cos(®c)cos(2^y) + 0,001cos(3®c)cos(2^y), рп = 0, р— =(1 - N )/v. (6)

Здесь e - погрешность метода Рунге - Кутта. Убедимся численно, что для представленных схем I - IV N = const. Так как N = (N^ + n , то это условие выполняется, если (n) = 0. Вычислим максимальное за период времени t = [600,800] значение (n): для схемы I оно равно 2,56 -10-14 ; II - 2,01-10-14; III - 4,95 -10-22; IV - 9,15 -10-21. Данные результаты

и

h

h

+

h

+

9

1

+

h

h

x

x

1

+

h

у

1

+

h

h

y

y

свидетельствует о том, что используемые схемы сохраняют N с высокой точностью.

В [2] проведен линейный анализ устойчивости однородного нетривиального стационарного равновесия системы (1)-(3) (Я*, п*, S*) = (1 -(Ы),0, (1 -(Ы))/у), О = [0,1]. Получено условие его устойчивости:

г(к2,(N)9SR,M(S* }Ss,x(s*) v)> 0

Vk 2 m 2.

(7)

где

72 (( 2,( N 8 , А* (= к3 ^ (5 *))+ + к 2 (я * +ф2 (*)+ 28К а(^2* )+

+ к^Я* +V )-(ы)Я*х( * )+8яу(+2Я*))+(я* +у)я\-

т - номер моды разложения в ряд Фурье малого пространственного возмущения стационарного режима (Я*, п*, 5*). Показано, что при потере устойчивости режима (Я*, п*, 5*) имеет место бифуркация Пуанка-ре-Андронова-Хопфа: в окрестности теряющего устойчивость однородного стационарного режима рождается периодический режим.

Проверим, как согласуются результаты вычислительных экспериментов с условием устойчивости (7).

В системе (1)-(3) с параметрами при изменении (ы)

первой возбуждается мода т = 1, т.е. условие (7) для параметров (4) впервые нарушается при т = 1. График значений функции Т2 для параметров (4) и т = 1 изображен на рис. 1: при {ы) <(N1, как и при (ы) > (ы)2, однородный режим устойчивый; при N1 <(Ы) <(Ы)2 - неустойчивый; (ы) 1 - 0,80057, (ы) 2 - 0,91403 являются бифуркационными (потеря устойчивости).

Рис. 1. График зависимости от функции Т2 для моды т = 1 и параметров (4)

В вычислительных экспериментах со схемой I и сеткой с 201 узлами бифуркационные значения

смещены примерно на 3-10"4: при (ы) 1 - 0,80070 режим теряет устойчивость; при 2 - 0,91438 становится устойчивым. На рис. 2 представлены проекции

решения системы (1) - (3) с параметрами (4) и близкими к бифуркационным 1 и (Ы^2 значениями

параметра (ы) на фазовую плоскость Я и (ЯЫ^, где ^Я) = (1/|п|){О Я(к - средняя плотность жертвы; (ЯЫ} = (1/|п|){О ЯЫ(х - среднее количество съеденных жертв.

0.1601

0.1598

0.1595

0.1592

А

X

0,199

0.1994

0.1998

0.2002

СЮ

Рис. 2. Проекция на фазовую плоскость (^ЯЫ^, (Я^) траектории решения системы (1)-(3)

Когда выполняется условие (7), стационарный режим (Я*, п*, 5*) является устойчивым (рис. 2, режимы А, Е). В области неустойчивости при значениях параметра (ы), близких к (ы) 1 и (ы) 2, наблюдаются

периодические режимы, соответствующие единственной неустойчивой первой моде (рис. 2, режимы Б, Г), возникающие в окрестности неустойчивого стационарного режима (Я*, п*, 5*) (рис. 2, режимы В, Д).

В случае двумерной прямоугольной области О = X х У = [0,1] х [0, ь] условие устойчивости стационарного режима (Я*, п*, 5*) к2 =п2 (т2 + р2/1?), где т, р жения в ряд Фурье малого пространственного воз

аналогично (7); номер моды разло-

0 25

0.15

У 4

0.2028 з

02026

с

02024 1

0.2022

1 0

У 4

мущения стационарного режима (Я*, п*, Б*) по пространственным переменным х е X и у е У . Вид наблюдаемых пространственных структур во многом зависит от соотношения длин сторон области О . Так, в прямоугольнике с длиной Ь = 1 в системе (1) - (3) с параметрами (4), согласно условию устойчивости (7), при N и 0,80057 первой теряет устойчивость квазиодномерная мода максимального пространственного масштаба (т = 0, р=1 или т=1, р=0). На рис. 3а изображено характерное распределение переменных модели системы (1) - (3), соответствующее единственной неустойчивой первой квазиодномерной моде т=1, р=0. Если Ь=4, то в системе (1) - (3) с параметрами (4), согласно условию устойчивости (7), при N и 0,79728

первой возбуждается истинно двумерная мода, для которой т=1, р = 2.

На рис. 3б изображено соответствующее распределение переменных модели, полученное для и 0,79750 , близкого к бифуркационному значению параметра . При увеличении параметра

наблюдается усложнение пространственно-неоднородной динамики (рис. 3в).

Выводы

Представлены схемы численного решения системы дифференциальных уравнений (1) - (3) для случаев одномерной и двумерной областей О, сохраняющие равновесие и свойство консервативности общей численности популяции хищников в модели трофической системы гарпактициды - диатомовые водоросли. Для схем с высокой точностью выполняется условие устойчивости ненулевого стационарного режима, аналитически полученного для исходной системы (1) - (3). Экспериментально показано, что модель (1) - (3) демонстрирует сложную пространственно-неоднородную динамику.

б

У 4

0.798

0.7978 3

0 7976

2

0.7974

0.7972 1

0.797

— У 4

0 22

0.2 3

0.13 2

0.16

1

I 0.14

■ 0

Mr = My, = 110,

Рис. 3. Пространственное распределение плотностей Я, N и сытости Б в модели (1) - (3) с параметрами (4) и начальными данными (6) для двумерной области: а - Ь=1, = 0,82,

г = 696,15 ; б - Ь=4, N = 0,7975, Мх = 100, Му = 300, г = 799,63;

в - Ь=4, N = 0,82, Мх = 100, Му = 300, г = 486,93

Литература

1. Моделирование уравнения потока популяционной плотности организмов с периодическими миграциями / Ю.В. Тютюнов [и др.] // Океанология. 2010. T. 52, № 1. С. 1-10.

2. Микромасштабная пятнистость распределения веслоногих рачков как результат трофически-обусловленных миграций / Ю.В. Тютюнов [и др.] // Биофизика. 2009. Т. 54, № 3. С. 508-514.

3. Самарский А.А. Введение в численные методы: учеб. пособие. М., 1987. C. 286.

4. Самарский А.А., Гулин А.В. Численные методы. М., 1989. C. 432.

5. Роуч П. Вычислительная гидродинамика. М., 1980. C. 612.

6. Morton K.W., Mayers D.F. Numerical Solution of Partial Differential Equations. Cambridge, 1994. Р. 227.

7. Soliton-like phenomena in one-dimensional cross-diffusion systems: a predator-prey pursuit and evasion example / M.A. Tsyganov [et al.] // Physica D. 2004. Vol. 197. P. 18-33.

8. Каханер Д., Моулер К., Неш С. Численные методы и программное обеспечение. М., 1998. C. 575.

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

22 июня 2009 г.

в