Алгоритм численного решения дифференциальных уравнений с~разрывной правой частью

Математические структуры и моделирование 2005, вып. 15, с. 46-54

УДК 519.62

АЛГОРИТМ ЧИСЛЕННОГО РЕШЕНИЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С РАЗРЫВНОЙ ПРАВОЙ ЧАСТЬЮ

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

Для систем обыкновенных дифференциальных уравнений X

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

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

• перенос биомассы в цепи питания организмов (растительность - травоядные - хищники);

• изменение экономических индексов в различных ситуациях;

• поведение индивидов в социальной среде.

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

1. Динамическая система с клетчатой структурой

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

X

Copyright © 2005 В.В. Коробицын, Ю.В. Фролова.

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

E-mail: korobits@univer.omsk.su, frolova@univer.omsk.su

Математические структуры и моделирование. 2005. Вып. 15.

47

с разрывной правой частью F. Решение таких систем не может быть представлено в классическом понимании как дифференцируемая кривая Х(£), скорость которой равна заданному векторному полю F. Иногда является естественным принять некоторое X(t) как решение системы (1), пренебрегая тем, что X(t) ф F(X(t)), например скользящий режим в [1].

Пусть — некоторая открытая область в Rn и F(X),X Є разрывное векторное поле в Df. Обозначим через Ср и Dp подмножества точек Г2', где F является непрерывной и разрывной соответственно. Поскольку нас не интересует поведение траектории динамической системы вблизи границы области, то мы будем рассматривать систему (1) в некоторой подобласти Q С ГУ, замыкание Q которой является компактом и находится внутри ГУ.

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

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

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

п

Cf = Ц) Сі.

і=1

Определение 2. Динамическая система X = F(t,X) имеет клетчатую структуру CS) если объединение замыканий клеток (Д, і = 1,..., гг совпадает с замыканием области определения Д функции F(t,X):

п

Cs: |jQ = a

і—1

Отметим, что для динамической системы с клетчатой структурой подмножество Dp точек, где функция F{t,X) разрывна, составлена из точек замыкания всех клеток:

п п

І—1 і—1

2. Алгоритм численного решения динамических систем с клетчатой структурой

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

48

В.В. Коробицын, Ю.В. Фролова Алгоритм численного решения...

Численная схема метода Рунге-Кутта 4(5) имеет вид [2]

Кг = F(t,X(t)),

К2 = F(t + h/2, X(t) + КгІг/2),

К3 = F(t + h/2, X(t) + K2h/2),

K4 = F(t + h,X(t) + K3h),

X(t + h) = X(t) + (Кг + 2 K2 + 2K3 + K4)h/6.

Здесь X(t) — известное значение функции в момент времения £, X(t + К) — неизвестное значение в момент t + /і, h — шаг интегрирования. Этот метод имеет четвертый порядок точности, является устойчивым и сходящимся.

Используя оценку Ричардсона, относительная погрешность численного решения определяется формулой

7 _ II Xh — X2h II 1

\\xh\\

здесь Xh и X2h — значения решения, найденные двумя шагами длины h и одним шагом длины 2h. Параметр к — порядок метода {к = 5). Для нахождения оптимального значения шага интегрирования воспользуемся оценкой

-0.9'ф'.

Эта оценка гарантирует приближение относительной погрешности к заданному значению точности є. Но шаг интегрирования не может меняться слишком быстро, поэтому его изменение ограничивается минимальным и максимальным значениями изменений

Щ 7TiCLx(^7TliTl(^S, Smin) і $max) ч

здесь smin = 0.1 и smax = 5.0 — экспериментальные значения параметров. Тогда новое значение шага интегрирования вычисляется по формуле

hnew h ' S\.

Этот алгоритм позволяет управлять погрешностью численного решения в пределах области непрерывности (Д функции F{t,X), но когда кривая решения пересекает границу между клетками, тогда этот алгоритм не может обеспечить приемлемое решение.

Пусть найденное решение X(t) лежит в клетке (Д, а значение в некоторый последующий момент времени X*(t + К) = X{t) + h • F(t,X(t)) принадлежит другой клетке СД j ф г. Полученное решение X*(t + h) показывает, что решение приблизилось к границе клетки. Мы предлагаем алгоритм склеивания решения на границе между клетками.

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

Шаг 1. (Идти к границе.) Найдем шаг интегрирования hi такой, что точка Х\ = Х(П), t\ — t + hi находится в клетке (Д и расстояние до точки решения на границе Хь Є Dp меньше, чем є • ||Хі||:

Математические структуры и моделирование. 2005. Вып. 15.

49

где є — заданная точность. Тогда мы говорим, что точка Х\ является приближением точки Хь-

Шаг 2.

Если (F(ti,Xi),n(Xb)) • (п(Хь), F(t2, Х2)) < 0, то найдена стационарная точка и программа завершается, иначе нужно перейти к шагу 3. (Здесь п(Хь) — вектор нормали к поверхности разрыва функции F(t,X) в точке Хь, (-, •) — скалярное произведение. Предложенное условие означает, что векторы F(t\, Х\) и F(t2,X2) направлены в сторону пересечения с поверхностью разрыва, а это означает, что решение не покинет поверхности разрыва).

Шаг 3. (Идти в клетку Cj.) Пусть следующий шаг интегрирования h2 будет таким, что точка X(t\ F h2j2) = Хь, тогда решение в клетке Cj находим с помощью схемы

Кг =

К2 = Р{іг + Н2,Хг + Іь2Кг),

X(t\ + h2) ~ Х2 = Х\ + h2 •--------.

(Здесь Xi~\~h2Ki представляет первое приближение точки Х2. Среднее значение скорости кривой в точке Хь задается вектором V = (Кг + К2)/2).

Далее докажем теорему об аппроксимации решения схемой (2). Введем некоторые обозначения. Пусть функция X(t) представляет точное решение динамической системы X = F(t,X). Рассмотрим решение около границы между клетками Сі и Cj. Будем считать, что решение в граничной точке Хь является непрерывным, но негладким. Кроме того, будем считать, что точка Хь не является стационарной. Обозначим Х\ = X(t\) Є Сі, Х2 = X(t2) Є Cj, ht = t2 — t\, a < ti < t2 < b. Тогда граничная точка определяется формулой Хь = X(t\ + Oht), в Є (0,1). Пусть точка Х% — приближенное решение, полученное с помощью численной схемы (2). Обозначим отклонение от точного решения АХ2 = Х2 — Х2.

Теорема 1. Пусть векторное поле F(t,X) ограничено: ||F(£,X)|| < М для всех I Gfi, a <t <Ъ, дифференцируемо по X в некоторой окрестности точки Х2; и производная в точке Х2 ограничена ||^(t2,X2)|| = L < ос. Тогда отклонение приближенного решения Х£ от точного значения Х2 определяется формулой ДХ* = R-ht + 0(h\), где R < 2М(в - 1/2).

Доказательство. Пусть Х7 = lim X(t\ + rht) и XT = lim X(t\ + т/ц). По-

0 т^в-о 0 т^в+О

скольку решение X(t) является гладким в клетках Сі, Cj, то значения X^ и X^ можно представить в виде

X- = Хг+ вИЖиХг) + 0

Х+ = Х2-{1- 0)htF{t2, Х2) + 0

Так как решение X(t) является непрерывным в точке Хь, то эти значения совпадают X^ = X^ = Хь. Приравняем правые части выражений (3) и (4),

50

В.В. Коробицын, Ю.В. Фролова Алгоритм численного решения...

получим уравнение для Х2:

Х2 = Хг + + (1 - 6)F(t2,X2))ht + 0(h2).

Вычислим приближенное значение Х2 = Х\ + Найдем отклонение

приближенного значения Х2 от точного решения Х2:

ДІ2 = Х2-Х2 = (1- e)ht(-F(t1,X1) + F(t2,X2)) + 0(h2).

Пусть — это приближенное значение точки Х2, полученное по схеме (2):

v F

А2 — Ai + tit-------2--------•

Вычислим отклонение точки Х% от точного значения Х2:

ах; = х2-х; = (e-^)htF(ti,Xi)+ht^i-e)F(t2,x2)-^F(t2,x2)^+o(h2).

Используя формулу Тейлора, представим функцию F{t) X) в окрестности точки Х2 в виде

F{t2,X2) - F{t2,X2) = LAX2 + 0(||ДХ2||2).

Это сделать возможно, поскольку в этой окрестности функция дифференцируема и ее производная ограничена. Кроме того, очевидно, что 0(| | АХ2| |2) = 0(Л|). Получаем

дх; = h,{(e - і) + (0 - 1Д )(F(i„X,) - F(t2,X2)) + O(ftr),

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

Ах; = R-ht + 0(h2t),

где R = {в —1/2) • (F(ti, Xi) — F(£2, Х2)). Поскольку функция F(t, X) ограничена числом М, то \ \F(ti,Xi) — F(£2,X2)|| < 2М, поэтому і? < 2М(в —1/2). Теорема 1 доказана.

Замечание 1. Если в условиях теоремы 1 положить, что шаг ht выбран таким образом, что в — 1/2 < ht. Тогда АХ% = 0(/г|).

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

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

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

• сохранение и загрузка файла данных системы. Формат файла текстовый, и файл может быть отредактирован при необходимости текстовым редактором;

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

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

Математические структуры и моделирование. 2005. Вып. 15.

51

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

dR

dt

dC

dt

dP

dt

P

R r 1 —

C

IP K/ zrcXrcR

XRcC

uppXppP

1 + IrcXrcR 1 + uppXpphppR ucpXcpP

ucpXcphcpC J

1 + hpcXpcR 1 + UppXpphppR uppXppeRpR + ucpXcpCcpC 1 + UppXpphppR + ucpXcphcpC

UcpXcphcpC

rrip

mc

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

J /Ж  А

(© <£> А

Рис. 1. Топология цепи питания адаптивных хищников.

Результаты моделирования приведены на рисунках 1 и 3. На рисунках представлены траектории решений в трехмерном фазовом пространстве (R,C,P), спроектированные на плоскость (R,C). Начальные значения переменных выбраны произвольно.

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

£ > £

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

R*

________еСР_________

^Rp{&RphcP — ecphRP)

Топология цепи питания не является фиксированной, она меняется от исключительно конкурентной цепи (Urp = 1; Ucp = 0 в модели (5)) до цепи всеядного

52

В.В. Коробицын, Ю.В. Фролова Алгоритм численного решения...

Таблица 1. Определения параметров и их значения, используемые при моделировании.

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

г Максимальная скорость роста растительности 0.3 час-1

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

ARC Скорость поиска растительности травоядным 0.037 л/

ARP Скорость поиска растительности хищником 0.025 л/

Аср Скорость поиска травоядного хищником 0.025 л/

URP Вероятность нападения хищников на растительность 1

ucp Вероятность нападения хищников на травоядных 1

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

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

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

crc Эффективность усвоения биомассы растительности для травоядных 0.6

crp Эффективность усвоения биомассы растительности для хищников 0.36

ecp Эффективность усвоения биомассы травоядного для хищников 0.6

mc Скорость отмирания травоядных 0.03 h-1

ТПр Скорость отмирания хищников 0.0275 h-1

хищника (urp = ucp = 1). Изменение стратегии происходит при переходе через переключательную плотность растительности R* (рис. 1(A)).

В случае, когда травоядные являются более питательными для хищников, чем растительность, неравенство (6) инвертируется. Когда плотность травоядных меньше переключательной плотности

= ________e_RP_______

Acp{ecphRp — eRPhcp) ’

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

4. Выводы и заключение

Результаты моделирования системы трехуровневой цепи питания в двух режимах приведены на рис. 1, 3. В фазовом пространстве системы имеется две клетки. В первом режиме клетки разделены плоскостью R = R*, а во втором — плоскостью С = 67*.

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

Математические структуры и моделирование. 2005. Вып. 15.

53

решения на границе клеток дает более точное решение, чем классический метод Рунге-Кутта, что продемонстрировано на графике решений системы двумя методами (рис. 4).

Рис. 2. Фазовое пространство системы в первом режиме. Параметры системы: ерр = 0.2, еср = 0.1, R* = 10, остальные указаны в таблице 1.

Рис. 3. Фазовое пространство системы во втором режиме. Параметры системы: ерр =0.1, еср = 0.3, С* = 5, остальные указаны в таблице 1.

54

Рис. 4. Кривые решения системы двумя методами: 1) классический метод Рунге-Кутты; 2) предложенный алгоритм. Параметры системы: ерр = 0.1, еср = 0.3, С* = 5, остальные указаны в таблице 1. Начальные значения: Я(0) = 15,С(0) = 1,Р(0) = 7.

Литература

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

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

3. Krivan, V. Optimal Intraguild Foraging and Population Stability // Theoretical Population Biology. 2000. V.58. P.79-94.

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