Журнал вычислительной математики и математической физики, 2020, T. 60, № 7, стр. 1170-1192
Неявные гнездовые методы Рунге–Кутты типа Гаусса и Лобатто с локальным и глобальным контролем точности для жестких обыкновенных дифференциальных уравнений
Г. Ю. Куликов *
Центр вычисл. и стохастической математики,
высший техн. институт, Лиссабонский университет
1049-001 Лиссабон, пр-т Ровишку Паиш, Португалия
* E-mail: gkulikov@math.ist.utl.pt
Поступила в редакцию 31.08.2019
После доработки 31.08.2019
Принята к публикации 10.03.2020
Аннотация
Исследуется задача эффективной оценки и контроля глобальной погрешности неявных гнездовых методов Рунге–Кутты типа Гаусса и Лобатто, используемых при интегрировании жестких обыкновенных дифференциальных уравнений (ОДУ). Жесткие задачи возникают во многих прикладных областях и их точное численное решение является важным вопросом вычислительной и прикладной математики. Эффективный метод оценки глобальной ошибки, недавно разработанный для встроенных методов Рунге–Кутты, может сильно завышать глобальную ошибку при применении к жестким ОДУ и, следовательно, снижать эффективность этих решателей. В настоящей работе исследована причина завышения этой ошибки и показано, как улучшить указанные методы для решения жестких задач. Такие модификации не только повышают эффективность численного интегрирования жестких ОДУ, но и делают встроенные пары неявных гнездовых методов Рунге–Кутты с автоматическим управлением масштабированными модифицированными локальными и глобальными ошибками численного решения превосходящими жесткие встроенные решатели ОДУ в пакете MATLAB с автоматическим управлением локальной ошибкой. Библ. 79. Фиг. 4.
1. ВВЕДЕНИЕ
Обыкновенные дифференциальные уравнения (ОДУ), возникающие в прикладной науке и технике, обычно представляются в виде следующей задачи Коши (ЗК):
(1.1)
$x'(t) = g(t,x(t)),\quad t \in [{{t}_{0}},{{t}_{{{\text{end}}}}}],\quad x({{t}_{0}}) = {{x}^{0}},$Точность численного решения ЗК (1.1) может иметь решающее значение для достоверности математического моделирования исследуемого явления. Например, в [1]–[15] показано, что применение адаптивных ОДУ-решателей с автоматическим управлением погрешности в непрерывно-дискретной калмановской фильтрации позволяет существенно улучшить качество оценивания стохастических систем в задачах слежения, навигации, управления, описания химических реакций, финансового моделирования и т.д., где широко используются стохастические динамические системы с непрерывным временем [16]–[24]. Поэтому настоящая статья посвящена решению этой важной задачи и представляет вычислительные средства для обработки жестких ОДУ (1.1). Эти адаптивные численные методы эффективны для практического использования, а также они превосходят жесткие встроенные решатели ОДУ в пакете MATLAB на ряде важных тестовых задач.
Далее будем считать, что ЗК (1.1) имеет единственное решение $x(t)$ на интервале интегрирования $[{{t}_{0}},{{t}_{{{\text{end}}}}}]$.
Традиционная вычислительная технология для решения уравнения (1.1) опирается на формулы Рунге–Кутты, которые привлекают внимание практиков уже более века. Такие методы эффективны для решения различных проблем, в том числе нежестких и жестких ОДУ, обратимых и гамильтоновых уравнений, задач, возникающих при полудискретизации уравнений в частных производных методом прямых и т.д. (см., например, [25]–[32]). При применении к ЗК (1.1) общая схема метода Рунге–Кутты (РК) имеет вид
где ${{x}_{0}} = {{x}^{0}}$, ${{t}_{{kj}}}: = {{t}_{k}} + {{c}_{j}}{{\tau }_{k}}$, ${{\tau }_{k}}$ – постоянный или переменный размер шага численного интегрирования и векторы ${{x}_{{ki}}}$ являются стадийными (промежуточными) величинами $l$-стадийной схемы РК (1.2). Заметим, что этот метод полностью определен его фиксированными действительными коэффициентами ${{a}_{{ij}}}$, ${{b}_{j}}$ и ${{c}_{j}}$, $i,j = 1,2, \ldots ,l$, которые могут быть удобно представлены в виде таблицы Бутчера где $A$ – квадратная матрица размера $l$, $b$ и $c$ суть $l$-мерные векторы. По форме матрицы $A$ различают подклассы формул РК, т.е. явные, диагонально неявные и полностью неявные схемы. Явные методы являются самыми “дешевыми” и наиболее удобными для численного решения прикладных задач. Однако этот подкласс неэффективен для решения многих жестких ОДУ, которые часто возникают на практике.
Таким образом, неявные схемы РК также необходимы. Более того, некоторые полностью неявные формулы (как, например, методы Гаусса и Лобатто) являются эффективными средствами для решения обратимых и гамильтоновых задач [27]–[31].
К сожалению, явление снижения порядка аппроксимации, которое сопровождает упомянутые выше методы Гаусса и Лобатто, а также их трудоемкость снижают практическую значимость этих методов РК, см., например, [25], [26], [29].
Наше исследование будет посвящено неявным гнездовым методам Рунге–Кутты NIRK (Nested Implicit Runge–Kutta) типа Гаусса и Лобатто, которые наследуют преимущества классических методов Гаусса и Лобатто, но не имеют их недостатков. Эти схемы NIRK впервые были построены и изучены в [33]–[37]. Они представляют собой специальное семейство эффективных моно-неявных формул РК (MIRK), предложенных в [38]–[41]. По определению (см., например, [35]–[37]), $l$-уровневый метод NIRK, примененный к ЗК (1.1), может быть записан в виде
(1.3a)
$x_{{kj}}^{2} = a_{{j1}}^{2}{{x}_{k}} + a_{{j2}}^{2}{{x}_{{k + 1}}} + {{\tau }_{k}}(d_{{j1}}^{2}g({{t}_{k}},{{x}_{k}}) + d_{{j2}}^{2}g({{t}_{{k + 1}}},{{x}_{{k + 1}}})),\quad j = 1,2,$(1.3б)
$\begin{gathered} x_{{kj}}^{i} = a_{{j1}}^{i}{{x}_{k}} + a_{{j2}}^{i}{{x}_{{k + 1}}} + {{\tau }_{k}}(d_{{j1}}^{i}g({{t}_{k}},{{x}_{k}}) + d_{{j2}}^{i}g({{t}_{{k + 1}}},{{x}_{{k + 1}}})) + \\ + \;{{\tau }_{k}}\sum\limits_{m = 1}^{i - 1} \,d_{{j,m + 2}}^{i}g(t_{{km}}^{{i - 1}},x_{{km}}^{{i - 1}}),\quad i = 3,4, \ldots ,l,\quad j = 1,2, \ldots ,i, \\ \end{gathered} $(1.3в)
${{x}_{{k + 1}}} = {{x}_{k}} + {{\tau }_{k}}\sum\limits_{i = 1}^l \,{{b}_{i}}g(t_{{ki}}^{l},x_{{ki}}^{l}),$Численные методы вида (1.3) мотивированы необходимостью дешевых и эффективных решателей, которые применимы к ЗК (1.1) различного рода. Существенно, что они наследуют главную особенность всех формул MIRK и не увеличивают размер решаемых ОДУ, т.е. нелинейные задачи, возникающие при дискретизации (1.3), всегда имеют ту же размерность $n$, что и исходная система ОДУ вне зависимости от порядка метода NIRK и количества его этапов (см., например, [33]–[37]).
Это особенно важно для эффективного интегрирования больших жестких задач, особенно тех, которые получены с помошью полудискретизации уравнений в частных производных (УрЧП) методом прямых, поскольку приведенные выше методы Гаусса и Лобатто, а также многие полностью неявные схемы РК расширяют размер обрабатываемых нелинейных уравнений до $l \times n$, где $l$ – число неявных стадий в используемой формуле РК (1.2).
Отметим, что другие эффективные неявные формулы РК, такие как однократно диагонально-неявные методы Рунге–Кутты (SDIRK) или однократно неявные схемы Рунге–Кутты (SIRK), имеют некоторые недостатки, ограничивающие область их применения (см., например, [35], [36]).
Кроме того, построенные формулы NIRK (1.3) типа Гаусса обладают многими привлекательными свойствами для практического использования. Например, это $A$-устойчивость, жесткая точность, симметричность и сопряженность симплектическими методам РК как минимум до порядка 6. Они также имеют достаточно высокие стадийные, псевдо-стадийные и классические порядки, а их некоторые стадийные величины, вычисленные на каждом шаге, достаточно точны и позволяют осуществлять плотную выдачу результатов интегрирования того же порядка, что и порядок сходимости основной формулы NIRK [35]–[37], [42], [43]. Именно поэтому методы NIRK (1.3) эффективны для решения различных ЗК вида (1.1), которые возникают в прикладных и инженерных исследования и охватывают нежесткие и жесткие ОДУ, а также обратимые и гамильтоновы задачи.
В [1]–[15] подробно изучено и описано применение методов NIRK для обработки математических моделей в гидромеханике и задачах нелинейной фильтрации Калмана.
В настоящей статье мы рассмотрим методы NIRK (1.3), удовлетворяющие условию
Потребуем также, чтобы стадийные величины $x_{{kj}}^{i}$ хорошо аппроксимировали точное решение $x(t)$ ЗК (1.1) в стадийных узлах $t_{{kj}}^{i}$, определяемых заданными коэффициентами В [35]–[37] доказано, что любой метод NIRK (1.3) с порядком сходимости $s \geqslant 1$ и удовлетворяющий условиям (1.4) может быть представлен в виде неявной формулы РК (1.2) с таблицей Бутчера видаВ настоящее время на практике принято применять адаптивные решатели для ОДУ, которые подразумевают, что методы РК (1.2), предназначенные для вычисления приближений $\{ {{x}_{k}}\} _{{k = 0}}^{K}$ к точному решению $x(t)$ ЗК (1.1), включают и механизм для автоматическего построения подходящей сетки
Это делает процедуру численного интегрирования более гибкой, эффективной и точной. Например, все встроенные в пакете MATLAB решатели ОДУ организованы именно таким образом, см. [44], [45].
С другой стороны, некоторое исключение из этого правила, т.е. когда методы с фиксированным шагом используются в реальных приложениях, также могут существовать [46]. Это касается случаев, когда гамильтоновы задачи интегрируются численно и где изменение размера шага может нарушать особые свойства геометрического численного интегрирования [27], [31]. В литературе также известно, что обсуждаемая реализация с переменным шагом интегрирования может создавать определенные неудобства для решения прикладных задач и требует непрерывного расширения реализуемой схемы РК для получения численных решений в заранее заданных точках на интервале интегрирования (см., например, [47]–[53]).
К счастью, последняя проблема решается автоматически во всех методах NIRK (1.3) из-за их специальной процедуры построения, которая позволяет применить дешевую процедуру плотной выдачи результатов численного интегрирования того же порядка аппоксимации, что и численное решение, полученное основной формулой NIRK.
Одной из важнейших характеристик любого приближенного решения является достигнутая точность численного интегрирования, т.е. величина ее глобальной ошибки $\{ \Delta {{x}_{k}}: = x({{t}_{k}}) - {{x}_{k}}\} _{{k = 0}}^{K}$ на автоматически генерируемой сетке $\{ {{t}_{k}}\} _{{k = 0}}^{K}$.
Однако большинство программных кодов для практического применения реализуют только тот или иной вариант управления локальной ошибкой численного интегрирования. Другими словами, программный код вычисляет оценку локальной погрешности схемы РК (или любой другой схемы) и сетка $\{ {{t}_{k}}\} _{{k = 0}}^{K}$ генерируется исходя их условия, что оценка локальной ошибки приближенного решения, вычисленного в узлах такой сетки, не может превышать заданной пользователем величины допуска Tol (см., например, описание часто используемых решателей ОДУ в пакете MATLAB [44], [45]).
Последний подход эффективно реализуется для методов РК (1.2) в виде их вложенных пар, т.е. приведенная выше таблица Бутчера дополнена коэффициентами ${{\tilde {b}}_{j}}$, $j = 1,2, \ldots ,l$, и принимает вид
(1.5)
$\begin{array}{*{20}{c}} c&\vline & A \\ \hline {}&\vline & {{{b}^{{\text{т}}}}} \\ \hline {}&\vline & {{{{\tilde {b}}}^{{\text{т}}}}} \end{array}\,{\kern 1pt} {\kern 1pt} .$Последняя строка таблицы (1.5) используется для вычисления второго (обычно более низкого порядка) численного решения
(1.6)
$x_{{k + 1}}^{{{\text{emb}}}} = {{x}_{k}} + {{\tau }_{k}}\sum\limits_{j = 1}^l \,{{\tilde {b}}_{j}}g({{t}_{k}} + {{c}_{j}}{{\tau }_{k}},{{x}_{{kj}}}),\quad k = 0,1, \ldots ,K - 1.$Тогда локальная ошибка оценивается с помощью простой формулы
на $(k + 1)$-м шаге вложенной пары РК. Далее, оценка локальной погрешности (1.7) масштабируется следующим образом:(1.8)
${\text{|}}l{{e}_{{k + 1}}}{{{\text{|}}}_{{sc}}}: = \mathop {max}\limits_{i = 1,2, \ldots ,n} \left\{ {{\text{|}}le_{{k + 1}}^{i}{\text{|/}}({\text{atol}} + \operatorname{rtol} \; \times \;{\text{|}}x_{{k + 1}}^{i}{\text{|}})} \right\},$Далее мы проверяем, что оценка локальной ошибки (1.8) удовлетворяет условию
Затем оптимальный размер $\tau _{k}^{ * }$ $(k + 1)$-го шага определяется по формуле(1.10)
$\tau _{k}^{ * }: = min\left\{ {1.5,0.8{\text{/}}\left| {l{{e}_{{k + 1}}}} \right|_{{sc}}^{{1/(p + 1)}}} \right\}{{\tau }_{k}},$Если условие (1.9) не выполняется, то локальная оценка погрешности (1.8) будет считаться слишком большой и, следовательно, численное решение ${{x}_{{к + 1}}}$, вычисленное на $(к + 1)$-м шаге с помощью встроенной пары РК (1.5), будет недостаточно точным и неприемлемым. Поэтому в этом случае все значения, вычисленные на текущем шаге, будут отклонены, размер шага ${{\tau }_{k}}$ будет уменьшен на основе формулы (1.10) и вычисления решения ${{x}_{{k + 1}}}$ и соответствующей оценки локальной ошибки ${\text{|}}l{{e}_{{k + 1}}}{{{\text{|}}}_{{sc}}}$ будут повторены той же парой (1.5), но для скорректированного размера шага ${{\tau }_{k}} = \tau _{k}^{*}$.
Если условие точности (1.9) выполняется для исходного (или скорректированного) размера шага ${{\tau }_{k}},\,$ то все найденные значения численного решения ${{x}_{{к + 1}}}$ и локальной ошибки $l{{e}_{{k + 1}}}$ будут приняты и $\tau _{k}^{*}$ будет рассматриваться как хороший прогноз следующего размера шага в адаптивной паре РК (1.5), т.е. ${{\tau }_{{k + 1}}} = \tau _{k}^{*}$ (см. подробное объяснение представленной процедуры управления размером шага в [25], [28], [29], [32]). Отметим, что формула (1.10) автоматически увеличит размер шага численного интегрирования, если текущий шаг принят и оценка ${{\left| {l{{e}_{{k + 1}}}} \right|}_{{sc}}}$ допущенной локальной ошибки много меньше 1.
Подчеркнем, что оценка локальной ошибки (1.6)–(1.8) является малозатратной, поскольку обе схемы РК во встроенной паре (1.5) используют одни и те же стадийные величины ${{x}_{{kj}}}$, $j = 1,2, \ldots ,l$. Более того, оптимистичный пользователь надеется, что контроль масштабированной локальной оценки погрешности (1.8) позволит получить решение $\{ {{x}_{k}}\} _{{k = 0}}^{K}$ с точностью, которая более или менее соответствует заданному допуску погрешности Tol. К сожалению, последнее не всегда имеет место для большинства численных методов и ОДУ. Основная причина этого заключается в том, что локальные и глобальные ошибки согласованных численных схем, которые включают все формулы РК, слабо связаны и могут отличаться значительно. Например, как показано в [54]–[57], отношение “глобальная ошибка”/“допуск” явного и неявного решателей ОДУ в пакете MATLAB может превышать ${{10}^{4}}$ даже для простых тестовых задач с достаточно высокими требованиями к точности вычислений.
Кроме того, в [58] дан пример жесткой задачи, для которого известный код RODAS дает для указанного отношения величину порядка ${{10}^{8}}$. В [56], [57] также показано, что эталонный решатель MATLAB для жестких задач ODE15s бесполезен для упомянутого примера. Таким образом, никакой стандартный локальный контроль точности не способен должным образом гарантировать глобальную погрешность используемых численных схем. Последнее привело к необходимости развития методов оценки и управления глобальной ошибкой в различных схемах численного интегрирования [59]–[65]. Современный прогресс в оценке глобальной ошибки в схемах РК в основном базируется на трех наиболее популярных методиках, т.е. на оценкe погрешности с помощью псевдозадачи Задунайского, на использовании специальных формул для вычисления оценки глобальной ошибки и при использовании двух численных формул разных порядков сходимости (подробнее см. [47]–[49], [55]–[57], [66]–[73]). Здесь мы также можем упомянуть оценку глобальной погрешности, основанную на глобальной экстраполяции, т.е. когда два независимых численных интегрирования выполняются одним и тем же адаптивным методом, но для разных значений локального допуска Tol. Однако последний метод подразумевает, что глобальная ошибка представима в виде асимптотического ряда по степеням параметра Tol, т.е. этот метод работает в РК-схемах со свойством пропорциональности допуска (см. точное определение этого свойства и его ограничения в [74]–[77]). Тем не менее все перечисленные методы оценки глобальной ошибки довольно дорогостоящи. Обычно необходимо повторное численное интегрирование ЗК (1.1) на всем интервале $[{{t}_{0}},{{t}_{{{\text{end}}}}}]$ для оценки глобальной ошибки, что, по крайней мере, вдвое увеличивает время выполнения таких расчетов с глобальным регулированием ошибки [63], [64]. Это являлось главным препятствием для широкого использования глобального контроля ошибки в используемых решателях ОДУ.
Недавно в [54] была представлена “дешевая” схема оценки глобальной ошибки, которая основана на простом суммировании оценок локальных ошибок (1.7). В этой работе также показано, что эта оценка ошибки может применяться в парах NIRK типа Гаусса (1.5). Конечно, представленный в [54] метод оценки глобальной ошибки не является асимптотически корректным, поскольку он только частично оценивает глобальную ошибку метода меньшего порядка в каждой встроенной паре, но он дешев и достаточно хорошо работает на практике. С другой стороны, численные тесты, проведенные с жесткими примерами в [54], показывают, что упомянутая глобальная оценка ошибки может сильно завышать истинную ошибку получаемого решения. Последнее обстоятельство значительно снижает эффективность опубликованных решателей NIRK для жестких ОДУ. Для решения этой проблемы представим пары NIRK типа Гаусса и Лобатто с автоматическим управлением масштабированными модифицированными локальной и глобальной ошибками. В разд. 4 будет исследована эффективность предложенных численных методов решения жестких ОДУ и показано их превосходство над жесткими решателями ОДУ в пакете MATLAB с автоматическим контролем локальных ошибок. Подчеркнем, что эффективность любого практического метода численного интегрирования зависит как от точности и устойчивости применяемых формул дискретизации, так и от качества реализованного способа оценки погрешности и генерации сетки. Вопросы, связанные с построением и изучением формул типа Гаусса и Лобатто для NIRK, уже рассматривались и были опубликованы в [35], [37]. В настоящей работе основное внимание уделяется автоматической генерации эффективных сеток, которые позволяют вычислять приближенные решения жестких ОДУ, соответствующие установленным пользователем требованиям к их точности в автоматическом режиме.
2. МЕТОДЫ NIRK ТИПА ГАУССА С АВТОМАТИЧЕСКИМ УПРАВЛЕНИЕМ ЛОКАЛЬНОЙ И ГЛОБАЛЬНОЙ ОШИБКАМИ
Здесь мы обсуждаем две встроенные пары методов NIRK типа Гаусса. Первая состоит из формул РК порядков 2 и 4, а вторая – из встроенных методов NIRK порядков 4 и 6. Оба численных метода используют комбинированный локальный и глобальный механизм контроля ошибки численного решения и предназначены для решения различных ОДУ с требуемой точностью.
2.1. Встроенная пара NIRK порядков 2 и 4
Настоящий адаптивный вычислительный метод разработан, исследован и протестирован численно в [33]–[35], [37], [54]. Его таблица Бутчера имеет вид
(2.1)
$\begin{array}{*{20}{c}} 0&\vline & 0&0&0&0 \\ {{{c}_{2}}}&\vline & {\frac{{6({{c}_{2}} + \theta ) - 5}}{{12}}}&{\frac{{1 - \theta }}{2}}&{\frac{{1 - \theta }}{2}}&{\frac{{6({{c}_{2}} + \theta ) - 7}}{{12}}} \\ {1 - {{c}_{2}}}&\vline & {\frac{{7 - 6({{c}_{2}} + \theta )}}{{12}}}&{\frac{\theta }{2}}&{\frac{\theta }{2}}&{\frac{{5 - 6({{c}_{2}} + \theta )}}{{12}}} \\ 1&\vline & 0&{\frac{1}{2}}&{\frac{1}{2}}&0 \\ \hline {}&\vline & 0&{\frac{1}{2}}&{\frac{1}{2}}&0 \\ \hline {}&\vline & {\frac{1}{2}}&0&0&{\frac{1}{2}} \end{array}.$Замечание 1. Параметр $\theta $ является свободным во встроенной паре NIRK (2.1) и может быть любым вещественным числом. Для любого $\theta $ метод более высокого порядка этой пары удовлетворяет упрощающим условиям $C(2)$ и $B(4)$ и имеет стадийный порядок 2 и классический порядок 4 (т.е. точно такие же, как и известный метод Гаусса порядка 4, представленный в [29, Table 5.1]). Кроме того, при $\theta = 1{\text{/}}2 + 2\sqrt 3 {\text{/}}9$ упомянутый метод NIRK дополнительно удовлетворяет упрощающему условию $C(3)$ и приобретает стадийный порядок 3 (см. доказательство в [33], [35], [37]). Напомним, что высокий стадийный порядок особенно важен для преодоления явления снижения порядка численного интегрирования, с которым можно столкнуться при решении очень жестких ОДУ [25], [29].
Замечание 2. Функция устойчивости формулы NIRK четвертого порядка в паре (2.1) является диагональной аппроксимацией Паде ${{R}_{{22}}}(z)$ независимо от $\theta $. Поэтому этот метод (с любым параметром $\theta $) эффективен для решения жестких ОДУ. Кроме того, он также подходит для решения очень жестких задач и дифференциально-алгебраических систем (см. дальнейшее объяснение в [25], [29]).
Замечание 3. Как следует из [28, Теорема 8.8], рассматриваемый метод является симметричным. Это свойство имеет решающее значение для его эффективности при численном интегрировании сингулярных краевых задач в [5]. Кроме того, это свойство делает указанный метод эффективным для решения обратимых ОДУ.
Замечание 4. Из [35] следует, что метод NIRK гауссовского типа и порядка 4 сопряжен с симплектическим методом, по крайней мере, до 6-го порядка включительно. Именно поэтому он эффективен для решения гамильтоновых задач на больших интервалах интегрирования, что подтверждено численным интегрированием задачи Кеплера (см. [34]).
Замечание 5. Встроенная пара NIRK (2.1) является неявной, поэтому ее практическая эффективность в значительной степени зависит от используемых итераций. Следующий упрощенный метод Ньютона рекомендован в [33], [37], [78]:
(2.2a)
$x_{{k1}}^{{2,\ell - 1}} = a_{{11}}^{2}{{\bar {x}}_{k}} + a_{{12}}^{2}x_{{k + 1}}^{{\ell - 1}} + {{\tau }_{k}}[d_{{11}}^{2}g({{t}_{k}},{{\bar {x}}_{k}}) + d_{{12}}^{2}g({{t}_{{k + 1}}},x_{{k + 1}}^{{\ell - 1}})],$(2.2б)
$x_{{k2}}^{{2,\ell - 1}} = a_{{21}}^{2}{{\bar {x}}_{k}} + a_{{22}}^{2}x_{{k + 1}}^{{\ell - 1}} + {{\tau }_{k}}[d_{{21}}^{2}g({{t}_{k}},{{\bar {x}}_{k}}) + d_{{22}}^{2}g({{t}_{{k + 1}}},x_{{k + 1}}^{{\ell - 1}})],$(2.2в)
${{Q}_{2}}({{\tau }_{k}}{{J}_{{k + 1}}})(x_{{k + 1}}^{\ell } - x_{{k + 1}}^{{\ell - 1}}) = - x_{{k + 1}}^{{\ell - 1}} + {{\bar {x}}_{k}} + {{\tau }_{k}}[{{b}_{1}}g(t_{{k1}}^{2},x_{{k1}}^{{2,\ell - 1}}) + {{b}_{2}}g(t_{{k2}}^{2},x_{{k2}}^{{2,\ell - 1}})],$Подчеркнем, что представленная схема Ньютона (2.2) использует тривиальный предиктор, означающий, что численное решение ${{\bar {x}}_{k}}$ с предыдущего шага метода NIRK используется в качестве его начального приближения. В этом случае ${{J}_{{k + 1}}}: = {{\partial }_{x}}g({{t}_{k}},{{\bar {x}}_{k}})$. Далее, в [78] доказано, что двух итераций в каждой точке сетки достаточно для обеспечения сходимости порядка 4 при ${{\tau }_{k}} \to 0$. Кроме того, векторы ${{\bar {x}}_{{k + 1}}}: = x_{{k + 1}}^{N}$, $x_{{k1}}^{{2,N}}$, $x_{{k2}}^{{2,N}}$ (заметим, что формулы (2.2а) и (2.2б) применяются еще один раз для вычисления стадийных величин, соответствующих искомому численному решению ${{\bar {x}}_{{k + 1}}}$) и вектор $\mathop {\bar {x}}\nolimits_k $, вычисленный на предыдущем шаге численного решения, могут быть использованы для построения полинома Эрмита третьей (или более высокой) степени, который можно применить для плотной выдачи результатов численного интегрирования с четвертым порядком точности. Кроме того, последний полином может быть использован для вычисления улучшенного начального приближения в итерациях (2.2). Тогда число итерационных шагов, требуемое для обеспечения сходимости четвертого порядка понижается до $N = 1$.
Также стоит отметить, что замена точного якобиана правой части формулы NIRK порядка 4 во встроенной паре (2.1) на ее аппроксимацию полиномом второй степени ${{Q}_{2}}({{\tau }_{k}}{{J}_{{k + 1}}})$ не влияет на устойчивость нашего метода и, следовательно, ее можно использовать для решения жестких ОДУ [33], [78]. Кроме того, эффективное решение линейной задачи (2.2в) осуществляется в виде последовательного решения двух линейных систем с одной и той же матрицей коэффициентов ${{I}_{n}} - {{\tau }_{k}}{{J}_{{k + 1}}}{\text{/}}4$.
Замечание 6. Схема оценки локальной ошибки (1.6), (1.7) для встроенной пары методов NIRK (2.1) принимает вид
(2.3)
$l{{e}_{{k + 1}}} = x_{{k + 1}}^{{{\text{emb}}}} - {{x}_{{k + 1}}} = \frac{{{{\tau }_{k}}}}{2}[g({{t}_{k}},{{\bar {x}}_{k}}) - g(t_{{k1}}^{2},x_{{k1}}^{{2,N}}) - g(t_{{k2}}^{2},x_{{k2}}^{{2,N}}) + g({{t}_{{k + 1}}},{{\bar {x}}_{{k + 1}}})],$(2.4)
${{R}_{{{\text{emb}}}}}(z) = \frac{{1 + 1{\text{/}}2z + 1{\text{/}}12{{z}^{2}} + 1{\text{/}}12{{z}^{3}}}}{{1 - 1{\text{/}}2z + 1{\text{/}}12{{z}^{2}}}},$Для устранения этого недостатка в цитируемых работах предлагается использовать модифицированную оценку локальной погрешности ${{\widetilde {le}}_{{k + 1}}}$, полученную в качестве решения линейной системы
где матрица коэффициентов вычисляется с помощью матричного полинома третьей степени ${{Q}_{3}}({{\tau }_{k}}{{J}_{{k + 1}}}): = {{({{I}_{n}} - {{\tau }_{k}}{{J}_{{k + 1}}}{\text{/}}4)}^{3}}$. Отметим, что эта модифицированная оценка локальной ошибки довольно дешева, потому что линейная задача (2.5) подразумевает последовательное решение трех линейных систем с одной и той же матрицей коэффициентов ${{I}_{n}} - {{\tau }_{k}}{{J}_{{k + 1}}}{\text{/}}4$. Кроме того, последняя матрица коэффициентов уже была вычислена и факторизована в итерациях (2.2). Положительным обстоятельством является то, что функция устойчивости, соответствующая модифицированной оценке локальной ошибки ${{\widetilde {le}}_{{k + 1}}}$, имеет вид(2.6)
${{R}_{{{\text{mod}}}}}(z) = \frac{{1 - 1{\text{/}}4z - 5{\text{/}}48{{z}^{2}} + 19{\text{/}}192{{z}^{3}} + 1{\text{/}}128{{z}^{4}} - 1{\text{/}}768{{z}^{5}}}}{{1 - 5{\text{/}}4z + 31{\text{/}}48{{z}^{2}} - 11{\text{/}}64{{z}^{3}} + 3{\text{/}}128{{z}^{4}} - 1{\text{/}}768{{z}^{5}}}},$Замечание 7. Оценка глобальной ошибки в паре NIRK (2.1) основана на простой идее из [54], т.е. истинная ошибка $\Delta {{x}_{{k + 1}}}$ встроенного решения $x_{{k + 1}}^{{{\text{emb}}}}$ в узле сетки ${{t}_{{k + 1}}}$ определяется по формуле
с начальной ошибкой интегрирования $\Delta {{x}_{0}}$, равной нулю, так как начальное значение решения ОДУ (1.1) известно точно. В формуле (2.7) учтен знак локальной ошибки, т.е. формулы для оценки локальной ошибки, использованные здесь и в [54], имеет разные знаки. Последний член в (2.7) находится по формуле (2.3). Затем вычисленные оценки локальной и глобальной погрешностей совместно регулируются с помощью достаточно сложного механизма управления размером шага численного интегрирования, описанным в [54, Алгоритм 3.2]. Цитируемая статья показывает разумную производительность пары NIRK (2.1) с комбинированным локальным и глобальным контролем точности, а также показывает значительную переоценку ошибки этого метода для жестких ОДУ.Причина завышения этой ошибки связана с формой функции устойчивости ${{R}_{{{\text{emb}}}}}(z)$ встроенного TR, заданной формулой (2.4). Кроме того, переоценка локальной ошибки, обсуждаемая в замечании 6, усиливается формулой суммирования (2.7). Вот почему глобальный контроль точности, разработанный в [54], неэффективен для решения жестких ОДУ. Чтобы устранить этот недостаток, необходимо снова обратиться к идее Shampine и заменить исходную оценку локальной ошибки $l{{e}_{{k + 1}}}$, использованную в формуле (2.7), на ее модифицированную версию ${{\widetilde {le}}_{{k + 1}}}$, полученную с помощью решения линейной задачи (2.5). Другими словами, теперь мы оцениваем модифицированную глобальную ошибку $\Delta \mathop {\tilde {x}}\nolimits_{k + 1} $ и масштабируем ее следующим образом:
(2.8)
$\Delta {{\tilde {x}}_{{k + 1}}} = \Delta {{\tilde {x}}_{k}} - {{\widetilde {le}}_{{k + 1}}},\quad {\text{|}}\Delta {{\tilde {x}}_{{k + 1}}}{{{\text{|}}}_{{sc}}}: = \mathop {max}\limits_{i = 1,2, \ldots ,n} \left\{ {{\text{|}}\Delta \tilde {x}_{{k + 1}}^{i}{\text{|/}}({\text{atol}} + {\text{rtol}}\; \times \;{\text{|}}\bar {x}_{{k + 1}}^{i}{\text{|}})} \right\},$Для обеспечения выполнения условия точности (2.9) на практике мы используем алгоритм автоматического управления размером шага численного интегрирования, представленный в [54, Алгоритм 3.2], но сформулированный теперь для масштабированных модифицированных оценок локальной и глобальной ошибок, вычисленных в замечаниях 6 и 7. Насколько нам известно, идея масштабированного глобального контроля точности численного решения является довольно новой и впервые обсуждается в явных и неявных равнозначных блочных схемах [55], [56]. Эффективность пары NIRK гауссовского типа порядков 2 и 4 с автоматическим управлением масштабированными модифицированными оценками локальной и глобальной ошибок будет продемонстрирована на жестких и крупномасштабных тестовых задачах и сравнена с жесткими встроенными решателями ОДУ пакета MATLAB в разд. 4.
2.2. Встроенная пара NIRK порядков 4 и 6
Этот адаптивный вычислительный метод разработан, изучен и протестирован численно в [35]–[37], [54]. Его таблица Бутчера имеет вид
(2.10)
$\begin{array}{*{20}{c}} 0&\vline & 0&0&0&0&0&0&0 \\ {{{c}_{2}}}&\vline & {\frac{{{{c}_{3}}}}{6}}&0&0&{\frac{{35 - 40{{c}_{3}}}}{{108}}}&{\frac{{56 - 64{{c}_{3}}}}{{108}}}&{\frac{{35 - 40{{c}_{3}}}}{{108}}}&{\frac{{{{c}_{3}} - 1}}{6}} \\ {{{c}_{3}}}&\vline & {\frac{{{{c}_{2}}}}{6}}&0&0&{\frac{{35 - 40{{c}_{2}}}}{{108}}}&{\frac{{56 - 64{{c}_{2}}}}{{108}}}&{\frac{{35 - 40{{c}_{2}}}}{{108}}}&{\frac{{{{c}_{2}} - 1}}{6}} \\ {{{c}_{4}}}&\vline & {\frac{{20{{c}_{6}} - 3}}{{200}}}&{\frac{{18{{c}_{3}} - 9}}{{100}} + \theta }&\theta &{\frac{{{{{v}}_{1}}}}{{360}} - \frac{{5\theta }}{9}}&{\frac{{{{{v}}_{1}}}}{{225}} - \frac{{8\theta }}{9}}&{\frac{{{{{v}}_{1}}}}{{360}} - \frac{{5\theta }}{9}}&{\frac{{3 - 20{{c}_{4}}}}{{200}}} \\ {{{c}_{5}}}&\vline & {\frac{1}{{32}}}&{\frac{{\sqrt {27} }}{{32}}}&{ - \frac{{\sqrt {27} }}{{32}}}&{\frac{5}{{36}}}&{\frac{2}{9}}&{\frac{5}{{36}}}&{ - \frac{1}{{32}}} \\ {{{c}_{6}}}&\vline & {\frac{{20{{c}_{4}} - 3}}{{200}}}&{ - \theta }&{\frac{{18{{c}_{2}} - 9}}{{100}} - \theta }&{\frac{{{{{v}}_{2}}}}{{360}} + \frac{{5\theta }}{9}}&{\frac{{{{{v}}_{2}}}}{{225}} + \frac{{8\theta }}{9}}&{\frac{{{{{v}}_{2}}}}{{360}} + \frac{{5\theta }}{9}}&{\frac{{3 - 20{{c}_{6}}}}{{200}}} \\ 1&\vline & 0&0&0&{\frac{5}{{18}}}&{\frac{4}{9}}&{\frac{5}{{18}}}&0 \\ \hline {}&\vline & 0&0&0&{\frac{5}{{18}}}&{\frac{4}{9}}&{\frac{5}{{18}}}&0 \\ \hline {}&\vline & {\frac{1}{6}}&0&0&0&{\frac{2}{3}}&0&{\frac{1}{6}} \end{array},$Ниже мы суммируем свойства метода (2.10) и особенности его практической реализации в виде следующих замечаний.
Замечание 8. Параметр $\theta $ является свободным во встроенной паре NIRK (2.10) и может быть любым вещественным числом. Для любого $\theta $ метод более высокого порядка этой пары удовлетворяет упрощающим условиям C(3) и B(6). В [35] также доказано, что он имеет стадийный порядок 3 и классический порядок 6 (т.е. точно такой же, как метод Гаусса порядка 6, представленный в [29, табл. 5.2 ]). Кроме того, при $\theta = (36{{c}_{6}} - 18{{c}_{3}} - 9){\text{/}}200$ упомянутый метод NIRK дополнительно удовлетворяет упрощающему условию $C(3,5)$, которое означает, что его стадийные величины третьего уровня удовлетворяют упрощающему условию $C(5)$ (см. подробности в [35]). Поэтому их можно использовать для плотной выдачи результатов численного интегрирования с шестым порядком точности.
Замечание 9. Функция устойчивости формулы NIRK шестого порядка в таблице (2.10) является диагональной аппроксимацией Паде ${{R}_{{33}}}(z)$ в независимости от $\theta $. Таким образом, этот метод NIRK (с любым $\theta $) эффективен для решения жестких ОДУ. Кроме того, он также жестко точен и, следовательно, подходит для решения очень жестких задач и дифференциально-алгебраических систем (см. дальнейшее объяснение в [25], [29]).
Замечание 10. Переставив строки в таблице (2.10) так, чтобы коэффициенты ${{c}_{i}}$ монотонно возрастали, из [28, Теорема 8.8] следует, что этот метод симметричен. Таким образом, он подходит для решения обратимых ЗК вида (1.1).
Замечание 11. В замечании 8 сказано, что формула NIRK порядка 6 во встроенной паре (2.10) по своим свойствам аналогична методу Гаусса порядка 6, но значительно “дешевле” последнего. Поэтому есть надежда, что она будет эффективна для решения гамильтоновых уравнений на больших интервалах интегрирования. Последнее было подтверждено численно интегрированием задачи Кеплера на интервале $[0,{{10}^{6}}]$ в [36].
Замечание 12. Встроенная пара NIRK (2.10) является неявной и, следовательно, ее эффективность в значительной степени зависит от используемой итерации. Следующий упрощенный метод Ньютона разработан в [36], [37] для реализации этого вычислительного инструмента:
(2.11а)
$x_{{kj}}^{{2,\ell - 1}} = a_{{j1}}^{2}{{\bar {x}}_{k}} + a_{{j2}}^{2}x_{{k + 1}}^{{\ell - 1}} + {{\tau }_{k}}[d_{{j1}}^{2}g({{t}_{k}},{{\bar {x}}_{k}}) + d_{{j2}}^{2}g({{t}_{{k + 1}}},x_{{k + 1}}^{{\ell - 1}})],\quad j = 1,2,$(2.11б)
$x_{{kj}}^{{3,\ell - 1}} = a_{{j1}}^{3}{{\bar {x}}_{k}} + a_{{j2}}^{3}x_{{k + 1}}^{{\ell - 1}} + {{\tau }_{k}}[d_{{j1}}^{3}g({{t}_{k}},{{\bar {x}}_{k}}) + d_{{j2}}^{3}g({{t}_{{k + 1}}},x_{{k + 1}}^{{\ell - 1}})] + {{\tau }_{k}}\sum\limits_{m = 1}^2 \,d_{{j,m + 2}}^{3}g(t_{{km}}^{2},x_{{km}}^{{2,\ell - 1}}),\quad j = 1,2,3,$(2.11в)
${{Q}_{3}}({{\tau }_{k}}{{J}_{{k + 1}}})(x_{{k + 1}}^{\ell } - x_{{k + 1}}^{{\ell - 1}}) = - x_{{k + 1}}^{{\ell - 1}} + {{\bar {x}}_{k}} + {{\tau }_{k}}\sum\limits_{i = 1}^3 \,{{b}_{i}}g(t_{{ki}}^{3},x_{{ki}}^{{3,\ell - 1}}),\quad \ell = 1,2, \ldots ,N,$Как и в приведенной выше встроенной паре NIRK порядков 2 и 4, упрощенная схема Ньютона (2.11) тоже применяет тривиальный предиктор, означающий, что численное решение $\mathop {\bar {x}}\nolimits_k $ с предыдущего шага метода NIRK используется в качестве его начального приближения. В этом случае ${{J}_{{k + 1}}}: = {{\partial }_{x}}g({{t}_{k}},\mathop {\bar {x}}\nolimits_k )$. Далее, в [36] доказано, что трех итераций в каждой точке сетки достаточно для обеспечения сходимости максимального шестого порядка при ${{\tau }_{k}} \to 0$. Кроме того, векторы ${{\bar {x}}_{{k + 1}}}: = x_{{k + 1}}^{N}$, $x_{{k1}}^{{3,N}}$, $x_{{k2}}^{{3,N}}$, $x_{{k3}}^{{3,N}}$ (заметим, что формулы (2.11а) и (2.11б) применяются еще один раз для вычисления стадийных величин, соответствующих искомому численному решению ${{\bar {x}}_{{k + 1}}}$) и вектор $\mathop {\bar {x}}\nolimits_k $, вычисленный на предыдущем шаге численного решения, могут быть использованы для построения полинома Эрмита пятой (или более высокой) степени, который можно применить для плотной выдачи результатов численного интегрирования с шестым порядком точности. Кроме того, последний полином может быть использован для вычисления улучшенного начального приближения в итерациях (2.11). Тогда число итерационных шагов, требуемое для обеспечения сходимости шестого порядка, понижается до $N = 1$ при ${{\tau }_{k}} \to 0$. Последнее вытекает из [78, Теорема 2.3].
Также стоит отметить, что замена точного якобиана правой части формулы NIRK шестого порядка во встроенной паре (2.10) на его аппроксимацию полиномом третьей степени ${{Q}_{3}}({{\tau }_{k}}{{J}_{{k + 1}}})$ приемлемо при решении жестких ОДУ [36]. В этом случае эффективное решение линейной задачи (2.11в) осуществляется в виде последовательного решения трех линейных систем с одной и той же матрицей коэффициентов ${{I}_{n}} - {{\tau }_{k}}{{J}_{{k + 1}}}{\text{/}}6$.
Замечание 13. Общая методика оценки локальной ошибки (1.6), (1.7), реализованная во встроенной паре (2.10), приобретает вид
(2.12)
$l{{e}_{{k + 1}}} = \frac{{{{\tau }_{k}}}}{3}\left[ {\frac{1}{2}g({{t}_{k}},{{{\bar {x}}}_{k}}) - \frac{5}{6}g(t_{{k1}}^{3},x_{{k1}}^{{3,N}}) + \frac{2}{3}g(t_{{k2}}^{3},x_{{k2}}^{{3,N}}) - \frac{5}{6}g(t_{{k3}}^{3},x_{{k1}}^{{3,N}}) + \frac{1}{2}g({{t}_{{k + 1}}},{{{\bar {x}}}_{{k + 1}}})} \right],$(2.13)
${{R}_{{{\text{emb}}}}}(z) = \frac{{1 + 1{\text{/}}2z + 1{\text{/}}10{{z}^{2}} + 1{\text{/}}120{{z}^{3}} + 1{\text{/}}2880{{z}^{5}}}}{{1 - 1{\text{/}}2z + 1{\text{/}}10{{z}^{2}} - 1{\text{/}}120{{z}^{3}}}},$Опять же, мы используем идею Shampine и вычисляем оценку модифицированной локальной ошибки ${{\widetilde {le}}_{{k + 1}}}$, решая линейную систему
матрица коэффициентов которой определяется матричным многочленом второй степени ${{Q}_{2}}({{\tau }_{k}}{{J}_{{k + 1}}}): = {{({{I}_{n}} - {{\tau }_{k}}{{J}_{{k + 1}}}{\text{/}}6)}^{2}}$. Подчеркнем, что модифицированная оценка локальной погрешности достаточно дешева, поскольку линейная задача (2.14) подразумевает последовательное решение двух линейных систем с одной и той же матрицей коэффициентов ${{I}_{n}} - {{\tau }_{k}}{{J}_{{k + 1}}}{\text{/6}}$. Кроме того, последняя матрица коэффициентов уже была вычислена и факторизована в итерациях (2.11). Теперь функция устойчивости, соответствующая модифицированной оценке локальной ошибки ${{\widetilde {le}}_{{k + 1}}}$, приобретает вид(2.15)
${{R}_{{{\text{mod}}}}}(z) = \frac{{1 + 1{\text{/}}6z - 7{\text{/}}180{{z}^{2}} - 1{\text{/}}90{{z}^{3}} + 1{\text{/}}1728{{z}^{5}}}}{{1 - 5{\text{/}}6z + 53{\text{/}}180{{z}^{2}} - 1{\text{/}}8{{z}^{3}} + 1{\text{/}}180{{z}^{4}} - 1{\text{/4}}320{{z}^{5}}}},$Из замечания 9 следует, что ${{R}_{{{\text{NIRK6}}}}}(z) \equiv {{R}_{{33}}}(z)$. Очевидно, что функция устойчивости (2.15) ограничена в левой комплексной полуплоскости и, следовательно, лучше соответствует управлению локальной ошибкой в жестких и очень жестких ОДУ (1.1). Затем управление локальной ошибкой осуществляется так же, как описано в разд. 1.
Замечание 14. Как уже говорилось в замечании 7, причина переоценки глобальной ошибки, наблюдаемая для встроенной пары (2.10) с управлением глобальной ошибкой на основе методики из [54], связана с формой функции устойчивости ${{R}_{{{\text{emb}}}}}(z)$ (2.13) встроенной квадратуры Лобатто порядка 4, используемой в формуле оценки погрешности (2.12). Кроме того, эта переоценка локальной ошибки усиливается формулой суммирования (2.7). Поэтому контроль глобальной ошибки, разработанный в [54], неэффективен для жестких ОДУ. Чтобы устранить этот недостаток, мы следуем идее Shampine и заменяем исходную оценку локальной ошибки $l{{e}_{{k + 1}}}$ в оценке глобальной погрешности (2.7) на значение модифицированной оценки локальной погрешности ${{\widetilde {le}}_{{k + 1}}}$, полученное с помощью решения линейной задачи (2.14). Другими словами, теперь мы оцениваем модифицированную глобальную ошибку по формуле
(2.16)
$\Delta {{\tilde {x}}_{{k + 1}}} = \Delta {{\tilde {x}}_{k}} - {{\widetilde {le}}_{{k + 1}}},$(2.17)
${\text{|}}\Delta {{\tilde {x}}_{{k + 1}}}{{{\text{|}}}_{{sc}}}: = \mathop {max}\limits_{i = 1,2, \ldots ,n} \left\{ {{\text{|}}\Delta \tilde {x}_{{k + 1}}^{i}{\text{|/}}\left( {{\text{atol}} + {\text{rtol}}\; \times \;{\text{|}}\bar {x}_{{k + 1}}^{i}{\text{|}}} \right)} \right\},$3. МЕТОДЫ NIRK ТИПА ЛОБАТТО С АВТОМАТИЧЕСКИМ УПРАВЛЕНИЕМ ЛОКАЛЬНОЙ И ГЛОБАЛЬНОЙ ОШИБКАМИ
Рассмотрим одну встроенную пару NIRK, основанную на двух формулах Лобатто порядков 2 и 4. Подробное описание этого инструмента численного интегрирования представлено в [37]. Приведенная пара NIRK типа Лобатто задается следующей таблицей Бутчера:
(3.1)
$\begin{array}{*{20}{c}} 0&\vline & 0&0&0 \\ {\frac{1}{2}}&\vline & {\frac{5}{{24}}}&{\frac{1}{3}}&{ - \frac{1}{{24}}} \\ 1&\vline & {\frac{1}{6}}&{\frac{2}{3}}&{\frac{1}{6}} \\ \hline {}&\vline & {\frac{1}{6}}&{\frac{2}{3}}&{\frac{1}{6}} \\ \hline {}&\vline & {\frac{1}{2}}&0&{\frac{1}{2}} \end{array}\;{\kern 1pt} .$Замечание 15. Сопоставив коффициенты таблицы метода (3.1) с коэффициентами метода Лобатто порядка 4, приведенными в [29, табл. 5.7 ], приходим к выводу, что формула NIRK четвертого порядка во встроенной паре (3.1) является формулой Лобатто IIIA порядка 4. Именно поэтому она удовлетворяет упрощающим условиям $C(3)$ и $B(4)$ и, следовательно, эта формула NIRK имеет стадийный порядок 3 и классический порядок 4.
Замечание 16. Хорошо известно, что функция устойчивости формулы РК более высокого порядка в паре NIRK (3.1) является диагональной аппроксимацией Паде ${{R}_{{22}}}(z)$ (см., например, [26, лемма 3.4.5]). Поэтому этот метод может быть использован для решения жестких ОДУ. Кроме того, он также жестко точен и поэтому подходит для решения очень жестких задач и дифференциально-алгебраических систем (см. дальнейшее объяснение в [25], [29]).
Замечание 17. Из [28, Теорема 8.8] следует, что рассматриваемый метод симметричен. Поэтому он может быть эффективен для решения обратимых ОДУ.
Замечание 18. В [35, Лемма 2] следует, что метод NIRK типа Лобатто и порядка 4 сопряжен с симплектическим методом, по крайней мере, до 6-го порядка включительно. Именно поэтому он подходит для решения гамильтоновых задач на больших интервалах интегрирования. Последнее было подтверждено численно интегрированием задачи Кеплера на интервале $[0,{{10}^{5}}]$ в [37].
Замечание 19. Встроенная пара NIRK (3.1) является неявной, поэтому ее практическая эффективность в значительной степени зависит от реализованной итерации. Кроме того, учитывая ее гнездовую структуру, для реализации численной формулы (3.1) в [37] рекомендован следующий упрощенный метод Ньютона:
(3.2a)
$x_{{k1}}^{{2,\ell - 1}} = a_{{11}}^{2}{{\bar {x}}_{k}} + a_{{12}}^{2}x_{{k + 1}}^{{\ell - 1}} + {{\tau }_{k}}[d_{{11}}^{2}g({{t}_{k}},{{\bar {x}}_{k}}) + d_{{12}}^{2}g({{t}_{{k + 1}}},x_{{k + 1}}^{{\ell - 1}})],$(3.2б)
${{Q}_{2}}({{\tau }_{k}}{{J}_{{k + 1}}})(x_{{k + 1}}^{\ell } - x_{{k + 1}}^{{\ell - 1}}) = - x_{{k + 1}}^{{\ell - 1}} + {{\bar {x}}_{k}} + {{\tau }_{k}}[{{b}_{1}}g({{t}_{k}},{{\bar {x}}_{k}}) + {{b}_{2}}g(t_{{k1}}^{2},x_{{k1}}^{{2,\ell - 1}}) + {{b}_{3}}g({{t}_{{k + 1}}},x_{{k + 1}}^{{\ell - 1}})],$Опять же, упрощенный метод Ньютона (3.2) применяет тривиальный предиктор, означающий, что численное решение ${{\bar {x}}_{k}}$ с предыдущего шага метода NIRK используется в качестве его начального приближения. В этом случае ${{J}_{{k + 1}}}: = {{\partial }_{x}}g({{t}_{k}},{{\bar {x}}_{k}})$. Кроме того, формулы NIRK типа Гаусса и Лобатто и порядка 4 используют одну и ту же аппроксимацию якобиана ${{Q}_{2}}({{\tau }_{k}}{{J}_{{k + 1}}})$. Поэтому замечание 5 доказывает, что двух итераций достаточно для обеспечения сходимости с максимальным порядком 4 при ${{\tau }_{k}} \to 0$. Далее, векторы ${{\bar {x}}_{{k + 1}}}: = x_{{k + 1}}^{N}$, $x_{{k1}}^{{2,N}}$ (напомним, что формула (3.2а) применяется еще один раз для вычисления стадийной величины, соответствующей искомому численному решению ${{\bar {x}}_{{k + 1}}}$) и вектор ${{\bar {x}}_{k}}$, вычисленный на предыдущем шаге численного решения, могут быть использованы для построения полинома Эрмита третьей (или более высокой) степени, который можно применить для плотной выдачи результатов численного интегрирования с четвертым порядком точности. Кроме того, последний полином может быть использован для вычисления улучшенного начального приближения в итерациях (3.2). Тогда число итерационных шагов, требуемое для обеспечения сходимости четвертого порядка, понижается до $N = 1$, что следует из [78, Теоремы 2.3].
Отметим, что замена точного якобиана правой части метода Лобатто IIIA на его аппроксимацию матричным полиномом второй степени ${{Q}_{2}}({{\tau }_{k}}{{J}_{{k + 1}}})$ подходит для решения жестких ОДУ, поскольку функция устойчивости итерации (3.2) с тривиальным предиктором и $N = 1$ задается формулой
где $z: = \lambda {{\tau }_{k}}$. Функция устойчивости (3.3) очевидно является $A$-устойчивой. Кроме того, напомним, что эффективное решение линейной задачи (3.2б) осуществляется в виде последовательного решения двух линейных систем с одной и той же матрицей коэффициентов ${{I}_{n}} - {{\tau }_{k}}{{J}_{{k + 1}}}{\text{/}}4$.Замечание 20. Общая схема оценки локальной ошибки (1.6), (1.7), реализованная во встроенной паре NIRK типа Лобатто (3.1), приводит к формуле оценки ошибки
(3.4)
$l{{e}_{{k + 1}}} = x_{{k + 1}}^{{{\text{emb}}}} - {{x}_{{k + 1}}} = \frac{{{{\tau }_{k}}}}{3}[g({{t}_{k}},{{\bar {x}}_{k}}) - 2g(t_{{k1}}^{2},x_{{k1}}^{{2,N}}) + g({{t}_{{k + 1}}},{{\bar {x}}_{{k + 1}}})],$(3.5)
${{R}_{{{\text{emb}}}}}(z) = \frac{{1 + 1{\text{/}}2z + 1{\text{/}}12{{z}^{2}} + 1{\text{/}}12{{z}^{3}}}}{{1 - 1{\text{/}}2z + 1{\text{/}}12{{z}^{2}}}}.$Очевидно, что функция устойчивости (3.5) неограничена в комплексной полуплоскости чисел $z$ с отрицательной вещественной частью. С другой стороны, эта функция устойчивости и функция устойчивости (2.4) встроенного TR в паре NIRK гауссовского типа (2.1) совпадают. Таким образом, модифицированная схема оценки локальной погрешности, представленная в замечании 6, также применима для преобразования оценки ошибки (3.4). Затем мы просто заменяем оценку локальной ошибки $l{{e}_{{k + 1}}}$ ее модифицированным значением $\mathop {\widetilde {le}}\nolimits_{k + 1} $, найденным путем решения линейной задачи (2.5). Замечание 6 объясняет, что это преобразование устраняет отмеченный недостаток и делает обсуждаемый контроль локальной ошибки подходящим для жестких ОДУ.
Замечание 21. Оценка и управление глобальной погрешностью во встроенной паре NIRK типа Лобатто (3.1) выполняются точно так же, как и в паре NIRK типа Гаусса (2.1) (см. замечание 7). Заметим, что все теоретические условия, лежащие в основе дешевого метода оценки глобальной ошибки, разработанного в [54], справедливы и для встроенной пары (3.1). Эффективность пары NIRK типа Лобатто порядков 2 и 4 с представленным выше управлением масштабированными модифицированными оценками локальной и глобальной ошибок будет продемонстрирована на жестких и крупномасштабных тестовых задачах в сравнении с жесткими встроенными решателями ОДУ пакета MATLAB в следующем разделе.
4. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ
Выполним ряд численных тестов для изучения трех пар NIRK, представленных в разд. 2 и 3, на жестких и крупномасштабных задачах. Сначала мы проверим их способность находить численные решения c заданной точностью. Требуемая точность вычислений в наших экспериментах определяется допуском на погрешность Tol. Другими словами, мы хотим, чтобы масштабированная глобальная погрешность найденного численного решения не превышала эту границу. Такая информация представлена на приведенных ниже фигурах для точности интегрирования. Для оценки эффективности наших адаптивных решателей ОДУ мы также изучим точность, достигнутую каждым методом, в зависимости от времени его работы. На практике нас интересуют инструменты численного интегрирования, которые обеспечивают более высокую точность интегрирования за более короткое время работы. Такие данные представлены на фигурах для эффективности численного решения. Отметим, что наши результаты (т.е. допущенные ошибки численного интегрирования) показаны в сравнении с используемыми допусками Tol и затраченным временем процессора в двойном логарифмическом масштабе.
На представленных фигурах использованы следующие сокращения для NIRK: пара NIRK типа Гаусса (2.1) с автоматическим управлением масштабированными модифицированными оценками локальной и глобальной ошибок отмечена как NIRK4(2)g; пара NIRK типа Гаусса (2.10) с автоматическим управлением масштабированными модифицированными оценками локальной и глобальной ошибок показана аббревиатурой NIRK6(4)g; пара NIRK типа Лобатто (3.1) с автоматическим управлением масштабированными модифицированными оценками локальной и глобальной ошибок указана как NIRK4(2)l. Мы сравним также точность и эффективность наших решателей NIRK с жесткими встроенными кодами ODE15s, ODE23s и ODE23tb пакета MATLAB. Для корректности сравнения все вычисления выполнены в MATLAB.
Дополнительно следует иметь в виду, что для повышения эффективности разработанных кодов NIRK все упрощенные итерации Ньютона, реализованные в решателях NIRK, используют интерполяционный полином второй степени, построенный по значениям численного решения, найденного в двух предыдущих узлах сетки. Этот полином применяется для более точного предсказания начального приближения. Поэтому первые два шага в каждом коде NIRK выполняются с тривиальным предиктором, но с увеличенным числом допустимых итераций. Кроме того, каждый упрощенный метод Ньютона использует переменное число итераций в узлах сетки, начиная с третьего шага в каждом решателе NIRK. Другими словами, мы должны решить две задачи. Первая из них – убедиться, что численные решения ${{\bar {x}}_{k}}$, полученные парами NIRK типа Гаусса и Лобатто (2.1) и (3.1), сходятся с порядком 4 (или численное решение ${{\bar {x}}_{k}}$ из пары NIRK типа Гаусса (2.10) сходится c порядком 6) при ${{\tau }_{k}} \to 0$. Вторая задача состоит в том, чтобы представить некоторые аргументы в пользу того, что итерационные ошибки, возникающие в примененных итерациях, незначительны по сравнению с ошибками дискретизации, имеющими место в использованных схемах NIRK. Для достижения обеих целей, начиная с третьего шага в каждом коде NIRK выполняется следующий расчет. Мы выполняем одну итерацию в формулах NIRK четвертого порядка (или три итерации в методе шестого порядка) без контроля точности. Это должно гарантировать, что найденное численное решение сходится с предопределенным теоретическим порядком. Затем мы начинаем мониторить невязку следующей итерации и так далее. Таким образом, вычислив очередную итерацию $x_{{k + 1}}^{\ell }$, мы проверяем ее точность с помощью критерия остановки
(4.1)
$\mathop {\max }\limits_{i = 1, \ldots ,n} \{ {\text{|}}{{(x_{{k + 1}}^{\ell })}^{i}} - {{(x_{{k + 1}}^{{\ell - 1}})}^{i}}{\text{|/}}(1\; + \;{\text{|}}{{(x_{{k + 1}}^{\ell })}^{i}}{\text{|}})\} \leqslant {{\varepsilon }_{l}}{\text{/}}10,$Наш первый жесткий пример довольно прост и взят из [46]. Он сформулирован в виде следующего ОДУ.
Задача I:
(4.3)
${\text{|}}\Delta x(t){{{\text{|}}}_{{sc}}}: = \mathop {\max }\limits_{k = 1,2, \ldots ,K} \mathop {\max }\limits_{i = 1,2, \ldots ,n} \{ {\text{|}}x_{{{\text{exact}}}}^{i}({{t}_{k}}) - x_{k}^{i}{\text{|/}}(1\; + \;{\text{|}}x_{{{\text{exact}}}}^{i}({{t}_{k}}){\text{|}})\} ,$В первом эксперименте мы установили допуски ${\text{Tol}}: = {{10}^{{ - 1}}},{{10}^{{ - 2}}}, \ldots ,{{10}^{{ - 10}}}$ на ошибку численного интегрирования во всех рассматриваемых кодах. Подчеркнем, что параметр точности ${{\varepsilon }_{g}}$, используемый при выборе размера шага в [54, Алгоритм 3.2], который реализован во всех парах NIRK, здесь обозначается как Tol. Решатели ОДУ ODE15s, ODE23s и ODE23tb в пакете MATLAB применяются с опциями AbsTol := Tol, RelTol := Tol и MaxStep := 10-1. Такое же ограничение на максимальный размер шага интегрирования установлено в наших парах NIRK с автоматическим управлением масштабированными модифицированными оценками локальной и глобальной ошибок.
Задача I решена численно как кодами NIRK NIRK4(2)g, NIRK6(4)g и NIRK4(2)l, так и с помощью встроенных в MATLAB решателей ODE15s, ODE23s и ODE23tb с вышеуказанными ограничениями на ошибку. Затем допущенная максимальная масштабированная глобальная ошибка вычислена по формуле (4.3) для каждого исследуемого решателя, и эти ошибки показаны в сравнении с используемыми допусками Tol и процессорным временем на фиг. 1. На левом графике представлено исследование точности этих решателей ОДУ, а на правой диаграмме показана их эффективность.
Фиг. 1.
Диаграммы точности и эффективности для жестких решателей ОДУ с локальным контролем ошибки из пакета MATLAB и для трех пар NIRK с масштабированным модифицированным локальным и глобальным контролем ошибки, протестированных на задаче I.
На фиг. 1а показано, что точность, достигнутая парами NIRK, лучше согласуется с нашими требованиями, представленными на каждом графике точности в виде жирной прямой линии, соответствующей условию “Масштабированная Глобальная Ошибка” = “Допуск” особенно при высоких требованиях к точности вычислений. В этой статье мы хотим получить численные решения, масштабированные глобальные ошибки которых (4.3) лежат ниже этой линии. В целом решатели NIRK удовлетворяют наш запрос в отличие от встроенных кодов пакета MATLAB. Это неудивительно, потому что локальный контроль ошибки, реализованный в ODE15s, ODE23s и ODE23tb, не может должным образом регулировать их глобальные ошибки. Наиболее точное численное интегрирование среди рассмотренных кодов MATLAB выполняется с помощью ODE23s. С другой стороны, из диаграммы эффективности фиг. 1 следует, что это наименее эффективный ОДУ-решатель для задачи I. Среди кодов NIRK наиболее эффективными являются NIRK4(2)g и NIRK4(2)l, которые на удивление так же эффективны, как ODE15s (этот решатель ОДУ переменного порядка в пакете MATLAB многими практиками признается эталоном для интегрирования жестких задач) и даже превосходят пару NIRK NIRK6(4)g более высокого порядка. Последнее объясняется простотой задачи I для высокоточного численного интегрирования и тем, что метод шестого порядка выполняет более крупные шаги по сравнению с методами четвертого порядка. Это требует большего числа повторных циклов интегрирования задачи I, начиная с начального узла сетки ${{t}_{0}}$, кодом NIRK6(4)g. Чтобы увидеть преимущества нашего кода NIRK более высокого порядка, нам потребуются более сложные тестовые примеры, рассмотренные ниже.
Следующей жесткой системой ОДУ, использованной в этой статье, является известный осциллятор Ван-дер-Поля [29, с. 144]. Он индуцирует устойчивые колебания и может быть сформулирован в виде следующего ОДУ.
Задача II:
(4.4)
${\text{|}}\Delta x({{t}_{{{\text{end}}}}}){{{\text{|}}}_{{{\text{sc}}}}}: = \mathop {\max }\limits_{i = 1,2, \ldots ,n} \{ {\text{|}}x_{{{\text{ref}}}}^{i}({{t}_{{{\text{end}}}}}) - x_{{{\text{end}}}}^{i}{\text{|/}}(1\; + \;{\text{|}}x_{{{\text{ref}}}}^{i}({{t}_{{{\text{end}}}}}){\text{|}})\} ,$Фиг. 2.
Диаграммы точности и эффективности для жестких решателей ОДУ с локальным контролем ошибки из пакета MATLAB и для трех пар NIRK с масштабированным модифицированным локальным и глобальным контролем ошибки, протестированных на задаче II.
Из диаграммы точности фиг. 2 следует, что все коды MATLAB не способны вычислить более или менее точные численные решения для осциллятора Ван-дер-Поля. Более того, мы не видим даже реальной сходимости этих кодов до тех пор, пока требования к точности численного решения не становится очень высокими, т.е. когда ${\text{Tol}} \leqslant {{10}^{{ - 8}}}$. Кроме того, среди решателей MATLAB для наиболее точного численного решения, полученного решателем ODE15s для ${\text{Tol}}: = {{10}^{{ - 10}}}$, соотношение “Истинная ошибка”/“Допуск на ошибку” составляет около ${{10}^{7}}$, что вряд ли приемлемо на практике. Напротив, все коды NIRK с масштабированным модифицированным локальным и глобальным контролем ошибки дают численные решения, соответствующие заданным требованиям к точности вычислений (см. фиг. 2а). Можно также наблюдать некоторую переоценку ошибки в решателях NIRK. Это результат разности на два порядка во встроенных методах NIRK, т.е. глобальная ошибка оценивается для метода, который является на два порядка менее точным, чем метод, используемый для вычисления итогового численного решения. Однако эта переоценка не приводит к существенному снижению эффективности кодов NIRK для задачи II. Последнее следует из диаграммы эффективности на фиг. 2. Теперь мы видим, что лучшим методом интегрирования осциллятора Ван-дер-Поля среди рассмотренных в данной работе является NIRK6(4)g, который обеспечивает требуемую точность численного решения за минимальное процессорное время.
На следующем жестком примере мы покажем, что разработанная стратегия оценки глобальной ошибки дешева, но не асимптотически корректна (см. объяснение в [54]). Это означает, что существуют жесткие ЗК, для которых коды NIRK могут автоматически не обеспечить требуемую точность численного интегрирования. Для этого рассмотрим наш наиболее сложный тестовый пример в [58], который представлен в виде следующего ОДУ.
Задача III:
(4.5)
${{x}_{1}}(t) = {{(t + 1)}^{2}},\quad {{x}_{2}}(t) = t + 1,\quad {{x}_{3}}(t) = exp\{ - 25{{(t - 1)}^{2}}\} .$Фиг. 3.
Диаграммы точности и эффективности для жестких решателей ОДУ с локальным контролем ошибки из пакета MATLAB и для трех пар NIRK с масштабированным модифицированным локальным и глобальным контролем ошибки, протестированных на задаче III.
График точности фиг. 3 показывает, что все рассматриваемые коды не в состоянии обеспечить требуемую точность численного интегрирования задачи III, за исключением NIRK6(4)g. Кроме того, решатели MATLAB ODE15s и ODE23tb практически бесполезны для этого тестового примера, поскольку их точность не зависит от уровня допуска на локальную ошибку, установленного в третьем вычислительном эксперименте, т.е. они находят численные решения с аналогичной ошибкой для всех значений порога ошибки Tol, примененных здесь. Остальные три кода ODE23s, NIRK4(2)g и NIRK4(2)l показывают более или менее одинаковые и лучшие результаты. Таким образом, диаграммы точности и эффективности на фиг. 3 показывают, что наиболее точным и эффективным кодом для решения задачи III снова является NIRK6(4)g.
Наконец, мы протeстируем наши адаптивные коды NIRK4(2)g, NIRK4(2)l и NIRK6(4)g в экстремальных условиях интегрирования крупномасштабных задач, возникающих при полудискретизации уравнений в частных производных (УрЧП). Наглядным примером для этого является известный двумерный брюсселатор с диффузией и периодическими граничными условиями [29, с. 151, 152]. Следуя цитированной книге, мы приходим к следующей задаче.
Задача IV:
Полудискретизация задачи IV означает, что мы берем сетку из 50 точек по каждой переменной. Пространственные производные второго порядка в правой части этого УрЧП заменяются центральными разделенными разностями. Это приводит к ОДУ размера 5000. В конечном итоге мы получаем достаточно большую тестовую задачу, удобную для исследования наших решателей NIRK в условиях, приближенных к решению реальных практических задач. Опять же, задача IV не имеет точного решения в замкнутом виде. Таким образом, мы используем эталонное решение, вычисленное с высокой точностью в конечной точке интервала интегрирования $[0,6]$ для оценки полученной масштабированной глобальной ошибки по формуле (4.4) с $n = 5000$. Требования к точности всех кодов NIRK определяются параметром ${\text{Tol}}: = {{10}^{{ - 1}}},{{10}^{{ - 2}}}, \ldots ,{{10}^{{ - 6}}}$, тогда как решатели MATLAB тестируются с теми же самыми допусками, что и в приведенных выше тестовых примерах. Глобальные ошибки всех решателей ОДУ, наблюдаемые в последнем эксперименте, показаны на фиг. 4.
Фиг. 4.
Диаграммы точности и эффективности для жестких решателей ОДУ с локальным контролем ошибки из пакета MATLAB и для трех пар NIRK с масштабированным модифицированным локальным и глобальным контролем ошибки, протестированных на задаче IV.
Прежде всего отметим, что задача IV достаточно проста и не содержит сильно нелинейных членов. Поэтому даже жесткие решатели MATLAB с только локальным контролем ошибок способны обеспечить приемлемые численные решения. Тем не менее все эти коды не в состоянии обеспечить требуемую точность численного интегрирования и занижают истинную погрешность примерно на два порядка при высоких требованиях к точности численного интегрирования. Напротив, NIRK4(2)g, NIRK4(2)l и NIRK6(4)g работают точно и возвращают численные решения, удовлетворяющие наперед заданным условиям точности, в автоматическом режиме (см. левый график фиг. 4). Однако такие достаточно малые масштабированные глобальные ошибки достигаются ценой повторных численных решений задачи IV, начиная с начальной точки интервала интегрирования. Эти дополнительные численные интегрирования требуют увеличения времени счета и, следовательно, снижают их эффективность, представленную на правой диаграмме фиг. 4. Последняя фигура также показывает, что наш код NIRK6(4)g превосходит оставшиеся два NIRK4(2)g и NIRK4(2)l по эффективности численного интегрирования задачи IV. Во-первых, его механизм оценки глобальной погрешности лучше оценивает истинную ошибку и приводит к ее меньшей переоценке, наблюдаемой по сравнению с другими решателями NIRK на диаграмме точности фиг. 4. Во-вторых, эта переоценка ошибки кодов NIRK4(2)g и NIRK4(2)l является причиной их меньшей эффективности, которая отмечена на фиг. 4б. Производительность методов NIRK четвертого порядка почти сравнима с производительностью наименее эффективного кода MATLAB ODE23s.
В конечном счете представленные численные тесты подтверждают успешность обновленных пар NIRK с автоматическим управлением масштабированными модифицированными оценками локальной и глобальной ошибок для численного интегрирования сложных жестких и крупномасштабных ЗК в практических расчетах, а также их превосходство над жесткими встроенными интеграторами MATLAB при решении ОДУ с наперед заданной точностью в автоматическом режиме.
В итоге адаптивный решатель NIRK6(4)g, по-видимому, является наиболее предпочтительным методом для решения жестких моделей ОДУ, возникающих часто в прикладной науке и технике. Первое успешное применение рассмотренных методов NIRK для обработки детерминистких и стохастических математических моделей в гидромеханике и нелинейной фильтрации Калмана описано в работах [1]–[15].
5. ЗАКЛЮЧЕНИЕ
В данной работе обсуждаются вычислительные методы для решения задач Коши вида (1.1), которые являются традиционным инструментом для описания математических динамических моделей с непрерывным временем в различных исследованиях. Точность их численного решения может иметь решающее значение для адекватности и согласованности используемого метода математического моделирования. В центре нашего исследования находятся адаптивные методы NIRK с автоматическим контролем глобальной ошибки, разработанные в [54]. К сожалению, приведенная схема оценки глобальной ошибки обнаруживает серьезную переоценку ошибки для жестких примеров, рассмотренных в [54], что снижает ее практическую эффективность. Именно поэтому в данной работе был исследован источник этой переоценки ошибки и найден способ решения указанной проблемы. Это привело к парам NIRK типа Гаусса и Лобатто с автоматическим управлением масштабированными модифицированными оценками локальной и глобальной ошибок, которые гораздо более эффективны для решения жестких ЗК, чем адаптивные решатели ОДУ, представленные в [54]. Кроме того, сложные тесты, выполненные в разд. 4, не только подтвердили высокую практическую эффективность наших вычислительных методов, но и выявили их превосходство над жесткими встроенными кодами ODE15s, ODE23s и ODE23tb пакета MATLAB для большинства используемых тестовых примеров. Все это создает твердую основу для дальнейшего изучения новых адаптивных методов NIRK в практических вычислениях. Например, решатели NIRK6(4)g, NIRK4(2)g и NIRK4(2)l, разработанные в этой статье, использованы для повышения эффективности непрерывно-дискретных обощенного, кубатурного и сигма-точечного методов фильтрации Калмана, применяемых для оценки состояния стохастических систем с непрерывным временем в моделях радиолокационного слежения как с гауссовским, так и с негауссовским шумом [3], [6], [7], [11], [14], [15], а также в электротехнических и химических инженерных исследованиях [2], [4], [8]–[10], [12], [13].
Список литературы
Kulikov G.Yu., Kulikova M.V. Accurate numerical implementation of the continuous-discrete extended Kalman filter // IEEE Trans. Automat. Contr. 2014. V. 59 (1). P. 273–279.
Kulikov G.Yu., Kulikova M.V. High-order accurate continuous-discrete extended Kalman filter for chemical engineering // Eur. J. Contr. 2015. V. 21. P. 14–26.
Kulikov G.Yu., Kulikova M.V. The accurate continuous-discrete extended Kalman filter for radar tracking // IEEE Trans. Signal Process. 2016. V. 64 (4). P. 948–958.
Kulikov G.Yu., Kulikova M.V. Estimating the state in stiff continuous-time stochastic systems within extended Kalman filtering // SIAM J. Sci. Comput. 2016. V. 38 (6). P. A3565–A3588.
Kulikov G.Yu., Lima P.M., Morgado M.L. Analysis and numerical approximation of singular boundary value problems with the $p$-Laplacian in fluid mechanics // J. Comput. Appl. Math. 2014. V. 262. P. 87–104.
Kulikov G.Yu., Kulikova M.V. Accurate cubature and extended Kalman filtering methods for estimating continuous-time nonlinear stochastic systems with discrete measurements // Appl. Numer. Math. 2017. V. 111. P. 260–275.
Kulikov G.Yu., Kulikova M.V. Accurate continuous-discrete unscented Kalman filtering for estimation of nonlinear continuous-time stochastic models in radar tracking // Signal Process. 2017. V. 139. P. 25–35.
Kulikov G.Yu., Kulikova M.V. Accurate state estimation of stiff continuous-time stochastic models in chemical and other engineering // Math. Comput. Simulation. 2017. V. 142. P. 62–81.
Kulikov G.Yu., Kulikova M.V. Square-root Kalman-like filters for estimation of stiff continuous-time stochastic systems with ill-conditioned measurements // IET Control Theory Appl. 2017. V. 11 (9). P. 1420–1425.
Kulikov G.Yu., Kulikova M.V. Moore-Penrose-pseudo-inverse-based Kalman-like filtering methods for estimation of stiff continuous-discrete stochastic systems with ill-conditioned measurements // IET Control Theory Appl. 2018. V. 12 (16). P. 2205–2212.
Kulikov G.Yu., Kulikova M.V. Estimation of maneuvering target in the presence of non-Gaussian noise: A coordinated turn case study // Signal Process. 2018. V. 145. P. 241–257.
Kulikov G.Yu., Kulikova M.V. Practical implementation of extended Kalman filtering in chemical systems with sparse measurements // Russian J. Numer. Anal. Math. Modelling. 2018. V. 33 (1). P. 41–53.
Kulikov G.Yu., Kulikova M.V. Numerical robustness of extended Kalman filtering based state estimation in ill-conditioned continuous-discrete nonlinear stochastic chemical systems // Int. J. Robust Nonlinear Control. 2019. V. 29 (5). P. 1377–1395.
Kulikov G.Yu., Kulikova M.V. Square-root accurate continuous-discrete extended-unscented Kalman filtering methods with embedded orthogonal and $J$-orthogonal $QR$ decompositions for estimation of nonlinear continuous-time stochastic models in radar tracking // Signal Process. 2020. V. 166. P. 107253.
Kulikov G.Yu., Kulikova M.V. NIRK-based Cholesky-factorized square-root accurate continuous-discrete unscented Kalman filters for state estimation in nonlinear continuous-time stochastic models with discrete measurements // Appl. Numer. Math. 2020. V. 147. P. 196–221.
Aggoun L., Elliot R.J. Measure Theory and Filtering: Introduction and Applications. Cambridge University Press, Cambridge, U.K., 2005.
Aström K.J. Introduction to Stochastic Control Theory. Academic Press, New York, 1970.
Bar-Shalom Y., Li X.-R., Kirubarajan T. Estimation with Applications to Tracking and Navigation. Wiley, New York, 2001.
Crassidis J.L., Junkins J.L. Optimal Estimation of Dynamic Systems, CRC Press LLC, New York, 2004.
Grewal M.S., Andrews A.P. Kalman Filtering: Theory and Practice. Prentice Hall, New Jersey, 2001.
Grewal M.S., Weill L.R., Andrews A.P. Global Positioning Systems. Inertional Navigation and Integration, Wiley, New York, 2001.
Jazwinski A.H. Stochastic Processes and Filtering Theory. Academic Press, New York, 1970.
Øksendal B. Stochastic Differential Equations: An Introduction with Applications. Springer, New York, 2003.
Rawlings J.B., Mayne D.Q. Model Predictive Control: Theory and Design. Bob Hill Publishing, LLC, Madison, Wisconsin, 2013.
Butcher J.C. Numerical Methods for Ordinary Differential Equations. John Wiley and Sons, Chichester, 2008.
Dekker K., Verwer M.P. Stability of Runge–Kutta Methods for Stiff Nonlinear Differential Equations. North-Holland, Amsterdam, 1984.
Hairer E., Lubich C., Wanner G. Geometric Numerical Integration: Structure Preserving Algorithms for Ordinary Differential Equations. Springer-Verlag, Berlin, 2002.
Hairer E., Nørsett S.P., Wanner G. Solving Ordinary Differential Equations I: Nonstiff Problems. Springer-Verlag, Berlin, 1993.
Hairer E., Wanner G. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Springer-Verlag, Berlin, 1996.
Jackiewicz Z. General Linear Methods for Ordinary Differential Equations. John Wiley and Sons, Hoboken, 2009.
Sanz-Serna J.M., Calvo M.P. Numerical Hamilton Problems. Chapman and Hall, London, 1994.
Shampine L.F. Numerical Solution of Ordinary Differential Equations. Chapman and Hall, New York, 1994.
Kulikov G.Yu., Shindin S.K. On a family of cheap symmetric one-step methods of order four, in: V.N. Alexandrov, et al. (Eds.), Computational Science – ICCS 2006, Proceedings, Part I, 6th International Conference, Reading, UK, May 28–31, 2006, Lecture Notes in Computer Science, vol. 3991, Springer-Verlag, Berlin, 2006. P. 781–785.
Kulikov G.Yu., Shindin S.K. Numerical tests with Gauss-type nested implicit Runge–Kutta formulas, in: Y. Shi, et al. (Eds.), Computational Science – ICCS 2007, Proceedings, Part I, 7th International Conference, Beijing, China, May 27–30, 2007, Lecture Notes in Computer Science, vol. 4487, Springer-Verlag, Berlin, 2007. P. 136–143.
Kulikov G.Yu., Shindin S.K. Adaptive nested implicit Runge–Kutta formulas of auss type // Appl. Numer. Math. 2009. V. 59 (3–4). P. 707–722.
Kulikov G.Yu. Automatic error control in the Gauss-type nested implicit Runge–Kutta formula of order 6 // Russian J. Numer. Anal. Math. Modelling. 2009. V. 24 (2). P. 123–144.
Kulikov G.Yu. Embedded symmetric nested implicit Runge–Kutta methods of Gauss and Lobatto types for solving stiff ordinary differential equations and Hamiltonian systems // Comput. Math. Math. Phys. 2015. V. 55 (6). P. 983–1003.
Cash J.R. A class of implicit Runge–Kutta methods for the numerical solution of stiff ordinary differential equations // J. ACM. 1975. V. 22 (4). P. 504–511.
Cash J.R. On a class of implicit Runge–Kutta procedures // IMA J. Numer. Anal. 1977. V. 19 (4). P. 455–470.
Cash J.R. On a note of the computational aspects of a class of implicit Runge–Kutta procedures // IMA J. Numer. Anal. 1977. V. 20 (4). P. 425–441.
Cash J.R., Singhal A. Mono-implicit Runge–Kutta formulae for the numerical integration of stiff differential systems // IMA J. Numer. Anal. 1982. V. 2. P. 211–227.
Skvortsov L.M. How to avoid accuracy and order reduction in Runge–Kutta methods as applied to stiff problems // Comput. Math. Math. Phys. 2017. V. 57 (7). P. 1124–1139.
Skvortsov L.M. Implicit Runge–Kutta methods with explicit internal stages // Comput. Math. Math. Phys. 2018. V. 58 (3). P. 307–321.
Shampine L.F., Reichelt M.W. The MATLAB ODE suite // SIAM J. Sci. Comput. 1997. V. 18. P. 1–22.
Higham D., Higham N. MATLAB Guide. SIAM, Philadelphia, 2005.
Kulikov G.Yu., Weiner R. Doubly quasi-consistent parallel explicit peer methods with built-in global error estimation // J. Comput. Appl. Math. 2010. V. 233 (9). P. 2351–2364.
Baker T.S., Dormand J.R., Gilmore J.P., Prince P.J. Continuous approximation with embedded Runge–Kutta methods // Appl. Numer. Math. 1996. V. 22. P. 51–62.
Dormand J.R., Prince P.J. Runge–Kutta triples // Comput. Math. Appl. 1986. V. 12A. P. 1007–1017.
Dormand J.R., Prince P.J. Practical Runge–Kutta processes // SIAM J. Sci. Statist. Comput. 1989. V. 10. P. 977–989.
Enright W.H. Analysis of error control strategies for continuous Runge–Kutta methods // SIAM J. Numer. Anal. 1989. V. 26. P. 588–599.
Enright W.H., Jackson K.R., Nørsett S.P., Thomsen P.G. Interpolants for Runge–Kutta formulas // ACM Trans. Math. Software. 1986. V. 12. P. 193–218.
Higham D.J. Highly continuous Runge–Kutta interpolants // ACM Trans. Math. Software. 1991. V. 17. P. 368–386.
Sharp P.W., Verner J.H. Generation of high-order interpolants for explicit Runge–Kutta pairs // ACM Trans. Math. Software. 1998. V. 24. P. 13–29.
Kulikov G.Yu. Cheap global error estimation in some Runge–Kutta pairs // IMA J. Numer. Anal. 2013. V. 33 (1). P. 136–163.
Weiner R., Kulikov G.Yu. Local and global error estimation and control within explicit two-step peer triples // J. Comput. Appl. Math. 2014. V. 262. P. 261–270.
Kulikov G.Yu., Weiner R. A singly diagonally implicit two-step peer triple with global error control for stiff ordinary differential equations // SIAM J. Sci. Comput. 2015. V. 37 (3). P. A1593–A1613.
Weiner R., Kulikov G.Yu., Beck S., Bruder J. New third- and fourth-order singly diagonally implicit two-step peer triples with local and global error controls for solving stiff ordinary differential equations // J. Comput. Appl. Math. 2017. V. 316. P. 380–391.
Kulikov G.Yu., Weiner R. Global error estimation and control in linearly-implicit parallel two-step peer-methods // J. Comput. Appl. Math. 2011. V. 236 (6). P. 1226–1239.
Aïd R., Levacher L. Numerical investigations on global error estimation for ordinary differential equations // J. Comput. Appl. Math. 1997. V. 82. P. 21–39.
Shampine L.F. Global error estimation for stiff ODEs // Lecture Notes in Mathematics. 1984. V. 1066. P. 159–168.
Shampine L.F. Error estimation and control for ODEs // J. Sci. Comput. 2005. V. 25. P. 3–16.
Shampine L.F., Watts H.A. Global error estimation for ordinary differential equations // ACM Trans. Math. Software. 1976. V. 2. P. 172–186.
Skeel R.D. Thirteen ways to estimate global error // Numer. Math. 1986. V. 48. P. 1–20.
Skeel R.D. Global error estimation and the backward differentiation formulas // Appl. Math. Comput. 1989. V. 31. P. 197–208.
Stetter H.J. Global error estimation in ODE-solvers, in: J. Hinze (Ed.), Numerical integration of differential equations and large linear systems, Proceedings, Bielefeld 1980. Lecture Notes in Mathematics, vol. 968, Springer-Verlag, Berlin, 1982. P. 269–279.
Dormand J.R., Duckers R.R., Prince P.J. Global error estimation with Runge–Kutta methods // IMA J. Numer. Anal. 1984. V. 4. P. 169–184.
Dormand J.R., Gilmore J.P., Prince P.J. Globally embedded Runge–Kutta schemes // Ann. Numer. Math. 1994. V. 1. P. 97–106.
Dormand J.R., Lockyer M.A., McGorrigan N.E., Prince P.J. Global error estimation with Runge–Kutta triples // Comput. Math. Appl. 1989. V. 18. P. 835–846.
Lang J., Verwer J.G. On global error estimation and control for initial value problems // SIAM J. Sci. Comput. 2007. V. 29. P. 1460–1475.
Macdougall T., Verner J.H. Global error estimators for order 7, 8 Runge–Kutta pairs // Numer. Algorithms. 2002. V. 31. P. 215–231.
Makazaga J., Murua A. New Runge–Kutta based schemes for ODEs with cheap global error estimation // BIT. 2003. V. 43. P. 595–610.
Shampine L.F., Baca L.S. Global error estimates for ODEs based on extrapolation methods // SIAM J. Sci. Stat. Comput. 1985. V. 6. P. 1–14.
Tirani R. A parallel algorithm for the estimation of the global error in Runge–Kutta methods // Numer. Algorithms. 2002. V. 31. P. 311–318.
Calvo M., Higham D.J., Montijano J.I., Rández L. Stepsize selection for tolerance proportionality in explicit Runge–Kutta codes // Adv. Comput. Math. 1997. V. 7. P. 361–382.
Calvo M., González-Pinto S., Montijano J.I. Global error estimation based on the tolerance proportionality for some adaptive Runge–Kutta codes // J. Comput. Appl. Math. 2008. V. 218. P. 329–341.
Higham D.J. Global error versus tolerance for explicit Runge–Kutta methods // IMA J. Numer. Anal. 1991. V. 11. P. 457–480.
Higham D.J. The tolerance proportionality of adaptive ODE solvers // J. Comput. Appl. Math. 1993. V. 45. P. 227–236.
Kulikov G.Yu., Merkulov A.I., Shindin S.K. Asymptotic error estimate for general Newton-type methods and its application to differential equations // Russian J. Numer. Anal. Math. Modelling. 2007. V. 22 (6). P. 567–590.
Kulikov G.Yu., Shindin S.K. One-leg variable-coefficient formulas for ordinary differential equations and local-global step size control // Numer. Algorithms. 2006. V. 43 (1). P. 99–121.
Дополнительные материалы отсутствуют.
Инструменты
Журнал вычислительной математики и математической физики
Пн.-пт.: 10.00-19.00