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

Математические структуры и моделирование 2007, вып. 17, с. 19-25

УДК 517.91

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

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

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

1. Введение

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

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

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

Омский государственный университет, The University of Kalmar (Швеция). E-mail: korobits@univer.omsk.su

при исследовании систем обыкновенных дифференциальных уравнений с разрывной правой частью.

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

2. Определение понятия решения разрывной системы

Рассмотрим поведение автономной системы, описываемой набором n параметров y = (y1 ,y2,... ,yn) в некотрой односвязной области D С Rn. Эта область разделена на две части Di, D2 некоторой поверхностью G, заданной алгебраическим уравнением g(y) = 0. Таким образом, заданы области Di = {y : g(y) < 0} и D2 = {y : g(y) > 0}. В этих областях поведение системы описывается решением системы обыкновенных дифференциальных уравнений:

dy = ш.у є Di, (1)

dy = My), y є D2. (2)

Функции fi, f2 удовлетворяют условию Липшица.

Решение внутри каждой области Di и D2 находится обычным способом и является непрерывно дифференцируемой кривой. Если рассматривать решение задачи в обеих областях одновременно, то возникает вопрос — что считать решением системы?

Определение 1. Решением объединенной системы (1)-(2) будем считать кривую y(t), которая совпадает с решением (1) в области D1 и решением (2) в D2, а также склеивается на границе G = {y : g(y) = 0}. Склеивание понимается как существование точки y* Є G и величины t* Є R, таких, что lim^y(t) = y*

и lim y(t) = y*. t^t*+oyK ' y

Определение 2. Множество Di будем называть клеткой системы обыкновенных дифференциальных уравнений dy = F (y), если Di является линейно связным подмножеством множества Dc непрерывности функции F(y).

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

n

Dc = U Di-

i=1

Определение 3. Считаем, что система обыкновенных дифференциальных уравнений dy = F(y) имеет клетчатую структуру Ds, если объединение замыканий

20

клеток Di, i = 1, ...,n совпадает с замыканием области определения D функции F (У):

n

D :у D = D.

i=1

Определение 4. Решением системы обыкновенных дифференциальных уравнений dy = F (у) с клетчатой структурой Ds будем считать непрерывную кривую y(t), удовлетворяющую дифференциальному уравнению в каждой клетке Di и склеенную на границе между ними. Склеивание понимается как существование точки у* Є G на поверхности разрыва функции F и величины t* Є R, таких, что lim y(t) = у* и lim y(t) = у*.

’ t^t*-0yK J y t^t*+0yK J y

Приведем пример линейной системы ОДУ с клетчатой структурой. Пример интересен наличием циклического решения, которое отсутствует в каждой отдельной клетке.

3. Рождение цикла в линейной системе с клетчатой структурой

Рассмотрим систему линейных обыкновенных дифференциальных уравнений с клетчатой структурой, составленной из двух клеток Di : у\ < а и D2 : у і > а на плоскости (у1,у2):

= а11(уі — С1) + а12(у2 — di), it = а21(у1 — с1) + а22(у2 — ^, (у1,у2) Є D1,

It = а22(у1 — с2) + а12(у2 — d2),

FT = а21(у1 - °2) + а11(у2 - d2) , (у1,у2) Є D2•

= а21 = 1. решение каждой из систем

у2 - d1, (5

у1 - с1, у1 < а,

у2 - d2, (6

у1 - с2, у1 >а

В случае, когда а11 = а22 = 0, а12

dyi =

dt

dm =

dt

dyi _

dt

dm

dt

имеет одну особую точку (седло). Когда с1 < а и с2 > а, то особые точки лежат в соответствующих клетках D1 и D2. Кроме того, когда d1 = d2, — особые точки лежат на прямой ортогональной границе между клетками, а их сепаратрисы пересекаются в некоторой точке на этой границе. В этом случае образуется замкнутая область D0, заключенная между сепаратрисами. Решение внутри этой области образует цикл и не выходит за ее пределы. Цикл составлен из двух гладких кривых, склеенных на границе между клетками (рис.1). Кривые решения описываются функциями

у1 (t) = A^ r + А2Є (t r) + ау1, у2^) = A1et-t* + A2e-(t-t*') + ау2,

21

где Ai = ((y* - ayi) + (y* - ay2))/2, A2 = ((y* - ayi) - (y* - ay2))/2. Значения t*, yl, y* задают начальные значения, а величины ay1 = c1}ay2 = d1 при yl < а и ayi = С2, ay2 = d2 при y* > a.

Рис. 1. Фазовый портрет системы (5)-(6)

Длина полупериода решения внутри D0 задается формулой

T = / |ayi - а\ + у/(ayl - a)2 - 4AiA2

2 = П\ 2Ai

4. Вычислительные эксперименты для системы с клетчатой структурой

Исследуем поведение численного решения на задаче с клетчатой структурой

(5)-(6). Очевидно, что в областях Di и D2 известные методы Рунге-Кутта без труда справятся с этой задачей, но что получается при переходе решения из Di в D2 и обратно? Именно на границе происходит скачок погрешности решения из-за резкого изменения правой части уравнения. И дело даже не в том, что сильно меняется правая часть, а в том, что малейшие отклонения от траектории приводят к большим отклонениям при переходе через границу G. Рассмотрим это на нескольких тестовых экспериментах.

Исследуем погрешность численного решения для системы (5)-(6) с заданными параметрами di = 0.5, ci = 0.2, d2 = 0.5, c2 = 0.8, а = 0.5 и начальными данными yi(0) = 0.499999, y2(0) = 0.3. Приведем результаты численных экспериментов для пяти явных методов решения систем обыкновенных дифференциальных уравнений:

22

1) явный метод Эйлера (Explicit Euler). Одностадийная явная схема с оценкой погрешности методом Ричардсона;

2) метод Рунге-Кутты 4-го порядка (Runge-Kutta 4). Классическая четырехстадийная схема с оценкой погрешности по Ричардсону и автоматическим управлением длиной шага интегрирования;

3) метод Дорманда-Принса 5-го порядка (Dormand-Prince 5(4)). Восьмистадийная схема 5-го порядка с оценкой погрешности по вложенной схеме

4-го порядка;

4) метод Рунге-Кутты 2-го порядка (Runge-Kutta 2). Явная двухстадийная схема, предложенная Рунге, с оценкой погрешности по Ричардсону;

5) метод Рунге-Кутты-Фельберга 4-го порядка (Felberg 4(5)). Наиболее распространенная вложенная схема четвертого порядка с 7-ю стадиями, использующая встроенную схему 5-го порядка для оценки локальной погрешности.

Первый эксперимент предназначен для исследования поведения схем на разрывах. Оценивалась относительная погрешность численного решения на точном решении после преодоления одного цикла длительностью T. В этом случае точное решения попадает в начальную точку, образуя цикл. Значения численного решения показывают насколько оно отклонилось от точного за один цикл. Все рассматриваемые методы решения имеют алгоритм оценки локальной погрешности и выбора длины шага интегрирования на основе заданной точности. Предполагается, что при уменьшении заданной точности погрешность численного решения также должна уменьшаться, обеспечивая тем самым сходимость численного метода. Однако для разрывных систем это не выполняется. Кроме того, заметно значительное уменьшение порядка точности метода. Результаты эксперимента приведены на рис. 2. График зависимости изображен в логарифмическом масштабе: по горизонтальной оси откладывается порядок заданной точности (-log10(loc_tolerance)), по вертикальной — порядок относительной погрешности решения (-log10(rel_error)), вычисляемой по формуле

rel_error

\\Vn - У(T)||

!Ы!

где y*(T) — точное решение в точке t = T, yn — численное решение, полученное после n шагов, соответствующих той же точке, \ \ • \ \ — евклидова норма.

Как можно видеть из диаграммы, метод Рунге-Кутты-Фельберга имеет большую погрешность, чем метод Эйлера. Метод Рунге-Кутты 2-го порядка ведет себя нестабильно, что отражается в изменении тенденции «роста-падения» погрешности решения. При достаточно высокой заданной точности метод Дорманда-Принса ведет себя лучше метода Эйлера. Наилучшим образом показал себя классический метод Рунге-Кутты 4-го порядка, который обеспечил заданную точность до 5-го порядка включительно. И это значительно лучше других методов.

23

-Iog10(loc_tolerance)

Рис. 2. Зависимость относительной погрешности решения от заданной локальной точности

Безусловно, для нахождения численного решения требуется разный объем вычислительных затрат для различных методов. Оценим вычислительные затраты каждого метода для получения решения с относительной погрешностью, не превышающей заданную. Результаты приведены в таблице 1. Вычислительные затраты оцениваются количеством вычислений правых частей уравнения, которых было достаточно для обеспечения заданной точности.

Таблица 1. Оценка вычислительных затрат для нахождения численного решения

rel error Explicit Euler RK 4 Dormand-Prince RK 2 Felberg

5•1Q-2 711 363 1627 69 9729

со 1 О т-Ч ю 7188 411 1627 2001 89979

1 о т-Ч ю 71736 819 15595 — —

Особенность численного решения задачи с разрывами проявляется в изменении длины периода циклического решения. Так, для наиболее удачной схемы RK4 численное решение постепенно смещается вдоль числовой оси, уменьшая тем самым длину периода решения. Результат смещения продемонстрирован на рис. 3. Начиная с одной точки (слева) решение постепенно смещается (центр) и входит в противофазу (справа).

В такой ситуации оценивать погрешность численного решения относительно точного традиционными методами (как разность между практикой и теорией) будет неразумно. Поэтому имеет смысл разделить погрешность численного решения на две составляющие: смещение по оси времени и отклонение от траектории точного циклического решения. Результаты оценки смещений численного решения относительно точного после прохождения одного периода (t = T) приведены в таблице 2. Обозначения, принятые в таблице: Tol — заданная точность, At — смещение по оси времени, Аг — отклонение от траектории точного решения. Как можно видеть из таблицы, сходимость численного решения к точному наблюдаются, причем отклонения от траектории значительно меньше,

чемр§му^льей:и!еііІ¥к]5іреиі]мнїйтов показали, что для преодоления разрывов в правой части системы обыкновенных дифференциальных уравнений лучше подходят

24

25

У У У2

t

Рис. 3. Смещение численного решения RK4 относительно точного по временной оси

Таблица 2. Оценка смещения решения по оси времени и отклонения от траектории

Tol At Ar

10-2 -0,02149378 0,0009338849

10-3 -0,00018583 0,0002725838

10-4 -0,00034308 0,0000122797

10-5 -0,00001584 0,0000016340

10-6 0,00000491 0,0000004262

10-7 0,00000491 0,0000000504

10-8 0,00000502 0,0000000062

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

ЛИТЕРАТУРА

1. Деккер К., Вервер Я. Устойчивость методов Рунге-Кутты для жестких нелинейных дифференциальных уравнений. М.: Мир, 1988. 334 с.

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

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

4. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М.: Мир, 1999. 685 с.