ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 6, с. 986-1007
УДК 519.622
ВЛОЖЕННЫЕ СИММЕТРИЧНЫЕ НЕЯВНЫЕ ГНЕЗДОВЫЕ МЕТОДЫ РУНГЕ-КУТТЫ ТИПОВ ГАУССА И ЛОБАТТО ДЛЯ РЕШЕНИЯ ЖЕСТКИХ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ И ГАМИЛЬТОНОВЫХ СИСТЕМ1
© 2015 г. Г. Ю. Куликов
(СЕМАТ, Institute Superior Тёетео, Universidade de Lisboa, Av. Rovisco Pais, 1049-001 Lisboa, Portugal)
e-mail: gkulikov@math.ist.utl.pt Поступила в редакцию 29.07.2014 г.
Предлагается методика построения неявных гнездовых методов Рунге—Кутты для решения обыкновенных дифференциальных уравнений, которые относятся к классу мононеявных формул этого типа. Их отличительной чертой является высокая практическая эффективность, которая вытекает из сохранения размерности исходной системы дифференциальных уравнений при применении, что невозможно в неявных многостадийных формулах Рунге— Кутты общего вида. С другой стороны, неявные гнездовые методы Рунге—Кутты наследуют все наиболее важные свойства общих формул этого вида, такие как ^-устойчивость, симметрия, симплектичность в определенном смысле. Кроме того, они могут иметь достаточно высокие стадийные и классические порядки, а также обеспечивать плотную выдачу результатов интегрирования той же точности, что и порядок основного метода, без больших дополнительных затрат. Таким образом, гнездовые методы эффективны для численного интегрирования дифференциальных уравнений самых разных видов, к которым относятся жесткие и нежесткие задачи, а также гамильтоновы системы и обратимые уравнения. Настоящая статья посвящена обобщению предложенных ранее гнездовых методов, основанных на квадратурных формулах Гаусса, на методы типа Лобатто. Более того, в статье представлена единая методика построения всех таких методов. Ее работоспособность продемонстрирована на вложенных примерах неявных гнездовых формул разного порядка. Все построенные методы снабжены механизмами для оценки локальной ошибки численного интегрирования и автоматической генерации переменных сеток на отрезке интегрирования за счет выбора оптимального шага. Такие вычислительные процедуры апробированы на тестовых задачах с известным решением и исследованы в сравнении со встроенными решателями системы матричных вычислений MATLAB. Библ. 73. Фиг. 3.
Ключевые слова: численные методы, решения обыкновенных дифференциальных уравнений, неявный гнездовой метод Рунге—Кутты, вложенные формулы, методы Рунге—Кутты типов Гаусса и Лобатто, автоматический контроль точности, жесткая задача Коши, гамильтоно-вы уравнения.
DOI: 10.7868/S0044466915030114
1. ВВЕДЕНИЕ
Минимизация вычислительных затрат, связанных с решением линейных систем, — одно из главных средств для повышения практической эффективности неявных методов Рунге—Кутты при решении жестких дифференциальных задач Коши большой размерности
X'(t) = g(t, x(t)) , t 6[t0, tend] , t0) = / , (1.1)
где правая часть g : Rm +1 —- Rm обладает необходимой гладкостью. Здесь и далее предполагается, что система дифференциальных уравнений (1.1) имеет единственное решение на всем отрезке
Работа выполнена при финансовой поддержке португальского фонда науки и технологии (Fundagao para a Ciencia e a Tecnologia) в рамках проекта Pest-OE/MAT/UI0822/2011 и программы Investigador FCT 2013.
интегрирования [/0, ?епй]. Следует отметить, что необходимость в решении таких задач возникает естественным образом в ряде прикладных исследований (см., например, [1]—[4]), а также при использовании более сложных вычислительных процедур, включая метод стрельбы для решения краевых задач, алгоритмы калмановской фильтрации и метод прямых для дискретизации задач математической физики (см., например, [5]—[14]). Поэтому разработке эффективных вычислительных формул и программ для решения уравнений (1.1) отводится важное место в отечественной и зарубежной вычислительной математике (см., например, [14]—[26]).
Как сказано выше, тематике удешевления вычислений при численном решении жестких дифференциальных уравнений большой размерности вида (1.1) уделялось и уделяется серьезное внимание в мировой литературе, а основной подход к реализации неявных методов Рунге-Кут-ты, снижающий указанные выше вычислительные издержки, заложен еще в [27] и [28]. Главная идея этих статей состоит в учете структуры матрицы коэффициентов А /-стадийного неявного метода Рунге-Кутты (РК-метода). Напомним, что все такие методы для решения задачи (1.1) могут быть записаны следующим образом:
= + ТС, (1.2а)
Хи; — Хг.
Kki
j — 1 l
+ Tk£aijg(tkj, xkj), i — 1, 2,..., l, (1.26)
— Xk + Tk£bjg(tkj, Xkj), k — 0, 1, ..., K- 1, (1.2в)
vk + 1 j — 1
где х0 = х0, а тк, как обычно, — размер к + 1-го шага численного интегрирования. Понятно, что любая РК-формула (1.2) однозначно задается вещественными коффициентами йц, Ъц и с, ¡,Ц = 1, 2, ..., /, которые принято оформлять в виде таблицы Бутчера
b1
в которой A обозначает матрицу размера l х l, a, b и c есть /-мерные векторы (см., например, [20], [26]). Неявность означает, что некоторые (или все) коэффициенты в верхней треугольной части матрицы A, включая главную диагональ, отличны от нуля при любой перестановке строк и столбцов, что эквивалентно изменению порядка уравнений и переименованию неизвестных в системе (1.26). В этом случае формулы (1.26) определяют систему нелинейных уравнений, которые необходимо решать на каждом шаге метода (1.2). Тогда если дифференциальная задача имеет большую размерность, то затраты на решение соответствующих нелинейных уравнений могут быть очень значительными, так как нелинейная система (1.2б) в общем виде имеет размер l х m, где l - число стадий в РК-формуле (1.2), а m - размер исходной дифференциальной задачи (1.1). Заметим, что возникающие нелинейные задачи требуют применения того или иного итерационного метода ньютоновского типа, а значит, приводят к большим линейным системам, решение которых по тем или иным причинам может быть неприемлемо на практике. В связи с этим и возникает необходимость в оптимизации решения таких линейных задач с целью максимального сокращения вычислительных издержек.
Подход, предложенный в [27], [28], приводит к замене решения одной большой линейной системы размера l х m на последовательное решение lлинейных задач размера m, что, очевидно, менее затратно. На этом принципе основано большое число эффективных неявных формул Рунге-Кутты, которые условно можно разделить на три группы: 1) однократно диагонально-неявные методы Рунге-Кутты, т.е. ОДНРК-методы (Singly diagonally implicit Runge-Kutta methods), которые требуют, чтобы матрица коэффициентов A была нижней треугольной с ненулевыми, но одинаковыми элементами на главной диагонали (сюда же относятся и ОДНРК-формулы с первой явной стадией); 2) однократно неявные методы Рунге-Кутты (Singly implicit Runge-Kutta methods), или ОНРК-методы, которые подразумевают, что их матрица коэффициентов A имеет единственное простое вещественное собственное значение, и 3) мононеявные методы Рунге-Кутты (Mono-implicit Runge-Kutta methods), или МНРК-методы, главное достоинство которых состоит в том, что все их стадийные величины (1.2б) могут быть последовательно подставлены в формулу (1.2в) продвижения на шаг. Последнее означает, что размер нелинейной задачи, решаемой на каждом шаге МНРК-метода, совпадает с размером дифференциального уравнения (1.1) вне зави-
симости от числа стадий и порядка такого метода. Наиболее серьезным недостатком МНРК-формул является тот факт, что их матрица Якоби, необходимость вычисления которой диктуется применением метода Ньютона в силу неявности этих РК-схем, представляет собой матричный многочлен, зависящий от матрицы Якоби правой части задачи Коши для системы (1.1). Степень этого многочлена определяется и не превосходит числа стадий в МНРК-формуле (см. процитированные ниже статьи Ката и Сингхал). Таким образом, полное вычисление матрицы Якоби МНРК-мето-дов может быть очень затратным на практике, сводящим на нет все преимущества этого класса неявных РК-формул. Поэтому эффективное использование МНРК-схем возможно только на базе упрощенных итераций Ньютона со специально подобранной матрицей для аппроксимации упомянутой выше матрицы Якоби численного метода, примененного для дискретизации дифференциальных уравнений (1.1). По этой причине эффективной реализации МНРК-формул уделено пристальное внимание в мировой литературе. Уже в [29]—[31] показано, что удачную аппроксимацию для матрицы Якоби МНРК-методов следует искать в виде степенной функции от некоторой подходящей матрицы. Понятно, что в этом случае затраты на вычисление матрицы Якоби сравнимы с затратами на вычисление такой матрицы для правой части системы дифференциальных уравнений (1.1), а одна итерация упрощенного метода Ньютона заключается в последовательном решении линейных систем с одной и той же матрицей коэффициентов, размер которой совпадает с размером исходной задачи (1.1).
По каждому из перечисленных выше направлений исследований известно большое число интересных статьей, среди которых следует отметить [29]—[56]. Мы не будем подробно останавливаться на анализе этих работ, так как подобный анализ с перечислением достоинств этих групп методов и присущих им недостатков выполнен, например, в [18], [21], [26], [53]—[55]. В конце этого раздела укажем только на некоторые новые идеи, которые появились сравнительно недавно и имеют непосредственное отношение к эффективному решению жестких дифференциальных уравнений вида (1.1).
В основу важного направления в построении эффективных численных методов для решения жестких задач Коши для системы (1.1) большой размерности легла идея, заключающаяся в применении явных вычислительных формул с расширенными областями устойчивости. Этот подход был предложен В.И. Лебедевым (см. [18]) и в настоящее время приобрел достаточно большую популярность (см. [22], [57]—[68]). Его суть состоит в том, что для некоторых практических жестких дифференциальных уравнений вида (1.1) условие ^-устойчивости применяемых численных методов, которое и приводит к их неявности, является избыточным и может быть с успехом заменено на боле
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.