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

УДК 517.91

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

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

The method for computing an intersection of solution for system of ordinary differential equations with continuous break surface is suggested. The feature of algorithm is an ability to compute the near points on different sides of the surface. These points are the approximation of accurate solution point. The convergence of the numerical method is proven.

1. Введение

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

Простым примером гибридной системы является линейная система управления, представленная в работе [16]. Решение такой системы может быть представлено замкнутыми циклами, скользящей траекторией и вибрирующим режимом. Все эти эффекты требуют внимательного изучения в различных прикладных областях. Для решения подобных систем используются разные модификации методов Рун го Купы. Так, в работе [18] векторное поле функции правой части переопределяется на поверхности разрыва. А решение находится численно с помощью встроенных средств математического пакета MATLAB: решателя ОДУ (ODE solver) и детектора событий (event detector). Статья [19] посвящена локализации событий для систем ОДУ. Авторы обеих работ предполагают, что функции f+ и f-, определяющие правую часть системы до и после разрыва, определены на всей области определения системы. Однако это не всегда реализуется. Часто правые части определены только в соответствующих областях

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

Омский государственный университет им. Ф. М. Достоевского. E-mail: [email protected]

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

2. Постановка задачи

Предположим, что в пространстве Е™ задана поверхность б = {х : х е К", <?(х) = 0}, где д : К™ —> К — гладкая функция, имеющая непрерывные частные производные. Поверхность С разделяет пространство К™ на две области Б+ = {х : х € Ега,д(х) > 0} и = {х : х € Ега,д(х) < 0}, Определим также замкнутые области Л+ = Б+ и С, Т>~ = Б~ и С.

Пусть в областях Л+ и Л- определены непрерывные функции ^(х) : Л+ —> М™, ^(х) : Л- —> М™ и функция ~x.it) : М. —> К™ является решением задачи Коши

(Ьс = Г Г^х), х е Л+,

(И \ ^(х), х е В~,

х(*0) = Хо,

где х,, € Л — заданная начальная точка, соответствующая

Задача состоит в том, чтобы численно определить координаты точки х* пересечения кривой х(£) с поверхностью С, если такое имеется.

Очевидным способом решения этой задачи представляется следующий: численно находим кривую решения задачи Коши до тех пор, пока она не пересечёт поверхность С. Тогда, имея две соседние точки кривой решения, находящиеся по разные стороны от поверхности С, уточняем положение точки пересечения с помощью интерполяции. Но проблема заключается в том, что правая часть задачи Коши состоит из двух функций, определённых в областях Л и Л-, которые пересекаются по С. Если решать задачу (1) традиционными методами без учёта разрыва, то на поверхности разрыва численное решение сильно отклоняется от точного.

Для решения этой проблемы предлагается метод, состоящий из трёх шагов (предполагается, что х0 € _0+):

1) с помощью метода Рун го Купы с управлением длиной шага находим точки кривой решения в области Л+ вплоть до точки, находящейся на расстоянии не более заданной точности То1 от поверхности С;

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

3) находим пересечение полинома Ньютона с поверхностью С.

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

шаг, чтобы приблизиться к С, или перейти к шагу 2, если расстояние до поверхности С мало (меньше То1).

Далее вместо кривой решения используется её приближение полиномом Ньютона, Полином заменяет кривую решения вблизи поверхности С в областях Т>+ и Т>~. Но строится полином по точкам решения в области Т>+, а это означает, что метод будет успешен только тогда, когда полином будет построен по точкам, находящимся в достаточной близости к поверхности С.

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

Доказательство сходимости представленного метода можно провести в несколько этапов:

1) доказать сходимость метода решения ОДУ, в частности, метода Рунге-Кутты четвёртого порядка. Решение ОДУ приближается конечными приращениями по схеме Рун го Купы. Сходимость обеспечивается уменьшением шага интегрирования;

2) доказать сходимость приближения полиномом Ньютона, построенного по заданным значениям точек кривой, к кривой функции в локальной области, Сходимость обеспечивается уменьшением шага между опорными точками;

3) доказать сходимость метода Ньютона для нахождения точки пересечения кривой, заданной гладкой функцией в явном виде, с поверхностью, заданной в виде уравнения д{х) = О,

Доказательство сходимости трёх этапов и согласованность их даёт доказательство сходимости всего метода,

3. Первый этап: вычисление опорных точек

Задача состоит в том, чтобы найти опорные точки для дальнейшей аппроксимации, Опорные точки являются приближением к точкам кривой решения ОДУ

<7т

— = /(ж), хЦо) = х0, х : Е ->• Ега, / : Ега ->• Ега (2)

Необходимо ВЫЧИСЛИТЬ ТОЧКИ Х\ ~ Х\(1), Х2 ~ Х2^), где £1 = £о + Т, ^2 = ^0 + 2т, Будем использовать для этого метод Рун го Купы 4-го порядка:

X! = ДХ4(/,ж0,^0,т), х2 = ЯК4(/,х1^1,т). (3)

Можно также вычислить решение с удвоенным шагом 2т.

х2 = КК4(/,х0,и, 2т).

Оценка погрешности метода Рун го Купы задаётся неравенством

\х2 - х(Ь)\ < С ■ тк+1,

где к - порядок метода, в нашем случае равен 4,

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

\х2 - х(Ь2)\ < С ■ тк+1, \х2 - х(Ь2)\ < С ■ (2т)к+1

и \х2 — х2 + х2 — х{Ь2)\ < \х2 — х2 \ + \х2 — ж(£2)|, то получаем

С ■ (2т)к+1 < \х2 - х2\ + С ■ тк+1.

Отсюда \х2 — х2\ > С ■ тк+1(2к+1 — 1), Следовательно,

С-г

к+1

<

\х2 - х2\

2к+1 - 1'

Тогда, подставив это выражение в (4), получим

Iх2 - х(г2) \ < с ■ тк+1 <

\х2 - х2\ 2к+1 - 1'

(4)

При точном значении х(^) = хо отклонения х2 от х^2) ограничено величиной 6 = Стк+1, которая стремится к нулю при т —> О,

4. Второй этап: построение полинома

Заданы три точки (^,^1), (Ь,х2), (£з,жз), причём Ь2 = + т, Ь3 = Ь2 + т, кро-

ме того, известны значения первой производной аппроксимируемой функции /ъ/2,/3 в этих точках (вектора ^(^1), ^(^2), ^(^з))- Необходимо построить интерполирующую функцию (многочлен Ньютона) по трём равноотстоящим узлам с кратностью 2, Таким образом, максимальная степень многочлена з < 6, Условия построения многочлена:

лгв(* 1) = хг ,N^1) = $1,м3{12) = х2,м'3(г2) = /2,ад3) = ж3,лг'(*3) = /3. Таблица разделённых разностей имеет вид

ь Х\ /1

и Х\ с1\

ь х2 /2

*2 х2 с12

и х3 /з

и х3

<11-ь

Т

Н-Л1

т

<^2-/2

г

/3~^2

/2-2^1+/!

<к — 2/г+^1 2 г2

/з-2^+/2

^2-4/2+5^1-2/1 4 т3

273-5^2+4/2-^1 4 т3

2/з-6^2+8/2-6^1+2/1 8 г4

Здесь используются следующие обозначения d\ = Х2 X1, d2 = ЖзгЖ2. Тогда полиномы 4-го и 5-го порядков будут иметь вил:

N4{t3 + 9)=х3 + 9R\ + 02Д2 + 02(0 + r)Rl + 02(0 + r)2i?2,

-/V5 (^з + 0) = Ni{t3 + 0) + 02(0 + т)2(0 + 2 t)Ri-После подстановки значений конечных разностей из таблицы получаем

N4(t3 + 9)=x3 + 9f3 + 92-^ + е2(9 + т) • /з~2Д2+/2 +

(5)

+02(0 + Г)2 • 2/3-5^+4/2-dx ^

Л^з + в) = N4(t3 + 0) + в2(в + т)2(0 + 2r)2/3~6d2+88^r6dl+2/l-

Ошибки аппроксимации задаются выражениями:

N4(t3 + 0) — x{t3 + 0) = 02(0 + т)2(0 + 2т) • ———, t\ < ^1 < t3 +

5!

^5(^3 + 0) — ж(^з + 0) = 02(0 + т)2(0 + 2т)2 • ———, t\ < ^2 < t3 +

6!

При 0 G [i3; t3 + т] получим

|^4(^з + 0) — ж(£3 + 9)\ < ^ , \N5(t3 + 9) — x(t3 + 9)\ < 2q > (6)

где = sup |ar(i)(i)|, ii < t < h + 3t.

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

\N4(t3 + 9) - x(t3 + 0)| ~ \N5(t) - N4(t)\.

Использование полинома N4 является более рациональным потому, что опорные узлы для аппроксимации получены на первом этапе численным методом 4-го порядка. Учитывая, что опорные точки были получены с помощью метода Рун го Купы 4-го порядка, оценим погрешность аппроксимации кривой решения полиномом Ньютона при заданном отклонении опорных точек.

Пусть значение в точке t\ аппроксимируемой функции задано точно, а значения в точках t2 и t3 известны с некоторой погрешностью, не превосходящей 6 = С ■ т5. Запишем случайное смещение для опорных точек х2 = х2 + Ах2, х3 = х3 + Ах3 с заданным ограничением \/S.x2\ < 5, |Джз| < 6. Тогда значения производных тоже будут получены С погрешностью /2 = f(x2), /3 = f(xз), А значит, аппроксимирующий полином будет отклонён от точной функции больше, чем полином, построенный по точным значениям опорных точек, В связи с этим существенным будет доказательство сходимости смещённого полинома к точному решению.

Лемма 1. Если опорные точки Xi (i = 2,3j возмущены не более чем на 8 (\х* — Xi| < 8) и производная аппроксимируемой функции ограничена величиной В в рассматриваем,ой области (\f'(x)\ < В), то многочлен N^(t), построенный по возмущённым узловым, точкам, будет равномерно сходиться к iV4(i) при 8 —> 0 с оценкой погреитости, \Щ{Ьз + 9) — -/V4(i3 + 9)\ < 8 + т8{21 + 12В) при 9 < т.

Доказательство. Для оценки погрешности полинома N£ необходимо вычислить погрешность каждого слагаемого из (5) и сложить их.

Оценим ошибку значения функции в узлах i = 2,3:

/* = f(x*) = f(xi) + f'(xi) ■ (x* - x^ + o(x* - Xi).

Отсюда получаем оценку — < B8. Будем считать /* = /j + A/j, |A/j| < В8.

Оценим отклонения d\ и d*2.

,* х\ - хх х2 + Ах2 - xi х2 - xi , Ах2 , ,

d1 = —------=---------------=-----------1----= а 1 + Adi.

т т т т

Величина d\ имеет погрешность Ad\ = ^г2-, которая оценивается по формуле

т

|Adi

Ах2

т

Для d2 получаем

х*3 - х*2 х3 + Ах3 — х2 — Ах2 хз - х2 Ах3 - Ах2 , . ,

do =-----------------------=----------------------------------------------------=-------------------------1---------------------------= d2 + A d2,

т

т

т

т

где

|А d

'2 —

Ахз - Ах2

т

т

Оценим конечные разности из (5).

е>2*

-/Т4

/з* ~ d*2 /3 + А/3 — d2 — Ad2 /3 — d2 А/3 — Ad2

Rl + A Rl

т т т

Отклонение АД2 оценивается неравенством

< А/з + Ad2

т т

(2 + В) 8

т

Аналогично оценим отклонения Aи AR\.

A А АЛ. А

I АДо1

L3| —

Ai?2 <

А/з Ad2 /2 В8 48 В8

г2 + 2 г2 + г2 т2 г2 г2

2А/3 1 5Ad2 1 4 А/2 1 Adi

4 г3 4 г3 4т3 4г3

(11 + 6Д)£ 4r3

Максимальная абсолютная погрешность полиномиальной интерполяции с заданными отклонениями опорных узлов оценивается неравенством

|-/V4(t3 + 9) — N4^3 + 0)| < 5 + т8{21 + 125),

где 5 — максимальное отклонение опорных точек; т5(21 + 125) — отклонение многочлена в точке, расположенной на расстоянии не более т от опорной.

Из этого неравенства следует сходимость метода приближения кривой многочленом Ньютона при 5 —> 0, т = const и \В\ < ос. Лемма 1 доказана, ■

Следствие 1. Оценка погрешности аппроксимации точного решения x(t) задачи (2) с помощью полипом,а (5), построенного по опорным, точкам, с помощью метода (3), удовлетворяет условию

\М*А(1з + в)-х(и + в)\<К-т\

где К = ff + <7(1 + т(21 + 12В)) при 9 < т.

Доказательство следует из оценки погрешности метода Рун го Купы (4) и леммы 1,

5. Третий этап: поиск пересечения полиномиальной кривой с гладкой поверхностью

Задача третьего этапа состоит в том, чтобы вычислить координаты х* € К™ ТОЧКИ пересечения кривой X = N4(^3 + 9) с поверхностью <?(х) = 0, где д : К™ —> К — гладкая функция, имеющая непрерывные частные производные. Введём функцию Н{9) := д(М4(£3 + 0)), Задача нахождения х* сводится к решению алгебраического уравнения

Цв) = 0. (7)

В общем случае решение может не существовать или не быть единственным. Поэтому предположим, что на интервале 0 < 9 т имеется единственное решение 9* уравнения (7), То есть Н(0) ■ Ь,{т) < 0 и имеется единственное число 9* такое, что 1г(9*) = 0,

Применим итерационный метод Ньютона для решения уравнения (7), получим следующую процедуру для поиска параметра 9*, соответствующего точке х*:

9 -9 (Я)

9i+1 - 0г - (8)

йв \“г>

В качестве начального приближения возьмём 9о = г/2,

Сходимость итерационной процедуры может быть обеспечена, если будут выполнены условия следующей теоремы.

Теорема 1 (о сходимости метода Ньютона), Если выполнены условия:

1)

< а\ < оо при 0 < 9 < т,

S)

ад)-ад)-(§)(02)-(0.-02)

< а216*2 — 0\\2 при 0 < 01,02 < т,

3) 0 < 6*о <Ъ при Ъ = min{r/2, }

то итерационный процесс Ньютона (8) сходится с оценкой погрешности

\ег-е*\<—{а1а2\е0-е*\)2\

а\а2

Условие теоремы переформулировано в принятых обозначениях. Доказательство теоремы приведено в [1], Мы докажем, что условия этой теоремы выполняются для нашей задачи. Тем самым мы докажем сходимость итерационной процедуры (8),

Лемма 2. Итерационная, процедура (8) сходится к решению уравнения (1), если выполнены условия:

1) Ъг > 0, где Ъг = inf Ш

0 <в<т сШ

dh I .

(1в I ’

2) Ь2 < оо, где Ь2 = sup ||^|;

0 <в<т

3) г <%.

Доказательство. Доказательство леммы сводится к проверке выполнения условий теоремы 1, откуда будет вытекать сходимость итерационной процедуры (8),

Первое условие теоремы будет выполнено, поскольку Ъ\ > 0, Достаточно взять а\ = 1/&1-

Так как Ь2 < оо, то достаточно взять а2 = 2Ь2 и условие 2 теоремы будет выполнено.

Поскольку т < то верно условие Следовательно, для во = г/2

будет выполнено условие 3 теоремы 1, В итоге все условия теоремы 1 выполнены, Лемма 2 доказана, ■

Из теоремы 1 также следует, что достаточным условием достижения требуемой ТОЧНОСТИ £ можно взять \вп — 0ra-i| < е, при этом будет выполнено условие

\0п-0*\ < £.

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

Теорема 2, Численное решение системы, (2), найденное по формулам (1), (5), (3), сходится к точному решению при т —> 0, если функция / имеет ограниченные производные в области определения, и кривая решения x(t) пересекает поверхность д(х) = 0.

6. Выводы

Теорема 2 даёт теоретическое обоснование сходимости предложенного метода

нахождения точки пересечения кривой решения с поверхностью разрыва.

Литература

1. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. М.: Бином. Лаборатория знаний, 2008.

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

3. Коробицын В. В., Маренич В. Б., Фролова Ю. В. Исследование поведения явных методов Рунге-Кутты при решении систем обыкновенных дифференциальных уравнений с разрывной правой частью // Математические структуры и моделирование. 2007. Вып. 17. С. 19-25.

4. Коробицын В. В., Фролова Ю. В. Алгоритм численного решения дифференциальных уравнений с разрывной правой частью // Математические структуры и моделирование. 2005. Вып. 15. С. 46-54.

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

6. Новиков Е. А., Каменщиков Л. П. Реализация полуявного (4,2)-метода решения жёстких систем на параллельной ЭВМ // Вычислительные технологии. 2004. Т. 9. Вестник КазНУ им. аль-Фараби. Сер.: Математика, механика, информатика. №3(42). (Совм. выпуск. Ч. 1). С. 235-241.

7. Фельдман Л. П. Сходимость и оценка погрешности параллельных одношаговых блочных методов моделирования динамических систем с сосредоточенными параметрами // Научные труды ДонНТУ. Сер. Информатика, моделирование и вычислительные методы (ИКВТ-2000). Вып. 15. Донецк: ДонНТУ, 2000. С. 34-39.

8. Фельдман Л., Назарова И. А. Параллельные алгоритмы численного решения задачи Коши для систем обыкновенных дифференциальных уравнений // Матем. моделирование. 2006. Т. 18. №9. С. 17-31.

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

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

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

12. Cash J. R., Karp А. Н. A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides // ACM Transactions on Mathematical Software. 1990. Vol. 16. No. 3. P. 201-222.

13. Gear C. W., Osterbv O. Solving ordinary differential equations with discontinuities // ACM Transactions on Mathematical Software. 1984. Vol. 10. P. 23-44.

14. Guglielmi N., Hairer E. Computing breaking points in implicit delay differential equations // Advances in Computational Mathematics. 2008. Vol. 29. No. 3. P. 229-247.

15. Jackson K., Norsett S. The potential for parallelism in Runge-Kutta methods. Part I: RK formulas in standard form // SIAM J. Numer. Anal. 1995. Vol. 32. P. 49-82.

16. Johansson K. H., Barabanov A. E., Astrom K. J. Limit cycles with chattering in relay feedback systems // IEEE Transations on Automatic Control. 2002. Vol. 247. P. 1414-1423.

17. Park T., Barton P. I. State event location in differential-algebraic models // ACM Transactions on Modeling and Computer Simulation. 1996. Vol. 6. No. 2. P. 137-165.

18. Piiroinen P. T., Kuznetsov Y. A. An event-driven method to simulate Filippov systems with accurate computing of sliding motions // ACM Transactions on Mathematical Software. 2008. Vol. 34. No. 13. P. 1-24.

19. Shampine L. F., Thompsons. Event location for ordinary differential equations // Computer and Mathematics with Application. 2000. Vol. 39. P. 43-54.