CC BY
109
Вычислительные технологии
Том 13, № 2, 2008
Алгоритм численного решения кусочно-сшитых систем
В. В. Коровицын, Ю.В. ФРОЛОВА Омский государственный университет, Россия e-mail: korobits@univer.omsk.su
В. Б. МАРЕНИЧ Кальмарский университет, Швеция e-mail: Valeri.Marenitch@hik.se
A PSS algorithm for numerical solving of the piecewise sewing systems is suggested. The algorithm realizes two modes of calculation: in areas of continuity and on the sewing surface. The Runge—Kutta method is used for the continuity mode. The Newton iteration method is used for the sewing mode. The developed PSS algorithm is compared with the classical methods. Usability of the PSS algorithm for piecewise sewing systems is shown.
Введение
Математические модели физических, химических, биологических и социальных систем зачастую описываются одним или несколькими дифференциальными уравнениями. Как правило, эти модели ведут себя адекватно лишь в некоторой ограниченной области фазового пространства. В природе подобные процессы имеют более сложную структуру, и поэтому приходится прибегать к использованию моделей, описываемых дифференциальными уравнениями с разрывной правой частью. Например, модель Кривана и Дей-ла [1] отражает поведение популяций в цепи питания, включающих всеядных особей, которые способны менять стратегию питания. А модель Мелина и Хюльгрена [2] описывает поведение резонансного преобразователя, фазовое пространство которого склеено из четырех областей.
Исследованию теории обыкновенных дифференциальных уранений с разрывной правой частью и их применению посвящено достаточно много работ [3-12]. Среди прочих систем выделяют класс кусочно-сшитых. В частности, понятие кусочно-сшитых систем на плоскости введено в работе Н.Н. Баутина и Е.А. Леонтовича [13].
Большинство практических задач решается численными методами. Исследованию аппроксимации, устойчивости и сходимости численных методов решения обыкновенных дифференциальных уравнений посвящено много публикаций. Однако работы по численным методам решения обыкновенных дифференциальных уравнений в основном рассматривают жесткие или нежесткие непрерывные задачи либо задачи с дифференциально-алгебраическими уравнениями [14-24]. Ряд работ представляет программное обеспечение для численного исследования систем обыкновенных дифференциальных
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.
уравнений [25-27]. Лишь немногие авторы, например, J.R. Cash, A.H. Karp в своей статье [28], уделяют внимание построению алгоритма численного решения с переменным порядком для задач с быстроменяющейся правой частью.
В данной работе предлагается алгоритм PSS для численного решения систем обыкновенных дифференциальных уравнений, имеющих разрыв на некотором подпространстве фазового пространства. Сравниваются погрешность численного решения и вычислительные затраты с традиционными численными методами решения обыкновенных дифференциальных уравнений.
1. Определение решения разрывной системы ОДУ
Рассмотрим автономную систему, описываемую набором n параметров y = (y1, y2,... , yn) в некоторой области D С Rn. Пусть эта область разделена на две части D1 и D2 гиперповерхностью G = {у 6 D : g(y) = 0}, где g(y) : D ^ R — непрерывная функция. Таким образом, заданы области D1 = {y 6 D : g(y) < 0} и D2 = {y 6 D : g(y) > 0}. Поведение системы описывается решением обыкновенного дифференциального уравнения в векторном виде
f = f(y), y 6 D, W
при заданном начальном условии
y(to) = yo, yo 6 D\G. (2)
Функция f задана следующим образом:
f fi (y) при y 6 Di U G, \ f2(y) при y 6 D2 U G,
где fi(y) и f2(y) — непрерывные функции в соответствующих областях определения. Согласно данному определению функции f (y) на гиперповерхности G принимают два разных значения f1(y*) и f2(y*), y* 6 G, которые в общем случае не совпадают.
Решением задачи (1)-(2) будем называть непрерывную функцию y(t), производная которой почти всюду совпадает с f (y). Как правило, в точке пересечения с гиперповерхностью G кривая решения непрерывна и не имеет производной. Но при стремлении к точке y* на G слева касательная к кривой решения совпадает с f1(y*), а справа — с
f2 (y*).
Будем рассматривать системы, для которых выполняется дополнительное условие: для любой точки y* 6 G векторы f1(y*) и f2(y*) лежат по одну сторону от касательной плоскости P к поверхности G в точке y*. В этом случае кривая решения задачи (1)-(2) составляется из траекторий фазового пространства в D1 и D2. Отсутствует так называемый скользящий режим. Предположим также, что начальная точка y0 6 D1. Тогда решение может быть представлено тремя способами:
— решение остается в области D1, y = y1(t) для всех t > t0;
— решение достигает границы D за конечное время T0, y = y1(t) для всех t0 < t < T0;
— решение достигает гиперповерхности G за конечное время t1 .
В третьем случае решение может пересекать гиперповерхность G конечное или бесконечное число раз. В этом случае решение можно представить в виде склеенных траекторий:
y(t) = y2i-1 (t) , t 6 (t2i-2 ,t2i-1), y2i-1(t) : R ^ D1,
f(y)
у(*) = Ы*), * £ (¿2г-1,^2г), У*(*): К ^ I = 1, 2,..., у(^) = У*, У* £ С, ^ = 1,2,...
Примером представленной системы могут служить кусочно-сшитые системы [13]. Одну из них приведем ниже. Она представляет интерес из-за наличия циклического решения, сшитого из обычных нециклических траекторий. Для исследования свойств численных алгоритмов подобные примеры имеют огромное значение. Аналитическое точное решение получить достаточно просто. При этом алгоритм можно исследовать на аппроксимацию к точному решению, а также на сходимость, учитывая, что решение всегда остается в некоторой ограниченной области.
2. Пример кусочно-сшитой системы с циклическим решением
Рассмотрим систему линейных обыкновенных дифференциальных уравнений, сшитую из двух систем на плоскости (у1,у2) вдоль прямой у1 = а:
(У1 (У2
(У1
«п(у1 - С1) + «12 (у2 - (1),
«21 (У1 - С1) + «22(У2 - (1), (У1, У2) £ АЬ
«22(У1 - С2) + «12(У2 - (2),
«21 (У1 - С2) + «11 (У2 - (2), (У1, У2) £ ^2,
(3)
(4)
Здесь А1 «12 = «21 :
= {(У1,У2) : У1 < а}, Аз = {(У1,У2) : У1 > а}. В случае, когда «и = «22 1 решение каждой из систем
(У1
У2 - (1,
У1 - С1, У1 < а, У2 - (2,
У1 - С2, У1 > а,
(5)
(6)
имеет одну особую точку (седло). Если с1 < а и с2 > а, то особые точки лежат в соответствующих клетках А1 и А2. Кроме того, если (1 = (2, то особые точки лежат на прямой ортогональной границе, а их сепаратрисы пересекаются в некоторых точках на этой границе. В этом случае образуется ограниченная область А0, заключенная между сепаратрисами. Решение внутри этой области образует цикл и не выходит за ее пределы. Цикл составлен из двух гладких кривых, склеенных на границе (рис. 1). Кривые решения описываются функциями
У1(*) = А^-* + + Аз,
У2(*) = А^-* + А2в-(4-*) + А4,
(7)
где А1 = ((у* - Аз) + (у* - А4))/2, А2 = ((у* - Аз) - (у* - А4))/2. Величины у*, у* задаются начальными значениями, А3 = с1, А4 = (1 при у* < а и А3 = с2, А4 = (2 при У* > а.
0
Рис. 1. Фазовый портрет системы (5)—(6)
Длина периода решения внутри D0 вычисляется по формуле
T = 2ln
A - а| + у/(A3 - а)2 -
2A-,
3. Вычислительный эксперимент с традиционными методами
Исследуем поведение численного решения на задаче (5)-(6). Очевидно, что в областях Di и D2 известные методы Рунге—Кутты без труда справятся с этой задачей. Но что получается при переходе решения из D1 в D2 и обратно? Именно на границе происходит скачок погрешности решения из-за резкого изменения правой части уравнения. И дело даже не в том, что сильно меняется правая часть, а в том, что малейшие отклонения от траектории приводят к большим отклонениям при переходе через границу G.
Исследуем погрешность численного решения для системы (5)-(6) с заданными параметрами d1 = 0.5, c1 = 0.2, d2 = 0.5, c2 = 0.8, а = 0.5 и начальными данными y1(0) = 0.49999999999, y2(0) = 0.3. Приведем результаты численных экспериментов для четырех явных одношаговых методов решения систем обыкновенных дифференциальных уравнений, реализованных согласно [18]:
1) явный метод Эйлера (Explicit Euler). Одностадийная явная схема с оценкой погрешности методом Ричардсона;
2) метод Рунге—Кутты четвертого порядка (Runge—Kutta 4). Классическая четы-рехстадийная схема с оценкой погрешности по Ричардсону и автоматическим управлением длиной шага интегрирования;
3) метод Дорманда—Принса пятого порядка (Dormand-Prince 5(4)). Восьмистадий-ная схема пятого порядка с оценкой погрешности по вложенной схеме четвертого порядка;
4) метод Рунге—Кутты—Фельберга четвертого порядка (Felberg 4(5)). Наиболее распространенная вложенная схема четвертого порядка с семью стадиями, использующая встроенную схему пятого порядка для оценки локальной погрешности.
Первый эксперимент предназначен для исследования поведения схем на разрывах. Оценивалась относительная погрешность численного решения в точном решении после
преодоления одного цикла длительностью T. В этом случае точное решение попадает в начальную точку, образуя цикл. Значение численного решения показывает, насколько оно отклонилось от точного за один цикл. Все рассматриваемые методы решения имеют алгоритм оценки локальной погрешности и выбора длины шага интегрирования на основе заданной точности. Предполагается, что при уменьшении заданной точности погрешность численного решения также должна уменьшаться, обеспечивая тем самым сходимость численного метода. Однако для разрывных систем это не всегда выполняется. Кроме того, заметно значительное уменьшение порядка точности некоторых методов. Результаты эксперимента приведены на рис. 2. График зависимости изображен в логарифмическом масштабе: по горизонтальной оси откладывается порядок заданной точности (— log 10(loc_tolerance)), по вертикальной — порядок относительной погрешности решения (— log 10(rel_error)), вычисляемой по формуле
rel
error
l|yn - y(T)||
где у(Т) — точное решение в точке £ = Т, уп шагов, соответствующих той же точке, || • ||
llynll
- численное решение, полученное после n евклидова норма.
4 5 6 7 8 9 -loglO(loc_tolerance)
-Explicit Euler -Runge-Kutta 4 ■Dormand-Prince5(4) -Feiberg 4(5)
Рис. 2. Зависимость относительной погрешности решения от заданной локальной точности
2 3 4 5 6 7 -loglO(rel_error)
■ Explicit Euler »Runge-Kutta 4
• Dormand-Prince 5(4)
• Felberg 4(5)
Рис. 3. Оценка вычислительных затрат для достижения заданной точности
Из диаграммы (рис. 2) видно, что метод Рунге—Кутты—Фельберга имеет большую погрешность, чем явный метод Эйлера. При достаточно высокой заданной точности метод Дорманда—Принса ведет себя чуть лучше метода Эйлера. Эксперимент показал, что эти методы не достигают заданной точности. Лучшим оказался классический метод Рунге—Кутты четвертого порядка с оценкой погрешности методом Ричардсона. Он обеспечил заданную точность.
Безусловно, для нахождения численного решения требуется разный объем вычислительных затрат для различных методов. Оценим вычислительные затраты каждого метода для получения решения с относительной погрешностью, не превышающей заданную (результаты даны на рис. 3). Вычислительные затраты оцениваются количеством вычислений правых частей уравнения, которых было достаточно для обеспечения заданной точности. На рисунке логарифм от количества вычислений правых частей уравнения обозначен (log10(fc)). Хуже всех оказался метод Рунге—Кутты—Фельберга, почти одинаковые результаты у методов Эйлера и Дорманда—Принса. Значительно лучше работает метод Рунге—Кутты. Для задачи с разрывной правой частью механизм оценки погрешности методом Ричардсона оказался более пригодным, чем методы вложенных схем.
Хотя метод Рунге—Кутты показал себя лучше других, тем не менее используемый механизм управления длиной шага интегрирования недостаточно успешен для задач с разрывной правой частью. На рис. 4 приведена гистограмма количества шагов интегрирования в зависимости от времени на отрезке [0,Т] при заданной локальной точности ¡ос_£о1етапсе = 10-8. Величина Т — длина периода решения, в нашем примере Т = 3.2188758251282. Решение при заданных начальных условиях начинается около границы и через половину периода Т/2 = 1.6094379125641 достигает границу в первый раз. В этот момент происходит повышение погрешности численного решения, и механизм управления длиной шага уменьшает его, что приводит к росту числа шагов на интервале [1.6; 1.8], представленного на диаграмме самым большим столбиком. Поэтому для улучшения метода интегрирования задач с разрывной правой частью необходимо уделить внимание модернизации механизма выбора шага интегрирования вблизи границы.
Результаты проведенных экпериментов показали, что для преодоления разрывов в
40--
35 30 25 20 15 10
5 -- ——II—II—1_ —II—II——
0 0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2
□ steps
Рис. 4. Гистограмма количества шагов численного решения методом Рунге—Кутты четвертого порядка
правой части системы обыкновенных дифференциальных уравнений лучше подходят методы с оценкой погрешности, основанные на экстраполяции Ричардсона, чем методы с вложенной схемой оценки погрешности. А самый удачный — классический метод Рунге—Кутты четвертого порядка. Тем не менее и этот метод нуждается в улучшении, когда речь идет о задачах с разрывной правой частью.
4. Описание алгоритма
Предлагаем алгоритм поиска решения кусочно-сшитой системы (Piecewise Sewing System) обыкновенных дифференциальных уравнений без скользящего режима.
Алгоритм PSS (поиск решения кусочно-сшитой системы). Обеспечивается поиск решения кусочно-сшитой системы обыкновенных дифференциальных уравнений на интервале времени [а,Ь] для начального значения y0. Решение сохраняется в массивах T и Y. Используется схема Рунге—Кутты четвертого порядка точности для интегрирования уравнения в области непрерывности. Для поиска решения на границе используется метод экстраполяции полиномом Ньютона.
Шаг 0. (Инициализация.) Установить значения y ^ y0, t ^ а, i ^ 1. Сохранить в массив T[0] ^ t, Y[0] = y. Снять все признаки ind_bound ^ false, ind_cham ^ false. Выбрать шаг интегрирования h.
Шаг 1. (Сделать шаг РК4.) Совершить шаг интегрирования h методом Рунге— Кутты четвертого порядка. (Получим приближенное значение решения y1 к y(t + h). Также будет установлен признак границы клетки ind_bound ^ true, если при вычислении правых частей обнаружен переход в другую область.)
Шаг 2. (Решение около границы?) Если ind_bound = true, то перейти на шаг 3, иначе перейти на шаг 6.
Шаг 3. (Решение переходит границу?) Если ind_cham = true, то перейти на шаг 4, иначе перейти на шаг 5.
Шаг 4. (Поиск решения вблизи границы.) Используя таблицу последних двух принятых точек решения T[k], Y[k] (k = i — 2,i — 1) и значений правых частей уравнения в этих точках, построить полином Ньютона N (в) к y(t — Oh). Применяя итерационный метод Ньютона, вычислить значение параметра в < 0 для достижения нужной точности приближения к границе. Определить значения полинома до границы y1 = N(Oi) и за границей y2 = N(O2), выбирая подходящие в1, в2. Сохранить результаты T[i] ^ t — O1h, Y[i] ^ yi, T[i + 1] ^ t — O2h, Y[i + 1] ^ y2, i ^ i + 2 и увеличить время t ^ t — O2h. Снять признаки ind_cham ^ false, ind_bound ^ false. Перейти на шаг 7.
Шаг 5. (Определение оптимального шага до границы.) Вычислить длину шага интегрирования h, обеспечивающую наибольшее приближение к границе области. Установить признак перехода в другую область ind_cham ^ true. Перейти на шаг 7.
Шаг 6. (Сделать два шага РК4.) Совершив два шага интегрирования h/2 методом Рунге—Кутты четвертого порядка, получить y4. Оценить погрешность решения и выбрать новый шаг интегрирования. Если погрешность в пределах заданной точности, то сохранить результат T[i] ^ t + h, Y[i] ^ y4, i ^ i + 1 и увеличить время t ^ t + h.
Шаг 7. (Закончить вычисления?) Если t < b, то вернуться к шагу 1, иначе закончить вычисления.
В основе алгоритма выступает классический метод Рунге—Кутты четвертого порядка. Различие проявляется на линии сшивания областей, там где правая часть уравнения терпит разрыв. При приближении к границе происходит переключение режима
гпд,_Ъоипд,. Приближение определяется попаданием в другую область при вычислении правых частей в схеме Рунге—Кутты. Если хотя бы для одной функции потребовалось значение из другой области, то считаем, что решение приблизилось к границе.
В режиме гпй_Ъоипй (шаг 5) определяется шаг интегрирования, который обеспечивает приближение к границе, при этом остается в той же области. Шаг вычисляется по формуле
Н = а
где уг — принятое значение решения в момент времени ¿¿, д(у) — функция границы, -Уд(у) — градиент к границе, f (у,Ь) — значение правой части в точке у в момент времени ¿, а — коэффициент, определяемый опытным путем. В знаменателе дроби стоит скалярное произведение градиента к границе и вектора правых частей системы, описывающего направление касательной к кривой решения. Таким образом, вся дробь оценивает время, за которое решение преодолеет расстояние до границы. Для того чтобы наверняка не пересечь границу, выбираем коэффициент а = 0.9.
При численном интегрировании методом Рунге—Кутты требуется вычисление правых частей в различных точках фазового пространства и, если одна из точек попадает в другую область, получаемое решение имеет большую погрешность. Эту погрешность приходится компенсировать значительным уменьшением шага интегрирования, как это происходит при традиционном методе Рунге—Кутты. Идея улучшения алгоритма состоит в том, что мы будем использовать найденные точки решения для прогноза поведения кривой решения при приближении к границе, не прибегая к значениям правой части в другой области (шаг 4). Воспользуемся интерполяционной формулой Ньютона [29] для интерполирования назад:
N (в) = у (и - вн) = Уг + (-вн) ( л + (-в) (л - d + (1 - еш - 2d + Л-1))),
где используются два предыдущих шага уг-1 ~ у (и - К), уг ~ у(£г) и значения правых частей в этих точках Шг-1, Ш; значение d = (Ш - Шг-1)/Н; параметр в определяет промежуточные значения. Но мы применяем эту формулу не для интерполяции, а для экстраполяции с отрицательными значениями параметра в. Тогда получим полином Ньютона, приближающий кривую решения не только внутри интервала \Ь г - Н, ¿¿], но и в интервале + Н]. Именно в этом интервале найдется пересечение кривой решения с поверхностью границы. Мы воспользовались формулой интерполирования назад, поскольку она обеспечивает более высокую точность для экстраполирования. Как показали эксперименты, увеличение степени полинома Ньютона не дают выигрыша в точности или скорости сходимости. Построенный полином используется для нахождения пересечения кривой решения с поверхностью границы. Для этого применяем итерационный метод Ньютона с формулой изменения параметра в:
е ~ в - в д^(е))
-,^д(уг)Ш (угЛ)Н"
В этой формуле используем фиксированное значение для скорости изменения кривой решения (знаменатель дроби). Это значение мы определяем в точке у2, такое приближение оказывается вполне допустимым. Коэффициент в выбирается экспериментально от 0.5 до 0.9 для обеспечения наибольшей скорости сходимости. Условием завершения
итерационном процедуры является выполнение выражения
Tol Ав < к
nh
где к — коэффициент от 0.01 до 0.1, То1 —заданная точность, п = т/(1-т), т = Д9/Д92. Величина ДО определяет изменение параметра 9, Дв2 — значение ДО на предыдущем шаге итерационной процедуры. Как показала практика, число итераций зависит от заданной точности и находится в пределах от 5 до 30 итераций для точности от 10-3 до 10-10.
5. Вычислительный эксперимент с алгоритмом PSS
Результаты проведенного эксперимента показаны на рис. 5. На диаграммах продемонстрирована зависимость относительной погрешности численного решения от заданной локальной точности и зависимость вычислительных затрат от погрешности решения. Из рисунка видно, что погрешности обоих методов (классический метод — Runge— Kutta 4, предложенный алгоритм — Algorithm PSS) одинаковы и соответствуют заданной локальной точности. Однако вычислительные затраты на выполнение алгоритмов разные. И они существенно меньше для предложенного алгоритма.
1 2 3 4 5 6 7 8 9 10 0123456789 10 -logl0(loc_tolerance) -logl0(rel_error)
«- Runge-Kutta 4 — Algorithm PSS
Рис. 5. Зависимость относительной погрешности решения и вычислительных затрат от заданной локальной точности для классического метода и алгоритма РЯЯ
□ Algorithm PSS ■ Runge-Kutta 4
Л М М М М М М М М М М М М М М М ГП
О 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2
Рис. 6. Гистограмма сравнения количества шагов численного решения методом Рунге—Кутты и алгоритма РЯЯ
35 30 25 20 15 10 5
Оценка времени решения для метода Рунге—Кутты и алгоритма PSS
Toi Runge— -Kutta 4 Algorithm PSS
100 циклов 1000 циклов 100 циклов 1000 циклов
10-4 13 126 3 26
10-5 17 136 5 42
10-6 22 215 6 62
10-7 28 268 10 98
10-8 33 343 14 153
10-9 44 443 23 241
Снижение уровня вычислительных затрат происходит благодаря вводу дополнительного механизма управления длиной шага при переходе через границу между областями. Резкое уменьшение числа шагов в районе границы продемонстрировано на гистограмме (рис. 6).
В предложенном алгоритме помимо самой схемы Рунге—Кутты требуется составлять полином Ньютона и находить точку пересечения этого полинома с поверхностью границы. Возникает вопрос: время, потраченное на дополнительные вычисления по поиску пересечения кривой решения с границей итерационным методом Ньютона, больше, чем время, потраченное на дополнительные шаги по методу Рунге—Кутты, или нет? Однако эксперименты показывают, что время вычисления полинома значительно меньше времени просчета мелких шагов около границы методом Рунге—Кутты с выбором шага по оценке Ричардсона (результаты экспериментов даны в таблице). Эксперименты проведены для задачи (5)-(6) в одной программе с разными алгоритмами решения. Но поскольку время на поиск решения на отрезке одного периода слишком мало, то для оценки использовалась задача на интевале 100 и 1000 периодов при разной заданной локальной точности Tol. Время дано в секундах при расчетах на компьютере с процессором Intel Pentium-M 1.8 ГГц с оперативной памятью 256 Мбайт.
Результаты экспериментов показывают целесообразность использования предложенного алгоритма PSS для решения кусочно-сшитых систем дифференциальных уравнений с разрывной правой частью.
Список литературы
[1] KRivan V., Diehl S. Adaptive omnivory and species coexistence in tri-trophic food webs // Theoretical Population Biology. 2005. N 67. P. 85—99.
[2] Melin J., Hültgren A. A limit cycle of a resonant converter // Engell H.G.S., Zaytoon J. Analysis and design of hybrid systems 2003. IFAC, Elsevier, 2003. P. 169-174.
[3] АйЗЕРМАН М.А., Пятницкий Е.С. Основы теории разрывных систем I, II // Автоматика и телемеханика. 1974. № 7. С. 33-47; № 8. С. 39-61.
[4] Дерр В.Я., КинЗЕВУЛАТОВ Д.М. Дифференциальные уравнения с обобщенными функциями, допускающими умножение на разрывные функции // Вестн. Удм. ун-та. Ижевск. 2005. № 1. С. 35-58.
[5] Зуев А.В. О периодических решениях обыкновенных дифференциальных уравнений с разрывной правой частью // Математические заметки. 2006. Т. 79, № 4. С. 560-570.
[6] Зуев Л.В., Филиппов В.В. О задаче Неймана для обыкновенного дифференциального уравнения с разрывной правой частью // Дифференциальные уравнения. 2005. Т. 41, № 6. С. 755-760.
[7] МАТРОСОВ И.В. О единственности справа решений невырожденных алгебро-дифферен-циальных уравнений с разрывами // Автоматика и телемеханика. 2007. № 1. С. 11-19.
[8] САМОйЛЕНКО А.М., Перестюк Н.А. Дифференциальные уравнения с импульсным воздействием. Киев: Вища школа, 1987. 288 с.
[9] Сесекин А.Н. О нелинейных дифференциальных уравнениях в классах функций ограниченной вариации // Дифференциальные уравнения. 1989. Т. 25, № 11. С. 1925-1932.
[10] Уткин В.И. Скользящие режимы в задачах оптимизации и управления. М.: Наука, 1981. 368 с.
[11] Филиппов А.Ф. Дифференциальные уравнения с разрывной правой частью. М.: Наука, 1985. 224 с.
[12] Sarrico C.O.R. The linear Cauchy problem for a class of differential equations with distributional coefficients // Portugalie Math. 1995. Vol. 52. P. 379-390.
[13] Баутин Н.Н., ЛЕОНТОвич Е.А. Методы и приемы качественного исследования динамических систем на плоскости. М.: Наука, 1990. 486 с.
[14] АйДА-ЗАДЕ К.Р. О численном решении систем дифференциальных уравнений с нелокальными условиями // Вычисл. технологии. 2004. Т. 9, № 1. С. 11-25.
[15] Деккер К., Вервер Я. Устойчивость методов Рунге—Кутты для жестких нелинейных дифференциальных уравнений. М.: Мир, 1988. 334 с.
[16] Новиков Е.А. Неоднородный метод второго порядка для жестких систем // Вычисл. технологии. 2005. Т. 10, № 2. С. 74-86.
[17] ОлемскоЙ И.В. Явный метод типа Рунге—Кутты пятого порядка // Вычисл. технологии. 2005. Т. 10, № 2. С. 87-105.
[18] Хайрер Э., НЁрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М.: Мир, 1990. 512 с.
[19] Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М.: Мир, 1999. 685 с.
[20] Deuflhard P., Hairer E., Zugck J. One-step and extrapolation methods for differential-algebraic systems // Numerische Mathematik. 1987. Vol. 51, N 5. P. 501-516.
[21] Gupta G.K., Sacks-Davis R., Tescher P.E. A review of recent developments in solving ODEs // ACM Computing Surveys. 1985. Vol. 17, N 1. P. 5-47.
[22] Hairer E., Lubich Ch., Roche M. Error of Runge—Kutta methods for stiff problems studied via differential algebraic equations // BIT Numerical Mathematics. 1988. Vol. 28, N 3. P. 678-700.
[23] Nedialkov N.S., Pryce J.D. Solving differential-algebraic equations by Taylor series (II): Computing the System Jacobian // BIT Numerical Mathematics. 2007. Vol. 47, N 1. P. 121-135.
[24] Park T., Barton P.I. State event location in differential-algebraic models // ACM Transactions on Modeling and Computer Simulation. 1996. Vol. 6, N 2. P. 137-165.
[25] Berzins M., Brankin R.W., Gladwell I. The stiff integrators in the NAG library // ACM SIGNUM Newsletter. 1988. Vol. 23, N 2. P. 16-23.
[26] Enright W.H., Pryce J.D. Two FORTRAN packages for assessing initial value methods // ACM Transactions on Mathematical Software. 1987. Vol. 13, N 1. P. 1-27.
[27] Kamel M.S., Enright W.H., Ma K.S. ODEXPERT: an expert system to select numerical solvers for initial value ODE systems // ACM Transactions on Mathematical Software. 1993. Vol. 19, N 1. P. 44-62.
[28] Cash J.R., Karp A.H. A variable order Runge—Kutta method for initial value problems with rapidly varying right-hand sides // ACM Transactions on Mathematical Software. 1990. Vol. 16, N 3. P. 201-222.
[29] Петровская Н.Б. Применение метода Ньютона в расчетных схемах высокого порядка для стационарных задач // Вычисл. технологии. 2005. Т. 10, № 2. С. 106-113.
Поступила в редакцию 29 октября 2007 г.
