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

Вычислительные технологии

Том 8, № 6, 2003

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

А. Н. РОГАЛЕВ Институт вычислительного моделирования СО РАН,

Красноярск, Россия e-mail: ran@krsk.info

An inclusion of solutions and global error estimation are proved for the guaranteed methods, based on approximations of the shift operator along a trajectory and conversion of the analytical formulas of solutions. Algorithms of evaluation of numerical values of guaranteed bounds are considered. The convergence of guaranteed boundaries to a set of exact solutions, considered in projections to coordinate axes, is proved.

1. Включение множеств решений систем ОДУ с интервальными начальными данными

В этой статье описывается класс методов построения гарантированных границ множеств решений с неточно заданными (интервальными) данными [1-4]. Используемые интервальные числа и векторы записываются как выделенные жирным шрифтом буквы латинского алфавита, например, X, A, с,... , здесь

X = (x^), x E IIIR, A = (aj), aj E IIR, с e IIIR.

Условимся также, что имена символьных переменных функций будут обозначаться буквами в каллиграфическом стиле A, Z и т.п. Необходимые в дальнейшем сведения из области интервальной математики содержатся в [1-4] (см. также приведенную в этих статьях библиографию) .

Рассматривается решение системы обыкновенных дифференциальных уравнений (ОДУ) с интервальными начальными данными

dt = f (t,y),v(t°) = v0 E y0, (1.1)

где Y0 E HR”-.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2003.

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

Пусть f (£, у) — функция, непрерывная по первому аргументу и удовлетворяющая условию Липшица по второму аргументу на любом отрезке вида [£0,£0 + Т^Т < Т = £м с константой Липшица, равной Ь. Обозначим множество точных решений задачи (1.1) как

У*(«) = {у(()|у(го) е Уо}. (1.2)

Чтобы оценить множество всех точных решений, соответствующее вариации начальных данных на интервале У0, разумно найти для него гарантированную оценку (включение), т. е. множество С, для которого выполняется

с Э У * = {у(£)|у(£0) е У0}. (1.3)

Эффективным описанием множества включения, связывающим утверждения о существовании множеств включения с возможностью их построения, является интервал (интервальный вектор) У, поскольку для его представления в пространстве Кп достаточно использовать 2п граней, что существенно меньше 2п вершин, требуемых в других алгорит-

мах.

При построении гарантированных границ множеств решений используются понятия символьных формул и символьных методов [3].

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

Определение 2. Пусть У(£, У0) = иу°^°У(£, У0) — символьная формула для множества точных решений, а Уь(£,у0) — значение гарантированной оценки, полученное на основе символьной формулы приближенного решения Уй(£, У0). Если выполнено

аь(уа1(У(1,У0)),Уь(£,У0)) < г(к,у0), (1.4)

где йь(•, •) — метрика Хаусдорфа в пространстве ЕЖп [3, 4], г(к, у0) — однородная функция порядка 5, такая, что г (Л,, у0) ^ 0 при к ^ 0 для любого начального значения у0 е У0, то назовем Ун £ -решением (Бе^решением) задачи (1.1).

Определение 3. Гарантированной (внешней интервальной) оценкой множества решений ОДУ с интервальными начальными данными при £ е [£0, £м] называется интервальная функция У(£,У0), такая, что для любого у0 е У0 и е [£0,£м] у(£) е У(£,У0).

Выполнение гарантированных методов, основанных на аппроксимации оператора сдвига вдоль траектории, разделено на два этапа.

На шаге предиктор гарантированного метода находится символьная формула приближенного решения, зависящая от символьной переменной у0, на ее основе — форма ^-решения.

Построение (запись) символьных формул приближенных решений как векторных функций

^П(У0) о5П-1 0) ◦ ...5^У0),

где вектор У0 = (^0, У0, • • • , У°) — вектор начальных значений, описано в статье автора [3].

Определение 4. Пусть существуют нормированное пространство У и элемент ун € У такие, что в У значение ун получено после подстановки в символьную формулу 5(£,у°) соответствующего числового значения, тогда ун = уа1 (5(£, у0)) будет называться числовым значением символьной формулы.

Назовем интервальным значением символьной формулы значение Ун, полученное после подстановки в формулу 5(£,у0) в интервальном пространстве БЖга интервальных величин и вычисления по некоторому алгоритму.

На втором этапе определяется гарантированная оценка глобальной ошибки приближенного решения. Оценка глобальной ошибки получается на основе символьной формулы 0(£, У0) решения уравнения для глобальной ошибки, выписанного в некоторой промежуточной точке интервала [у(^г), у(^г+1)]. Для нахождения этой оценки предварительно строится “брус” точных решений, а затем определяется интервальная оценка как символьная формула решения этого уравнения. Таким образом, мы избегаем накопления суммы интервалов, включающих глобальную ошибку, т. е. устраняем влияние так называемого '^арр1^-эффекта, проявляющегося практически для всех известных интервальных и двусторонних методов решения ОДУ с неточно заданными данными. Далее к множеству (значению многозначной или интервальной функции) Ян(£, У0) добавляется величина, задающая гарантированные границы глобальной ошибки. В итоге будет получена требуемая гарантированная оценка множества точных решений

У(*,У0)=уа1(£(*, У0) У Я(*, У0)).

у °ет°

2. Гарантированные оценки глобальной ошибки и схема выполнения гарантированных методов

Пусть в процессе численного интегрирования значение индекса I фиксировано и символьные переменные у 1+к определяются из соотношений

Уг = (й), I - к < г < к,

к т

у+г), (2.1)

г=0 г=0

к < I < N.

Здесь Уг — символьные переменные начальных значений; уг+к — значения символьного решения в точке £1+к, ^ — произвольные точки из дискретного множества точек . Назовем этот способ символьной формой к-шагового (при к =1 одношагового) разностного метода.

Пусть у(£) — точное решение задачи (1.1), = £п — точное решение, начинающееся в

точке Уп-1, У (£п) = Уп — символьная формула приближенного решения гарантированного метода (2.1), точка Уп определяется формулой (2.1).

Локальной ошибкой символьного метода (2.1) называется величина

£п = УВ1 (Уп) - £п,

разность

Яп = Ув1 (Уп) - у(Г)

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

Для того чтобы величина локальной ошибки совпадала со значением погрешности дискретизации или невязки метода, нормируем [5, 6] формулу локальной ошибки

ег = уа1 (У(Г)) — , 0 < г < к,

в/(*п+г, Уп+г)^ , (2.2)

0 < п < М — к,

^(1) = ^

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

Для гарантированных символьных методов на первом этапе идет накопление смешанных символьно-числовых формул без применения численного счета по какому-либо приближенному методу. Кроме этого, решение систем ОДУ с неточными (интервальными) параметрами закладывает факт возмущения начальных данных и правых частей в саму

постановку задачи. Значит, исследование устойчивости методов решения систем ОДУ с

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

Применяются два способа переноса вдоль траектории решения погрешности, совершаемой на каждом шаге символьного метода. На основе этого в гарантированных методах строится оценка, гарантированно включающая все глобальные ошибки, соответствующие вариации начальных данных системы ОДУ. Далее эта оценка добавляется ко множеству приближенных решений (Б-решений) и формирует гарантированную оценку.

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

Рассмотрим первый способ оценки глобальной ошибки с учетом построенной символьной формулы, которая позволяет вычислить числовое и интервальное значения приближенного решения в любой точке [£0; £].

Х>с ‘ |{,1

г=0

(2.3)

е + = ^/(1)^ (агуа1 (У(^п+г))

г=0

Теорема 1. Если символьная формула дискретного метода (2.1) после подстановки числового значения у0 € У0

— имеет степень к > 2,

— вектор-функция /(£,у) является достаточно гладкой функцией,

то глобальная ошибка Д(£) метода(2.1) является решением следующего дифференциального уравнения:

лэд =д/(г,у)| _ „(/), Б(/) (24)

^ ду 1у=у(*)+0(*)хй(*) ( ) ( ), ( . )

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

Б (() = ^УТ — /<(' У<())- (2*5)

Доказательство. Подставив в исходную систему ОДУ (1.1) символьную формулу приближенного решения У (£), получим

^ = / (*, У (*))+ Б ((). (2.6)

функция Б (£) может быть эффективно вычислена как невязка приближенного решения по формуле (2.5). Вычитая выражения (1.1) из (2.6), получим равенство

~ УМ) = /(*,у) - /У(«))+ Б(«), (2.7)

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

^= СМВД+ Б М, (2.8)

где

гт- д/^ I — (29)

^(^ ду |У=У(*)+^(*)хД(*)- (29)

Здесь вектор 0(£) составлен из п компонент, являющихся числами на интервале [0,1], в

силу этого вектор у(£) + 0(£) х Д(£) состоит из компонент, принимающих некоторые значения точного решения, принадлежащие множеству (“брусу”) точных решений на интервале

[£к,^к+1] при варьировании начальных данных на интервале У0. В

В начальный момент £° известно точное включение начальных данных, поэтому можно записать формулу решения

£

Я(«) = упк,()5 К )С (2.10)

£0

где П(£, Л,) — это фундаментальная матрица системы.

В известных методах оценки глобальной ошибки [8, 9] разложения, подобные (2.4), используются как асимптотическая формула с вычислением значений разложений в узлах сетки.

Возможно считать, что система возмущений (2.4) представляет собой линейную неоднородную систему ОДУ с коэффициентами, заданными в некоторых промежуточных точках интервалов дискретной сетки, которые мы считаем далее символьными величинами. В книге [10] описываются два способа линеаризации нелинейной системы ОДУ, названные локальной и глобальной линеаризацией системы. Данный подход более применим для глобальной линеаризации системы, но он не приводит к системе ОДУ с переменными коэффициентами. Остаточные члены разложения вычислены в некоторых промежуточных точках. После этого для системы ОДУ, полученной после разложения и являющейся линейной относительно вектора возмущений, с коэффициентами, вычисленными в указанных промежуточных точках, выписывается формула точного решения. В этой формуле коэффициенты, содержащие промежуточные точки, описываются как символьные величины (однако символьные величины не являются функциями, зависящими от £). В силу этого метод позволяет получить формулу решения, в которую включены неизвестные значения точек разложения, обозначенные как символьные величины. Далее на этапе предиктор в методе определяется включение всех точных решений (достаточно грубое). После этого, подставляя в формулу вместо символьных величин интервалы, включающие все значения точного решения, мы получим гарантированную оценку глобальной ошибки.

Для удобства записи введем обозначение

д^(*, У) |

дУ 1У=У(£)+в(£)хП(£) ’

тогда определение формулы, описывающей поведение возмущений, можно получить, записав решение задачи

^ = ДМУ (()+ Б ((), (2.11)

У (О = У0.

Для этой системы решение может быть представлено следующим образом:

£

У (£) = ехр((А(£ — 0)У0 + у'ехр((Д(£ — т )))Б (т)^т. (2.12)

£0

Идея использовать формулу точного решения при численном решении линейных систем ОДУ появилась давно в численном анализе, например [5, 6]. В данной статье этот подход видоизменяется и применяется к нахождению символьной формулы решения уравнения относительно глобальной ошибки. Преобразуем уравнение (2.12) в предположении, что свободный член Б (£) задан в виде

Б (т) = (а ехр(в1т ),а ехр(в2т),..., а ехр(вй т))т, (2.13)

вектор е = (0,... , 0,1, 0,... , 0)т имеет все компоненты, равные нулю, кроме г-й, которая равна 1 (г — орт). Подставив (2.13) в (2.12) и воспользовавшись соотношением

Е о ехр(вгТ) = ехр(вг£Т),

где Е — оператор сдвига по узлам Е= уга+1, имеем

ь I

У(£) = ехр(Л(^ — ^0)У0) + ^ / ехр(Л(^ — Т))а» ехр(вг)^Т ■ в»

Ь Т

= ехр(Л(£ — т))У0 + ехр(Л(£ — т)+ в»Ет)^т ■ в». (2.14)

Положим Т = £ — £0, тогда

ср(ЛТ)у0 + ^ а [ <

»=1 0

м т

= ехр(ЛТ)У0 + ^ а» ехр(в»(£0 + Т)) / ехр(Л — в»Е)(Т — т))^т ■ в». (2.15)

»=1 0

Вводя обозначения

Ф»(Т)= I ехрА(т-£) ^ = I ехр(Л — в»Е)(Т — £)^£

У(£0 + Т) = ехр(ЛТ)у0 + } у а» у ехр((Л(Т — т) + в»Е£ + в»Е£0)^т ■ в» =

0

т т

и

&(Т) = Ф»(Т )в»

получаем следующее выражение для у(£) :

м

У (£) = У (£0 + т) = ехр(ЛТ )У0 + а» ехр(вг(*0 + Т ))&(Т). (2.16)

а»

»=1

Матрицы ехр(ЛТ), Ф»(Т) и векторы ^»(Т) удовлетворяют рекуррентным соотношениям [5, 6]:

ехр(ЛТ) = ехр(ЛТ/2) ехр(ЛТ/2),

Ф»(Т) = Фг(Т/2)(Е + ехр(ЛТ/2)ехр(—АТ/2)), (2.17)

£г(Т) = (Е + ехр(^4Т/2))^г(Т/2),

где

и, кроме того,

ехр(ЛТ/2) = ехр(ЛТ/2) ехр(—в»Т/2)

“ (ЛТ)к к!

(Л — вг)Й-1ТЙ

ехр(ЛТ) = ^ —— , (2.18)

к! ’ й=0

Ф.(Т) = £

к!

Й=1

Решение линейной системы ОДУ (2.11) с символьными начальными данными может быть найдено, если известны матрица ехр(ЛТ) и векторы ^г(£). А эти матрица и векторы

могут быть вычислены на основе последовательного применения рекуррентных соотношений (2.17) и преобразования ехр(Ак), £(к) при достаточно малом к = Т/2М с помощью частичных сумм рядов (2.18):

и £г(к) = ФДк)^. Далее решение задачи (2.11) вычисляется по формуле (2.16,) которую можно представить в виде

Второй способ гарантированного оценивания глобальной ошибки основывается на разложениях приближенных решений вблизи полученных значений. Для гарантированных методов применяются формулы разложения значения приближенного решения, в которых остаточные члены не отбрасываются, они включаются в множество (“коридор”) всех точных решений. Глобальная ошибка численного метода у(£га) — — это разность между

точным и приближенным решениями в каждой точке £”".

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

которое было бы равномерным при 0 < nh < const. В него необходимо еще добавить отвечающие возмущению члены [11]

которые будут компенсировать нерегулярность вблизи £0. Если возмущения затухают экс-

асимптотическое разложение в £га.

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

Разность точного и приближенного решений — у(£га) можно представить в виде

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

Формула (2.20) является основой для построения гарантированной границы ошибки. Чтобы компенсировать влияние '^арр1^-эффекта, суммирование локальных ошибок и значений решений не проводится. Вначале выписываются формулы локальной ошибки,

k=0

M

Y(tn+1) = Y(tn + T) = ехр(ЛТ)Y(tn) + Y, аг exp(A(tn + T))Gi(T). (2.19)

= y(tn) + dp(tn)hp + dp+l(r)hp+1 + ...,

- y(tn) = (dp(tn) + <)hp + (dp+l(tn) + <+>р+1 + ...,

поненциально вблизи t0 при n ^ 0, то на некотором расстоянии от t0 они не влияют на

R(tn) = (y(tn) - yn-l(tn)) + (y„-l(tn) - y„-2(tn)) +

+ (yn-2(tn) - yn-3(tn)) + ■ ■ ■ + (y1(tn) - y0(tn)).

(2.20)

Здесь у-(£) — точное решение системы (1.1), начинающееся в точке уа! (У- *)> точка эта

затем независимо друг от друга подсчитываются их интервальные расширения и на завершающей стадии объединяются значения Б-решений с включением глобальной ошибки.

Выражение для величины локальной ошибки имеет вид

к«+1 / к к \

р+ = (з + 1)1а'(1) ^ * + а — (з + 1) ^ * в")у( + )(^га+к-1 + ’

ст'(1> = (Е !<=1-

В этой формуле рга+к — вектор глобальной ошибки в точке £га+к; у(,5+1) (£га+к-1 + £) — векторная функция производной правой части, вычисленная в разных точках, зависящих от номера компоненты; — коэффициенты разностной схемы. Производная точно-

го решения у(5+1)(^га + С) оценивается некоторой интервальной функцией, гарантированно содержащей производные всех точных решений, лежащих на интервале [^,^+1]. Эффективный способ, позволяющий сделать это, состоит в том, чтобы после определения символьной формулы производной найти ее интервальное расширение. Для этого необходимо, чтобы зависимость решения от вектора начальных данных была бы дифференцируемой функцией (з + 1)-го порядка.

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

Гарантированный метод оценки множества решений выполняется согласно следующим шагам.

Шаг 1. Инициализируется модуль оценки множества решений, свои значения получают переменные, идентифицирующие систему: размерность, правая часть, начальные данные.

Шаг 2. Модуль, в котором определяется символьная формула Б-решения, возвращает эту символьную формулу (вектор с символьными компонентами) Б(£к, У0). Каждая компонента символьного вектора Бг(£к, У0), * = 1,... , п, определяется заново в каждой точке £к. Эта символьная формула является основой, на которой строится интервальное расширение (£, У0) по всем у0 € У0.

Шаг 3. Символьная формула (вектор с символьными компонентами) У, аппроксимирующая оператор сдвига вдоль траектории, позволяет определить прообразы эктремальных значений, принадлежащие вектору начальных значений У0 для всех узлов сетки £к.

Шаг 4. Для нахождения прообразов применяется одна из форм решения экстремальной задачи. Обозначим найденные точки й = 1,..., 2п. Они служат основой для определения множеств (“брусов”) точных решений, включающих все точные решения из окрестности точек у°.

Шаг 5. Для каждого начального значения, являющегося прообразом граничных значений Б-решений (Бе^решений), всего таких точек 2п, производится выбор произвольных интервалов, содержащих их и обозначенных где индекс й принимает значения 1 или

2 в зависимости от того, верхняя или нижняя граница Б-решения была выбрана, а индекс * принимает значения 1, 2,... , п, т. е. размерность по пространству. Выбор интервала Увгив а г, в котором начинается “брус” точных решений, производится на основе применения принципа сжимающих отображений в форме алгоритмов типа гомотопии. Конструктивно этот алгоритм заключается в поиске “пестрых” симплексов с изменением известных алгоритмов Ивза и Сайгала [12] и настройке их на случай интервальных чисел и многозначных

функций. Подробно эта часть гарантированного метода оценки множеств решений будет описана в другой статье. Результаты работы (графики гарантированных оценок и численные значения) этого алгоритма гомотопии, включая построение начальных “брусов”, представлены в настоящей статье ив [4].

Шаг 6. Начиная с полученных точек Y°Brus d j, d = 1, 2 i = 1,... ,n, производится построение S-решений (расширенных за счет выбора начальных “брусов” и привязанных к экстремальным значениям S-решений). Как упоминалось ранее, эти S-решения впоследствии называются “брусами”.

Шаг 7. Инициализация включений завершается снова применением алгоритма типа гомотопии. Он строит 2п-интервальные функции, включающие множества точных решений. Эти функции включают граничные гиперплоскости для множества решений системы ОДУ и позволяют сформировать гарантированные границы множества.

Шаг 8. Этот шаг завершающего определения гарантированной оценки глобальной ошибки реализуется двумя способами.

1. Определяется символьная формула решения задачи (2.11) согласно формуле (2.16), которую можно представить в виде

M

Y (tn+i) = Y (tn + T) = exp A(t)T a exP A(tn + T )G (T), (2.21)

i=1

матрица exp(AT) и векторы Gi(t) могут быть вычислены на основе последовательного применения рекуррентных соотношений (2.17) и вычисления exp(Ah), G(h) при достаточно малом h = T/2n с помощью частичных сумм рядов (2.18).

2. Производится суммирование оценок локальной ошибки вдоль траектории, заданной символьным решением, причем эти ошибки привязаны к каждому экстремальному значению S-решения. Формула локальной ошибки имеет вид (2.21), где i = 1,... , п — индексы компонент точного решения.

Выражение Yi(s+1) (tn + С) — это символьная формула производной i-компоненты точного решения, tn + С — некоторая промежуточная точка интервала [tn, tn+1]. На основе этой формулы определяется интервальное расширение производной. При этом интервальное расширение вычисляется для той величины интервала, которая найдена для каждой из точек прообраза YNw, Y^,, на основе принципа неподвижной точки.

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

Шаг 9. В интервальной арифметике к границам S-решений добавляется оценка глобальной ошибки (таких операций будет 2п):

Y(t, Y0) = val (s(t, Y°) U R(t, Y°)

V Y°eY°

Величина Y(t, Y°) — это гарантированная оценка множества точных решений с учетом всех видов ошибок.

3. Вычисление числовых значений гарантированных границ множеств решений на основе символьных формул

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

Теорема 2. Пусть интервальная функция Е(і, У) получена как интервальное расширение правой части системы ОДУ (1.1) и удовлетворяет условиям

• Е ограничена, т. е. существует К > 0, такое, что || Е(і, У) ||< К для любого вектора (і, уаі (У)) Є Б, Б Є К”;

• если последовательности точек уа1 (У^) Є К”- и з^, принадлежащие интервальной функции Е(і, У), сходятся

уа1(У^) ^ г, (3.1)

^ Є Е(і, г), (3.2)

з^ (^, уа1 (У^)) ^ з, (3.3)

то для предела этой последовательности следует включение з Є Е(і, У) для любого

і и любого уаі (У) Є К”,

• символьная формула Б (і, У0) строится на основе символьного линейного многошагового метода, удовлетворяющего условиям согласованности

к к т

= о, ^ іа = ^ ві, (3.4)

і=0 і=0 і=0

• для стартовых значений символьного многошагового метода выполняется включение у7 Є У 7, і = 0,1,... , к — 1, здесь у7 — значения точного решения в первых к точках,

тогда для любого у0, принадлежащего интервалу (интервальному вектору) У0, символьное решение, построенное по методу (2.1), сходится к точному решению задачи (1.1).

Доказательство. Доказательство сходимости основано на символьном представлении формулы приближенного решения Б (і, У0) = У (і, У0) и формулы глобальной ошибки ^(і0, У0). Пусть [і0, і0 + Т] — компактный интервал, содержащий более одной точки. Для любого у0 Є У0 в методе (2.1) могут быть вычислены последовательности числовых значений

уаі (У(ік)) = уа1(Ук(і0 + кЬ)), к = 1, 2,... (3.5)

В силу свойств аппроксимации (3.4) символьного метода (2.1)

|уа1 (У(ік+1)) — уа1(У(ік))| < ЬР,

где Р = Р(Е) > 0. Следовательно, функции У, определяющиеся выбором начальных данных У0 и имеющие последовательности числовых значений уа1 (У(£к)), будут равномерно непрерывными и равномерно ограниченными на интервале [£0,£0 + Т]. Согласно теореме Арцела — Асколи [13] возможно выбрать равномерно сходящуюся подпоследовательность уа1 (У(^)) при ^ 0, ^ то, такую, что ^^=0 кЛ = £к, и имеющую предел

У. Функции, соответствующие значениям уа1 (У(£к)), удовлетворяют условию Липшица равномерно на интервале [£°,Т]. Предельная функция у(£) = уа1 (У(£)) также будет лип-шицевой.

Чтобы показать, что у(£) удовлетворяет задаче (1.1), выбираем £,£ + в Е [£0,Т] в > 0.

Пусть для любого 0 < Л = < в числа п = п^ и т = т^ таковы, что для = £° + пЛ и

^п+йш = £«. + ктЛ выполняются неравенства

|Г - *| < Л, (3.6)

|Г+йш - ^ - в| < кЛ. (3.7)

Поскольку к=0 ам = 0, в силу (3.5) мы имеем

к й

ё^ = Л^АЛГ+г,уа1 (Ук(Г*)). (3.8)

г=0 Г г=0

Действуя аналогично, получим

ткЛ

м=0 ^=0 г=к

к—1 т—1 I

= — V У V (£п+*+^, уа1 (Ук(£п+*+^+М)). (3.9)

тк '

^=0 ^=0 г=0

Так как уа1 (Ук) равномерно сходятся к у, переходя к пределу, получим справа в (3.9)

в (ё ■- Щ' “

Из ограниченности уа1 (^), свойства (3.1) и равномерной сходимости ук ^ у имеем

^(^™+*+^к+^ У«+г+^к+^) е р(^™+*+^к+^ Y(t”+*+vk+M))

Из свойства монотонности по включению

к—1 т— 1

А = — V V ^(Г+^к+^, уа1 (Ук(£п+^к+м))) Е тк

^=0 V=0

Е Г(-£п+^к+м, уа1 (Ук(£п+^к+м))) (3.11)

для любого /. Поскольку

кк

, ви >

^ ^ ] Ди, Ди > 0,

^=0 ^=0

Графики гарантированных границ (а) и оценок (б) множества решений системы, носящей название “брюсселятор” и описывающей модели химической реакции, в которой концентрация реагентов меняется во времени. Проекция гарантированных оценок на плоскость Ь, у 1 (а) и Ь,У2 (б). Интервальный вектор начальных данных равен [1, 2] х [1.5, 2.5].

получим включение

1 к

^2вгАг е ¥(г, У(г)). (3.12)

^=°

На основе замкнутости множества уа1 (Т(г, Ук(г))), к = 1, 2,... при в ^ 0, что эквивалентно к ^ ж и дифференцируемости у(Ь), имеем

уа1 (тР)е ¥(г,у(г)). (3.13)

Поскольку включение (3.13) справедливо для любого у° е У°, строя объединенное расширение функции у (г) = У(г, У°) и выражения, задающего глобальную ошибку Щ(Ь°, У°),

получим гарантированную оценку множества точных решений

У (г, У°) = уа1 (г, У°) и Щг, У°)

V у°ет°

■

Интересной тестовой задачей является математическая модель реакции, в которой концентрации компонент меняются во времени по периодическому закону са = А + а вт^Ь) (см. рисунок). Тогда концентрации промежуточных продуктов у\,у2 в одном реакторе описываются уравнениями

ТУ1 ! , 2 ,Лу2 2

~Ж = 1 + уу2 - 4у- И = 3у - у У2'

Эта система называется моделью “брюсселятор”, ее решения очень чувствительны к изменению начальных данных. Выберем в качестве начальных данных интервальные векторы У° = [0.5,1.5] х [1.5, 2.5]. Получаемые гарантированные оценки множеств решений дают хорошее качественное совпадение с картиной поведения решений системы “брюсселятор”, а также точное включение множества точных решений.

Список литературы

[1] ROGALEV A.N. The interval methods for ordinary differential equations and shift along trajectories // IMACS/GAMM Intern. Symp. on Sci. Computing, Computer Arithmetic and Validated Numerics, Sept. 10-12, 1997. Ecole Normale Superieure de Lyon (France). P. X9-X13.

[2] РОГАЛЕВ А.Н. Использование границ глобальной ошибки в гарантированных оценках решений обыкновенных дифференциальных уравнений // Вычисл. технологии. 2002. Т. 7. (совместный вып. с Вестником КазНУ). Ч. 4. С. 88-95.

[3] РОГАЛЕВ А.Н. Гарантированные методы решения систем обыкновенных дифференциальных уравнений на основе преобразования символьных формул // Вычисл. технологии. 2003. Т. 8, № 5. С. 102-116.

[4] РОГАЛЕВ А.Н. Границы множеств решений систем обыкновенных дифференциальных уравнений с интервальными начальными данными // Вычисл. технологии. 2003. (В печати).

[5] БАХВАЛОВ Н.С., Жидков Н.П., КОБЕЛЬКОВ Г.М. Численные методы. М.: Лаборатория базовых знаний, 2001. 632 с.

[6] ИВАНОВ В.В. Методы вычислений на ЭВМ. Киев: Наук. думка, 1986. 583 с.

[7] SCHRODER, J. Fehlerabschatzung mit Rechenanlagen bei gewonlichen Differential-gleichungen erster Ordnung // Numerische Mathematik. 1961. B. 3. S. 39-61.

[8] ESTEP D. A posterori error bounds and global error control for approximation of ordinary differential equations // SIAM J. on Numer. Analysis. 1995. Vol. 32, N 1. P. 1-49.

[9] SkEEL R. Thirteen ways to estimate global error // Numerische Mathematik. 1986. Vol. 48, N 1. P. 1-20.

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

[11] ШТЕТТЕР Х. Анализ методов дискретизации для обыкновенных дифференциальных уравнений. М.: Мир, 1978. 461 с.

[12] EAVES R.C., Saigal R. Homotopies for computation of fixed points on unbounded regions // Mathematical Programming. 1972. Vol. 3, N 2. P. 225-237.

[13] ТРЕНОГИН В.А. Функциональный анализ. М.: Наука, 1980. 495 с.

[14] БАХВАЛОВ Н.С. К оценке ошибки при численном интегрировании дифференциальных уравнений экстраполяционным методом Адамса // Докл. АН СССР. 1955. Т. 104, № 5. С. 683-686.

[15] ЛОЗИНСКИЙ С.М. Недостаточные и избыточные методы численного интегрирования обыкновенных дифференциальных уравнений // Вестник Ленингр. ун-та. Сер. Математика, механика, астрономия. 1967. № 7. С. 74-86.

[16] ФИЛИППОВ А.Ф. Получение на ЭВМ строгих оценок решений дифференциальных уравнений // Журн. вычисл. математики и мат. физики. 1991. Т. 31, № 7. С. 9941005.

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

[18] BERZ M., MAKINO K. Verified integration of ODE’s and flows using diifferential algebraic methods on high-order Taylor models // Reliable Computing. 1998. Vol. 4, N 4. P. 361-369.

[19] NEDIALKOV N., JACKSON K., Pryce J. An Effective High-order interval method for validating existence and uniqueness of the solution of an IVP for an ODE // Reliable Computing. 2001. Vol. 7, N 6. P. 449-465.

[20] NEUMAIER A. The wrapping effect, ellipsoidal arithmetic, stability and confidence regions // Computing. 1993. Suppl. 9. P. 12-28.

Поступила в редакцию 25 февраля 2003 г., в переработанном виде — 12 мая 2003 г.