ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2010, том 50, № 6, с. 1047-1059
УДК 519.62
ДИАГОНАЛЬНО-НЕЯВНЫЕ МЕТОДЫ РУНГЕ-КУТТЫ ДЛЯ ДИФФЕРЕНЦИАЛЬНО-АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ ИНДЕКСОВ 2 И 3
© 2010 г. Л. М. Скворцов
(105005 Москва, 2-я Бауманская, 5, МГТУ) e-mail: lm_skvo@rambler.ru Поступила в редакцию 07.10.2009 г.
Рассматриваются диагонально-неявные методы Рунге—Кутты, удовлетворяющие дополнительным условиям порядка. Эти условия позволяют повысить точность решения дифференциально-алгебраических уравнений индексов 2 и 3. На тестовых задачах показано преимущество предложенных методов по сравнению с известными методами. Библ. 16. Фиг. 4. Табл. 7.
Ключевые слова: дифференциально-алгебраические уравнения, индекс дифференцирования, диагонально-неявные методы Рунге—Кутты, методы SDIRK, ESDIRK.
1. ВВЕДЕНИЕ
Рассмотрим систему дифференциально-алгебраических уравнений (ДАУ) вида
y' = f(y, z), y( to) = Уо, (1.1а)
0 = g(y, z), z(to) = zo. (1.1б)
Предположим, что начальные условия y0, z0 согласованны, т.е. удовлетворяют алгебраическому уравнению (1.1б). Согласно определению из [1], индекс дифференцирования системы (1.1) есть наименьшее число аналитических дифференцирований, требующихся для того, чтобы из уравнений (1.1) посредством алгебраических преобразований можно было бы получить систему обыкновенных дифференциальных уравнений (ОДУ) в явном виде. Системы высших индексов (2 и выше) имеют вырожденную матрицу dg/dz. В этом случае невозможно разрешить алгебраическую подсистему (1.1б) относительно вектора z, поэтому алгебраические и дифференциальные уравнения решают совместно, используя неявный метод (т.е. применяется метод s-вложения, см. [1]).
Один шаг неявного метода Рунге—Кутты для решения задачи (1.1) выполняется согласно формулам
s s
y1 = yo + h£ЬХ, Y = yo + h£a^, Y' = f(Y, Z),
i = 1 J = 1
ss
zi = zo + h£bZ', Zi = zo + h£ajZj, o = g(Yi, Zi), i = 1, 2, ..., s,
i = i ; = i
где ay, bi, i, j = 1, 2, ..., s, — коэффициенты метода. Эти формулы задают систему алгебраических уравнений относительно векторов Y', Z', i = 1, 2, ..., s (векторы стадийных значений Yi, Ziнетруд-но исключить).
Метод Рунге—Кутты удобно представить в виде таблицы Бутчера
c1 a11 . a1s
c A
С s aS1 . ass ьт
b . ■ bs
где c = Ae, e = [1, ..., 1]т. При решении жестких и дифференциально-алгебраических уравнений преимущество имеют неявные жестко точные методы, у которых asi = b, i = 1, 2, ..., s. Среди неявных методов наиболее просто реализуются диагонально-неявные (DIRK — Diagonally Implicit Runge—Kutta), матрица A которых имеет нижнюю треугольную форму. Обычно также требуют, чтобы все ненулевые диагональные элементы матрицы A были равны между собой, что позволяет выполнять не более одного LU-разложения на шаге интегрирования.
Будем рассматривать два типа таких методов: жестко точные методы SDIRK (Singly DIRK) с таблицей Бутчера вида
С1 Y
С, - 1 a, -1,1 • ■ Y
1 bi . ■ bs -1 Y
b1 - ■ b, -1 Y
(такие методы предлагались в [1]—[3]) и жестко точные методы ESDIRK (Explicit first stage Singly DIRK) с таблицей Бутчера вида
0 0
c2 a21 Y
С, -1 as -1,1 as -1,2 • Y
1 b1 b2 • ■ b, -1 Y
b1 b2 • ■ b, -1 Y
0 0 0
c a A
b1 bт
(1.2)
(такие методы предлагались в [4]—[9]). В дальнейшем будем опускать последнюю строку таблицы, поскольку она повторяет предпоследнюю строку. Отметим, что функция устойчивости Я(1) жестко точного метода SDIRK удовлетворяет условию Я(<х>) = 0, которое является необходимым условием Х(а)-устойчивости. Для методов ESDIRK также будем требовать выполнения этого условия.
Согласно классической теории, глобальная ошибка при численном решении ОДУ с достаточно малым размером шага ведет себя как 0(Нр), гдер — порядок метода. Однако при решении ДАУ, особенно высших индексов, проявляется феномен снижения порядка, когда реальный порядок метода оказывается ниже классического порядка. Вопросы точности методов Рунге—Кутты при решении систем индексов 2 и 3 рассматривались в [1], [10]—[12], где были получены оценки порядков сходимости различных компонент решения. При выполнении некоторых условий эти оценки могут быть улучшены, однако отсутствует общая теория, позволяющая вывести дополнительные условия заданного порядка.
Будем различать порядок дифференциальных компонент рА и порядок алгебраических компонент ра. Применительно к системе (1.1) это означает, что глобальные ошибки соответствующих компонент при интегрировании с постоянным шагом на заданном интервале имеют оценки
У (К) - Уп = О(НР), г (1п) - гп = О(НРа),
где у^„), г(?п) — точное, а уп, гп — численное решение в конце интервала.
Для исследования сходимости методов Рунге—Кутты рассмотрим простейшие ДАУ, моделирующие поведение различных составляющих ошибки. Эти уравнения имеют известные точные решения, а их численные решения можно получить в виде выражений, содержащих коэффициенты метода. Сравнивая точные и численные решения, можно получить условия, необходимые для достижения заданных порядков дифференциальных и алгебраических компонент. Подчеркнем, что полученные таким образом условия являются необходимыми, но у нас нет доказательств того, что они также и достаточные. Мы рассмотрим методы, удовлетворяющие этим условиям. Результаты решения тестовых задач показали, что эти методы обладают повышенной точностью при решении ДАУ индексов 2 и 3 по сравнению с другими методами такого же типа.
Таблица 1
Дерево Уравнение Хj(t) ~хц(, h)
• х1 = 1 t h Ьте
S х21 = Х1 t2/2 h2bтc
V , 2 х31 = х1 t3/3 h3Ьтc2
> х32 = Х21 t3/6 h3ЬтAc
Все предложенные методы являются Х(а)-устойчивыми и могут быть использованы также и для решения жестких задач.
2. МОДЕЛЬНЫЕ УРАВНЕНИЯ
Определение порядка методов Рунге-Кутты сводится к сравнению разложений в ряд Тейлора точного и численного решений. При выводе условий порядка используют соответствие между элементарными дифференциалами и корневыми деревьями (см. [13]). Эти же деревья будем использовать для формирования модельных уравнений, при этом каждой вершине дерева соответствует определенная переменная.
Принцип построения модельных дифференциальных уравнений изложен в [13, упражнение II.2.4] и в [14, табл. 314.(1)]. В [7] приведены аналогичные жесткие уравнения и ДАУ индекса 1.
Дереву-точке поставим в соответствие переменную, описываемую уравнением х1 = 1. Корневой вершине дерева Ту (i — порядок дерева, т.е. число его вершин, j — номер дерева среди различных деревьев порядка i) поставим в соответствие переменную Ху. Уравнения зададим рекуррентно согласно формуле хотец = Пхсыновья , где "сыновья" — деревья, полученные в результате удаления
корневой вершины дерева-"отца". Начальные значения всех переменных примем нулевыми. При таком подходе дерево можно рассматривать как сигнальный граф, ребра которого ориентированы по направлению к корневой вершине и осуществляют передачу сигналов, а вершины выполняют интегрирование произведения входных сигналов. Предположив, что вершины, не имеющие входов, вырабатывают сигнал x1 = t, получим на выходе корневой вершины дерева Ту переменную Ху.
Полученные уравнения и их решения (xy(t) — точное, x¡j (h) — численное на первом шаге) для деревьев до 3-го порядка включительно приведены в табл. 1. Здесь и далее предполагается покомпонентное выполнение операции возведения вектора в степень, а операцию покомпонентного умножения векторов будем обозначать точкой. При t = 0 разложение в ряд Тейлора переменной Ху содержит только один ненулевой элементарный дифференциал, соответствующий дереву Ту. Обозначим ошибку на первом шаге через 5^Ху) = Ху(И) — x¡j (h), тогда условия, обеспечивающие порядокp метода, запишутся в виде 5^Ху) = 0 при i < p.
На (n + 1)-м шаге численные решения запишутся в виде
*i( tn +1) = tn + h, X21 (tn +1) = X21 (tn) + h b %,
X31(tn + 1) = X31(tn ) + h ЬтХ2, X32 (tn + 1) = X32( tn ) + h ьт X21,
где векторы стадийных значений этого шага имеют вид
X1 = e tn + h c, X21 = ex21( tn) + hAX1,
X31 = e -X31 (tn) + h AX2, X32 = еХХз2 (tn) + hAX21.
Обозначим ошибку переменной Ху после выполнения n шагов через Дп(Ху) = Xy(tn) — х^ (tn). Вычитая из точных решений численные решения, получим выражения для глобальных ошибок в виде
,2
Д« + i(*2i) = Д„(Х21) + \(1 - 2b'с),
h 3
Д« +1 (Х31) = Д«(Х31) + h2tn( 1 - 2bTс) + | (1 - 3bTc2),
Д« +1 (Х32) = Д«(Х32) + ,Д«(Х21) + h2«(1 - 2bTc) + h-3(1 - 6bTAc).
В общем случае переменная Ху интегрируется точно методом порядкаp, если i < p. Если i = p + 1, то глобальная ошибка выражается формулой Дп + 1(Ху) = Дп(Ху) + Sh(xy), где локальная ошибка Sh(xy) пропорциональна hp +1. При интегрировании на заданном интервале с постоянным размером шага глобальная ошибка накапливается и пропорциональна hp, что полностью согласуется с классической теорией.
Определим теперь операцию численного дифференцирования с помощью неявного метода Рунге—Кутты. Система ДАУ индекса 2, осуществляющая дифференцирование переменной ф(0, запишется в виде
У = z, 0 = y - ф(t), Уо = ф(to), Zo = ф'(to), где г — искомая переменная. Используя неявный метод Рунге—Кутты, получаем
У« +1 = У« + hb'Z, Y = ey« + hAZ,
Z« +1 = z« + h Ьт Z', Z = ez« + hAZ', (2.1)
Y = Ф = [ф(t« + C1h), ...,ф(t« + csh)]т.
Исключим векторы стадийных значений. Тогда для методов с обратимой матрицей A, в том числе и SDIRK, получим
У« +1 = аоУ« + ьта-1ф , Z« +1 = ao Z« + аh Лу« + b4-2 ф,
-1 т -2 (2.2)
а0 = 1 - ЬтА e, а1 = -ЬтА е.
Для методов ESDIRK (1.2), а также аналогичных жестко точных методов, имеющих обратимую матрицу А, получим
У« +1 = аоУ« + b ТА-1Ф, z« +1 = ао z« + а1 h + hl ЬтА-2 <Ф, )
ао = 1 - ЬтА-2с, а1 = -ЬтА-3с, <Ф = [ф(t« + C2h),..., ф(t« + csh)]т.
Отметим, что а0 = R(o>), а1 = lim z (R(z) — а0). Формулы вида (2.1)—(2.3) использовались нами
z
при выводе выражений для глобальных ошибок модельных уравнений, при этом вместо переменной ф(0 задавались произведения переменных, входящих в модельные уравнения.
При выводе условий порядка для ДАУ используются деревья с вершинами двух видов (точки и кружки) [1]. Порядком такого дерева называется число вершин-точек минус число вершин-кружков
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.