ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 6, с. 978-984
УДК 519.62
НЕЯВНЫЙ МЕТОД ПЯТОГО ПОРЯДКА ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ ДИФФЕРЕНЦИАЛЬНО-АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
© 2015 г. Л. М. Скворцов
(105005 Москва, 2-я Бауманская, 5, МГТУ) e-mail: lm_skvo@rambler.ru Поступила в редакцию 29.11.2013 г. Переработанный вариант 12.01.2015 г.
Для численного решения дифференциальных и дифференциально-алгебраических уравнений предлагается использовать неявный двухшаговый метод Рунге—Кутты пятого порядка. Расположение узлов нового метода позволяет получить оценки старших производных в начальной и конечной точках шага интегрирования. Поэтому такой метод можно рассматривать как конечно-разностный аналог метода Обрешкова. Результаты численных экспериментов, некоторые из которых приведены в статье, показывают, что метод сохраняет порядок при решении жестких уравнений и уравнений индексов 2 и 3. Это является основным его преимуществом по сравнению с известными методами. Библ. 9. Табл. 3.
Ключевые слова: неявные методы, жесткие уравнения, дифференциально-алгебраические уравнения, индекс дифференцирования, феномен снижения порядка.
DOI: 10.7868/S0044466915060125
1. ВВЕДЕНИЕ
Численное решение жестких обыкновенных дифференциальных уравнений (ОДУ) и дифференциально-алгебраических уравнений (ДАУ) часто сопровождается явлением, которое получило известность как "феномен снижения порядка" (см. [1]). Это явление проявляется в том, что реальный порядок сходимости оказывается ниже порядка, полученного исходя из классических представлений. Снижение порядка тем заметнее, чем больше разность между классическим и стадийным порядком (стадийный порядок — наименьший порядок на всех стадиях метода, включая формулу интегрирования). Поэтому для расчетов с повышенной точностью желательно применять методы, имеющие достаточно высокий стадийный порядок. Уменьшить эффект снижения порядка позволяет также использование жестко точных методов, у которых последняя стадия совпадает с формулой интегрирования.
Стадийный порядок одношаговых методов Рунге—Кутты не может быть больше числа стадий. Оптимальные (имеющие наивысший порядок при заданном числе стадий) методы Рунге—Кутты высоких порядков имеют большую разность между классическим и стадийным порядком, что может привести к снижению реального порядка. Среди жестко точных Z-устойчивых методов оптимальными являются методы Радо IIA, которые имеют классический порядок 25 - 1 при стадийном порядке, равном числу стадий 5. Снижение порядка наиболее заметно при решении ДАУ высших индексов (определение индекса дифференцирования ДАУ приведено в [1]). Например, порядок сходимости алгебраических переменных при решении ДАУ индекса 3 методами Радо IIA равен 5 - 1 (см. теорему 6.1 из [2]). Тем не менее решатели RADAU5 и RADAU, реализующие эти методы, показывают превосходные результаты при решении многих тестовых задач. Можно ожидать, что методы, стадийный порядок которых превышает число стадий и равен классическому порядку, окажутся еще более эффективными.
Стадийный порядок может быть больше числа стадий в многошаговых методах либо в методах со старшими производными, к которым относятся методы Обрешкова (см. [3]). В предлагаемом методе сочетаются оба эти подхода. За основу взят метод Обрешкова 5-го порядка
Уп+1 = Уп + ЦЬу + a{y'n+i) + h 2(b2J« + a2y'n+i) + hay'+i, (1)
b = 2/5, b2 = 120, a1 = 3/5, a2 = -3/20, a3 = 1/60.
Его функция устойчивости
*(z) = 60 + 24z + 3z2 3 (2)
60 - 36z + 9z - z
является аппроксимацией Паде (2, 3) экспоненты и совпадает с функцией устойчивости трехста-дийного метода Радо IIA 5-го порядка. Недостатком метода (1) является необходимость вычисления старших производных, поэтому методы такого типа редко применяют на практике. Для устранения этого недостатка аппроксимируем старшие производные через конечные разности. В результате получим двухшаговый трехстадийный метод Рунге—Кутты 5-го порядка, стадийный порядок которого совпадает с классическим порядком.
2. ПОСТРОЕНИЕ МЕТОДА
Для аппроксимации 2-й и 3-й производных в начале и в конце интервала [[n, tn+1] выберем абсциссы метода в виде c1 = 1 -а, c2 = 1 + а, c3 = 1, где значение а достаточно мало (аналогично методам, предложенным в [4]). Формулы шага интегрирования при решении задачи Коши
y' = f (t, y), y (to) = Уо
запишем в виде
Yt = yn + h
Z aYl +Z bY
j=1 j=1
, Y; = f(tn + ch,Yi), t = 1,2,3, yn+1 = Y3, (3)
где У у — векторы стадийных значений, полученные на предыдущем шаге. Такой метод относится к классу двухшаговых методов Рунге—Кутты (см. [5]).
Чтобы ограничиться оценкой 2-й производной при t = tn, примем
Ьп = -Ьа> 1 = 1,2,3. (4)
Из условий 5-го стадийного порядка получим
3 3 к
Уа^к'1 + УЬк = ^, 1 = 1,2,3, к = 1,...,5,
~ ~ 1 к (5)
с1 =-а/м>, с2 = а/м>, с"3 = 0,
где w — отношение размера текущего шага к размеру предыдущего шага. Из уравнений (4), (5) находим коэффициенты ау, Ьу (явные формулы для этих коэффициентов довольно громоздки, поэтому здесь не приводятся).
Поскольку метод двухшаговый, то на первом шаге следует использовать подходящий одноша-говый метод, содержащий абсциссы 1 - а, 1 + а и 1. Такие методы были предложены в [4].
3. СВОЙСТВА МЕТОДА
Исследуем устойчивость метода. Обозначим через А и В матрицы коэффициентов ау и Ьу. Решая уравнение Далквиста у' = Ау, получаем У = Н(г)У, где г = Кк, Н(г) = (I - гА)-1 (Е + гВ), I — единичная матрица, Е = [1,1,1]т[0,0,1] . Устойчивость метода определяется корнями характеристического многочлена Р(п,г) = |п1 -Н(г). Метод устойчив, если абсолютные значения этих корней не превышают 1, а среди корней, равных по модулю 1, нет кратных (см. определение У.1.1 из [1]). В нашем случае матрица Н (г) имеет ранг 2, поэтому характеристический многочлен может иметь не более двух ненулевых корней и запишется в виде Р(п,г) = п(П - Р (г)п + Р0 (г)).
При постоянном размере шага (w = 1) имеем
Po(z) = No(z)/D(z), Pi(z) = Ni(z)/D(z), No(z) = a4(8z + 7z2 + z3),
N1(z) = 60 + 120a2 + (24 + 24a2 + 16a4)z + (3 - 6a2 - 8a4)z2 - (2a2 + 8a4)z3,
D(z) = 60 + 120a2 - (36 + 96a2 - 8a4)z + (9 + 30a2 - 15a4)z2 - (1 + 4a2 - 5a4)z3.
Такой метод ^-устойчив при |a| < 0.498 и ^(89°)-устойчив при |а| < 0.553. При а = 0.5 он имеет 6-й порядок. При w ф 1 (т.е. при убывании либо возрастании размера шага) ^-устойчивость может нарушаться, но если значение w не очень велико, то сохраняется ^(ф)-устойчивость с углом ф, близким к 90°. Например, если а = 0.1, то ф = 87° при w < 3, ф = 88.3° при w < 1 и ф = 90° при 1 < w < 2.55.
Если а ^ 0, то p0(z) ^ 0, p1 (z) ^ R(z), где R(z) — функция устойчивости (2). Это означает, что при решении линейной автономной системы ОДУ новый метод с малым значением а практически равноценен по точности и устойчивости методам 5-го порядка (1) и Радо IIA. Формулу интегрирования можно записать в виде
Уп+1 = Y = Уп + h
IL_IL+ö y; гщ + Y + ßY' + ß2 Y Y-
2a a 2 a/ w
где aL = a31+a32 + a33, a2 = a(a32- a31), a3 = a (a31 + a32)/2, ßL = b33, ß2 = -2ab31/w. При a ^ 0 имеем
а = ai + O(a2), ßt = bt + O(a2), = hy+i + O(a2)),
2a
72' 2123' + 71' = Ь2(у:+1 + 0(а2)), = Н(у: + 0(а2)),
а 2 а/ м
где аI, Ь — коэффициенты из (1). Таким образом, при а ^ 0 абсциссы с1, с2 и с3 стягиваются в один кратный узел, что эквивалентно использованию старших производных, а в пределе получаем метод Обрешкова (1).
Перейдем теперь к решению системы ДАУ
У' = /(/,У,z), у(?о) = Уо, (6а)
0 = ^(/, у, z), z (?о) = Zо, (6б)
где / и g — достаточное число раз дифференцируемые векторные функции, размерности которых совпадают с размерностями векторов у и г. Предположим, что начальные условия согласованы (см. [1]). Воспользуемся методом е-вложения из [1], основанном на представлении уравнения (6б) в виде sz' = g(?,у,z) при б —> 0. Рассмотрим решение этой системы с помощью метода Обрешкова (1). Обозначив через gn, g'n и ¿П значения функции gи ее производных при ? = ¿П, получим
е - Zn) = + a^g„+l) + Ь\b2gn + а2^п+1) + ^п'+Р
тогда +1 = е-1gn+1, zn'+1 = е-1gn+1, zn'+1 = е~£П'+1. Если численное решение системы (6) существует, то значения zn+1, zn'+1 и ¿п+\ должны иметь конечные пределы при б ^ 0, а это возможно только когда gn+1 = gn+1 = gn'+1 = 0. Таким образом, один шаг численного решения уравнений (6) методом (1) запишется в виде
Уп+1 = Уп + К ь1/п + а1/п+1) + ь ХЬ/ + а2/п+1) + Ь а3/п+1,
(7а)
0 = gn+l, 0 = gn+l, 0 = gn,+l. (7б)
Размерность полученной системы нелинейных алгебраических уравнений (7) больше размерности исходной системы (6), но дополнительные уравнения 0 = gn+1 и 0 = gП,+1 необходимы для нахождения векторов zn+1 и zn'+1, которые входят в выражения для /п'+1, /п'+1, gn+1 и gn'+1. Например, /' = / + /У/ + /^', где /, /У и / — вектор и матрицы соответствующих частных производных.
Обсудим формулы (7). Известно, что жестко точные методы обеспечивают точное выполнение алгебраического соотношения (6б). Благодаря этому они позволяют решать ДАУ индекса 1 без снижения порядка (см. теорему У1.3.3 из [1]). Наше предположение состоит в том, что выполнение соотношений (7б) в методе со старшими производными позволяет решать без снижения порядка также и ДАУ индексов 2 и 3. Это предположение основывается на проведенном в [4] исследовании и подтверждается при решении тестовых задач. Но строгого доказательства у нас нет.
Рассмотрим теперь решение ДАУ (6) предложенным методом. В этом случае формулы шага интегрирования запишутся в виде
y; = f ( + ch y, Z),
yn+1 = Y3, ¿n+1 = Z3 ■
(8)
3 3
У = у« + К X аУк +Х ЬкУк
_ у=1 у=1 _
0 = я (п + сК, У, ^), ¡ = 1,2,3,
Преимущество этих формул по сравнению с (7) — отсутствие старших производных, которые фактически заменяются конечными разностями. Численные эксперименты показали, что при малых значениях а методы (7) и (8) показывают близкие результаты без какого-либо снижения порядка.
При уменьшении а увеличиваются абсолютные значения коэффициентов (наибольшие по
модулю коэффициенты а3 ~ -а~2/30). Как следствие, возрастают вычислительные ошибки, что может потребовать повышенной разрядности вычислений. Поэтому был проведен ряд экспериментов с целью нахож
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.