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

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

УДК 519.62

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

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

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

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

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

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

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

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

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

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

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

X(t) = F(t,X) (1)

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Кг), (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(h2), (3)

Х+ = Х2-{1- 0)htF{t2, Х2) + 0(h2). (4)

Так как решение 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(tuX2) + F(t2,X2)

А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

(5)

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

J® /Ж ® А

(© <£> А

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

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

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

£ > £ (6)

Оптимальная стратегия питания хищников состоит в поглощении растительности (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.