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

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

Том 9, № 1, 2004

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

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

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

Guaranteed methods are applied for the estimation of solution sets of ODE systems with interval initial data. ODE systems with permanently operating perturbations (problem of practical stability), controlled ODE systems with the unknown data (extreme states of these systems and reachable sets) and some classes of the “determined stochastic”

ODE problems are considered. Bounds of solution sets are calculated and their graphs are presented.

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

ddt = f (t,y), y(t0) = y0 e Y0 e HR", (1)

где Y0 — это интервальный вектор (прямоугольный параллелепипед в пространстве R").

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

Результаты применения методов, основанных на аппроксимации оператора сдвига вдоль траектории, дают хорошие оценки множеств решений, отклонение которых от точных решений стремится к нулю, при уменьшении шага сетки. В статье приводятся примеры оценки множеств решений систем ОДУ с интервальными параметрами в задачах, типичных для данной ситуации нахождения гарантированных оценок множеств решений. Было бы интересно сравнить численные значения гарантированных границ множеств решений с аналогичными расчетами, полученными другими методами. К сожалению, подобные результаты практически не включаются в публикуемые работы. Реально проанализировать численные значения возможно лишь, рассматривая статью [6] и пользуясь программной системой AWA [7]. В других работах последних лет [8-17] не приводятся примеры расчета ОДУ с начальными интервальными данными, отличными от нуля, кроме того, в них

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

печатается лишь ширина ошибки интервального решения, но отсутствует иная информация о самом решении. В принципиальном плане ни один из этих методов не использует в процессе решения иной методики, отличающейся от QR-вращения интервалов. Применение программной системы AWA, в которой реализован авторский вариант метода Лонера QR-вращения, не позволило довести до конца расчеты ни примера 1, ни примера 2, приведенных ниже. Гарантированные границы множества решений системы “брюсселятор” сравнивались с двусторонними оценками решений из [6].

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

Рассматривается система ОДУ

-У^ = 2 (У4(() K.W. ^д(г) = 10 tyt(t)y’t(t),

(2)

-»>(t) ^ -y4(t) ,w

ST = 2ty4(t)- ~ = -2t(y3 (t) - 1)1

для которой поставлена задача с неточно заданными начальными данными (принадлежащими некоторому интервалу).

Считая, что начальные данные — некоторые символьные переменные

yi(0) = Y0. y2(0) = Y0 уз(0) = Y0, y4(0) = Y0, (3)

выписываем точное решение

Y /.) = Yi0 exp(- cos(t2) + cos(t2)Y° + Y4 sin(t2))

1 cos h(-1 + Y°) + sin h(-1 + Y3) ’

Y2(t) = (Y0)5 exp(-5 cos(t2) + 5 cos(t2)Y° + 5Y40 sin(t2)) —

-(Y0)5 cos h(-5 + 5Y0) - (Y))5 sin h(-5 + 5Y30) + Y20 cos h(-5 + 5Y3°) + (4)

+Y0 sin h(-5 + 5Y0) / [cos h(-5 + 5y°) + sin h(-5 + 5y°)],

Y3(t) = - cos(t2) + cos(t2)Y3 + Y° sin(t2) + 1,

Y4(t) = sin(t2) - sin(t2)Y0 + Y4 cos(t2).

Полученные формулы точного решения могут быть эффективно использованы для описания границ множества точных решений системы (2) при изменении начальных данных системы в пределах заданных интервалов. Для этого решается несложная задача нахождения экстремальных значений отображения, заданного (4) при t = 2. Однако слишком трудоемко продолжать решения подобных экстремальных задач при всех значениях t. Более того, подобный подход бесполезен для подавляющего большинства практических задач, поскольку точное решение для них неизвестно. Полученные значения Y1 (2), Y2 (2), Y3(2), Y4(2) позволили сравнить числовые значения найденных интервалов с границами множества точных решений при t = 2. Приведенная ниже табл. 1 содержит числовые значения гарантированных границ множеств решений системы ОДУ из четырех уравнений (2) (вторая строка) и числовые значения границ множества точных решений (третья строка).

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

Таблица 1

Числовые значения гарантированных границ множеств решений системы ОДУ из четырех уравнений (2)

t Yi Y2 Y3 Y4

t0 = 0 t200 = 2.0 t = 2.0 [0.8; 1.2] [0.2817; 0.6292] [0.2744; 0.6178] [0.8; 1.2] [-8.0723;31.0065] [-6.9251;30.1800] [0.8; 1.2] [0.0893; 0.4011] [0.0012; 0.3952] [0.8; 1.2] [-0.8144; [-0.6330; -0.5048] -0.4291]

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

Время, затраченное на вычисление гарантированных границ множеств решений системы (2), для t <Е [0; 2] составляет 5.3 с на компьютере с процессором Intel Pentium 200 MMX.

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

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

Приведенные ниже графики проекций гарантированной оценки (множества значений): рис. 1, 3, 5, 7 отображают поведение границ гарантированной оценки множества решений ОДУ с интервальными начальными данными [0.8; 1.2] х [0.8; 1.2] х [0.8; 1.2] (в проекциях

Рис. 1. Проекция гарантированных границ Рис. 2. Графики точных решений системы

множества решений системы (2) на плос- ОДУ (2), начинающихся в параллелепипеде

кость £,у2. Y0-проекции на плоскость £,У2-

Рис. 3. Проекция гарантированных границ множества решений системы (2) на плоскость *,У3-

Рис. 4. Графики точных решений системы ОДУ (2), начинающихся в параллелепипеде Y0-проекции на плоскость уз.

Рис. 5. Проекция гарантированных границ множества решений системы (2) на плоскость *,У4.

Рис. 6. Графики точных решений системы ОДУ (2), начинающихся в параллелепипеде Y0-проекции на плоскость у4.

Рис. 7. Проекция гарантированных границ множества решений системы (2) на плоскость У1,У2,У3.

Рис. 8. Графики точных решений системы ОДУ (2), начинающихся в параллелепипеде Y0-проекции на плоскость У1, У2, уз.

на координатные плоскости) при изменении £ на интервале от 0 до 2.

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

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

На рис. 9, 10 изображены проекции гарантированной оценки множества решений на плоскости У1,Уз, У4 и У2,У3,У4-

П р и м е р 2. Гарантированные оценки решений системы ОДУ, описывающей модель “брюсселятор”.

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

(5)

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

Рис. 9. Проекция гарантированных границ множества решений на плоскость у\, уз, у4.

Рис. 10. Проекция гарантированных границ множества решений на плоскость У2,Уз,У4.

Рис. 11. Графики гарантированных границ множества решений системы, носящей название “брюсселятор” и описывающей модели химической реакции, в которой концентрация реагентов меняется во времени: проекции границ на плоскость Ь,у\ (а) и Ь,у2 (б); интервальный вектор начальных данных равен [—0, 001; 0, 001] х [—0, 001; 0, 001].

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

поведения решений системы “ брюсселятор”, а также точное включение множества точных решений. В качестве сравнения полученных гарантированных оценок привлекались результаты численных расчетов, приведенные в работе [6]. Двусторонние оценки в этой работе получены с помощью некоторой разновидности метода рядов Тейлора и достаточно точного разложения и оценивания сверху и снизу правой части исходной системы. Эта методика делает результаты, полученные и напечатанные в [6], надежными. Однако это не защищает метод, изложенный в статье [6], от влияния так называемого '№гарр1^-эффекта [12, 18]. На интервале от 0 до 5 гарантированные границы, полученные в символьном методе, и оценки работы [6] совпадают, затем в результатах [6] обнаруживается катастрофический рост границ оценки под влиянием '№гарр1^-эффекта [7]. Гарантированные границы не подвержены влиянию этого эффекта. На рис. 11, 12 отображены результаты расчетов на интервале от 0 до 20. На этих рисунках пучок траекторий начинается с интервала начальных данных, имеющего ненулевую ширину. Пример свидетельствует

Таблица 2

Числовые значения гарантированных границ множеств решений системы ОДУ (5) (модель “брюсселятор” )

1 Ух У 2

£и =0 [-0.01; 0.02] [-0.01; 0.02]

5 2. 0. Ю 2 [0.1563498; 0.1672700] [0.0573972; 0.1011260]

£50 = 0.5 [0.217113; 0.2214096] [0.2008253; 0.2491820]

£100 = 1.00 [0.2512616; 0.2525505] [0.5485139; 0.5986268]

со о о 0 О [0.28576071; 0.28703354] [1.97658802; 2.024520273]

о о 0 О [0.3065726828; 0.308223171] [2.660495993; 2.706604232]

сл о о 5. 0 О [0.33440087; 0.33669275] [3.31337170; 3.35689718]

£600 = 6.00 [0.37510527; 0.37868067] [3.91977277; 3.95915558]

£900 = 9.00 [3.18777105; 3.48042329] [1.23234796;1.70489849]

£1000 = 10.00 [1.6412143; 1.7527928] [1.4171911; 1.4823045]

£1900 = 19.00 [0.3838418; 0.3852392] [3.2064150;3.2212531]

ю о о о 2 0 о [0.377646211; 0.378265766] [3.776662503; 3.789646006]

Таблица 3

Двусторонние оценки решений системы (5), представленные в [6]

t Yi Y2

t = 0.25 [0.158118; 0.158118] [0.0688847; 0.0688847]

t = 0.5 [0.217156; 0.217156] [0.211771; 0.211779]

t = 1.0 [0.251224; 0.251224] [0.558443; 0.558443]

t = 3.0 [2.85991; 0.285885]; [1.98496; 1.98496]

t = 4.0 [0.306753; 0.306975] [2.66802; 2.66819]

t = 5.0 [0.32805; 0.341509] [3.31503; 3.32521]

t = 6.0 [-0.0940371; 0.713255] [3.6932; 4.302599119]

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

Время, затраченное на вычисление гарантированных оценок множеств решений системы (5) для t 6 [0, 20], составляет 11.3 с на компьютере с процессором Intel Pentium 200 MMX.

Анализируя числовые значения из табл. 2 и 3, легко увидеть, что в начале вычислений автор статьи [6] полагает известными точные значения решения системы ОДУ (5), что практически невозможно. Нужно стремиться к тому, чтобы получить строгие включения точных решений. Это само по себе сложно, но выполнимо. Далее, начиная счет с интервалов, ширина которых равна величине ошибки округления машинной арифметической операции (^ 10-7, ~ 10-8), к моменту времени t = 6 ширина интервала становится равной 0.7, т. е. возрастает в 107 раз. Влияние wrapping-эффекта в задачах построения двусторонних и интервальных оценок и отсутствие роста границ в гарантированных методах демонстрируют табл. 2 и 3.

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

[1] РОГАЛЕВ А.Н., Шокин Ю.И. Исследование и оценка решений обыкновенных дифференциальных уравнений интервально-символьными методами // Вычисл. технологии. 1999. Т. 4, № 4. C. 51-76.

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

[3] РОГАЛЕВ А.Н. Исследование практической устойчивости при постоянно действующих возмущениях // Там же. Ч. 5. C. 148-150.

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

[5] РОГАЛЕВ А.Н. Включение множеств решений и гарантированные оценки глобальной ошибки // Там же. 2003. Т. 8, № 6. C. 80-94.

[б] Некрасов С.А. Двусторонние методы численного интегрирования начальных и краевых задач // Журн. вычисл. математики и мат. физики. 199Б. Т. ЗБ, № 1G. C. 14Тб-149З.

[Т] Lohner R.J. Program AWA (Anfangs-Wert-Aufgabe). 1994.

Available at ftp://iamk4515.mathematik.uni-karlsruhe.de/pub/awa

[8] Anguelov R. Wrapping effect of the initial value problems for ODE: Applications // Reliable Computing. 1999. Vol. Б, N 2. P. 14З-1б4.

[9] 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. 3б1-3б9.

[1G] Berz M., Makino K. and Hoefkens J. Verified Integration of Dynamics in the Solar System // Nonlinear Analysis: Theory, Methods and Appl. 2GG1. Vol. 4Т. P. 179-190.

[11] Nedialkov N., Jackson K., Corliss G. Validated solutions of initial value problems for ordinary differential equations. Technical Report, 1998. 4G p.

Available at http://www.cs.toromto.edu/NA/reports.html

[12] Nedialkov N., Jackson K., Corliss G. Validated solutions of міпі^і value problems for ordinary differential equations // Appl. Math. and Comp. 1999. Vol. ЮБ, N 1. P. 21-б8.

[13] 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. Б, N Б. P. 289-31G.

[14] NEDIALKOV N., Jackson K. Some recent advances in validated methods for IVP’s for ODE’s // Ibid. P. 289-31G.

[1Б] NEDIALKOV N., Jackson K. The design and implementations of an object-oriented validated ODE solver. Technical Report. 2GG1. 2б p.

Available at http://www.cs.toromto.edu/NA/reports.html

[1б] 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. 2GG1. Vol. Т, N б. P. 449-4бБ.

[1Т] 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, 2GG1. P. 219-2б4. Available at http://www.cs.toromto.edu/NA/reports.html/ ned.scanGG.ps.Z

[18] Lohner R. On the ubiquity of the wrapping effect in the computation of error bounds // SCAN-2GGG, Interval-2GGG. GAMM-IMACS Intern. Symp. on Sci. Computing, Computer Arithmetic and Validated Numerics, Sept. 18-22, 2GGG. Univ. Karlsrue, Institut fur Angewandte Mathematic (Germany). P. Зб.

[19] Калмыков С.А., Шокин Ю.И., Юлдашев З.Х. Методы интервального анализа. Новосибирск: Наука, 198б. 221 с.

[20] RoGALEv A.N. On optimal interval bounds for reachable set ot dynamic systems // IMACS-GAMM Intern. Conf. on Interval Methods and Computer Aided Proofs, Sept. 29 — Oct. 2, 1996. Wurzburg Univ. (Germany). Editorial Board of the Intern. J. Reliable Computing. P. 94-96.

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