Журнал вычислительной математики и математической физики, 2020, T. 60, № 7, стр. 1111-1125

1. ВВЕДЕНИЕ

Ограничение на размер шага не позволяет использовать традиционные явные методы Рунге–Кутты и Адамса для эффективного решения жестких задач. Поэтому для решения таких задач обычно применяют неявные методы. В то же время неявные методы требуют вычисления матрицы Якоби и решения системы нелинейных алгебраических уравнений, что усложняет их реализацию и затрудняет использование при наличии сложных нелинейных зависимостей. Поэтому наряду с неявными методами разрабатывают и применяют специальные явные методы, позволяющие эффективно решать жесткие задачи.

Класс жестких задач настолько разнообразен, что существует множество определений таких задач, но среди них нет общепринятого. Авторы известной книги по численному решению жестких задач [1] Э. Хайрер и Г. Ваннер считают, что наиболее практичным определением понятия “жесткий” является самое раннее, данное в 1952 г. Кертиссем и Хиршфельдером [2]: “Жесткие уравнения – это уравнения, для которых определенные неявные методы, в частности ФДН, дают лучший результат, обычно несравненно более хороший, чем явные методы”. В настоящее время известны специальные явные методы, эффективные для многих жестких задач, поэтому под явными методами в этом определении следует понимать классические явные методы Рунге–Кутты и Адамса. Отметим, что термины “жесткие уравнения” или “жесткая задача” применимы к конкретной задаче Коши, т.е. при заданных начальных условиях и интервале интегрирования, поскольку задача, жесткая на большом интервале, может быть нежесткой на меньшем интервале или при других начальных условиях.

Будем решать задачу Коши для системы ОДУ

где t – независимая переменная, y – вектор переменных состояния, ${\mathbf{f}}(t,{\mathbf{y}})$ – нелинейная векторная функция. Предположим, что задача (1.1) имеет единственное решение ${\mathbf{y}}(t)$, а функция f – гладкая в любой точке решения на интервале интегрирования. Тогда на всем интервале определена и непрерывна матрица Якоби системы (1.1) Требование гладкости правой части не всегда согласуется с реальными задачами, правые части которых могут содержать разрывные характеристики типа реле или переключателей. При наличии таких элементов будем предполагать, что число переключений конечно, а весь интервал интегрирования можно разбить на несколько интервалов, на каждом из которых функция f остается гладкой. Тогда решения “склеиваются”, т.е. на каждом последующем интервале в качестве начального условия принимается решение, полученное в конце предыдущего интервала.

Для дальнейшего изложения нам нужно определить меру жесткости задачи Коши (1.1). В качестве такой меры можно использовать минимальное число шагов, необходимое для устойчивого решения задачи простейшим явным методом. Рассмотрим сначала линейную задачу

где матрица J имеет собственные значения ${{\lambda }_{1}},\; \ldots ,\;{{\lambda }_{n}}$ с отрицательными действительными частями. Для такой задачи число шагов пропорционально интервалу интегрирования T и обратно пропорционально минимальной постоянной времени ${{\tau }_{{\min }}}$. Поэтому для линейных задач вычислительные затраты можно оценить числом Значение ${{M}_{{{\text{stf}}}}}$, округленное до большего целого, – минимальное число шагов явного метода Эйлера, обеспечивающее устойчивое решение всех задач вида (1.2), удовлетворяющих условиям $\operatorname{Re} \left( {{{\lambda }_{i}}} \right) < 0$, $\left| {\operatorname{Im} \left( {{{\lambda }_{i}}} \right)} \right| < \left| {\operatorname{Re} \left( {{{\lambda }_{i}}} \right)} \right|$.

В случае нелинейной задачи (1.1) спектр матрицы Якоби зависит от t, поэтому вычислительные затраты явных методов можно оценить с помощью интеграла

Величину ${{M}_{{{\text{stf}}}}}$ назовем мерой жесткости задачи Коши. В формуле (1.3) оцениваются вычислительные затраты, вызванные только жесткостью задачи. Реальное число шагов может превышать эти оценки вследствие негладкости правой части, наличия больших собственных значений вблизи мнимой оси или в правой полуплоскости (колебательные и плохо обусловленные задачи), а также по другим причинам.

Аналогично можно определить меры колебательности и неустойчивости задачи Коши. Величину

назовем мерой колебательности, а величину назовем мерой неустойчивости задачи Коши. Значения этих мер, а также число переменных n и величина интервала интегрирования T для восьми жестких тестов из [1] приведены в табл. 1. Отметим, что задача BEAM имеет чисто мнимый спектр матрицы Якоби. Поэтому формально ее можно отнести к колебательным задачам, но она проявляет свойство жесткости, поскольку эффективно может быть решена только неявными методами. Этот факт еще раз подтверждает трудность формального определения класса жестких задач.

Наряду с неявными методами, для решения жестких задач с вещественным жестким спектром успешно применяют специальные явные методы. Принципы построения таких методов рассмотрим на примере явного двухстадийного метода Рунге–Кутты 1-го порядка с функцией устойчивости $R(z) = 1 + z + d{{z}^{2}}$. На фиг. 1 приведены области устойчивости такого метода при различных значениях d. При $d > 1{\text{/}}8$ область устойчивости односвязная, а ее длина равна $1{\text{/}}d$. Уменьшая d, можно увеличить длину области устойчивости. При $d < 1{\text{/}}8$ область устойчивости перестает быть односвязной и состоит из двух разделенных областей, наиболее удаленная из которых позволяет обеспечить стабилизацию метода в жесткой части спектра.

Приведенные области позволяют сформулировать два способа построения явных методов для жестких задач. Первый способ основан на максимальном расширении области устойчивости. Построение многочленов устойчивости одношаговых методов такого типа обычно выполняют исходя из условия чебышёвского альтернанса, поэтому их часто называют методами Рунге–Кутты–Чебышёва. Методы с расширенными областями устойчивости рассматривались в [1], [3]–[7] и многих других работах.

Идея второго способа заключается в получении на основе предварительных стадий оценок наибольших по модулю собственных значений матрицы Якоби, которые используются в заключительной формуле для стабилизации расчетной схемы в полученных точках жесткого спектра. Как отмечалось в [7]–[9], эта идея (обычно в замаскированном виде) использовалась во многих явных нелинейных методах (например, в методах из [10]–[14]). В [10] впервые были четко выделены этапы получения оценок собственных значений и последующего использования этих оценок для стабилизации численного решения. Несмотря на большое число работ, посвященных явным нелинейным методам, они так и не получили широкого практического применения. В [14] на конкретных примерах было показано, что предложенные ранее нелинейные методы пригодны для решения с низкой точностью лишь весьма узкого класса жестких задач. Основной причиной низкой эффективности этих методов является снижение реального порядка при решении жестких задач. Например, двухстадийный метод имеет 2-й порядок при решении нежестких задач, но при решении жестких задач реальный порядок снижается до нулевого (см. [15], [16]). Под реальным порядком мы понимаем порядок, который оценивается при реальных размерах шага h, позволяющих обеспечить заданную точность, тогда как классический порядок оценивается при $h \to 0$.

Более эффективные явные нелинейные одношаговые и многошаговые методы, названные адаптивными, были предложены в [15], [16]. Рассмотренные в этих работах одношаговые методы позволяют обеспечить 1-й порядок для жестких задач при числе стадий $s = 3$ и 2-й порядок при $s > 3$. Общий принцип построения адаптивных одношаговых методов изложен в [7]–[9], а в [8] предложен метод, имеющий для жестких задач 3-й порядок. Столь же эффективные для жестких задач одношаговые методы более высоких порядков построить не удалось, но аналогичные многошаговые методы могут иметь практически любой порядок, и на их основе в [17] построен явный адаптивный метод переменного порядка. Явные адаптивные методы реализованы в программных комплексах МВТУ [18], [19] и SimInTech [20], а также в виде программ в системе MathCad, которые доступны в дополнительных материалах к книге [7] по адресу: http://3v-services.com/books/978-5-97060-636-0. В настоящей статье исследуются явные адаптивные одношаговые методы, которые заметно проще многошаговых методов и в то же время позволяют эффективно решать многие жесткие задачи при умеренных требованиях к точности.

Оценим вычислительные затраты явных методов в зависимости от жесткости при решении задач с вещественным жестким спектром. Определим вычислительные затраты как число вычислений правой части на всем интервале Nf. Тогда Nf = N × s , где N – число выполненных шагов, s – число стадий на одном шаге. Для методов с односвязной областью устойчивости (к ним относятся классические методы и методы с расширенными областями), исходя из условия устойчивости численного решения, получаем

где $l(s)$ – длина области устойчивости вдоль вещественной оси. При числе стадий, равном порядку метода, для метода 2-го порядка имеем $c = 1$, для метода 3-го порядка $c = 1.19$ и для метода 4-го порядка $c = 1.44$. Таким образом, вычислительные затраты классических методов пропорциональны жесткости задачи ${{M}_{{{\text{stf}}}}}$.

Оценим теперь затраты чебышёвских методов. Для них размер шага h не ограничен условием устойчивости, а выбирается исходя из заданной точности, и может быть таким же большим, как в неявных методах. Но при этом увеличивается число стадий s, которое может достигать нескольких сотен и определяется из условия устойчивости $l(s) > h\mathop {\max }\limits_i \left( { - {{\lambda }_{i}}} \right)$, где $l(s) \approx K{{s}^{2}}$ (для методов 2-го порядка $K \approx 0.8$, а для методов 3-го порядка $K \approx 0.48$). Тогда оценку затрат получаем в виде

Если задача не очень жесткая, то число шагов N зависит от характера решения и почти не зависит от жесткости, поэтому вычислительные затраты пропорциональны $\sqrt {{{M}_{{{\text{stf}}}}}} $.

Методы с расширенными областями устойчивости применяют для решения задач с распределенным вещественным спектром матрицы Якоби, к которым относятся, прежде всего, задачи, полученные в результате дискретизации уравнений в частных производных параболического типа (например, задача BRUSS из [1]). В таких задачах жесткие составляющие, соответствующие различным собственным значениям, сильно связаны между собой. Другой тип задач с вещественным спектром – когда жесткие составляющие слабо связаны между собой. Это проявляется в том, что для получения с необходимой точностью оценок всех жестких собственных значений матрицы Якоби достаточно выполнить несколько итераций покомпонентного степенного метода. К таким задачам относятся задачи VDPOL, ROBER, OREGO, HIRES и CUSP из табл. 1. Именно при решении таких задач явные адаптивные методы наиболее эффективны и при расчетах с двойной точностью позволяют эффективно решать задачи с жесткостью до ${{M}_{{{\text{stf}}}}} = {{10}^{{155}}}$ (см. результаты решения задачи Капса в [7], [9]). Никакие другие явные методы не могут решать такие задачи.

2. ПОСТРОЕНИЕ ЯВНЫХ АДАПТИВНЫХ МЕТОДОВ

В общем случае явные адаптивные одношаговые методы строятся на основе явных стадий Рунге–Кутты, выполняемых по формулам

Пусть ${{a}_{{i,i - 1}}} \ne 0$, $i = 2,...,s$. Определим векторы ${{{\mathbf{u}}}_{1}},...,{{{\mathbf{u}}}_{s}}$ в виде где $[{{\beta }_{{ij}}}] = [{{{\mathbf{A}}}^{0}}{\mathbf{e}},{{{\mathbf{A}}}^{1}}{\mathbf{e}}, \ldots ,{{{\mathbf{A}}}^{{s - 1}}}{\mathbf{e}}]$, ${\mathbf{A}} = [{{a}_{{ij}}}]$, ${\mathbf{e}} = {{[1, \ldots ,1]}^{{\text{T}}}}$. Тогда при решении линейной системы (1.2) будут выполняться соотношения которые позволяют воспользоваться степенным методом для получения вектора покомпонентных оценок больших по модулю собственных значений матрицы $h{\mathbf{J}}$ в виде В этой и последующих формулах предполагаем покомпонентное выполнение операций с векторами.

Шаг интегрирования выполняем согласно формуле

где d_s – вектор настраиваемых параметров. Учитывая (2.2), эту же формулу можно записать в виде где ${{{\mathbf{d}}}_{{s - 1}}} = {{\mathbf{e}} \mathord{\left/ {\vphantom {{\mathbf{e}} {\left( {s - 1} \right)\,!}}} \right. \kern-0em} {\left( {s - 1} \right)\,!}} + {{{\mathbf{d}}}_{s}}{\mathbf{\tilde {z}}}$. В дальнейшем будем использовать формулу (2.4), которая содержит меньше членов и менее чувствительна к вычислительным ошибкам. Вектор ${{{\mathbf{d}}}_{{s - 1}}}$ задаем из условия, чтобы численное решение уравнения $y{\text{'}} = \lambda y$ выражалось формулой ${{y}_{1}} = Q\left( z \right){{y}_{0}}$, где $z = h\lambda $, $Q(z)$ – скалярная функция устойчивости. Тогда для вычисления очередной компоненты вектора ${{{\mathbf{d}}}_{{s - 1}}}$ (обозначим ее ${{d}_{{s - 1}}}$) можно использовать формулы

В случае линейной системы (1.2) получаем

где есть матричная функция устойчивости, ${{{\mathbf{D}}}_{{s - 1}}} = {\text{diag}}\left( {{{{\mathbf{d}}}_{{s - 1}}}} \right)$.

Скалярную функцию устойчивости $Q\left( z \right)$ задаем из условий заданного порядка сходимости нежестких компонент, затухания жестких компонент, расхождения неустойчивых компонент и ограниченности коэффициента ${{d}_{{s - 1}}}$ (см. [7]–[9]). Анализ численных экспериментов показал, что для больших по модулю значений z порядок согласованности с экспонентой неважен. Главное – обеспечить качественно правильное поведение (быстрое затухание либо расхождение) соответствующих компонент. Например, при $s = 3$ можно задать

а при $s = 4$ можно задать

Вычисление коэффициентов вектора ${{{\mathbf{d}}}_{{s - 1}}}$ организуем таким образом, чтобы избежать операций с большими числами и деления на 0. Рассмотрим, например, вычисление ${{d}_{2}}$ при $s = 3$ и $Q(z)$ в виде (2.6). Если очередная компонента нежесткая, т.е. если $\left| {{{u}_{3}}} \right| \leqslant 1.6\left| {{{u}_{2}}} \right|$, то вычисляем $\tilde {z} = {{{{u}_{3}}} \mathord{\left/ {\vphantom {{{{u}_{3}}} {{{u}_{2}}}}} \right. \kern-0em} {{{u}_{2}}}}$, ${{d}_{2}} = {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2} + {{\tilde {z}} \mathord{\left/ {\vphantom {{\tilde {z}} 6}} \right. \kern-0em} 6}$ (при ${{u}_{2}} = 0$ принимаем $\tilde {z} = 0$), а иначе вычисляем

Согласно (2.1), первая стадия следующего шага выполняется по формуле

В [8], [9] предложена модификация адаптивных методов, при которой первая стадия на всех шагах, кроме первого, выполняется согласно формуле Из (2.2), (2.5) следует, что формулу (2.9) можно записать также в виде Использование формулы (2.10) вместо (2.8) позволяет сэкономить на каждом шаге одно вычисление правой части, а также уменьшает влияние на устойчивость неточности оценок собственных чисел. Поэтому стадия в виде (2.10) была названа стабилизированной первой стадией.

Другой способ повышения устойчивости заключается в выполнении дополнительной коррекции жестких компонент. Например, при $s = 3$ коррекцию жесткой компоненты можно выполнять по формуле

(индекс компоненты опускаем). В результате появляется дополнительный корень многочлена устойчивости, близкий к оценке $\tilde {z}$ для этой компоненты, что приводит к расширению области устойчивости в окрестности $\tilde {z}$. Отметим, что коррекцию (2.11) можно применять итерационно, а если применить ее к методу Эйлера ${{{\mathbf{y}}}_{1}} = {{{\mathbf{y}}}_{0}} + h{{{\mathbf{f}}}_{0}}$, то получим формулу адаптивного метода (2.4).

Таким образом, можно выделить три типа адаптивных методов: обычный, со стабилизированной 1-й стадией, с коррекцией. Рассмотрим особенности каждого из этих типов, их преимущества и недостатки.

3. УСТОЙЧИВОСТЬ АДАПТИВНЫХ МЕТОДОВ

В [15], [16] были предложены явные адаптивные методы, стадии которых выполняются согласно формулам

Исходя из классической теории, оптимальное значение β в этих методах равно ${2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-0em} 3}$. В этом случае при $s \geqslant 3$ метод может иметь 3-й порядок. Однако анализ решения простейших модельных уравнений показал, что для жестких задач реальный порядок методов, основанных на стадиях вида (3.1), не превышает $s - 2$ и не может быть выше 2-го, а оптимальное значение β равно 1 (см. [7], [8]). Функция устойчивости i-й стадии выражается формулой Чтобы ограничить рост внутренних функций устойчивости, следует задать достаточно малое значение α. Мы рекомендуем задавать этот параметр в виде где $\bar {\alpha } = {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}$ или $\bar {\alpha } = {1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-0em} 3}$, ${{h}_{{{\text{old}}}}}$ – размер предыдущего шага, ${{\tilde {z}}_{i}}$ – покомпонентные оценки вида (2.2), полученные на предыдущем шаге.

Исследуем устойчивость адаптивных методов, основанных на стадиях (3.1) при $s = 3$, в зависимости от неточности вычисления оценок ${{\tilde {z}}_{i}}$. Решаем уравнение Далквиста $y{\text{'}} = \lambda y$. Пусть $\bar {z} = h\lambda $ и оценка этой величины получена с относительной ошибкой ε, тогда $\tilde {z} = \bar {z}(1 + \varepsilon )$. Предположим также, что

Для обычного метода, исходя из условия $Q(\tilde {z}) = 0$, получаем ${{d}_{2}} = - {{\tilde {z}}^{{ - 1}}} - {{\tilde {z}}^{{ - 2}}}$. Функция устойчивости такого метода при $z = \bar {z}$ имеет вид и при выполнении (3.3) получаем условие устойчивости в виде ограничения на значения $\bar {z}$ в виде $\left| {\bar {z}} \right| < {{\left| \varepsilon \right|}^{{ - 1}}}$. Отметим, что если вместо (2.4) использовать формулу (2.3), то интервал допустимых при данной ошибке значений $\bar {z}$ заметно сузится: $\left| {\bar {z}} \right| < \sqrt 2 {{\left| \varepsilon \right|}^{{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}$.

Рассмотрим теперь метод с дополнительной коррекцией. В этом случае функция устойчивости скорректированного метода будет имееть вид

откуда, учитывая (3.3), получаем ${{R}_{{{\text{cor}}}}}(\bar {z}) \approx {{\varepsilon }^{2}}\bar {z}$ и условие устойчивости $\left| {\bar {z}} \right| < {{\left| \varepsilon \right|}^{{ - 2}}}$.

Для метода со стабилизированной первой стадией формулы одного шага решения уравнения Далквиста можно записать в виде

Эти формулы задают разностную схему с характеристическим многочленом При выполнении (3.3) оба корня этого многочлена по модулю меньше 1, что обеспечивает устойчивость разностной схемы. Таким образом, неточность оценки собственного значения практически не накладывает ограничения на значение $\bar {z}$. Отметим, что если вместо (2.10) использовать эквивалентную при точном вычислении формулу (2.9), то условие устойчивости будет таким же, как у обычного метода.

Аналогичным образом получены ограничения вида $\left| {\bar {z}} \right| < {{\left| {\bar {z}} \right|}_{{\max }}}$ для методов, основанных на четырех стадиях вида (3.1), и для методов, основанных на шести стадиях, приведенных в [7], [8]. Значения ${{\left| {\bar {z}} \right|}_{{\max }}}$ приведены в табл. 2, где s – число стадий, на основе которых построен адаптивный метод, $\tilde {p}$ – реальный порядок метода при решении жестких задач, nf – число вычислений правой части на одном шаге. Обозначения методов приняты в виде ARK (Adaptive Runge–Kutta) – обычный метод, ARKc – с дополнительной коррекцией, ARKs – со стабилизированной 1-й стадией.

Приведем пример, подтверждающий полученные оценки. Будем решать задачу

с точным решением ${{y}_{1}}(t) = \sin t$, ${{y}_{2}}(t) = \cos t$. Используем методы со стадиями (3.1) при $s = 3$, $\beta = 1$. При точном выполнении всех вычислений получим точные оценки собственного значения. Внесем ошибки в эти оценки в виде Ошибку решения вычисляем по формуле Зависимости ошибки от жесткости задачи μ для методов рассмотренных трех типов при $h = {\pi \mathord{\left/ {\vphantom {\pi {30}}} \right. \kern-0em} {30}}$ и $\varepsilon = {{10}^{{ - 3}}}$ приведены на фиг. 2. Здесь и далее в обозначении методов первая цифра – классический порядок, а вторая – реальный порядок для жестких задач. Если эти порядки совпадают, то оставляем только одну цифру.

Как и следовало ожидать из проведенного выше анализа, расхождение решения методом ARK21 наблюдается при $h\mu > {{\varepsilon }^{{ - 1}}}$, а методом ARK21с – при $h\mu > {{\varepsilon }^{{ - 2}}}$. Внутри интервалов устойчивости ошибки этих двух методов уменьшаются при увеличении μ, что объясняется видом уравнений и свойством жесткой точности при $\beta = 1$ (см. анализ ошибки решения уравнения Протеро–Робинсона в [7], [8]). Метод ARK21s, в отличие от ARK21, не расходится при $h\mu > {{\varepsilon }^{{ - 1}}}$, но ошибка перестает уменьшаться с ростом μ и сохраняет свое значение вплоть до $\mu = {{10}^{{155}}}$, и только при $\mu = 1.4 \times {{10}^{{155}}}$ выдается сообщение о переполнении.

Рассмотрим теперь методы со стадиями (3.1) при $s = 4$, $\beta = 1$, которые сохраняют 2-й порядок при решении жестких задач и обозначены через ARK2, ARK2c и ARK2s. Результаты решения задачи (3.4) этими методами при $\varepsilon = {{10}^{{ - 6}}}$ приведены на фиг. 3. В отношении устойчивости эти результаты полностью согласуются с табл. 2. Ошибки методов ARK2 и ARK2c, как и аналогичных методов ARK21 и ARK21c, уменьшаются при увеличении μ, а ошибка метода ARK2s сначала уменьшается, а при $h\mu \approx {{\varepsilon }^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}}$ начинает возрастать. Таким образом, можно заключить, что неточность оценки собственного значения может привести к расхождению численного решения, а для методов со стабилизированной 1-й стадией – также и к снижению точности при сохранении устойчивости.

4. ТОЧНОСТЬ АДАПТИВНЫХ МЕТОДОВ

Для исследования точности используем три задачи, точные решения которых приводим ниже. Первая из них – задача Капса (см. [21]):

Вторая задача – линейная неавтономная: Третья задача имеет вид

Жесткость задач определяется значением μ, которое (при больших μ) практически совпадает с мерой жесткости ${{M}_{{{\text{stf}}}}}$. Мы решали эти задачи с шагом $h = {1 \mathord{\left/ {\vphantom {1 {30}}} \right. \kern-0em} {30}}$ адаптивными методами, а также методом Розенброка 2-го порядка, который обозначим через Ros2 (этот метод реализован в решателе ode23s системы MATLAB, см. [22]). Ошибки решения в зависимости от μ приведены в табл. 3 для задачи (4.1), в табл. 4 для задачи (4.2) и в табл. 5 для задачи (4.3).

При решении задачи (4.1) все методы показали приемлемые и вполне объяснимые результаты. При повышении жесткости методы ARK21, ARK21c и ARK21s заметно снижают точность, что объясняется снижением реального порядка до 1-го. Снижение порядка подтверждается расчетами с уменьшением размера шага. Для всех методов при $\mu > {{10}^{4}}$ дальнейшее повышение жесткости уже не приводит к заметному увеличению ошибки, которая мало изменяется вплоть до значения μ, при котором решение расходится.

Задача (4.2) оказалась трудной для метода ARK21s, который показывает неудовлетворительную точность при $\mu \geqslant {{10}^{4}}$, хотя и сохраняет устойчивость вплоть до значения $\mu = {{10}^{{155}}}$. Но при этом сходимость обеспечивается только при малых размерах шага. Например, при $\mu = {{10}^{6}}$ и уменьшении шага в 10 раз ошибка мало изменилась и равна $0.195$. И только при $h < {{10}^{{ - 3}}}$ ошибка пропорциональна ${{h}^{2}}$. Приведенный пример показывает, что при решении некоторых жестких задач метод со стабилизированной 1-й стадией может быть менее точным по сравнению с другими адаптивными методами.

Задача (4.3) – пример жесткой задачи, трудной для решения неявными методами. При больших значениях μ все неявные методы, которые мы использовали для решения этой задачи, требовали больших вычислительных затрат и не обеспечивали требуемой точности. Причиной низкой эффективности неявных методов является медленная сходимость простых ньютоновских итераций, которые используются в этих методах. Явные адаптивные методы были более эффективными. Несколько неожиданно методы ARK21 и ARK21c имели 2‑й порядок при больших значениях μ, хотя при решении других жестких задач реальный порядок снижается до 1-го. Метод ARK2 не позволил получить решение с данным шагом при $\mu > {{10}^{4}}$, а ARK2c – при $\mu > 2 \times {{10}^{5}}$.

При решении всех трех задач методы ARK и ARKc показали близкие результаты, но при этом методы ARKc позволяют решать более жесткие задачи. Например, задача (4.2) была успешно решена методом ARK21 вплоть до значения $\mu = {{10}^{{18}}}$, а методом ARK21c – вплоть до $\mu = {{10}^{{30}}}$.

5. ПОСТРОЕНИЕ МЕТОДОВ ARK32 И ARK32C

Методы, построенные на основе стадий (3.1) при $s = 3$, не могут иметь порядок выше 1-го для жестких задач, поэтому они эффективны только при низкой задаваемой точности. Более эффективные методы 2-го порядка получим, задав $s = 4$. Для эффективного решения жестких задач следует задать $\beta = 1$, но для нежестких задач оптимальное значение β равно 2/3 (тогда метод имеет 3-й порядок). Поэтому имеет смысл настраивать в процессе решения не только заключительную формулу интегрирования, но и параметр β, что позволяет наиболее эффективно решать не только жесткие, но и нежесткие задачи.

Характерной особенностью многих жестких задач является наличие в решении коротких пограничных участков (слоев) с быстрым изменением переменных. При прохождении этих участков все значения $h{{\lambda }_{i}}$ невелики, поэтому задача внутри такого участка не является жесткой. Таким образом, все решение можно разбить на “медленные” жесткие и “быстрые” нежесткие участки. На нежестких участках изменения переменных наиболее значительны, поэтому они вносят значительный вклад в ошибку численного решения. Это является еще одной причиной для построения методов, имеющих повышенный порядок сходимости при решении нежестких задач.

Построим такой метод на основе стадий (3.1) при $s = 4$. Параметры α и β задаем исходя из полученных на предыдущем шаге оценок ${{\tilde {z}}_{i}}$, при этом $\beta = 1 - \alpha $, а α принимаем в виде (3.2), где $\bar {\alpha } = {1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-0em} 3}$. Такие значения обеспечивают эффективное решение как жестких, так и нежестких задач. Шаг интегрирования выполняем по формуле

где ${{{\mathbf{u}}}_{i}} = {{({{{\mathbf{F}}}_{i}} - {{{\mathbf{F}}}_{{i - 1}}})} \mathord{\left/ {\vphantom {{({{{\mathbf{F}}}_{i}} - {{{\mathbf{F}}}_{{i - 1}}})} {(\beta {{\alpha }^{{i - 2}}})}}} \right. \kern-0em} {(\beta {{\alpha }^{{i - 2}}})}}$, а коэффициенты вектора ${{{\mathbf{d}}}_{3}}$ вычисляем по формулам (2.2), (2.5), (2.7).

При реализации алгоритмов с автоматическим выбором размера шага используют вложенную пару формул (методов) вычисления ${{{\mathbf{y}}}_{1}}$ и ${{{\mathbf{\hat {y}}}}_{1}}$, позволяющую получить оценку локальной ошибки как норму вектора ${{{\mathbf{y}}}_{1}} - {{{\mathbf{\hat {y}}}}_{1}}$. В [7], [8] при $\beta = 1$, $s = 4$ использовалась вложенная формула вида

где компоненты вектора ${{{\mathbf{\hat {d}}}}_{2}}$ вычисляем согласно (2.5) с использованием скалярной функции устойчивости вложенного метода $\hat {Q}(z)$. Но если $\beta < 1$ и внутри интервала $\left[ {{{t}_{0}}\, + \beta h,\;{{t}_{0}}\, + h} \right]$ происходит быстрое изменение правой части, то на значениях (5.1) и (5.2) это никак не сказывается, вследствие чего снижается точность решения (см. пример в разд. 2.4 книги [7]). Поэтому при $\beta < 1$ оценка ошибки на основе (5.2) ненадежна, и для построения формулы вычисления ${{{\mathbf{\hat {y}}}}_{1}}$ следует использовать также значение ${{{\mathbf{f}}}_{1}} = {\mathbf{f}}({{t}_{0}} + h,{{{\mathbf{y}}}_{1}})$. Это же значение является результатом выполнения первой стадии следующего шага, поэтому такие вложенные пары называют FSAL (First Same As Last, см. [23]).

Для FSAL-пары принимаем

Вычисление очередной компоненты векторов ${{{\mathbf{\hat {d}}}}_{2}}$, ${{{\mathbf{\hat {d}}}}_{3}}$ и ${{{\mathbf{\hat {d}}}}_{4}}$ выполняем, используя оценку $\tilde {z}$ по этой компоненте: Эти значения выбраны с учетом условия 2-го порядка формулы (5.3) для нежестких компонент и 1-го порядка для жестких компонент. Кроме этого, функция устойчивости $\hat {R}(z) = 1 + z + {{\hat {d}}_{2}}{{z}^{2}} + {{\hat {d}}_{3}}{{z}^{3}} + {{\hat {d}}_{4}}{{d}_{3}}{{z}^{4}}$ удовлетворяет условию $\hat {R}( - 1{\text{/}}\gamma ) \approx 0$, а параметр g остается свободным (мы принимаем $g = {1 \mathord{\left/ {\vphantom {1 8}} \right. \kern-0em} 8}$). Метод, основанный на формулах (5.1), (5.3), (5.4), обозначим через ARK32.

Рассмотрим теперь метод с коррекцией жестких компонент, которую выполняем после вычислений по формулам (5.1), (5.3), (5.4). Если для i-й компоненты ${{\tilde {z}}_{i}} < - 4.5$, то пересчитываем эту компоненту согласно формулам

где ${{y}_{1}}$, ${{u}_{3}}$, ${{v}_{4}}$ суть i-е компоненты векторов ${{{\mathbf{y}}}_{1}}$, ${{{\mathbf{u}}}_{3}}$, ${{{\mathbf{v}}}_{4}}$. Метод с коррекцией жестких компонент обозначим через ARK32c. При решении некоторых жестких задач этот метод значительно более эффектен, чем ARK32.

Чтобы адаптивные методы были эффективными, необходимо ограничить рост внутренних функций устойчивости. В методах, основанных на стадиях (3.1), это осуществляется с помощью настройки параметра α. Однако реальный порядок таких методов для жестких задач не может быть выше 2-го. Построить эффективный метод, имеющий при решении жестких задач 3-й порядок, удалось только на основе явного метода, имеющего 2-й псевдостадийный порядок (см. [7]). Формулы явного адаптивного метода 3-го порядка со стабилизированной 1-й стадией приведены в [7], [8]. Обозначим этот метод через ARK3s.

6. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ

Для экспериментов с переменным размером шага были выбраны 5 задач из [1], основные характеристики которых приведены в табл. 1. Для управления размером шага используем стандартную процедуру (см. [1], [7]). Допустимую ошибку обозначим через Tol. Для всех задач задаем допустимую относительную ошибку $Rtol = Tol$, а допустимую абсолютную ошибку принимаем в виде $Atol = Tol$ для задач VDPOL и OREGO, $Atol = {{10}^{{ - 6}}}{\kern 1pt} \times Tol$ для ROBER, $Atol = {{10}^{{ - 4}}}{\kern 1pt} \times Tol$ для HIRES и $Atol = {{10}^{{ - 2}}}{\kern 1pt} \times Tol$ для CUSP. Вычислительные затраты оцениваем числом вычислений правой части Nf, а в качестве оценки точности решения, как и в [24], используем значение

где ${{y}_{i}}$ – точное, а ${{\tilde {y}}_{i}}$ – численное решение i-й компоненты в конце интервала. Для решения, наряду с рассмотренными в [7], [8] методами ARK21s, ARK2s и ARK3s, применялись новые методы ARK2с, ARK32 и ARK32с. Метод ARK2с отличается от ARK32с только значением $\beta = 1$.

Результаты для трех значений Tol приведены в табл. 6. Обсудим сначала результаты решения задачи VDPOL, которая задается уравнениями

при $\mu = {{10}^{6}}$. Метод ARK21s эффективен только при низкой задаваемой точности, но он позволяет решать очень жесткие задачи, например, задача (6.1) была успешно решена при $\mu = {{10}^{{46}}}$. Метод ARK2с заметно уступает методу ARK32с по точности, что объясняется его меньшим порядком на нежестких участках. При решении других задач методы ARK21s и ARK2с также оказались менее эффективными по сравнению с другими методами, поэтому мы исключаем их из дальнейшего рассмотрения. Некоторые результаты методов ARK21, ARK21s и ARK21с приведены в [7] (эти методы можно рекомендовать для решения очень жестких задач с низкой точностью).

Обратим теперь внимание на результаты методов ARK2s, ARK3s и ARK32 при решении задачи ROBER. При $Tol = {{10}^{{ - 2}}}$ вычислительные затраты Nf намного больше, чем при меньших значениях Tol. Объясняется это тем, что при увеличении размера шага возникает неустойчивость численного решения, которая приводит к резким колебаниям величины шага и увеличению числа шагов. Аналогичное явление проявляется и при решении задачи CUSP методами ARK3s и ARK32 (хотя задача CUSP значительно менее жесткая, чем задачи VDPOL и OREGO, при решении которых этот эффект отсутствует).

Для более наглядного сравнения методов мы объединили результаты решения всех пяти задач и приводим на фиг. 4 зависимости общих затрат на их решение от усредненного для этих задач значения scd при $Tol = {{10}^{{ - i}}}$, $i = 2,\; \ldots ,7$. При $Tol = {{10}^{{ - 4}}}$ все методы показывают примерно одинаковые результаты. При повышении точности ($Tol < {{10}^{{ - 4}}}$) преимущество имеет метод более высокого порядка ARK3s, а при $Tol > {{10}^{{ - 4}}}$ метод ARK32c превосходит другие методы.

Анализ приведенных результатов позволяет сделать вывод, что наиболее эффективным явным адаптивным методом при решении жестких задач с умеренной точностью является ARK32c, а для расчетов с повышенной точностью можно использовать метод ARK3s.