Численный алгоритм склеивания решения на поверхности разрыва для систем с клеточной структурой Текст научной статьи по специальности «Математика»

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

Вестн. Ом. ун-та. 2010. №4. С. 156-163.

УДК 519.62

В.В. Коробицын, В.Б. Маренич, Ю.В. Фролова

Омский государственный университет им. Ф. М. Достоевского

ЧИСЛЕННЫЙ АЛГОРИТМ СКЛЕИВАНИЯ РЕШЕНИЯ

НА ПОВЕРХНОСТИ РАЗРЫВА

ДЛЯ СИСТЕМ С КЛЕТОЧНОЙ СТРУКТУРОЙ

Для систем обыкновенных дифференциальных уравнений предлагается алгоритм для нахождения точки пересечения траектории решения с поверхностью разрыва правой части. Вводится понятие динамической системы с клеточной структурой. Доказывается теорема об аппроксимации решения системы схемой склеивания решений на границе между клетками. Приводятся примеры систем с клеточной структурой и их фазовые портреты, построенные с помощью программы МЕР2, разработанной авторами статьи.

Ключевые слова: численный алгоритм, обыкновенное дифференциальное уравнение, система с клеточной структурой.

Введение

Математические модели, описывающие динамическое поведение системы во времени с изменением качественного поведения, зачастую используют аппарат дифференциальных уравнений с разрывной правой частью. Качественное изменение поведения системы обусловлено изменением самой структуры системы, ее подстройки к внешнему воздействию или сочетанию внутренних состояний. Например, поведение некоторых экологических симбиозов построено на изменении поведения одной или нескольких популяций. Так, некоторые особи способны быть либо травоядными, либо хищниками в зависимости от содержания окружающего корма. Такое качественное изменение поведения приводит к выбору разрывных функций при описании скорости изменения биомассы популяций. Не удается подобрать непрерывные функции для описания изменения скорости компонент системы. А это означает, что в правой части системы дифференциальных уравнений появляются разрывные функции.

Непрерывные модели значительно легче исследовать, но системы с разрывными функциями, как правило, имеют совершенно иное качественное поведение. А это требует использовать новые подходы при построении и анализе моделей. В системах с разрывами возникает много трудностей: теоретическое обоснование разрешимости задачи, построение специальных численных методов для расчета на компьютере, согласование теоретических и вычислительных результатов. Классической монографией по аналитическому исследованию дифференциальных уравнений с разрывной правой частью стала монография А.Ф. Филиппова [1]. В его работе приводятся обзор характерных изменений в поведении разрывных систем, теоретическое обоснова-

© В.В. Коробицын, В.Б. Маренич, Ю.В. Фролова, 2010

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

Развитие компьютерной техники привело к тому, что исследование моделей не обходится без вычислительного эксперимента. Любая модель проходит апробацию на компьютере. А раз так, то для каждой модели должны быть разработаны численные методы решения. Это касается и исследования систем дифференциальных уравнений с разрывной правой частью. Некоторые аспекты разработки таких методов рассматриваются в работах [2-4].

Разрабатываемый авторами проект МЕР2 (Modelling Evolution Processes 2) направлен на создание программного обеспечения для анализа динамических систем с разрывной правой частью. Модели процессов различной природы описываются подобными динамическими системами: перенос биомассы в цепи питания организмов (растительность - травоядные

- хищники); изменение экономических индексов в различных ситуациях; поведение индивидов в социальной среде; функционирование электрических схем с нелинейными элементами; сложные кинематические механизмы.

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

Динамические системы с клеточной структурой

Будем рассматривать системы вида

динамические

— = F(X) dt

(1)

с разрывной правой частью Т7. Решение таких систем не может быть представлено в классическом понимании как дифференцируемая кривая Х@, скорость которой равна заданному векторному полю К Иногда является естественным принять

некоторое Х(1) как решение системы (1), пренебрегая тем, что с1Х I ф Р(Х), например, в скользящем режиме.

Пусть О! - некоторая открытая область в 7?" и ЩХ), X е О.' разрывное векторное поле в О.'. Обозначим через С/. и £>/■ подмножества точек О', где Т7 является непрерывной и разрывной соответственно. Поскольку нас не интересует поведение траектории динамической системы вблизи границы области, то мы будем рассматривать систему (1) в некоторой подобласти О. с О.', замыкание О. которой является компактом и находится внутри О.'.

Если динамическая система составлена из конечного числа «хороших» систем, т. е. множество Д|. состоит из конечного числа кусочно-гладких гиперповерхностей, тогда будем говорить о динамической системе с клеточной структурой.

Определение 1. Множество С будем называть клеткой динамической системы (1), если С является линейно-связным подмножеством множества Ср.

Для всех кусочно-гладких систем множество С/, может быть представлено как объединение конечного числа клеток

О, I = 1, ..., п.

Определение 2. Динамическая система (1) имеет клеточную структуру Се,

если объединение замыканий клеток Сг, 1=1,..., п совпадает с замыканием области определения О. функции Щ1,Х).

Отметим, что для динамической системы с клеточной структурой подмножество £>/■ точек, где функция Щ1,Х) разрывна, составлена из точек замыкания всех

П П

клеток: /),

;=1 ;=1

Определение 3. Функцию Х(Ц, а < t < Ъ будем называть решением динамической системы (1) в строго обобщенном смысле, если для любого гладкого отображения ¥(1):[а,Ь]—> Е", обращающегося в нуль на концах У (а) = У(Ъ) = О, имеем ъ ъ

[ (Г ^ ),Х(Г ))Л = -{ ),Р(Х^ )))Л.

а а

Определение 4. Векторное поле Р(Х) допускает («^-направление (или направленное векторное поле) V в области О., если для любого Хое О. существует единичный вектор У(Хо) (направление) такой, что для любого X, удовлетворяющего условию

\Х-Хо\ < а, выполняется неравенство (Р(Х), ЩХо)) > Р-

Отметим, что если Т7 является непрерывным в точке Хо и не равно нулю, то оно допускает направление в Хо, для которого мы можем принять само направление ЩХо). Если Т7 - равномерно непрерывное поле на некотором замкнутом множестве Б с О., тогда Т7 само является направленным векторным полем для Т7 для некоторых (а,Р), и, кроме того, оптимальное значение р зависит от а, так что Да)—» 1 при а—>О.

Таким образом, выбор направления ]/(Хо) отличного от ЩХо) является естественным только в случае, когда Хо е Ор, т. е. в точке разрыва функции Т7. Для таких точек определим направление сНг(.Р,Х;)

- множество единичных векторов всех возможных пределов ЩХп) при Хп—>Хо. Если £(Хо) - единичная сфера Э1-1 является множеством всех единичных направлений в точке Хо, то 6Ь(Р,Х()) с £(Хо). Обозначим через У’(Хо) выпуклую оболочку направлений б1г(Р,Хо) в Б (Ха). Тогда направленное векторное поле для Т7 существует тогда и только тогда, когда выпуклая оболочка У’(Хо) полностью содержится в некоторой открытой полусфере 3+"~! в £(Хо). Если V является центром ЭТОЙ полусферы в+н-1, т. е. для любого вектора Ш из 8+"~1 угол между Ши Vбудет меньше л/2, тогда Vбудет естественным представлением для У(Хо).

Определение 5. Непрерывная функция Х(1), а < t < Ъ, удовлетворяющая условию Липшица называется решением динамической системы (1) в слабом обобщенном смысле, если эта функция всюду является решением дифференциального

включения Х+(1) еУ(Х(1)), гдеХ+(/) -множество правых производных непрерывной липшицевой функции Х(Ц.

Будет правильным называть решение в слабом обобщенном смысле «скользящим режимом», если оно не является решением в строго обобщенном смысле, следуя работе [1].

Легко видеть, что для динамической системы с непрерывной правой частью Т7 все решения Х(1) в слабом (или строгом) обобщенном смысле будут на самом деле решениями в обычном смысле. Это оправдывает предложенные определения 3 и 5 для решений в строгом или слабом

обобщенном смысле для систем с разрывной правой частью.

Общая схема алгоритма численного решения динамических систем с клеточной структурой

Поиск численного решения динамической системы с клеточной структурой состоит из двух этапов. Во-первых, для любой клетки С,, 1 = 1,..., п численное решение может быть найдено классическими методами Рунге-Кутта. Во-вторых, решение должно быть склеено на множестве Ор, где динамическая система не имеет гладкого решения. Для нахождения численного решения в области непрерывности функции Т7 можно воспользоваться методом, предложенным Дорманом и Принсом [7]. Этот метод имеет 5-й порядок точности и обладает механизмом оценки погрешности решения. Для нахождения точки пересечения решения с поверхностью разрыва функции Т7 можно воспользоваться непрерывным расширением метода Дормана и Принса [8, с. 191]. Воспользовавшись этой схемой можно без дополнительных вычислений правой части Т7 найти точку пересечения непрерывного решения с поверхностью разрыва с заданной точностью, применив итерационный метод решения нелинейного алгебраического уравнения.

Пусть найденное решение Х(Ц лежит в клетке С;, а значение в некоторый последующий момент времени Х*{Ь + /г) принадлежит другой клетке С/, 7 ф I. Полученное значение Х*^ + /г) принять за приемлемое решение в общем случае нельзя, поскольку на границе между клетками величина погрешности может значительно вырасти. Мы предлагаем следующий алгоритм склеивания решения на границе между клетками.

Алгоритм (поиска решения на границе между клетками)

Шаг 1 (идти к границе). Применим метод дихотомии для нахождения значения параметра в такого, что точка Х1 = Х(и), ^= 1+вН будет находиться в клетке С; и расстояние до границы будет меньше, чем е | |Х1| |, т. е. | |Х1- Хъ || < е-1 | Х1 | |, где е - заданная локальная точность, Хъ е Ор- некоторая точка на границе клеток. Тогда будем говорить, что точка Х\ является точкой приближенного решения на границе.

Шаг 2 (скользящий режим на границе). Если выполняется условие (.Р^Х]), п(Хь)) • {п(Хь), ЩХ2)) > О, то перейти на шаг 3, иначе реализуется скользящий режим на поверхности разрыва. Воспользуемся проекцией ртогЩ!) функции Т7 на поверхность разрыва Л>/. для решения задачи Коши на поверхности разрыва. Находим решение до тех пор, пока кривая решения не сойдет с поверхности разрыва, т. е. не нарушится указанное выше условие. (Здесь п(Хь) - вектор нормали к поверхности разрыва функции Р(Х) в точке Хъ. Предложенное условие означает, что вектора ЩХ1) и ЩХ) направлены в сторону пересечения с поверхностью разрыва, а это означает, что решение не покинет поверхности разрыва).

Шаг 3 (идти в клетку С,). Пусть {Щ -итерационная последовательность значений параметров в при выполнении шага

1, О,- - параметр остановки итерационной процедуры. Тогда точка Х^ + О,--/Н) е С,. Выберем шаг интегрирования /г2 = (вг-1-вг)Н и найдем решение в клетке С/ с помощью схемы:

Кг = ЩХг),

К = Щ1+ Ьг К1), (2)

Х(и+ /г2) и Х2 = XI + (Кх+ Щ • /г2/2.

Здесь Х\ + Н2К\ представляет первое приближение точки Х2. Среднее значение скорости кривой в точке Хъ представляется вектором V = (К] + Ко)/2, которое соответствует естественному выбору вектора V, приведенного выше.

Далее докажем теорему об аппроксимации решения системы (1) схемой (2) на границе клеток. Введем некоторые обозначения. Пусть функция X (I) представляет точное решение динамической системы (1). Рассмотрим решение вблизи границы между клетками С; и С,. Будем считать, что решение Х{!) в граничной точке Хъ является непрерывным, но негладким. Кроме того, будем считать, что точка Хъ не является стационарной (| | ЩХъ) | | > О). Обозначим Хх = Х(^)е Ci,

Х2 = Х(12 ) & С, Ы= £2 - и, а < и< ^ < Ъ. Тогда граничная точка определяется формулой Хь = Х(^ + вкг), ве (0,1). Пусть точка Х2* - приближенное решение, полученное с помощью численной схемы (2). Обозначим отклонение от точного решения ДХ2* = Х2 - Х2*.

Теорема. Пусть векторное поле Р(Х) ограничено: | |Р(Х)\\< М для всех Хе А, a<t< Ь, дифференцируемо по X в некоторой окрестности точки Х2, и производная

в точке Х2 ограничена

— Ь < 00.

Тогда отклонение приближенного решения Х2* от точного значения Х2 определяется формулой ДХ2* = И • К, + О(Нг), где Я<2М(<9-1/2).

Доказательство.

Пусть Хь = Нт Х(1] +тк{)

и

XI = Пт X(/, + г/?(). Поскольку решение

т->д+ О

Х(Ц является гладким в клетках С;, С/, то значения Хь" и Хь+ можно представить в виде

х-ь=хх

■вк^іХ^ + О^),

(3)

х;=х2-(1 -в)ь1Р(х2)+о(^). (4)

Так как решение Х(/) является непрерывным в точке Хь, то эти значения совпадают Хъ~ = Хъ+ = Хъ. Приравняем правые части выражений (3) и (4), получим уравнение для Х2:

Х2=Х1+(вР(Х1) +

+(1 -в)Р(Х2Щ+Оф.

Для приближенного решения этого уравнения применим схему (2) и оценим погрешность найденного решения. Вычислим приближенное значение для Х2 по

формуле X= X! + к( 1< ( X /) (аргумент функции во второй формуле схемы (2)). Найдем отклонение приближенного значения Х2 от точного решения Х2:

Ах 2 =(1- 0 А (~Р(Х1) + Р(Х2 )) + 0(И? ).

Пусть Х2* - это приближенное значение точки Х2, полученное в итоге применения всех формул схемы (2) и вычислим отклонение точки Х2* от точного значения

х2:

АХ*2=(в-\/2ЩР(Х1) +

+к{((1-в)Р(Х2)-Р(Х2)/2)+0(И?).

Погрешность вычисления функции Р в точке X2 можно оценить, используя формулу Тейлора. Для этого представим функцию Р(Х) в окрестности точки Х2 в виде

Р(X2 ) — Р(X2 ) = ЬАХ2 + О АХ

м2

-2

Это сделать можно, поскольку в этой окрестности функция дифференцируема и ее производная ограничена. Кроме того,

очевидно, что ОЦАУоЦ у = 0(}7,'). В итоге

получаем оценку погрешности численного решения Х-2'.

ДХ2* = А, ({в -1/2) + {в - 1)Щ/2)(Р(Х1)~

-Р(Х2)) + 0(И;).

Теперь мы можем записать

ах *2 = д-л, +о(и/

где И= (0-1/2) • {Р[Х\) - ЩХ2)). Поскольку функция Р(Х) ограничена числом М, то | | ^Щ) - Щз) | | < 2М, поэтому И < 2М (0 --1/2). Теорема доказана.

Замечание. Если в условиях теоремы предположить, что шаг Уц удовлетворяет условию 6-1/2<1ъ. Тогда схема (2) аппроксимирует решение со вторым порядком точности АХ/- = О (Кг).

Программа МБР2 для исследования динамических систем с клеточной структурой

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

• ввод системы с клеточной структурой в графическом режиме;

• сохранение и загрузка файла данных системы;

• вывод графиков решений системы;

• построение проекций фазового пространства на плоскость.

Рождение цикла в линейной системе с клеточной структурой. Рассмотрим систему линейных обыкновенных дифференциальных уравнений с клеточной структурой, составленной из двух клеток Б г. х< хь и Д>: х> хо на плоскости (х,у):

= а]](х-с]) + а]2(у-с1]),

= а21(х-с1) + а22(у^1), (x,y)eD1 = a22(x-C2) + a12(y-d2),

(5)

dx

dt

4у dt

dx

dt

dУ = a2l(x-C2) + all(y-d2), (x,y)&D2. dt

В случае, когда ап = аз2 = О, 012 = 021 = 1 решение каждой из систем имеет одну особую точку (седло). Когда С1 < хо и 02 > хо, особые точки лежат в соответствующих клетках Р\ ж ГЬ,- Кроме того, когда с1\ = с12, особые точки лежат на прямой ортогональной границе между клетками, а их сепаратрисы пересекаются на границе клеток. В этом случае образуется замкнутая область Ц), заключенная между сепаратрисами. Решение внутри этой области образует цикл и не выходит за ее пределы. Цикл составлен из двух гладких кривых, склеенных на границе между клетками (рис. 1).

Рис. 1. Фазовый портрет линейной клеточной системы (5)

Как можно видеть на рис. 1, траектории, заключенные между сепаратрисами на границе клеток Б\ и ГЬ,, образуют цикл. Проведенный вычислительный эксперимент показал устойчивость численной схемы на этом цикле - решение не уходило с циклической траектории. Однако отмечалось смещение вдоль оси времени, которое выражалось в том, что численное решение имело более короткий временной период, чем точное решение. Это можно объяснить погрешностью возникающей при вычислении точки пересечения границы клеток. Каждый раз происходит переход из одной клетки в другую с некоторой величиной перескока, причем всегда в большую сторону. Тем не менее получен-

ный фазовый портрет (рис. 1) вполне соответствует точному решению системы (5).

Модель трехуровневой цепи питания. Для демонстрации возможностей программы МЕР2 представим результаты исследования модели трехуровневой цепи питания [10; 11]. Эта модель описывает изменение биомассы растительности (Р), травоядных (С) и хищников (Р}. Хищники могут питаться как травоядными, так и растительностью. Динамика популяций описывается системой дифференциальных уравнений:

—=к

А

(

R

лксс

ч I К) 1 + кксАкс11

__________икр^црР____________

1 ы рр ^ррНррК ыСрАСРЬСРС

— = с л

екс^ксК 1 + ИксАкс11

иср^срР

1 + иррХрркррК + иср1срксрс

~тс

*Р = Р

А

иНР^НРеНрВ- + иСР^СРеСР^

1 + ЫррХрркррЯ + иСрХСркСрС

(6)

-тР

Определения и значения параметров модели представлены в табл.

Определения параметров и их значения, исследуемой модели (6)

Параметр Определение Значение, единицы

г Максимальная скорость роста растительности 0,3 час'1

К Вместимость растительности в окружающей среде 10 мг С/л

^РС Скорость поиска растительности травоядным 0,037 л/(мг • час)

Арр Скорость поиска растительности хищником 0,025 л/(мг • час)

Лер Скорость поиска травоядного хищником 0,025 л/(мг • час)

ирр Вероятность нападения хищников на растительность 1 (безразмерна)

иср Вероятность нападения хищников на травоядных 1 (безразмерна)

hрс Время, необходимое травоядным для усвоения растительности 3 час

hрр Время, необходимое хищникам для усвоения растительности 4 час

hср Время, необходимое хищникам для усвоения травоядных 4 час

ерс Эффективность усвоения биомассы растительности для травоядных 0,6 (безразмерна)

ерр Эффективность усвоения биомассы растительности для хищников 0,36 (безразмерна)

еср Эффективность усвоения биомассы травоядного для хищников 0,6 (безразмерна)

тс Скорость отмирания травоядных 0,03 час'1

тр Скорость отмирания хищников 0,0275 час'1

Результаты моделирования системы (6) приведены на рис. 2 и 3. Представлены траектории решений в трехмерном фазовом пространстве (И, С, Р), спроецированные на плоскости (Р, С) ж (Р, С). Начальные значения траекторий выбраны произвольно из допустимой области определения.

Согласно модели В. Кривана и С. Дейла хищник выбирает свою диету по следующему правилу. Возможно два варианта питания: 1) растительность более продуктивна для питания хищников, чем травоядные; 2) травоядные более продуктивны для питания хищников, чем растительность. В первом случае имеет место неравенство

>■

к к • *' о о * 1 г

(7)

"СР

Оптимальная стратегия питания хищников состоит в поглощении растительности (к/.:/'= 1), травоядные же включены в диету хищников только тогда, когда плотность растительности меньше переключательной плотности Р.

Топология цепи питания не является фиксированной, она меняется от исключительно конкурентной цепи (иыр= 1; иср = 0 в модели (6)) до цепи всеядного хищника (ирр=иср= 1). Изменение стратегии происходит при переходе через переключательную плотность растительности Р* (рис. 4 (А)).

В случае, когда травоядные являются более питательными для хищников, чем растительность, неравенство (7) инвертируется. Когда плотность травоядных меньше переключательной плотности С хищники выбирают стратегию всеядного (ирр=иср= 1), а когда больше переключательной плотности, то они потребляют только травоядных (ирр= 0; иср= 1). Таким образом, топология цепи питания не является фиксированной и переключается между цепью питания «настоящего» хищника и цепью питания всеядного хищника, зависящей от плотности травоядных (рис. 4(В)).

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

Также было исследовано изменение поведения численного решения с учетом и без учета разрыва правой части системы. Траектории численного решения системы

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

Рис. 2. Проекции фазового портрета модели (6) в первом режиме на плоскости (Я?, С) и (Я?, Р). Параметры системы: еКР = 0,2, еСр= 0,1, Я?* = 10, остальные указаны в табл.

Р к

Рис. 3. Проекции фазового портрета модели (6) во втором режиме на плоскости (Р, С) и (Я?, С). Параметры системы: еяр= 0,1, еср = 0,3, С* = 5, остальные указаны в табл.

Рис. 4. Схема смены стратегий поведения адаптивных хищников

методом Дормана-Принса: 1 - без учета разрыва, 2 - с использованием предложенного алгоритма. Параметры системы: еяр= 0,1, еср= 0,3, С* = 5, остальные указаны в табл.

Начальные значения: Я?(0) = 15, С(0) = 1, Р{0) = 7

Заключение

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

Полученные фазовые портреты подтверждают целесообразность введения клеточной структуры в динамической системе (6). Решения системы имеют выраженные особенности на границе между клетками. Заранее определенное геометрическое расположение клеток значительно уменьшает затраты на исследование системы. В частности, предложенный численный алгоритм поиска решения на границе клеток дает более точное решение, чем традиционные методы Рунге-Кутты (в частности, метод Дормана-Принса), что продемонстрировано на графике решений системы двумя методами (рис. 5).

ЛИТЕРАТУРА

[1] Филиппов А. Ф. Дифференциальные уравнения с разрывной правой частью. М. : Наука,

1985.

[2] Gear C. W, Osterby O. Solving ordinary differential equations with discontinuities // ACM Trans. Math. Software. 1984. Vol. 10. P. 23-44.

[3] Piiroinen P. T., Kuznetsov Yu. A. An event-driven method to simulate Filippov systems with accurate computing of sliding motions // ACM Trans. Math. Software. 2008. Vol. 34. № 13. P. 1-24.

[4] Shampine L. F., Thompson S. Event location for ordinary differential equations // Computer and Mathematics with Application. 2000. Vol. 39. P. 43-54.

[5] Коробии^т B. S., Фролова Ю. В., Маренич В. Б. Алгоритм численного решения кусочно-сшитых систем // Вычисл. технологии. 2008. Т. 13. № 2. С. 70-81.

[6] Коробицын В. В., Фролова Ю. В. Алгоритм вычисления скользящего режима для системы с гладкой границей разрыва // Вычисл. технологии. 2010. Т. 15. № 2. С. 56-72.

[7] Dormand J. R., Prince P. J. A family of embedded Runge-Kutta formulae // J. Comp. Appl. Math. 1980. Vol. 6. P. 19-26.

[8] Хайрер Э., Hepcemm С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М. : Мир, 1990.

[9] Коробицын В.В. МЕР2 - программа для моделирования эволюционных процессов. URL: http://users.univer.omsk.su/~korobits/mep2 (дата обращения: 01.06.2010).

[10] Krivan V. Optimal Intraguild Foraging and Population Stability // Theoretical Population Biology. 2000. Vol. 58. P. 79-94.

[11] Krivan V., Diehl S. Adaptive omnivory and species coexistance in tri-trophic food webs // Theoretical Population Biology. 2005. Vol. 67. P. 85-99.