ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2012, том 52, № 10, с. 1801-1811
УДК 519.622.1
КОЛЛОКАЦИОННЫЕ МЕТОДЫ РУНГЕ-КУТТЫ ДЛЯ ДИФФЕРЕНЦИАЛЬНО-АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ ИНДЕКСОВ 2 И 3
© 2012 г. Л. М. Скворцов
(105005Москва, ул. 2-я Бауманская, 5, МГТУим. Н.Э. Баумана) e-mail: lm_skvo@rambler.ru Поступила в редакцию 13.09.2011 г.
Рассматриваются жестко точные коллокационные методы Рунге—Кутты с явной первой стадией. Выбор параметров этих методов осуществляется исходя из условия минимизации ошибок решения дифференциально-алгебраических уравнений индексов 2 и 3. В результате построены методы, превосходящие по точности известные методы Рунге—Кутты при решении таких уравнений. Библ. 22. Табл. 5.
Ключевые слова: дифференциально-алгебраические уравнения, индекс дифференцирования, неявные методы Рунге—Кутты, феномен снижения порядка.
Будем решать систему дифференциально-алгебраических уравнений (ДАУ), заданную в полунеявной форме:
при согласованных начальных условиях (см. [1]). Неавтономную систему всегда можно привести к виду (1.1), добавив уравнение V = 1. Предположим, что функции Г и g дифференцируемы необходимое число раз, а размерности вектора г и функции g совпадают. Согласно определению из [1], индекс дифференцирования системы (1.1) есть наименьшее число аналитических дифференцирований, необходимых для того, чтобы из уравнений (1.1) посредством алгебраических преобразований получить систему обыкновенных дифференциальных уравнений в явном виде.
Введем обозначение gy = дg (у, г )/ду, gz = дg (у, г )/дг. Продифференцировав (1.1б), получим
0 = gyf (у, г) + gz z'. Если матрица gz обратима в окрестности решения, то г' = 1^^ (у, г) и система (1.1) имеет индекс 1. Наиболее трудны для численного решения системы высших индексов (2 и выше), которые имеют вырожденную матрицу gz. В этом случае алгебраическая подсистема (1. 1б) не может быть решена относительно вектора г, поэтому алгебраические и дифференциальные уравнения решают совместно, используя неявный метод. Многие динамические объекты описываются уравнениями высших индексов, причем наиболее часто встречаются ДАУ индексов 2 и 3. Например, механические системы со связями описываются уравнениями индекса 3 (см. [1], [2]).
В системах индекса 1 можно исключить вектор г, решая аналитически либо численно уравнение (1.1б) относительно г и подставляя полученное решение в (1.1а). В системах высших индексов присутствуют алгебраические переменные, которые невозможно исключить подобным образом. Если матрица gz вырожденна и имеет постоянный ранг, то путем исключения некоторых алгебраических переменных можно привести уравнения (1.1) к виду
1. ВВЕДЕНИЕ
у' = f (y, z), y (to) = y 0 0 = g (y, z), z (to) = z o,
(1.16)
(1.1a)
y' = f (y, z), y (to) = y 0 o = g (y), z (to) = z o
(1.26)
(1.2a)
1801
1802
СКВОРЦОВ
(детали см. в [1]). Продифференцировав (1.2б), получим
0 = (у, г). (1.3)
Если матрица gyfz обратима, то уравнения (1.2а), (1.3) образуют систему индекса 1, а продифференцировав (1.3), можно получить выражение для г'. Таким образом, если матрица gyfz обратима в окрестности решения, то система (1.2) имеет индекс 2. В этом случае начальные условия будут согласованными, если они удовлетворяют уравнениям (1.2б) и (1.3).
Уравнения индекса 3 обычно представляют в виде
у' = f (у, г), у (*0) = у 0, (1.4а)
г' = к (у, г, и), г ) = г о, (1.4б)
0 = g (у), и (/о) = и о, (1.4в)
предполагая, что матрица gyfzku обратима в окрестности решения. Начальные условия этой системы будут согласованными, если они удовлетворяют уравнению (1.4в), а также уравнениям, полученным в результате однократного и двукратного дифференцирования (1.4в). Представления уравнений в виде (1.2) и (1.4) позволяют выделить переменные, имеющие при численном решении разные порядки сходимости. Это переменные индекса 1 (у-компонента), переменные индекса 2 (^-компонента) и переменные индекса 3 (и-компонента).
Один шаг неявного метода Рунге—Кутты для решения задачи (1.1) выполняется согласно формулам
у п+1 = у п + ¿X ¿/V/, У, = у „ + иX ацХ), у; = f (У,, ъ,),
,=1 м
5 5
гв+1 = гп + ¿XЬЪ, Ъ/ = гп + ¿Xа^Ъ), 0 = g(У,,Ъ,), / = 1,2, ..., 5,
/=1 )=1
которые задают систему алгебраических уравнений. Коэффициенты метода а,), Ь, ,,) = 1,2,..., 5, удобно представить в виде таблицы Бутчера:
С1 a11 • ■ ak
c A
cs as1 • ■ ass bт
b • ■ bs
где c = Ae, e = [1, ..., 1]т.
При решении ДАУ высших индексов реальный порядок метода может быть значительно ниже классического порядка, т.е. в полной мере проявляется феномен снижения порядка. Основные результаты по сходимости методов Рунге—Кутты при решении ДАУ приведены в [1], [3]—[5]. Используя жестко точные методы, можно в ряде случаев повысить порядок сходимости у-компо-ненты, но порядки сходимости других компонент остаются низкими. Обычно порядок z-компо-ненты не превышает стадийного порядка q, а порядок и-компоненты не превышает q - 1. Поэтому диагонально-неявные методы Рунге—Кутты (DIRK), которые наиболее просто реализуются, но имеют q < 2, мало подходят для решения ДАУ высших индексов. Среди них однократно диагонально-неявные методы (SDIRK) из [1], [6], [7] вообще непригодны для решения уравнений индекса 3, а методы c явной первой стадией (ESDIRK) из [8]—[12] обеспечивают только 1-й порядок сходимости и-компоненты. Путем использования дополнительных условий удалось повысить реальный порядок методов ESDIRK при решении ДАУ индекса 3 до 2-го (см. [13], [14]), но такие методы обеспечивают эффективное решение только при низких требованиях к точности. К тому же разный порядок для разных компонент затрудняет автоматический выбор размера шага интегрирования.
Таким образом, для эффективного решения ДАУ высших индексов следует использовать жестко точные методы, имеющие достаточно высокий стадийный порядок. Условия, обеспечивающие стадийный порядок q, запишутся в виде
Ack-1 = ck/k, bV-1 = 1/k, k = 1,2, ..., q (1.5)
(здесь и далее предполагается покомпонентное возведение вектора в степень). При заданном числе стадий s максимально возможный стадийный порядок q = s имеют коллокационные методы, коэффициенты которых однозначно определяются по вектору узлов c из соотношений (1.5).
Будем рассматривать жестко точные коллокационные методы с явной первой стадией, тогда первая строка матрицы A состоит из нулей, а последняя строка совпадает с bт (из этого следует, что c1 = 0 и cs = 1). Эти методы не требуют вычислений на первой стадии, поскольку она совпадает с последней стадией предыдущего шага. Среди методов такого типа наивысший классический порядок, равный 2s - 2, имеют методы Лобатто IIIA. Однако их функция устойчивости R (z) такова, что |R (да) = 1, поэтому методы Лобатто IIIA не очень хороши для решения жестких и дифференциально-алгебраических уравнений.
В [15], [16] исследовались жестко точные коллокационные методы с явной первой стадией, названные SAFERK (Stiffly Accurate, First Explicit, Runge-Kutta). При заданном числе стадий s они имеют один свободный параметр cs_1, а остальные узлы задаются из условия обеспечения
максимального порядка. При c** 1 < cs-1 < 1, где c** 1 — значение cs_1 для метода Лобатто IIIA, имеем 0 < |R(да) < 1. На ряде жестких задач в [15] было показано, что эти методы превосходят по точности методы Радо IIA и Лобатто IIIA при одинаковом числе неявных стадий.
Мы исследовали методы SAFERK с целью найти оптимальное значение cs _1 для ДАУ высших индексов. Оказалось, что для таких уравнений преимущество имеют методы с близким к 1 значением cs-1 (для жестких задач в [15] использовалось значение, близкое к c*_1). В результате более детального исследования были выделены две группы коллокационных методов, имеющих преимущество при решении ДАУ индексов 2 и 3. При заданном s каждая из этих групп образует од-нопараметрическое семейство со свободным параметром а.
Методы первой группы имеют значения узлов
c = 0, cs-1 = 1 -а, cs = 1, (1.6)
а остальные узлы задаются из условия обеспечения порядка 2s - 3. Эти методы рекомендуется использовать для решения ДАУ индекса 2. Методы второй группы имеют значения
c1 = 0, cs-2 = 1 -a, cs-1 = 1 + а, cs = 1, (1.7)
а остальные узлы задаются из условия обеспечения порядка 2s - 4. Эти методы рекомендуются для ДАУ индекса 3. Численные эксперименты показали, что при небольших значениях а предложенные методы превосходят известные методы Рунге—Кутты при решении ДАУ индексов 2 и 3.
Отметим, что при малом а расположение узлов (1.6) позволяет получить хорошую оценку второй производной решения на последней стадии, а расположение узлов (1.7) позволяет получить также и оценку третьей производной. В пределе при а ^ 0 получаем кратный узел, что соответствует использованию старших производных (коллокационные методы со старшими производными рассматривались в [17]—[20]). Ниже будет показано, что такие методы усиливают свойство жесткой точности применительно к уравнениям высших индексов. Практически это означает отсутствие снижения порядка при решении ДАУ индексов 2 и 3.
2. МОДЕЛЬНЫЕ УРАВНЕНИЯ
Феномен снижения порядка при решении жестких задач удалось объяснить c помощью модельного уравнения, которое предложено в [21] (см. также [1]). Минимизация ошибок решения этого и других простых уравнений позволила построить явные и диагонально-неявные методы Рунге—Кутты, обладающие повышенной точностью при решении жестких задач (см. [10]—[12], [22]). Модельные уравнения для исследования ошибок решения ДАУ высших индексов были предложены в [13], [14]. В этих работах были получены выражения для глобальных ошибок при решении некоторых модельных уравнений. Минимизация этих ошибок позволила повысить реальный порядок методов DIRK при решении ДАУ индексов 2 и 3. Таким образом, простые мо-
1804 СКВОРЦОВ
дельные уравнения оказались удобным инструментом для конструирования методов Рунге—Кутты повышенной точности.
При численном решении ДАУ высших индексов самыми неточными являются алгебраические переменные: ^-компонента в уравнениях (1.2) и и-компонента в (1.4). Исследуем точность этих компонент на примере системы
0 = у -ф(t), уо =ф(t
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.