CC BY
61
Вычислительные технологии
Том 8, № 5, 2003
ГАРАНТИРОВАННЫЕ МЕТОДЫ РЕШЕНИЯ
СИСТЕМ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ НА ОСНОВЕ ПРЕОБРАЗОВАНИЯ СИМВОЛЬНЫХ
ФОРМУЛ
А. Н. РОГАЛЕВ Институт вычислительного моделирования СО РАН,
Красноярск, Россия e-mail: [email protected]
We consider the systems of ordinary differential equations containing parameters which are not defined exactly, but specified in a certain interval. It is required to define a set of solutions which covers all solutions arising due to the variations of the parameters. The guaranteed method is based on the fact that the solution of a system can be expressed at any point as the analytical function depending on the initial data. Then the interval extension including all sets of solutions is constructed for this analytical function. These methods ensure the inclusions of sets of solutions of ODE systems in the obtained interval bounds.
1. Постановка задачи оценки множеств решений систем ОДУ
Рассмотрим решение системы обыкновенных дифференциальных уравнений (ОДУ) с интервальными начальными данными
где V0 £ БЖ" — интервальный вектор или прямоугольный параллелепипед в К"". Практически все известные численные методы решения дифференциальных уравнений мало пригодны для решения задач с такими данными. Кроме того, для всех интервальных и двусторонних методов решения ОДУ с неточно заданными данными [1-6] имеет место влияние так называемого 'отарр1^-эффекта [7, 8], проявляющегося в экспоненциальном росте границ интервалов. Применение гарантированных символьных методов позволяет получать верхние и нижние границы решений, в которые, безусловно, включаются все решения, соответствующие изменениям параметров, например коэффициентов и начальных
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2003.
(1.1)
данных в пределах заданных интервалов. Эти границы не подвержены влиянию эффекта, т.е. дают достаточно близкие включения. Описанные методы эффективно применяются для оценки решений систем ОДУ при постоянно действующих возмущениях (задачи практической устойчивости), оценки состояний управляемых систем при неточно заданных данных (экстремальные состояния, множества достижимости) и для некоторых классов задач детерминированной стохастики (поведение решений которых сильно изменяется при малых возмущениях входных данных, ошибках реализации модели). Получение результатов и изучение свойств гарантированных методов оценки множеств решений систем ОДУ могут быть сведены к решению задач с неточно заданными начальными данными, что не нарушит общности наших рассмотрений. В задаче (1.1) выполняются условия, обеспечивающие существование и единственность решений при любых начальных данных из заданного интервала.
Пусть f (¿, у) — функция, непрерывная по первому аргументу и удовлетворяющая условию Липшица по второму аргументу на любом отрезке вида [¿о, ¿о+7\], Т1 < Т с константой Липшица, равной Ь.
Обозначим множество точных решений задачи (1.1) как
^(*) = (у(£) I у(*0) е V0}. (1.2)
Для последующего описания алгоритма оценки множества решений поставленной задачи нам потребуются определения и свойства многозначной и интервальной функций.
Определение 1. Отображение Г(г) : Z ^ У, ставящее в соответствие каждому элементу г множества Z некоторое подмножество Г(г) множества У, называется многозначным отображением (функцией).
Многозначное отображение можно трактовать как однозначное отображение Z в 2¥, т.е. множество всех подмножеств У.
Определение 2. Для двух многозначных функций Г^г) : Z ^ У, Г2(г) : Z ^ У естественным образом вводится операция включения: Г1 С Г2, если Г^г) С Г2(г) для всех
г е ^
Определение 3. Сечением многозначного отображения Г называется такое однозначное отображение f : Z ^ У, что f (г) е Г(г) для всех г е Z.
Доказательство многих характеристик многозначных отображений основано на следующем свойстве многозначных функций.
Определение 4. Многозначное отображение является полунепрерывным сверху, если для каждого открытого множества и е У множество Г+(и) = {г е Z : Г(г) С и, Г(г) = 0} открыто в Z.
Для численного нахождения значений гарантированной оценки эффективным является использование подмножества многозначных функций — интервальнозначных функций.
Множество всех интервалов [7] с вещественнозначными границами ( вещественных отрезков числовой оси [а, а] ) обычно обозначают как ЕЖ, множество интервальных векторов — ЕЖ", множество интервальных матриц — ЕЕЯ"*"". Сами интервальные числа и векторы записываются как выделенные жирным шрифтом буквы латинского алфавита, например, X, А, с,... здесь
X = (х,), х е ЕЖ, А = (а^), а^- е Ш, с е ЕЕ!.
Условимся также, что имена символьных переменных и функций в этой статье будут обозначаться буквами в каллиграфическом стиле А, 2 и т.п. Арифметические операции, которые вводятся над интервалами, удовлетворяют свойству монотонности относительно включения [7] и в общем виде записываются так:
а * Ь = {х * у| х € а,у € Ь}, (1.3)
здесь символом * обозначена любая арифметическая операция. Над интервальными величинами определены также другие операции, например:
|а| = тах(| а |, | а |), т(а) = (а + а)/2, адг^(а) = (а — а).
Достаточно обстоятельное описание интервальных величин, интервальных методов и их свойств приводится в монографиях [7, 9, 10], а также в ряде статей. В этой работе кратко приведены необходимые сведения, а базовые свойства и определения можно найти в указанных выше работах.
Определение 5. Введем интервальное 'расширение Е функции f как отображение со значениями в 111", такое, что f (¿) € Е(1) для любого £ € Т. Часто полагают, что в самой формуле функции f заложен способ построения интервального расширения (естественное интервальное расширение, центрированная форма и т.п.).
В качестве расстояния в пространстве интервалов БЖга (п > 1) используется метрика Хаусдорфа.
На интервальных пространствах БЖ, БЖга метрика Хаусдорфа определяется по следующим формулам:
¿ь = ¿ьаиз(а, Ь) = тах{|а — Ь|, |а — Ь|}, (1.4)
¿ь = ¿ьаиз(А, В) = тахтах{|ак — Ь*|, |а* — Ьк|}- (1.5)
к
Пространство интервальнозначных функций
) = {д(г)1 7(г) < д(г) < 7(г)}
состоит из функций, каждая из которых является совокупностью (множеством) непрерывных вещественных функций #(£), заданных на интервале [¿о, ¿к] и удовлетворяющих двойному неравенству, содержащемуся в определении этого множества. Для этого пространства метрика
dh(f, g) = тахтах{ |7(г) — д(г)| , \/(г) — д(г)|} .
Пространство интервальных функций обозначим через С/([¿о,¿к])-
Гарантированные методы включения множеств решений должны обладать хорошими аппроксимационными свойствами в плане близости построенного интервального решения и множества решений исходной системы ОДУ, а также свойствами устойчивости относительно возмущения структуры ОДУ. Это связано с тем, что построение практически всех методов оценки множеств решений основано на выполнении неравенств относительно границ множеств (интервальных расширений) значений функций правых частей, которые будут справедливы не только для исходных решений, но и для многих близких функций. Разработка двусторонних и интервальных методов решения систем ОДУ [5, 7, 11, 12] первоначально была направлена на получение гарантированных оценок решения с учетом
влияния всех видов ошибок, в том числе ошибок выполнения арифметических операций. Практическое применение этих методов показало, что задача построения интервальной оценки единственного решения системы ОДУ преобразуется к задаче оценки совокупности таких решений (пучка траекторий).
Использование символьных формул, задающих аппроксимацию оператора сдвига вдоль траектории, является очень эффективным и полезным инструментом в области решения этих задач. Это подтверждают покоординатно сходящиеся оценки множеств решений [13], построенные для многих систем ОДУ с неточно заданными параметрами с использованием этих методов.
Определение 6. Оператором сдвига по траектории (оператором Коши) [14] системы обыкновенных дифференциальных уравнений
^ = I(г,у), у е К" (1.6)
называется оператор Б, зависящий от параметров в и т, такой, что Б (в, т) : К" ^ К"", и сопоставляющий значению всякого решения у (г) системы в точке г = в значение этого же решения в точке г = т :
Б(в,т) у(т) = у(в).
Для линейной системы ОДУ с постоянными коэффициентами оператор сдвига вдоль траектории задается формулой
Б (п,т) = ехр(п — т )А. (1.7)
Из этой формулы видно, что оператор Коши зависит только от разности параметров п — т:
Б (п + г,т + г) = Б (п,т).
Это равенство — следствие автономности системы ОДУ. Для нелинейной системы (1.1), удовлетворяющей условиям теоремы существования и единственности решения задачи Ко-ши, также справедливы подобные равенства с должными оговорками относительно области определения входящих в нее операторов.
Система ОДУ с неточно заданными данными описывается как ансамбль систем. Для этого правые системы части (1.1) должны быть непрерывно дифференцируемы как по г, так и по пространственной переменной у в некоторой области Е, а решения системы должны оставаться в области Е при их продолжении как вправо, так и влево по г. Тогда ансамбль системы ОДУ (1.1) — это множество идентичных систем вида (1.1), отличающихся друг от друга лишь начальными состояниями. Ансамблю системы (1.1) будет соответствовать в области Е при каждом г ансамбль изображающих точек. Каждая из изображающих точек у0 е Е4о, двигаясь по траекториям системы (1.1) согласно преобразованию у (г) = Бу0, переместится соответственно в точку у(г) е Е4. Таким образом, оператор сдвига Б(г, г0) переводит область Е4о в область Е4 С К" за время г — г0. Траектории ансамбля систем (1.1) распределены (вовсе не обязательно по закону распределения случайной величины) определенным образом в ансамбле, что зависит от свойств решений этой системы.
2. Описание метода сдвига вдоль траектории, основанного на преобразовании символьных формул
Чтобы оценить множество всех точных решений, соответствующее вариации начальных данных на интервале V0, разумно найти для него гарантированную оценку (включение), т.е. множество О, для которого выполняется
с Э У* = {у(1)|у(1°) € V0}. (2.1)
Эффективным описанием множества включения, позволяющим связать утверждение о существовании множеств включения с возможностью их построения, является интервал (интервальный вектор) V, поскольку для его представления в пространстве Кп достаточно использовать 2п граней, что существенно меньше 2п вершин, использующихся во многих других алгоритмах.
Пусть на интервале [¿°,£М], где рассматривается задача (1.1), введена сетка узлов
1° <11 <...<1М, К = шах{£г+1 - ?}.
г
Значения вектора численного решения и вектора правой части определяются в узлах разностной сетки уГ = у(1т) = (уГ ,УГ, • • • ,У2), Г = / (П = (/Г, /Г, ■■■, /Г).
В основе построения гарантированных символьных методов лежит символьная формула аппроксимации оператора сдвига вдоль траектории.
Определение 7. Символьная формула (аналитическое выражение) — запись имен переменных и совокупности действий, которые нужно проделать в определенном порядке над значениями этих переменных, чтобы получить значение функции. В силу этого символьный метод (аналитический метод) — запись численного метода как метода преобразования символьной информации (символьных формул) на языке математического анализа. В дальнейшем при записи символьных формул, аппроксимирующих оператор сдвига вдоль траектории, допускается включение в них числовых констант с отложенным выполнением арифметических действий над ними.
В этом символьный метод отличается от численного алгоритма, основанного на исполнении конечной последовательности действий над конечным множеством чисел. Чтобы строить символьные формулы, аппроксимирующие оператор сдвига вдоль траектории ОДУ и позволяющие получать достаточно точные включения множеств решений (например, интервальные расширения), необходимо получить формулу, хорошо приближающую точное решение, и выполнить алгоритм преобразования этой формулы. Следует также учесть, что требование устойчивости разностных методов как непрерывной зависимости их решений от возмущений входных данных трансформируется в случае реализации интервальных символьных методов.
Определение 8. Пусть Нг — последовательность нормированных пространств; Тг,1 = 1,... ,п — 1, — последовательность символьных формул непрерывных отображений Гг, таких, что Гг определены на прямом произведении Н1 х Н2 х • • • х Нг, отображают это произведение в пространство Нг+1 и задают зависимость между значениями решения в каждой точке области определения и начальными значениями. Тогда результат
последовательного исполнения преобразований формул
у1 = т ЧЛЛ У, У ') = 5^у°),
У2 = Уу1, У2) = 52(у°) О 51(У°),
Уг = Г(г°,...,гг,у°,у1,...,Уг) = 5г(У°) о5 1-1 (у°) о5 1(у°), (2.2)
ут = тт(г°,...,гт, у°, у1,..., ут) = 5 т(у°) о5 т-1(у°) о-о5 1(у°)
будет называться символьной формулой метода сдвига вдоль траектории решения системы ОДУ.
Очевидно, что для большинства методов Нг = Кг,1 = 1, 2, 3,... ,т — 1, и Нт = Лт, где Ят есть т-мерное евклидово пространство и Тг есть формула отображения, основанного на элементарных арифметических операциях.
Алгоритм исполнения этого метода особенно просто выглядит в случае применения явных разностных схем. Применяя последовательные подстановки и приведение подобных членов в формуле (2.2), можно перейти к выражению, зависящему только от у°. В общем случае в методе предлагается модель вычислений (преобразований и вычислений) символьных формул, основанная на поэтапном статичном хранении информации и преобразовании ее в завершающей стадии метода. Таким образом, формула будет представлять рекурсивную структуру, размер которой изменяется. Для записи такой формулы в компьютере используются линейные динамические структуры [15].
В силу этого модель вычислений (преобразований и вычислений) символьных формул осуществляется без явного выписывания суперпозиций компонент формулы, определяемых на каждом шаге. Связь между этими компонентами определяется посредством задания механизма адресации. Ссылки на адреса различных уровней хранятся в стековой памяти в виде дерева. Генерация кода вычислений по формуле (2.2) осуществляется в процессе обхода этого дерева, начиная с вершин.
В обозначенном выше алгоритме получения символьных формул
Уп = Тп(£°,£ь... ,£п, уу1,..., Уп) = 5п(у°) О 5п-1(У°) О ... О 51(у°)
используется следующая методика обработки последовательности символьных формул. Пусть ф(у°) есть однозначное отображение единичного интервала из Л1 на гиперкуб из Кп, которое каждой точке £ € Л1 сопоставляет некоторую точку у = ф(£). При помощи такого отображения можно построить алгоритм исполнения, который для каждой точки £ € Л1 позволяет определить формулу отображения Т (У°, У°,..., УП) и процесс ее сборки по адресам. Для этого предлагается использовать в качестве отображения ф(у°) — непрерывное, однозначное отображение единичного интервала на п-мерный куб, известное также под названием кривой Пеано, заполняющей пространство. Фактически кривая Пеано представляет собой непрерывную, нигде не дифференцируемую кривую, которая проходит через все точки единичного гиперкуба в пространстве Лп. Изобразить кривую Пеано нельзя, возможно лишь дать последовательность кривых [15], которая в пределе сходится к ней. Каждая такая кривая называется приближением кривой Пеано и имеет номер, определяющий ее номер в последовательности кривых.
Таким образом, т-е приближение можно рассматривать как некоторую аппроксимацию т-й функции в рекурсивной формуле (2.2). Это соответствие задано отображением
элементов конечного множества отрезков из единичного интервала и элементами конечного множества гиперкубов, входящих в Rn.
Такое соответствие строится следующим образом. Гиперкуб разбивается на 2n гиперкубов, называемых квантами первого разбиения. Длина ребра каждого такого кванта равна 1/2, кванты нумеруются индексом i\ от 2n — 1 так, чтобы номера квантов, имеющих общую грань, отличались на 1. Обозначать кванты первого разбиения будем как r(ii),ii = 0,..., 2n — 1.
Каждый квант первого разбиения разбивается таким же образом, как гиперкуб в Rn, на кванты второго разбиения с длиной ребра 1/4, которые нумеруются по тому же закону, что и кванты первого разбиения. При этом нулевой квант второго разбиения, входящий в r(i1), должен иметь общую грань с (2n — 1)-квантом второго разбиения, входящим в r(i1 — 1). Кванты второго разбиения обозначим через r(i1, i2), где i2 — номер кванта второго разбиения, входящего в квант r(ii).
Аналогично получаем кванты любого разбиения с номером m
r(il, i2, . . . ,im),
0 < ij < 2n. Способ соединения между собой квантов, номера которых различаются на
1 (в лексикографическом порядке), лежит в основе алгоритма, связывающего адреса, по которым хранятся компоненты символьной формулы (2.2).
Итогом работы описанного выше алгоритма будет возможность в любой точке tk построить символьную формулу решения (2.2) и вычислять на основе этой формулы значения решений.
Определение 9. Пусть существуют нормированное пространство Y и элемент yh G Y, такие, что значение yh, получено после подстановки в символьную формулу Y(t, Y0) = S(t, Y0) соответствующего числового значения в Y, тогда yh = val (S(t,y0)) будет называться числовым значением символьной формулы.
Определение 10. Символьная формула метода сдвига вдоль траектории сходится в точке tk, если после подстановки в нее числовых значений из области определения символьного метода справедливо соотношение
| y(tk,y0) — val (Yk(y0)) H 0 при h ^ 0, (2.3)
где y0,y0 G Y0.
Определение 11. Назовем Yh интервальным значением символьной формулы, если в формулу S(t, y0) подставлены интервальные величины и произведены вычисления по некоторому алгоритму в интервальном пространстве IIIRn.
Если для полученного интервального значения символьной формулы будет выполняться соотношение, описывающее меру близости этого значения и множества значений всех точных решений, то назовем его S-решением (Set-решением).
Определение 12. Пусть Y(t, Y0) = UyoeYoY(t, y0) — символьная формула для множества точных решений, а Yh(t,y0) — значение оценки, полученное на основе символьной формулы приближенного решения Yh(t, Y0). Если выполнено
dh(val(Y (t, Y0)), Yh (t, Y0)) < r(h,y0), (2.4)
где r(h,y0) — однородная функция порядка s, такая, что r(h,y0) ^ 0 при h ^ 0 для любого начального значения y0 G Y0, то назовем Yh — S-решением (Set-решением) задачи (1.1).
Это определение S-решения будет применяться в алгоритме для вычисления гарантированной внешней оценки Yh ,in D Y * и гарантированной внутренней оценки y h ,out С Y* множества решений задачи (1.1).
Пусть в процессе численного интегрирования значение индекса l фиксировано и символьные переменные yl+k определяются из соотношений
Yr = Sr(h), l — k < r < k,
k m
Y, аУ1+г = У+г), (2.5)
i=о г=о
k < l < M.
Здесь Уг — символьные переменные начальных значений; yl+k — символьные решения в точке tl+k, где tl — точки из дискретного множества точек (сетки). Назовем этот способ представления символьной формой k-шагового (при k = 1 одношагового) дискретного метода.
Такие символьные формулы позволят получить S-решения (1.2) задачи (1.1) с начальными данными в виде интервала. Для доказательства этого свойства приходится привлекать условия, обеспечивающие сходимость выражений, полученных при замене символьных переменных на числовые значения из интервалов начальных данных к решению.
Включение множеств решений достигается за счет свойства аппроксимации дифференциального оператора и надежного оценивания глобальной ошибки.
Исследуем разрешимость системы нелинейных уравнений, полученных из (2.5) при замене символьных переменных на числовые значения, выбираемые из интервала начальных данных.
Если используется полный набор начальных условий
yk+l = G (tk+l),l = 0, l,...,k — 1, (2.6)
где G(tl) — заданная непрерывная функция, такая, что val (G(t0)) = y0, то задача (2.5) является корректно поставленной в том и только том случае, когда разностный оператор Em=0 aival(Yk+i) устойчив [16, 17].
В случае гарантированных методов, построенных на основе символьных формул, в ходе преобразований вычисления не производятся. Поэтому условия, налагаемые требованиями устойчивости, ослабляются и преобразуются к ограничениям, налагаемым вычислениями в конечной точке интервала интегрирования.
Неопределенные коэффициенты аг, вг формулы (2.5) и начальные условия (2.6) следует выбирать так, чтобы они удовлетворяли следующим условиям:
— уравнение (2.5) является аппроксимацией уравнения (1.1), именно справедливо тождество
km
J^val (У (tn+i)) = h^piF (tn+i, val (Yn+i)) + 0(hr ф(К)), i=0 i=0
где У (t) — символьная формула решения задачи (2.5), (2.6), при этом ф(К) — функция типа модуля непрерывности, натуральное число r должно выбираться возможно наибольшим при заданном наборе параметров метода;
— задача (2.5), (2.6) однозначно разрешима, это означает существование таких чисел t > t0 и h0 > 0, что каждому значению h из интервала 0 < h < h0 сопоставляется определенное решение yk(У0) задачи (2.5),(2.6), k = 1, 2,... , M = (t — t0)/h;
— каждая последовательность решений val Yv(Y0), v = 1, 2,... , задачи (2.5), (2.6), соответствующая последовательности шагов hv ^ 0, v =1, 2,..., компактна;
— из сходимости последовательности значений формул val (Yk(Y0)), k = 1, 2,..., в точке t следует факт существования решения y(t), для которого || y(tk) — val (Yk(y0))
0 при hv ^ 0 равномерно относительно допустимых значений k и любого начального значения y0 G Y0.
Если вектор-функция f в параллелепипеде П имеет непрерывные частные производные до порядка s — 1 включительно, то необходимые и достаточные условия аппроксимации уравнения (1.1) уравнением (2.5) с порядком 0(hr0(h)) записываются в виде
k k m
= 0, ^a^ = v^A^-\ v = 1, 2,...,r. (2.7)
i=0 i=0 i=1
Пусть в уравнении (2.5) коэффициенты ai,ei выбираются так, чтобы оно аппроксимировало уравнение (1.1) с максимальным порядком при условии, что
k—1
Е A < A|, (2.8)
i=0, i=j
где
Ai = «0 + ■ ■ ■ + «i (2.9)
и j принимает значения 0,1,... , k — 1.
Учитывая первое из условий (2.7), перепишем уравнение (2.5 ), используя разности первого порядка A(Y1+i) = Y1+i+1 — Y1+i, в виде
k— 1 m
— Е AiA(Y1+i) = h Е AFl+\ k < l < M.
i=0 i=0
Разделив это уравнение на величину —Ai, получим
k—1 r
A(Yl) = — £ A*A(Yz—j+i) + j+i,
i=0,i=j i=0
где A* = A*/A¿, в* = —ei/A¿. Начальные условия будем выбирать в виде
Yl, l = k — 1,..., 1, y0 = y(t0), (2.10)
где Yl — такие n-мерные векторы, что график вектор-функции (2.10) принадлежит па-раллелипипеду П = {|t — t0| < D, || y — y0 ||< C } и
A(val (Yl)) = [val (Y1) — 0(h), val (Y1) + 0(h)], (2.11)
h ^ 0, O(h) ^ 0, l = k — 1, k — 2,..., 0.
При доказательстве сходимости символьных формул, построенных по символьному методу (2.5), используем оценки, подобные тем оценкам, которые применялись в доказательстве сходимости линейных многошаговых методов [11, 16].
Чтобы доказать разрешимость задачи (2.5), (2.10), определим метод последовательных приближений, все переменные в котором — векторы n-измерений. Нулевое приближение выберем в следующем виде:
(Y1 )(0), l = k - 1,k - 2,..., 0, (2.12)
так, чтобы график вектор-функции (2.12) принадлежал параллелепипеду П и
Д((У)(0)) = 0(h), l = 1, 2,... , k - 1, h ^ 0. (2.13)
Последующие приближения будем строить по формулам
val((Y')(s)) = y1, l = k - 1,..., 0,
k— 1 m
A(val((Yj)(s+1))) = - £ (Ai)*A(val((Y1—j+i)(s))) + h J>*(f1—j+T,
i=0,i=j i=0
I = k, k + 1,...,M - 1,
val((Y1)(s+1)) = val((Yn)(s+1)), l>M, s = 0,1,..., (2.14)
(f1 )s = f (t1, (y')(s)).
При соответствующем выборе чисел t > t0 и h0 последовательные приближения определяются для всякого неотрицательного целого s. Доказательство этого факта проводится путем суммирования формул последовательного приближения и применения условия Липшица.
Действительно, при s = 0 из (2.14) получается оценка
r
II A(val (Yk)(1)) ||< h(AC + J] |вЛМ) = P, (2.15)
i=0
где C — константа, определяемая равенствами (2.11), (2.13), M = maxten |f (t,y)| и A = EI=0 W |A*| < 1. Затем при s = 1 оценка
II A(val((Yk)(2))) ||< (A + 1)P,
и, наконец, методом индукции получим
II Aval (Y(s)) ||< (As—1 + +As—2 + ■ ■ ■ + 1)P. (2.16)
Для упрощения записи вводится обозначение
(Vk )(s) = A(val((Yk )(s))).
На s-м шаге описываемого процесса выводятся оценки
II A((V'))(s) ||< hR(A + t - t0)Rs, (2.17)
II ((V')(s)) ||< (t - t0)R(A + t - t0)Rs. (2.18)
Выбирая t > t0 так, чтобы выполнялись условия
n
(t - t0)(AC + J] |в*|M)(1 - A)—1 < B, (2.19)
i=0
t - t0 < H, jh0 < H, h0 < t - t0,
и условие
С = (Л + t - t0)R< 1, (2.20)
из (2.18) получим оценку
(V¿)(s) < (t -t0)£s.
Из этой оценки вытекает, что последовательные приближения (2.14) сходятся к некоторой вектор-функции yk при любом k = 1, 2,... , M. Переходя в (2.14) к пределу при s ^ то, убедимся, что эта вектор-функция является решением задачи (2.5) ,(2.11). Так как система неравенств (2.19) ,(2.20) относительно t > t0, h0 > 0 разрешима, вопрос о существовании метода вида (2.5),(2.6) с указанными условиями разрешимости сводится к построению формул вида (2.14), для которых выполняется условие (2.8).
Используем обозначения yk = val (Yk(Y0)) — семейство решений задачи (2.5), (2.11). Будем полагать, что символьная формула приближенного решения Y(t) выбирается как кусочно-линейное восполнение сеточных функций val (Yk)(t) при каждом t. Тогда производная Y (t) определится по формуле
Y'(t) = А(У: 1), tk—1 < t < tk, k = 1, 2,... , M - 1. (2.21)
h
Теорема 1. Если задача (2.5),(2.11) разрешима и каждая последовательность Yk ее решений, отвечающая последовательности шагов hv, сходящейся в нуле, компактна и выполняются условия
= 0, (2.22)
i=0
k— 1 m
Y, Ai = E A = 0, (2.23)
i=0 i=0
то эта задача сходится.
Доказательство. Пользуясь условием (2.22), преобразуем равенство (2.5) к виду
km
Ai(Yk)'(t + ih) = £ в/(t + ih, (Yk)(0)(t + ih)), (2.24)
i=0 i=0
где tk < t < tk+1, k = 0,1,...,M — 1, и проинтегрируем последнее равенство от t0 до t. Тогда получим
m—1 n „ t
Ai[Yk (t + ih) -Y0] = 5] вг / / (С + ih, Yk (С + ih))d£.
i=0 i=0 ^to
Перейдем в последнем тождестве к пределу при h = hv ^ 0, v ^ +то вдоль последовательности шагов hv, v = 1, 2,... , для которых последовательность соответствующих решений val (Yk(hv)) задачи (2.5), (2.10) сходится в пространстве к некоторой функции y(t). В результате получаем равенство
k—1 m r.±
0
-]>>i[Y(t)-Y°] = J]А/ /(С,У(С))dC,
i=0 i=0
k
из которого в силу условия (2.23) имеем
y(t) - y0 = /f (£, val (Y(£)M,
Jt°
что и требовалось доказать. □
3. Схема построения гарантированных границ множеств решений систем ОДУ
Определение 13. Гарантированной (внешней интервальной) оценкой множества решений ОДУ с интервальными начальными данными при t G [t0,tM] называется интервальная функция Y(t, Y0), такая, что для любого y0 G Y0 и Vt G [t0,tM] решение y(t,y0) G Y(t, Y0) или
(Y¿(t))¿=i,n = ([Y¿)iow, Yt,Up])í=li„ D ([Y*low, Y*up])i=i;ra = (Y*(t))¿=i,ra = Y*(t). (3.1)
Выполнение гарантированных методов, основанных на аппроксимации оператора сдвига вдоль траектории, разделено на два этапа — предиктор и корректор. На первом этапе (предиктор), происходит построение (запись) символьных формул приближенных решений как векторных функций
Sn(Y0) n-1 (Y0) ◦ ...S 1(Y0),
где Y0(t) = (Y0, Y°, • • • , Y^ — вектор начальных значений, рассматриваемых как символьные величины. Затем вычисляется область значений этой формулы, т.е. S-решение Sh(t, Y0) = val (Sh(t, h, Y0)). Подробно этот этап описан в разд. 1.
На втором этапе (корректор) определяется гарантированная оценка глобальной ошибки приближенного решения val (R(t, Y0). Оценка глобальной ошибки получается как включение решения уравнения для глобальной ошибки, выписанного в некоторой промежуточной точке интервала [y(t¿), y(ti+1)]. Для нахождения этой оценки строится предварительная интервальная оценка ("брус") точных решений, а затем определяется уточненная интервальная оценка как формула решения уравнения относительно глобальной ошибки. Подробное описание этих шагов гарантированного метода дано в работах [13, 18]. Таким образом, мы избегаем накопления суммы интервалов, включающих глобальную ошибку, т.е. устраняем влияние так называемого wrapping-эффекта. Далее к множеству (значению многозначной или интервальной функции) Sh(t, Y0) в точке tk добавляется величина, задающая гарантированные границы глобальной ошибки. В итоге будет получена требуемая гарантированная оценка множества точных решений
Y(t, Y0) = val ÍS (t, Y0) U R(t, Y0)) .
V Y°eY° J
Результаты применения методов, основанных на аппроксимации оператора сдвига вдоль траектории, дают хорошие оценки множеств решений, отклонение которых от точных решений стремится к нулю при уменьшении шага сетки и в определенных случаях позволяет получать новые результаты.
Пример 1. Гарантированные оценки решений нелинейной системы ОДУ четвертого порядка.
Рассматривается система ОДУ
(У1рО (И
2¿уд*) у!(*),
10 *У4(*)У1(*)5
(Уз(*0
= 2*У4(*),
(У4 (*)
= -2*(Уз(*) - 1),
(3.2)
для которой поставлена задача с неточно заданными начальными данными (принадлежащими некоторому интервалу).
Результаты вычислений и графики гарантированной оценки множества решений и экземпляров точных решений подтверждают, что гарантированные оценки для каждого координатного измерения включают экстремальные значения множеств точных решений. Оценка глобальной ошибки имеет порядок константы, умноженной на величину Нк, где к — порядок символьной формулы, аппроксимирующей оператор сдвига вдоль траектории. В действительности каждая траектория решения лежит в пространстве Л4, и в силу выполнения условий теоремы существования и единственности решений эти траектории нигде не пересекаются. Однако при проектировании этих графиков на координатные плоскости проекции траекторий могут налагаться друг на друга, что не означает их пересечения.
Рис. 1. Проекция гарантированной оценки множества решений на плоскость Ь,у 1.
Рис. 2. Графики экземпляров точных решений описываемой системы ОДУ, начинающихся в параллелепипеде V0 — проекция на плоскость у1.
Список литературы
[1] Добронец Б.С., ШАйдуров В.В. Двусторонние численные методы. Новосибирск: Наука, 1990. 208 с.
[2] Nedialkov N., Jackson K., Corliss G. Validated solutions of initial value problems for ordinary differential equations // Appl. Math. and Comput. 1999. Vol. 105, N 1. P. 2168.
[3] 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. 361369.
[4] Nedialkov N., Jackson K. An interval Hermite — Obreschkoff method for computing rogorous bounds on the solution of an initial value problem for an ordinary differential equation // Reliable Computing. 1999. Vol. 5, N 5. P. 289-310.
[5] 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.
[6] Nedialkov N., Jackson K. A New Perspective on the Wrapping Effect // Interval Methods for Initial Value Problems for Ordinary Differential Equations. Perspectives on Enclosure Methods / Ed. by Kulisch U., Lohner R., Facius A. B.: Springer-Verlag, 2001. P. 219-264. (http://www.cs.toromto.edu/NA/reports.html/ ned.scan00.ps.Z)
[7] Moore R.E. Interval Analisis. Englewood Cliffs: Prentice Hall, 1966.
[8] Lohner R. On the ubiquity of the wrapping effect in the computation of error bounds // SCAN-2000, Interval 2000. GAM / IMACS Intern. Symp. on Sci. Computing, Computer Arithmetic and Validated Numerics, Sept. 18-22, 2000. Univ. Karlsrue, Institut fur Angewandte Mathematic (Germany). P. 36.
[9] Алефельд Г., Херцбергер Ю. Введение в интервальные вычисления. М.: Мир, 1987. 356 с.
[10] Калмыков С.А., Шокин Ю.И., Юлдашев З.Х. Методы интервального анализа. Новосибирск: Наука. 1986. 221 с.
[11] Горбунов А.Д., ШАхов Ю.А. О приближенном решении задачи Коши с наперед заданным числом знаков // Журн. вычисл. математики и мат. физики. 1963. Т. 3, № 2. С. 239-257.
[12] Филиппов А.Ф. Получение на ЭВМ строгих оценок решений дифференциальных уравнений // Журн. вычисл. математики и мат. физики. 1991. Т. 31, № 7. C. 9941005.
[13] Рогалев А.Н. Границы множеств решений систем обыкновенных дифференциальных уравнений с интервальными начальными данными // Вычисл. технологии. 2003 (в печати).
[14] НЕМЫЦКИй В.В., СТЕПАНОВ В.В. Качественная теория дифференциальных уравнений. М.: Гос. изд-во научно-техн. лит-ры, 1949. 550 с.
[15] ВИРТ Н. Алгоритмы + структуры данных = программы. М.: Мир, 1985. 407 с.
[16] DAHLQUIST G. Convergence and stability on the numerical integration of ordinary differential equations // Math. Scand. 1956. Vol. 4, N 1. P. 33-53.
[17] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Лаборатория базовых знаний, 2001. 632 с.
[18] РогАЛЕв А.Н. Включение множеств решений дифференциальных уравнений и гарантированные оценки глобальной ошибки // Вычисл. технологии. 2003 (в печати).
[19] Новиков В.А., РогАЛЕв А.Н. Построение сходящихся верхних и нижних оценок решений систем обыкновенных дифференциальных уравнений // Журн. вычисл. математики и мат. физики. 1993. Т. 33, № 2. C. 219-231.
[20] РогАЛЕв А.Н. Нахождение оптимальных гарантированных оценок множеств решений систем ОДУ с интервальными данными // Вычисл. технологии. 1995. Т. 4, № 13. С. 58-64.
[21] РогАЛЕв А.Н., Шокин Ю.И. Исследование и оценка решений обыкновенных дифференциальных уравнений интервально-символьными методами // Вычисл. технологии. 1999. Т. 4, № 4. C. 51-76.
[22] РогАЛЕв А.Н. Динамические разностные схемы для построения интервальных методов решения ОДУ с неточно заданными начальными данными // Вычисл. технологии. 2001. Т. 6, ч. 2 (спецвыпуск). С. 510-518 (http://www.ict.nsc.ru/ws/NikNik/).
[23] РогАЛЕв А.Н. Использование границ глобальной ошибки в гарантированных оценках решений обыкновенных дифференциальных уравнений // Вычисл. технологии. 2002. Т. 7, № 4. С. 88-95.
[24] ANGUELOV R. Wrapping effect of the initial value problems for ODE: Applications // Reliable Computing. 1999. Vol. 5, N 2. P. 143-164.
[25] LOHNER R.J. Program AWA (Anfangs-Wert-Aufgabe). 1994. (ftp://iamk4515.mathematik.uni-karlsruhe.de/pub/awa)
[26] Shampine L.F., GORDON M.K. Computer Solution of Ordinary Differential Equations: The Initial Value Problems. San Frantisco, 1975.
Поступила в редакцию 25 февраля 2003 г., в переработанном виде —12 мая 2003 г.
