CC BY
38
Вычислительные технологии
Том 16, № 4, 2011
Об определении точки пересечения кривой решения задачи Коши с поверхностью разрыва
В.В. Коробицын Омский государственный университет им. Ф.М. Достоевского, Россия
e-mail: [email protected]
Вычисление точки пересечения кривой решения с поверхностью разрыва является важным этапом решения систем обыкновенных дифференциальных уравнений с разрывной правой частью. Рассмотрен вопрос о существовании решения данного этапа. Предложен алгоритм численного определения точки пересечения. Доказана сходимость метода к точному решению. Приведены результаты вычислительных экспериментов для линейной и нелинейной поверхностей разрыва.
Ключевые слова: обыкновенные дифференциальные уравнения с разрывной правой частью, численный метод, экстраполяция, метод Ньютона.
Введение
Поведение систем управления, скорость изменения которых зависит от переключательной функции, часто моделируют с помощью гибридных систем, записывающихся в виде систем обыкновенных дифференциальных уравнений (ОДУ) с разрывной правой частью. Теоретическое исследование таких систем детально проведено в работе [1]. Поведение решения данных систем более сложное, чем непрерывных. Для численного исследования гибридных систем требуются новые методы решения, отличные от традиционных методов Рунге — Кутты, Необходимо учитывать качественное изменение поведения решения в области разрыва правой части системы, что отмечалось еще в [2].
Простым примером гибридной системы является линейная система управления [3]. Решение этой системы может быть представлено замкнутыми циклами, скользящей траекторией и вибрирующим режимом. Все эффекты требуют внимательного изучения в различных прикладных областях. Для решения подобных систем используются разные модификации методов Рунге —Кутты. Так, в [4] векторное поле функции правой части переопределяется на поверхности разрыва, а решение находится численно с помощью встроенных средств математического пакета MATLAB: решателя ОДУ (ODE solver) и детектора событий (event detector). Ряд публикаций (см., например, [5, 6]) посвящен локализации событий для систем ОДУ, т. е. определению места попадания траектории решения на поверхность разрыва правой части системы. Авторы этих работ предполагают, что функции f+ и fзадающие правую часть системы до и после разрыва, определены на всей фазовой области системы. Однако это не всегда реализуется. Часто правые части определены только в соответствующих областях системы по разные стороны от поверхности разрыва. В этом случае использование предлагаемых методов невозможно.
Представленный в настоящей работе метод позволяет находить точки пересечения кривой решения системы с поверхностью разрыва, используя определение функций
только в областях непрерывности, С помощью интерполяционного многочлена Ньютона строится продолжение решения. Кроме того, доказана теорема о существовании точки пересечения кривой решения ОДУ с поверхностью разрыва и о сходимости численного решения к точному. Метод используется как часть более сложного алгоритма [7, 8], предназначенного для численного решения систем ОДУ с разрывной правой частью,
1. О существовании решения
Предположим, что в пространстве М задана поверхность С = {х € Кга| д(х) = 0}, где д : М — М — непрерывная функция, имеющая непрерывные частные производные первого и второго порядка. Поверхность С разделяет проетранетво М на две области: Б+ = {х € Ега|д(х) > 0} и Б- = {х € Ега|д(х) < 0}, Определим замкнутые области V+ = Б+ и С, V- = Б- и С.
Пусть в областях V+ и V- определены непрерывные функции ¥ +(х) : V+ — Мга, ¥- (х) : V- — Мга, Функция х(£) : М — является решением задачи Коши
( ¥+ (х), х € Б+,
х' = ¥(х), где ¥(х) = 1 (¥+ (х) + Г(х))/2, х € С, (1)
I Г(х), х € Б-,
х(£о) = хо.
Необходимо найти х* и ¿* такие, что
х* = х(£*) и д(х*) = 0. (2)
Положим, что х0 € Б+ — заданная начальная точка, соответствующая ¿0. Отметим, что кривая х(£) может пересечь поверхность С несколько раз. Кроме того, при реали-
С
множеетво [¿', ¿''] такое, что д(х(£)) = 0 для всех £ € [¿', £"], Поэтому задача уточняется — необходимо найти первое пересечение кривой х(£) с поверхностью С, которое соответствует значению ¿* такому, что д(х(£*)) = 0 и для всех ¿0 < £ < ¿* выполняется условие д(х(£)) < 0 при д(х(£0)) < 0 (д(х(£)) > 0 при д(х(£0)) > 0),
Введем функцию Л,(£) := д(х(£)), описывающую приближение точки решения задачи (1) к поверхности С. Поскольку х(£) — решение задачи (1) и д(х) — непрерывная функция, то функция Л,(£) непрерывная и имеет непрерывную производную вплоть до пересечения с С, т.е. та отрезке [¿0,£*]. Существование решения задачи отыскания точки пересечения гарантируется, если выполнены условия следующего утверждения.
Утверждение 1. Кривая х(£) имеет пересечение (¿*, х*) с поверхностью д(х) = 0 но отрезке [£0, £ если Л.(£0) ■ /¿(¿^ < 0.
Доказательство. Поскольку функция х(£) является решением задачи Коши (1), то она является непрерывной. Функция д(х) непрерывна по условию. Следовательно, функция Л,(£) является непрерывной как суперпозиция непрерывных функций. Поскольку Л.(£0) ■ ^(¿1) < 0, то по теореме Больцано — Коши эта функция обращается в нуль на отрезке [¿о, ¿1] - П
Утверждение 2. Если, функция Л,(£) строго монотонна на, отрезке [¿0,^] м Л.(£0) ■ ^(¿1) < 0, то на, отрезке [£0,£1] найдется единственное пересечение (¿*, х*) кривой х(£) с поверхностью д(х) = 0.
Доказательство. Поскольку функция h(t) непрерывна (см, доказательство утверждения 1) и строго монотонна, то по теореме о существовании обратной функции она имеет единственное значение t* = h-1(0), которое и будет единственным пересечением кривой x(í) с поверхностью д(х) = О, □
Рассмотрим некоторую область U С Rn, содержащую точку xo и некоторую часть границы G. Ответ на вопрос о существовании пересечения кривой решения x(t) задачи (1) с поверхностью G в области U дает следующее утверждение.
Утверждение 3. Если скалярное произведение векторов n = g'(x) и f = f (x) для, всех x G U строго меньше нуля при x0 G D+ (строго больше нуля при x0 G D~), то x( t) G U
Доказательство. Строгое неравенство нулю скалярного произведения указывает,
f(x) U
U
U
Однако отсутствие циклов также гарантирует условие строгого сохранения знака
nf
x( t) G
вернуться в какую-нибудь предыдущую точку,
x(t) G
границы области U. □
Утверждение 4. Пусть на отрезке [t0,t0 + т] функция h(t) непрерывна и имеет непрерывную производную первого порядка. Предположим,, что h(t0) > 0 и h'(t) <
— R< 0. Тогда, при т = ^ функция h(t) принимает значение нуля, на отрезке
R
[t0,t0 + т].
h(t)
изводные первого порядка на отрезке [t0, t0 + т], то по теореме Лагранжа существует в G [t0,t0 + т] такое, что
h(t0 + г) - hjt0) = h, т
Следовательно, h(t0 + т) = h(t0) + Л/(в)т, А поскольку h'(t) < — R для всех t G [t0, t0 + т], то
hito + т)< h(t0) -Rt = h(t0) - R ■ Щ^- = 0.
R
Таким образом, h(t) знак на отрезке [t0,t0 + т]. Но эта функция непрерывна,
поэтому на данном отрезке она принимает значение нуля (см, утверждение 1), □
x(t) x(t)
шением задачи (1), то ее производная x' = f (x). Вторая производная определяется как производная сложной функции x'' = f'f (x), здесь f' — матрица частных производных fx Для доказательства существования решения задачи (1) (2) необходимо подмноже-U x0
f0 = f(x0) с заданным углом раствора а. Обозначим конуе в Rn как
C(xo, fo, а) = I x G
(х - Хр) • f0
0 < 71-и и „ и < cos a
||x - Xo||||to||
Здесь в числителе — скалярное произведение векторов в Rn. Подпространство S, ортогональное вектору f0, ^^^^^етим как S = {x G Rn|x • f0 = 0}.
Любой вектор в можно разложить на сумму двух, один из которых параллелен вектору 10, а второй ортогонален ему, т. е, лежит в S. Для разложения будем использовать операцию проецирования. Обозначим проекцию вектора х € на направление 10:
х ■ 10 10
РГГ0Х
ilfoll НУ
Проекцию вектора х на подпространетво S обозначим как
рг5х = х - ргГох.
Получим разложение х = ргГох + рг5х, причем ргГох ■ рг5х = 0,
Введем обозначения подмножеств С+ = С(х0,10, а) П С0 = С(х0,10, а) П С, С- = С(х0,10, а) П Б-.
Существование пересечения кривой х(£) с поверхностью С обеспечивается при выполнении условий следующей теоремы.
Теорема 1. Пусть х0 € Б+, векторное поле 1 (х) имеет непрерывные частные производные первого порядка в Б+. Если существуют положительные числа К, М и
л д(х0) Мт ^ +
К такие, что при т = —, а = a,тctg—^- для всех хбС+ вы/полнены, условия
1) 1Ко(1 (х))||> К
2) ||рг5(ГС(х))|| < М,
3) д'(х) ■ 1 (х) < — Д
то решение х(£) задачи, (1) пересекает поверхность С внутри кон уса С(х0,10, а) на отрезке [¿0, ¿0 + т].
Доказательство. Вычислим проекцию вектора 1 (х) на подпростран етво S для вектора х(£0 + т) € С+, удовлетворяющего (1), Внутри С + векторное поле 1 (х) имеет первые частные производные, поэтому согласно формуле конечных приращений получаем
1 (х(£0 + т)) = 10 + т? 1 (х(£0 + 01)), где 01 € [0, т].
Тогда проекция па подпространство S будет
рг51 (х(£0 + т)) = рг510 + т рг5 (ГС (х(£0 + 01))).
Так как для произвольного х € С + выполнено перавенство ||рг5 (1 1 (х)) || < М и вектор 1 ортогонален то ||рг51 (х(£0 + т))|| < Мт,
Получим разложение вектора (х(£0 + т) — х0) по 10 и Поскольку функция х(£) имеет непрерывные частные производные первого порядка, то
х(£0 + т) = х0 + т1 (х(£0 + 02)), где 02 € [0, т].
Поскольку ||ргГо(1 (х))|| > К, то
||рг,о(х(£0 + т) — х0)|| = т||рг,о 1 (х(£0 + 02»| > Кт (3)
и, как показано выше, ||рг51 (х(£0 + т))|| < Мт, в силу чего
||рг5(х(£0 + т) — х0)| = т||рг51 (х(£0 + 02)) | < Мт2. (4)
Угол в между векторами (х(£0 + т) — х0) и 10 вычисляется по формуле
/3 = агс!§ 11Рг^(х(^ + г)-хо)И
||ргГо(х(*о + г) -х0)|Г
Согласно полученным оценкам (3) и (4) получим
Мт2 Мт
р < аг^ —— = arctg —— = а. Кт К
Отсюда следует, что траектория х(£) на интервале [¿0,£0 + т] не выйдет за пределы конуса С,
Из условия 3 теоремы 1 и утверждения 3 следует, что траектория решения задачи (1)-(2) либо пересечет поверхность С внутри кон уса С, либо дойдет до границы конуса. Однако выше доказано, что на интервале [¿0,£0 + т] траектория решения х(£) не выйдет за пределы конуса. Из условия 3 теоремы 1 и утверждения 4 следует, что на интервале [¿0,£0 + т] функция Л,(£) = д(х(£)) принимает значение нуля. Следовательно, траектория х(£) пересечет поверх ноеть С внутри кон уса С+ на интерв але [¿0,£0 + т ], Теорема 1 доказана, □
Если задача (1) удовлетворяет теореме 1, то задача о поиске пересечения траектории решения с поверхностью С имеет решение. Приближенное значение точки пересечения будем искать с помощью численного метода, предлагаемого ниже,
2. Описание численного метода
Численный метод решения задачи (I) (2) представляется следующим: определяем точ-
х0
рееечет поверхность С, Затем, имея две соседние точки кривой решения, находящиеся по разные стороны от поверхности С, уточняем положение точки пересечения с помощью интерполяции [9]. Однако в этом случае правая часть задачи Коши включает две функции ¥ + и ¥ _, заданные в областях "Р+ и V-, которые пересекаютея по С, Поскольку функции определены не во всей области, то невозможно найти корректное решение задачи Коши традиционными методами решения ОДУ, Для решения этой проблемы предлагается следующий алгоритм.
Алгоритм А (приближенное вычисление решения задачи (1)-(2)). Предполагается, что х0 € Алгоритм позволяет найти две точки по разные стороны от поверхности С, приближающие точку пересечения х*, по нижеприведенной схеме:
1) вычисляется приближенное значение шага интегрирования до поверхности С по формуле
д'Ы-ЦкоУ и
2) выполняется к шагов методом Рунге — Кутты длины т/к го точки х0, В итоге получим точки хь х2,..., хк;
3) по точкам х0, хь..., хк строится ин терполяционный многочлен Ньют она ^(¿),
Г к + 1
приближающий кривую решения хт при £ € ¿о ¿о Н--;—'
к
4) итерационным методом Ньютона строится последовательность точек
= + ь = , V (6)
= 0, x0 = xfc, tk = to + т, i = 1, 2,
обеспечивающую приближение к точному решению x*; 5) выдается результат xl е, xl,
Алгоритм начинается с 1ценкГшага т до поверхности G но формуле (5). Коэффициент 0 < а < 1 вводится для обеспечения надежности оценки с целью предотвращения длины шага больше t* — t0, Дадее в теореме 2 приведено условие выбора параметра а. Используя вычисленное т, производится k шагов методом Рунге — Кутты порядка p с постоянной длиной шага интегрирования т/k. Это дает опорные точки для построения
k
многочлена. Однако эффективное значение степени полинома определяется порядком
p
степени s > p. Учитывая, что в опорных точках известны значения не только функции x(t), то и производных x'(t) = f (x), получаем k = |~(s —1)/2], Так, для схемы четвертого порядка необходимы три опорные точки x0, xe, x2,
Далее вместо кривой решения x(t) используется ее приближение мн ого членом Ns(t), который заменяет кривую решения вблизи поверхности G в областях D+ и D_, то строится только по точкам из области Т>+. Предполагается, что на отрезке to + т, to + —^—r
k
многочлен обращается в нуль. Только на этом отрезке экстраполяция многочленом Ньютона будет удовлетворять требуемой точности,
G
цедуры (6), Коэффициент b выбирается таким образом, чтобы на каждой итерации
G
вая приближение к точке пересечения с обеих сторон от этой поверхности. Как правило, достаточно принять b = 1.1. В знаменателе итерационной формулы присутствует производная от сложной функции, которая представляется в виде скалярного произведения двух векторов g'(x*_e)NS(tk + Выражения для производных могут быть вычисле-
ны аналитически или численно. Условие завершения итераций ||xl — xl_e| < To/, где To/ — заданная локальная точность,
2.1. Вычисление шага до границы
Определение шага до границы осуществляется на основе оценки скорости убывания функции h(t) := g(x(t)). Необходимо найти т такое, чтобы, сделав шаг т методом Рунге—Кутты, остаться внутри области D+, Однако при этом до границы остается расстояние, преодолеваемое не более чем за т/k. Таким образом, т должно удовлетворять неравенству
А < < (7)
k + 1 t* —10 w
t*
Утверждение 5. Предположим, что h(t0) > 0 и R = min h'(0) < 0. Тогда
öe[io;i*j
t* e (t0 + T,t0 +
ah(t0) k
для t = — —-—- при --< а < 1.
h'(tQ) F k +1
Доказательство. Так как функция к(£) непрерывна и имеет непрерывные производные первого порядка, то будет верна формула конечных приращений
+ г) - /¿(¿о) = ^ т
л г т т-г як(со)
для некоторого 9 € [¿о,£о + Л- Поскольку т = — —-—- и на а задано условие, то значение
л'(¿о)
функции к в точке ¿0 + т будет удовлетворять неравенству
Если из этой точки сделать шаг длиной т/к, то функция к поменяет знак, поскольку
, ( к + 1 \ ,, , -тД (¿0 - Г)Я к ¿о + —г"г - + г) > —> 1 ;
k J ' - k k +1
Отсюда следует, что на интервале (to + г; to + —r ) Ф^'нкция hit) поменяет знак, а
V к )
в силу непрерывности примет на данном интервале значение нуля, □
2.2. Вычисление опорных точек
Задача этапа состоит в том, чтобы начиная с известной точки xo найти для дальнейшей аппроксимации опорные точки, являющиеся точками решения задачи (1), Предположим, что точка x0 б D+, Тогда опорные точки будут определяться решением ОДУ с непрерывной правой частью f +:
x
' = f+(x), x(to) = xo, x : R ^ Rn, f+ : D+ ^ Rn. (8)
1
Необходимо вычислить точки X! « х(^), х2 ~ х(£2), • • •, хд. « х(^), где ti = to + —г,
к
1 = 1, 2,...,к, а т — величина шага интегрирования, определяемая по формуле (5), Предполагается, что величина т должна быть таков а, что хь х2,..., хк € и суще-
т
ствует величина 0 < 9 < — такая, что + 9)) = О,
к
Для численного решения системы (8) будем использовать метод Рунге — Кутты (ПК) порядка р:
х* = ЕКр , г = 1,2 ... ,к, (9)
здесь процедура НКр(Л, В, С, Д) определяет значение точки решения в момент времени Сметодом Рунге — Кутты порядка р начиная с точки В и используя правую часть Л, Оценка погрешности метода Рунге —Кутты задается неравенством [10]
Х ( ¿0 + у ) - Xj
/т\Р+1
£*с(-к) ■ <11Г)
где С те зависит от 1, т., к.
Метод Рунге —Кутты обладает сходимостью, поэтому при т — 0 отклонение хк от х(£0 + т) будет стремиться к пулю.
2.3. Построение интерполирующего многочлена
Задача состоит в том, чтобы построить интерполирующую функцию (многочлен Ньютона) по к + 1 равноотстоящим узлам с кратностью 2, поскольку известны значения не только функции, но и ее производной. Максимальная степень многочлена 5 не превышает 2к +1, Заметим, что векторный многочлен строится отдельно по компонентам, поэтому мы имеем п многочленов с разными независимыми коэффициентами. Многочлен должен удовлетворять следующим условиям:
N,(¿0 = хг, ) = {(хг), г = 0,1, 2,..., к.
В общем виде многочлен Ньютона можно записать в виде
N,(4 + в) = + + 02х(4, 4, 4_1) + в2 (в + х(4, 4, 4-1,4-1) + ...
(Т \ 2 / СТ \ ^
0 + ■ ■ ■ [в + — ^ х(4,4,4-Ь4-Ь...,4_с),
5 — 1
где с = —-— , д = з — 2с, х(4,4,4-1) = (^ — (хд; — ~х.к-{)/т)/т — разделенная разность
второго порядка, х(4, 4, 4^,4^) = — 2(хк — хк-1)/т + ^^/т2 — разделенная разность третьего порядка, х(4,4, 4_ъ4_ъ..., 4_с) _ разделенная разность порядка 5, Погрешность аппроксимации многочленом определяется неравенством
(с!)2т т ||х(4 + б1) - N^(4 + б1)|| < \ \-n.lt для 0 < в < (И)
(5 + 1)!к« к
здесь М3+1 = тах ||^(С)||; ¿о < ( < 4 + т- Из (11) следует сходимость интерполирую-
к
щего многочлена ^(¿) к кривой х(£).
2.4. Поиск пересечения полиномиальной кривой с гладкой поверхностью
Задача заключается в вычислении координат точки х* € Кга пересечения кривой х = N«(4 + б) с поверхноетыо #(х) = 0 и сводится к решению алгебраического уравнения
#(N«(4 + 0)) = О. (12)
Решение задачи (12) находим с помощью итерационной процедуры (6), сходимость которой будет обеспечена, если выполнены условия теоремы Канторовича о сходимости метода Ньютона [11].
Утверждение 6. Если задача (1)-(2) удовлетворяет условиям теоремы 1, то итерационная процедура (6) сходится к решению уравнения (12).
Доказательство. Поскольку функция ^(¿) — многочлен, построенный по опорным точкам, являющийся приближением к решению задачи (1), то ^(¿) будет гладким продолжением кривой решения х(£) из области в область Гладкость обеспечивается ограниченностью функции := #(^(£)) и ее первых и вторых производных
на отрезке
т
4,4 +
Кроме того, если выполнено условие 3 теоремы 1, то производная Л/(£) будет отделена от нуля. Этих условий достаточно для удовлетворения теореме Канторовича, из которой следует сходимость итерационной процедуры (6) к решению уравнения (12), □
2.5. Сходимость метода
Сходимость алгоритма А к точному решению задачи (1)-(2) обеспечивается, если выполнены условия следующей теоремы.
Теорема 2. Предположим, что на, некотором, отрезке ¿2] задача, (1)-(2) удовлетворяет условиям теорем,ы, 1. Тогда, алгоритм А сходится к точному решению при ¿0 — ¿*, если, выполнены следующие условия:
р
2) порядок многочлена ^(¿) удовлетворяет условию з > р;
3) количество вычисляемых опорных точек к = |~(з — 1)/2];
4) параметр а удовлетворяет условию ^ ^ < а < 1.
Доказательство. Примем последовательность начальных точек ¿0,1, ¿0,2,..., ¿0,г,..., которая стремится к точке ¿* слева. Пусть для каждой ¿0,г задано точное значение х0,г, Рассмотрим сходимость точек х^ к точному решению х* задачи (1)-(2).
При уменьшении расстояния тг = ¿* — ¿0,г получаем снижение погрешности вычисления опорных точек х^г, х2,г,..., хк,г со скоростью, имеющей для метода Рунге —Кутты порядка р величину 0(тгр),
Каждая кривая семейства построенная по соответствующим опорным точ-
к + 1
кам х | •. х •_>..:.....х/,. •. имеет отклонение на отрезке
¿0,г, ¿о,г Н--^—Ч
порядка O(rf).
При выполнении условий 2 и 3 теоремы 2 обеспечивается согласованность методов Рунге —Кутты и Ньютона по порядку погрешности. Это дает равномерную сходимость семейства кривых ^,г(£) к соответствующему участку крив ой решения х(£).
При выполнении условия 4 согласно утверждению 5 получаем гарантию существования решения задачи (12) на интервале + Ъ, ¿о,г + ^ ^ ^ т^ , а согласно утверждению
6 обеспечивается сходимость итерационной процедуры (6) к решению задачи (12).
Равномерная сходимость семейства кривых к х(£) и сходимость итерационной
процедуры (6) дают сходимость численного решения, полученного с помощью алгоритма А, к точному решению задачи (1)-(2). □ Теорема 2 дает теоретическое обоснование применимости алгоритма А для решения задачи (1)-(2). Далее рассмотрим эффективность применения алгоритма А на практических примерах.
3. Вычислительные эксперименты с системами на плоскости
В вычислительных экспериментах использовались схема Рунге — Кутты — Фельберга четвертого порядка и полином Ньютона пятого порядка. Условием завершения итерационной процедуры (6) было достижение предела точности типа double (To/ = 2 • 10-15),
3.1. Система с линейной границей
Рассмотрим исследуемую ранее в [7] систему линейных ОДУ, сшитую из двух систем на плоскости (y 1, y2) вдоль прямой y1 = а:
у1 = У2 - db y2 = У1 - Съ У1 < а
у1 = У2 - ^ У2 = У! - С2, У! >
(14)
На прямой у1 = а система терпит разрыв правой части уравнения. Аналитическое решение системы в области у1 < а задается функциями
У1 (£) = Ле4 + Ле- + С1, У2(*) = А1б* - ^е- + ¿1,
(15)
здесь А = ((у? - С1) + (у?2 - ¿1))/2 А2 = ((у? - С1) - (у20 - ^))/2, у1(0) = у?, у2(0) = У?-
Выберем в качестве точки в момент времени £ = 0 точку пересечения кривой решения с прямой у1 = а. В эксперименте использованы следующие значения параметров: с1 = 0.2, = 0.5, А1 = 0.25, А2 = 0.05 а = 0.5 и точки пересечения с координатами у° = 0.5 и у?2 е (0.5; 0.8).
Выделенный круг па фазовом портрете системы (рис. 1) показывает интересующую пас область вблизи точки пересечения кривой решения с поверхностью разрыва. При-
у1 < а
у1 = а
точка Пересе чепия.
Исследуем поведение численного алгоритма определения точки пересечения кривой решения с поверхностью разрыва. Найдем вид зависимости погрешности вычисления точки пересечения от длины шага т.
Зададим величину т и определим значения точного решения по формулам (15) при £ = -т. Полученные значения (у1(-т),у2(-т)) используем как начальные данные для численного алгоритма. В качестве входного параметра алгоритму необходимо задать значение параметра а, который используется в оценке расстояния до границы в формуле (5). Исследуем влияние а на погрешность вычисления точки пересечения. Относительная погрешность найденной точки оценивается по формуле
Р
УЫ - У г)2 + (УГ^Ж2
л/ (у?)2 + Ш2
Уг п
Рис. 1. Фазовый портрет системы (13) (14); в круге подробнее показана область вблизи точки пересечения поверхности разрыва
Рис. 2. Графики зависимости относительной погрешности Р от длины шага т в логарифмических шкалах при разных значениях параметра а: 0.9 (i), 0.67 (2), 0.58 (3), 0.57 (4) (а), при изменении координаты y2 точки пересечения и а = 0.9 (б)
Результаты вычислительного эксперимента показали зависимость относительной погрешности P от длины шага т при различных значениях параметра а (рис. 2). При фиксированной точке пересечения у0 = 0.5, у0 = 0.7 и изменяемой величине шага т для значений а в пределах от 0.67 до 0.9 алгоритм обладает сходимостью к точному решению при т < 1 (см. на рис. 2, а, кривые 1, 2 при — lgт > 0), причем скорость сходит
предел точности 10-16 при т < 0.02 (— lg т > 1.7) обусловлен разрядностью мантиссы
а
личивается и происходит постепенный переход (см. рис. 2, а, кривая 3) от зависимости шестого порядка (кривая 2) к зависимости первого порядка (кривая 4).
Полученный результат в полис объясним. В численном алгоритме применен метод Рунге —Кутты для вычисления опорных точек с шагом h/2 и h, Эти опорные точки используются дня построения многочлена Ныотопа, который приближает продолжение кривой решения. Согласно известному результату о точности приближения многочленом Ньютона, заявляемая точность приближения выполняется только в пределах од-h/2
2
сходимость метода при — < а < 1.
3
При изменении точки пересечения у0 = 0.5, у0 £ (0.5; 0.8), заданных значениях шага т и а = 0.9 также была получена сходимость метода к точному значению при
т
сти метода показана па графике и оценивает порядок метода величиной 5.8031. 3.2. Система с нелинейной границей
В качестве системы дня второго эксперимента была выбрана модель резонансного преобразователя. Модели подобных устройств активно изучаются дня построения эффективных преобразователей электрической энергии, применяемых при разработке различных электронных устройств. Преобразователи резонансного типа имеют преимущества
перед импульсными и управляемыми преобразователями, такие как низкие потери при переключении с высокой частотой и удобство использования дня фильтрации электромагнитных помех. Модель преобразователя с обратной связью и постоянным напряжением нагрузки Ц0 описывается следующей системой [12]:
, _ £2
, Х1 + ДЖ2 - Щ х2 =----' (16)
здесь х1 — напряжение на емкости я2 — сила тока на индуктивности г^, Щ — управляющий параметр, принимающий значения щ1 = Е — [70, и2 = Ц0 — Е, и3 = — Ц0 и = и0 в соответствующих областях П^, Области определяются системой неравенств
П1 Х2 > 0, 2 1 я2 < г2
П2 Х2 < 0, 2 1 ья2 г г,
Пз Х2 > 0, 2 1 ья2 » г г,
П4 Х2 < 0, 2 1 ья2 » г г,
где гг > 0 — заданная величина силы тока. Данная система имеет две границы, которые разбивают плоскость па четыре области. Фазовый портрет системы (рис, 3) имеет симметрию, поэтому поведения решений в верхней и нижней полуплоскостях совпадают.
Исследовалось поведение алгоритма при определении точки пересечения кривой решения с частью границы X + я2 — г^ = 0 при х1 > 0, я2 > 0, В данном случае кривая разрыва нелинейна, что позволило оцепить поведение алгоритма при нелинейной границе, В эксперименте использовались следующие значения параметров элементов схемы: Я = 0.2 Ом, Ь = 31 ■ 10-6 Гн, С = 2 ■ 10-6 Ф, Е = 500 В, Ц0 = 100 В и гг = 50 А.
Как и в предыдущем эксперименте, начальные точки траектории решения выбирались согласно заданному шагу до границы т. Случайным образом принималась координата в интервале (0; 50), вычислялось значение ж2> на границе. Затем фиксировалось значение шага т в интервале (0; 10-6) и вычислялись значения начальных точек траектории (я1(—т),я2(—т)), используя функции точного решения. На фазовой плоскости
Рис. 3. Фазовый портрет системы (16)
50 40 30 20 10 0
-50 -40 -30 -20 -10 0 10 20 30 40 Х1
Рис. 4. Начальные точки траектории (о) и точки пересечения (□) кривой решения с поверхностью разрыва (штрих)
-1 %Р~
987650 2 4 6 8 ТхЮ"7
Рис. 5. Экспериментальные значения вычисленной погрешности Р в зависимости от длины шага т
(рис. 4) эти точки показаны знаком о. Дмее запускмся алгоритм с шагом т при заданном параметре а = 0.9, который определял точки пересечения кривой решения с границей (знак □ на рис. 4). Вычисленные точки пересечения сравнивались с точкой (ж0,ж0) и вычислялась погрешность Р. Полученные значения Р (рис. 5) показали, что погрешность найденных точек не превышает величины 10-7,
Таким образом, представленный алгоритм обладает сходимостью, теоретическое обоснование которой дает теорема 2. Результаты вычислительных экспериментов подтверждают сходимость метода па линейной и нелинейной поверхностях разрыва правой части системы. Алгоритм может быть использован для вычисления точек пересечения траектории решения систем дифференциальных уравнений па границе разрыва, что является одним из этапов численного решения таких систем. Две точки, найденные по разные от границы стороны, позволяют провести анализ о продолжении решения.
Список литературы
[1] Филиппов А.Ф. Дифференциальные уравнения с разрывной правой частью. М.: Наука, 1985. 224 с.
[2] Gear C.W., Osterby О. Solving ordinary differential equations with discontinuities // ACM Trans, on Math. Software. 1984. Vol. 10. P. 23-44.
[3] Johansson K.H., Barabanov A.E., Astrom K.J. Limit cycles with chattering in relay feedback systems // IEEE IVans. on Automatic Control. 2002. Vol. 247. P. 1414-1423.
[4] Phroinen P.Т., Kuznetsov Y.A. An event-driven method to simulate Filippov systems with accurate computing of sliding motions // ACM Trans, on Math. Software. 2008. Vol. 34, No. 13. P. 1-24.
[5] Shampine L.F., Thompson S. Event location for ordinary differential equations // Computer and Math, with Appl. 2000. Vol. 39. P. 43-54.
[6] Park Т., Barton P.I. State event location in differential-algebraic models // ACM Trans, on Modeling and Computer Simulation. 1996. Vol. 6, No. 2. P. 137-165.
[7] Коровицын В.В., фролова Ю.В., Маренич В.Б. Алгоритм численного решения кусочно-сшитых систем // Вычисл. технологии. 2008. Т. 13, № 2. С. 70-81.
[8] Коровицын В.В., Фролова Ю.В. Алгоритм вычисления скользящего режима для системы с гладкой границей разрыва // Там же. 2010. Т. 15, № 2. С. 56-72.
[9] Хайрер Э., нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М.: Мир, 1990. 512 с.
[10] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Бином. Лаборатория знаний, 2008. 632 с.
[11] Канторович Л.В., Акилов Г.П. Функциональный анализ. М.: Наука, 1984. 752 с.
[12] Melin J., hljltgren A. A limit cycle of a resonant converter // Analysis and Design of Hybrid Systems 2003 / Eds. H.G.S. Engell, J. Zavtoon. Brittany, Prance: IFAC Publ., Els., 2003.
Поступила в редакцию 24 декабря 2010 г.
