Журнал вычислительной математики и математической физики, 2023, T. 63, № 4, стр. 557-572

Оценка области абсолютной устойчивости численной схемы решения жестких задач Коши методом продолжения решения по параметру

Е. Б. Кузнецов^1, *, С. С. Леонов^1, 2, **, Е. Д. Цапко^1, ***

¹ МАИ (национальный исследовательский университет)
125993 Москва, Волоколамское ш., 4, Россия

² РУДН
117198 Москва, ул. Миклухо-Маклая, 6, Россия

^* E-mail: kuznetsov@mai.ru
^** E-mail: powerandglory@yandex.ru
^*** E-mail: zapkokaty@gmail.com

Поступила в редакцию 11.08.2022
После доработки 01.09.2022
Принята к публикации 15.12.2022

EDN: IPJBXS
DOI: 10.31857/S0044466923040129

Полный текст (PDF)

Аннотация

При моделировании физических и технологических процессов исследователи часто сталкиваются с решением жестких начальных задач. Нахождение их точного аналитического решения в большинстве случаев затруднительно. В то же время применение численных схем для их решения не всегда позволяет получить достаточно точное решение за приемлемое расчетное время. Более того, для некоторого класса задач численные схемы решения оказываются непригодными из-за недостаточной устойчивости. В статье рассматриваются численные методы на основе продолжения решения по аргументам различного вида, которые позволяют увеличить устойчивость явных численных схем. Наиболее часто используемый наилучший аргумент оказывается малоприменим для решения задач, скорость роста интегральных кривых которых является сверхстепенной или близка к экспоненциальной. Авторами ранее была предложена модификация наилучшего аргумента, которая позволила сгладить указанные недостатки. В настоящей работе получена оценка области абсолютной устойчивости явной схемы метода Эйлера при решении задач, преобразованных к модифицированному наилучшему аргументу специального вида, и уточнено доказательство аналогичной оценки для начальных задач, преобразованных к наилучшему аргументу. Проведена апробация полученных теоретических оценок и дан анализ применения модифицированного наилучшего аргумента продолжения решения на примере тестовой начальной задачи. Библ. 41. Фиг. 2. Табл. 1.

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

1. ВВЕДЕНИЕ

Жесткие начальные и краевые задачи возникают во многих областях науки и техники: аэро- и гидродинамика, химическая кинетика, теория горения и теория управления (см. [1]). Л. Прандтль был одним из первых, кто обнаружил эффект пограничного слоя, рассматривая движение вязкой жидкости с малым трением, описываемое системой уравнений Навье–Стокса (см. [2]). К подобным задачам относится приближенное решение уравнений с малым параметром при старшей производной. Они получили название сингулярно возмущенных уравнений и также относятся к классу жестких задач. Методы решения подобных задач получили развитие в работах А.Н. Тихонова [3] и А.Б. Васильевой [4]. В дальнейшем, в работах В.Ф. Бутузова, А.Б. Васильевой, Н.Н. Нефедова и их учеников (см., например, [5–7]) были рассмотрены задачи не только с пограничным слоем, но и с внутренними слоями. Последние получили название задач с контрастными структурами.

Вычислительные проблемы численного решения жестких задач были отмечены уже в конце 40-х годов прошлого века (см. [8]). Позже было разработано множество специализированных подходов, которые подробно описаны в [9] и монографии Э. Хайрера и Г. Ваннера [10]. Вопросы устойчивости методов семейства Рунге–Кутты для жестких начальных задач подробно рассмотрены в монографии К. Деккера и Я. Вервера [11]. Отметим основные и новые работы, связанные с численным решением жестких задач. Одним из первых применил формулы дифференцирования назад к решению жестких систем американский исследователь Гир (см. [12]), что привело к развитию целого ряда неявных схем. Созданию эффективных неявных схем посвящены циклы работ Н.Н. Калиткина и соавт. [13, 14], Г.Ю. Куликова и соавт. [15–20] и Л.М. Скворцова [21–24]. В статьях Н.Н. Калиткина исследуются особенности численного интегрирования нелинейных обыкновенных дифференциальных и дифференциально-алгебраических уравнений и предложен класс обратных (диагонально-неявных) схем Рунге–Кутты, а также неявные схемы для задач химической кинетики. В работах Г.Ю. Куликова разработан новый класс неявных (диагональных и симметричных) гнездовых схем численного решения жестких задач, в том числе, с глобальным контролем точности. Показана эффективность гнездовых методов на ряде тестовых задач. Помимо этого, исследована устойчивость симметричных формул Рунге–Кутты при решении гамильтоновых задач, а также при квадратичной экстраполяции численного решения. Предложенные Л.М. Скворцовым эффективные реализации неявных схем и модификации диагонально-неявных схем показали свою применимость к решению жестких начальных задач как для дифференциальных, так и для дифференциально-алгебраических уравнений. Указанные работы показывают, что неявные схемы дают высокую точность и устойчивость при решении жестких задач, но по быстродействию они намного уступают явным схемам.

Применению явных методов к решению жестких задач посвящена монография Е.А. Новикова [25] и его учебное пособие в соавторстве с Ю.В. Шорниковым [26]. К этой же тематике относятся статьи этого же автора [27–29] и Л.М. Скворцова [30], [31]. В статьях Е.А. Новикова приводятся результаты сразу по нескольким областям исследования: настройка (адаптация) численной схемы на основе критерия L-устойчивости по схеме Ческино, конструирование явных схем с расширенными областями устойчивости, решение жестких задач явными схемами с малой и умеренной точностью. Явным адаптивным (настраиваемым на конкретную задачу) схемам семейства Рунге–Кутты посвящены указанные выше работы Л.М. Скворцова. В работах В.И. Лебедева [32], [33] рассматривается применение многочленов Чебышева к решению жестких задач явными методами. Отмеченные статьи ярко указывают на недостатки применения явных схем к решению жестких задач. Большинство явных схем дают только малую или умеренную точности. При расширении области устойчивости в большинстве случаев явная схема усложняется, сокращается ее быстродействие. То же самое относится и к адаптивным явным методам, в которых возникают также трудности подбора структуры численной схемы. Тем не менее явные схемы предпочтительнее неявных, особенно для задач большой размерности.

Для устранения указанного недостатка и возможности использования явных схем, которые значительно проще алгоритмически и менее требовательны к вычислительным ресурсам, для решения жестких задач в монографии В.И. Шалашилина и Е.Б. Кузнецова [34] было предложено использовать метод продолжения решения по наилучшему аргументу. Ими же было доказано, что данный метод позволяет расширить область абсолютной устойчивости для явной схемы метода Эйлера. Отметим также значение наилучшей параметризации (позволяющей в ряде случаев улучшить ситуацию с применением методов с фиксированным размером шага интегрирования) для адаптивности квазисогласованных методов Нордсика и явных численных схем со свойством двойной квазисогласованности при решении умеренно жестких задач (см. [35], [36]).

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

2. МЕТОД НАИЛУЧШЕЙ ПАРАМЕТРИЗАЦИИ

Рассмотрим задачу

(1)

$\frac{{dy}}{{dt}} = f(t,y),\quad t \in [{{t}_{0}},T],$

с начальным условием

(2)

$y({{t}_{0}}) = {{y}^{{(0)}}},$

где ${{t}_{0}}$ – начальное значение аргумента задачи, ${{y}^{{(0)}}}$ – значение искомой функции $y(t)$ в начальной точке, $T$ – заданное значение правой границы отрезка изменения аргумента задачи, $f(t,y)$ – функция правой части, которая, в общем случае, может иметь особенности, такие как предельные особые точки или точки бифуркации, а также содержать множитель, принимающий большие абсолютные значения.

Применение метода продолжения решения по наилучшему аргументу, называемому иначе методом наилучшей параметризации или методом длины дуги, для решения задач Коши основывается на переходе к новому аргументу.

2.1. Наилучший аргумент

Наилучший аргумент $\lambda $, рассматриваемый В.И. Шалашилиным и Е.Б. Кузнецовым в монографии [34], отсчитывается по касательной вдоль интегральной кривой рассматриваемой задачи Коши, квадрат дифференциала дуги интегральной кривой имеет для задачи (1)–(2) вид

(3)

$d{{\lambda }^{2}} = d{{y}^{2}} + d{{t}^{2}}.$

Используя аргумент $\lambda $, задаваемый соотношением (3), применительно к задаче (1)–(2), получим преобразованную систему относительно нового аргумента $\lambda $ следующего вида:

(4)

$\begin{gathered} \frac{{dy}}{{d\lambda }} = \frac{{f(t,y)}}{{\sqrt {1 + {{f}^{2}}(t,y)} }},\quad \frac{{dt}}{{d\lambda }} = \frac{1}{{\sqrt {1 + {{f}^{2}}(t,y)} }}, \\ y(0) = {{y}^{{(0)}}},\quad t(0) = {{t}_{0}}. \\ \end{gathered} $

В [34, 38] было изучено влияние такого преобразования на область абсолютной устойчивости явного метода Эйлера. Доказано, что переход к наилучшему аргументу увеличивает область абсолютной устойчивости явного метода Эйлера для задачи (4), что уменьшает вычислительные трудности, характерные для жестких систем.

2.2. Экспоненциальный наилучший аргумент

В [37] предложен новый подход к решению жестких и сверхжестких начальных задач. Он основан на использовании вместо наилучшего аргумента его модификации следующего вида:

(5)

$d{{\mu }^{2}} = d{{y}^{2}} + \exp ( - 2\gamma t){\kern 1pt} d{{t}^{2}},$

где $\gamma $ – заданный параметр, выбираемый из условия его превосходства над показателем жесткости решаемой задачи.

Задача (1), преобразованная к новому аргументу $\mu $ задаваемым соотношением вида (5), примет вид

(6)

$\begin{gathered} \frac{{dy}}{{d\mu }} = \frac{{f(t,y){\kern 1pt} \exp (\gamma t)}}{{\sqrt {1 + \exp (2\gamma t){\kern 1pt} {{f}^{2}}(t,y)} }},\quad \frac{{dt}}{{d\mu }} = \frac{{\exp (\gamma t)}}{{\sqrt {1 + \exp (2\gamma t){\kern 1pt} {{f}^{2}}(t,y)} }}, \\ y(0) = {{y}^{{(0)}}},\quad t(0) = {{t}_{0}}. \\ \end{gathered} $

В [37] на тестовых примерах была показана численная эффективность разработанного подхода в сравнении с традиционными методами и наилучшей параметризацией, когда скорость роста интегральных кривых рассматриваемой задачи является экспоненциальной. В настоящей статье будет дана оценка области абсолютной устойчивости явного метода Эйлера при решении преобразованной экспоненциальным модифицированным наилучшим аргументом задачи (6). Будет доказано и продемонстрировано на тестовой задаче, что применение нового подхода увеличивает область абсолютной устойчивости явного метода Эйлера и дает возможность более эффективно использовать явный метод Эйлера при решении жестких начальных задач.

3. УСТОЙЧИВОСТЬ ЯВНОЙ СХЕМЫ МЕТОДА ЭЙЛЕРА

Исследование области устойчивости и спектральных характеристик дифференциальных операторов, как было предложено в [39], рассматривают на примере задачи, получившей название задачи Далквиста:

(7)

$\frac{{dy}}{{dt}} = ay,\quad y(0) = {{y}_{0}},$

где $a$ – некоторая, в общем случае комплексная, константа. В стандартной задаче Далквиста константа $a$ предполагается действительной и отрицательной. Но рассматривая возможность комплексных значений константы $a$, можно области абсолютной устойчивости придать ясный геометрический смысл, как это показано ниже в замечании 2.

Задача (7) моделирует локальное поведение решения дифференциального уравнения (1) в том смысле, что в окрестности любой точки $({{t}_{0}},{{y}_{0}})$ решение уравнения (1) ведет себя так, как решение линеаризованного уравнения

$\frac{{dY}}{{dt}} = {{f}_{y}}({{t}_{0}},{{y}_{0}})Y$

в окрестности нуля. Тогда, если $a$ – собственное значение линеаризованной задачи, то по поведению решений разностного метода на уравнении (7) можно предсказать их поведение на произвольном дифференциальном уравнении.

Теорема 1. Разностная схема явного метода Эйлера для задачи (7) является абсолютно устойчивой, если шаг интегрирования удовлетворяет неравенству

(8)

$\left| {{{h}_{t}}} \right| \leqslant \frac{2}{{\left| a \right|}}.$

Доказательство. Применяя явную схему метода Эйлера к задаче (7) в окрестности точки $M({{t}_{m}},{{y}_{m}})$, получим разностное уравнение вида

(9)

${{y}_{{m + 1}}} = {{y}_{m}} + a{{h}_{t}}{{y}_{m}} = (1 + a{{h}_{t}}){{y}_{m}},$

где ${{h}_{t}}$ – шаг интегрирования.

По определению разностный метод абсолютно устойчив, если все корни характеристического уравнения не превосходят по модулю единицу (см. [11], [40]). Для схемы (9) собственное значение $\xi = 1 + a{{h}_{t}}$.

При положительных значениях ${{h}_{t}}$ неравенство $\left| {1 + a{{h}_{t}}} \right| \leqslant 1$ может быть справедливо только при $a < 0$. В этом случае условие устойчивости схемы (9) примет вид

$0 < {{h}_{t}} \leqslant \frac{2}{{\left| a \right|}}.$

В случае отрицательных значений ${{h}_{t}}$ неравенство $\left| {1 + a{{h}_{t}}} \right| \leqslant 1$ может быть справедливо только при $a > 0$. В этом случае условие устойчивости схемы (9) перепишется в виде

$ - \frac{2}{{\left| a \right|}} \leqslant {{h}_{t}} < 0.$

Объединяя два этих неравенства, получим условие (8).

В тривиальных случаях при ${{h}_{t}} = 0$ или $a = 0$ абсолютная устойчивость также имеет место. В обоих случаях максимальное собственное значение ${{\xi }_{{\max }}} = 1$. Однако при практических расчетах случай ${{h}_{t}} = 0$ исключается.

3.1. Область устойчивости метода Эйлера для задачи, преобразованной к наилучшему аргументу

Задачу (7) можно преобразовать к наилучшему аргументу $\lambda $ вида (3). Преобразованная задача примет вид

(10)

$\begin{gathered} \frac{{dy}}{{d\lambda }} = \frac{{ay}}{{\sqrt {1 + {{a}^{2}}{{y}^{2}}} }},\quad \frac{{dt}}{{d\lambda }} = \frac{1}{{\sqrt {1 + {{a}^{2}}{{y}^{2}}} }}, \\ y(0) = {{y}^{{(0)}}},\quad t(0) = {{t}_{0}}. \\ \end{gathered} $

Условие абсолютной устойчивости для нее получено в статье Е.Б. Кузнецова и В.И. Шалашилина [38]. Откуда следует

Теорема 2. Разностная схема явного метода Эйлера для преобразованной к наилучшему аргументу $\lambda $ вида (3) задачи (10) является абсолютно устойчивой, если шаг интегрирования удовлетворяет неравенству

(11)

$\left| {{{h}_{\lambda }}} \right| \leqslant \frac{{2{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}{{\left| a \right|}}.$

Здесь ${{y}_{m}}$ – решение задачи (10), полученное на $m$-м шаге явным методом Эйлера.

Доказательство. Приведем уточненное доказательство теоремы. Обозначим правые части уравнений задачи (10) как

$\begin{array}{*{20}{l}} {{{f}_{1}}(y) = \frac{{ay}}{{\sqrt {1 + {{a}^{2}}{{y}^{2}}} }},\quad {{f}_{2}}(y) = \frac{1}{{\sqrt {1 + {{a}^{2}}{{y}^{2}}} }}.} \end{array}$

В общем случае, правые части системы (10) могут зависеть от аргумента исходной задачи $t$, но в данном случае такой зависимости нет, поэтому представим функции ${{f}_{1}}(y)$ и ${{f}_{2}}(y)$, используя разложение по формуле Тейлора до первой степени включительно в окрестности $y = {{y}_{m}}$. Рассмотрим сначала первое уравнение задачи (10), разложив его правую часть

${{f}_{1}}(y) = \frac{{a{{y}_{m}}}}{{\sqrt {1 + {{a}^{2}}y_{m}^{2}} }} + \frac{a}{{{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}\left( {y - {{y}_{m}}} \right) + O\left( {{{{\left( {y - {{y}_{m}}} \right)}}^{2}}} \right).$

Рассмотрим линеаризованное первое уравнение вида

(12)

$\frac{{dy}}{{d\lambda }} = \frac{{a{{y}_{m}}}}{{\sqrt {1 + {{a}^{2}}y_{m}^{2}} }} + \frac{a}{{{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}\left( {y - {{y}_{m}}} \right).$

Введем замену

(13)

$y = \bar {y} + y_{*}^{{}},$

где $y_{*}^{{}}$ – координата точки покоя линеаризованной системы, получаемая из решения уравнения

$\frac{{a{{y}_{m}}}}{{\sqrt {1 + {{a}^{2}}y_{m}^{2}} }} + \frac{a}{{{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}\left( {y - {{y}_{m}}} \right) = 0$

и равная

$y_{*}^{{}} = - {{a}^{2}}y_{m}^{3}.$

Используя замену (13), преобразуем уравнение (12) к виду

$\frac{{d\bar {y}}}{{d\lambda }} = \frac{a}{{{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}{\kern 1pt} \bar {y}.$

Линеаризуем второе уравнение задачи (10):

${{f}_{2}}(y) = \frac{1}{{\sqrt {1 + {{a}^{2}}{{y}^{2}}} }} = \frac{1}{{ay}}\frac{{ay}}{{\sqrt {1 + {{a}^{2}}{{y}^{2}}} }} = \frac{1}{{ay}}{\kern 1pt} \left( {\frac{a}{{{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}{\kern 1pt} \left( {y - y_{*}^{{}}} \right) + O\left( {{{{\left( {y - y_{*}^{{}}} \right)}}^{2}}} \right)} \right).$

Разложим первый множитель по формуле Тейлора до членов первого порядка в окрестности точки $y_{*}^{{}}$:

$\frac{1}{{ay}} = \frac{1}{{ay_{*}^{{}}}} - \frac{1}{{ay_{*}^{2}}}\left( {y - y_{*}^{{}}} \right) + O\left( {{{{\left( {y - y_{*}^{{}}} \right)}}^{2}}} \right).$

Подставляя это разложение в выражение для ${{f}_{2}}(y)$, получим

${{f}_{2}}(y) = - \frac{1}{{{{a}^{2}}y_{m}^{3}{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}{\kern 1pt} (y - {{y}_{*}}) + O\left( {{{{(y - {{y}_{*}})}}^{2}}} \right).$

Используя полученное разложение, перейдем к линеаризованной системе для задачи (10), преобразованной к переменной $\bar {y} = y + {{a}^{2}}y_{m}^{3}$:

(14)

$\frac{{d\bar {y}}}{{d\lambda }} = \frac{a}{{{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}{\kern 1pt} \bar {y},\quad \frac{{dt}}{{d\lambda }} = - \frac{1}{{{{a}^{2}}y_{m}^{3}{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}{\kern 1pt} \bar {y}.$

Явная схема метода Эйлера для системы (14) имеет вид:

(15)

${{\bar {y}}_{{m + 1}}} = {{\bar {y}}_{m}} + {{h}_{\lambda }}\frac{a}{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}{\kern 1pt} {{\bar {y}}_{m}},\quad {{t}_{{m + 1}}} = {{t}_{m}} - {{h}_{\lambda }}\frac{1}{{{{a}^{2}}y_{m}^{3}{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}{\kern 1pt} {{\bar {y}}_{m}}.$

Поскольку линеаризованная система (14) получена использованием линейной замены и правая часть системы уравнений задачи (10) ограничена, то условие абсолютной устойчивости разностной схемы (15) будет эквивалентно условию абсолютной устойчивости для линеаризованной системы исходной задачи (10).

Найдем собственные значения $\xi $, рассмотрев характеристическое уравнение

$\left| {A - \xi E} \right| = \left| {\begin{array}{*{20}{c}} {1 + {{h}_{\lambda }}\frac{a}{{{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}} - \xi }&0 \\ { - {{h}_{\lambda }}\frac{1}{{{{a}^{2}}y_{m}^{3}{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}}&{1 - \xi } \end{array}} \right| = \left( {\xi - 1} \right)\left( {\xi - \left( {1 + {{h}_{\lambda }}\frac{a}{{{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}} \right)} \right) = 0.$

Получим собственные значения

(16)

${{\xi }_{1}} = 1,\quad {{\xi }_{2}} = 1 + {{h}_{\lambda }}\frac{a}{{{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}.$

Поскольку ${{\xi }_{1}} = 1$, разностная схема (15) будет абсолютно устойчивой, если второе собственное значение ${{\xi }_{2}}$ не превосходит по модулю единицу, т.е.

(17)

$\left| {1 + {{h}_{\lambda }}\frac{a}{{{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}} \right| \leqslant 1.$

При значениях шага ${{h}_{\lambda }} > 0$ абсолютная устойчивость схемы (15) возможна только при значениях параметра $a < 0$. Разрешая неравенство (17) относительно шага ${{h}_{\lambda }}$, полагая $a = - \left| a \right|$, получим

$0 < {{h}_{\lambda }} \leqslant \frac{{2{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}{{\left| a \right|}}.$

При отрицательных значениях шага ${{h}_{\lambda }}$ абсолютная устойчивость схемы (15) возможна только при значениях параметра $a > 0$. Разрешая неравенство (17) относительно шага ${{h}_{\lambda }}$, полагая $a = \left| a \right|$, получим

$ - \frac{{2{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}{{\left| a \right|}} \leqslant {{h}_{\lambda }} < 0.$

Объединяя два этих неравенства, получим условие (11).

В тривиальных случаях при ${{h}_{\lambda }} = 0$ или $a = 0$ абсолютная устойчивость также имеет место. В обоих случаях максимальное собственное значение ${{\xi }_{{\max }}} = {{\xi }_{1}} = {{\xi }_{2}} = 1$. Однако при практических расчетах случай ${{h}_{\lambda }} = 0$ исключается.

Замечание 1. Приводимое выше доказательство теоремы, сформулированной в работе Е.Б. Кузнецова и В.И. Шалашилина [34], является новым, обобщая доказательство, приводимое в статье [34]. Е.Б. Кузнецовым и В.И. Шалашилиным принималось, что второе уравнение системы задачи (10) вносит пренебрежимо малый вклад в характер устойчивости приближенного решения, что является верным только в окрестности предельной особой точки. Вследствие чего второе уравнение исключалось из рассмотрения. В приведенном здесь доказательстве этот малый недостаток устранен.

Замечание 2. Известно, что областью устойчивости явного метода Эйлера для исходной задачи (7) является круг единичного радиуса с центром в точке $( - 1,0)$. В [34] доказано, что для преобразованной задачи (10) областью устойчивости будет круг радиуса ${{\rho }^{{3/2}}}$ с центром в точке $({{\rho }^{{3/2}}},0)$, где $\rho = \left| {1 + {{a}^{2}}y_{m}^{2}} \right|$. При комплексных значениях $a$ выражение для $\rho $ можно переписать в виде

$\rho = \sqrt {{{{(1 + ({{\alpha }^{2}}{{\beta }^{2}})y_{m}^{2})}}^{2}} + 4{{\alpha }^{2}}{{\beta }^{2}}y_{m}^{4}} ,$

где $\alpha $ и $\beta $ – действительная и мнимая части параметра $a$. Доказательство данного утверждения здесь опускается.

Замечание 3. Само доказательство теоремы 2 дает множество интересных следствий. Одно из них заключается в том, что одно из собственных значений ${{\xi }_{1}} = 1$. Это означает, что даже при выборе значения шага из условия (11) максимальное значение ${{\xi }_{{{\text{max}}}}} = 1$, т.е. разностная схема находится в пограничном состоянии между абсолютной устойчивостью и ее потерей. Даже для линеаризованной системы абсолютная устойчивость может быть потеряна при наложении ошибок округления. При линеаризации системы уравнений задачи (10) отбрасываются нелинейные члены, учет которых также может привести к потере абсолютной устойчивости. Данное наблюдение позволяет объяснить тот факт, что для ряда жестких начальных задач для автономных систем дифференциальных уравнений преобразование к наилучшему аргументу либо не дает вычислительных преимуществ, либо эти преимущества минимальны.

Замечание 4. Представленное выше замечание показывает, что для анализа области абсолютной устойчивости разностных схем решения начальных задач, преобразованных к наилучшему аргументу, недостаточно рассмотрения тестовой задачи Далквиста вида (7). Более целесообразно исследовать область абсолютной устойчивости для обобщенной неавтономной задачи Далквиста вида

(18)

$\frac{{dy}}{{dt}} = ay + bt,\quad y(0) = {{y}_{0}}.$

Решение данной задачи является предметом дальнейших исследований.

3.2. Область устойчивости метода Эйлера для задачи, преобразованной к модифицированному наилучшему аргументу

Применение экспоненциального наилучшего аргумента вида (5) призвано, с одной стороны, упростить вид преобразованной системы и уменьшить трудозатраты, возникающие при решении преобразованных начальных задач. С другой стороны, экспоненциальный наилучший аргумент может позволить расширить область абсолютной устойчивости и снять ограничения, присущие наилучшему аргументу, отмеченные в замечании 3.

Рассмотрим преобразованную к экспоненциальному наилучшему аргументу $\mu $ вида (5) за-дачу (7):

(19)

$\begin{gathered} \frac{{dy}}{{d\mu }} = \frac{{ay\exp \left( {\gamma t} \right)}}{{\sqrt {1 + {{a}^{2}}{{y}^{2}}\exp \left( {2\gamma t} \right)} }},\quad y(0) = {{y}_{0}}, \\ \frac{{dt}}{{d\mu }} = \frac{{\exp \left( {\gamma t} \right)}}{{\sqrt {1 + {{a}^{2}}{{y}^{2}}\exp \left( {2\gamma t} \right)} }},\quad t(0) = 0. \\ \end{gathered} $

Теорема 3. Разностная схема явного метода Эйлера для преобразованной к экспоненциальному наилучшему аргументу $\mu $ вида (5) задачи (19) является абсолютной устойчивой для значений параметра $\gamma $, удовлетворяющих условию

(20)

$a{\kern 1pt} \gamma \leqslant 0,$

если шаг интегрирования удовлетворяет неравенству

(21)

$\left| {{{h}_{\mu }}} \right| \leqslant \frac{{4{{{\left[ {1 + {{a}^{2}}y_{m}^{2}\exp \left( {2\gamma {{t}_{m}}} \right)} \right]}}^{{3/2}}}}}{{\left| {{{D}_{{\max }}}} \right|\exp \left( {\gamma {{t}_{m}}} \right)}},$

где ${{y}_{m}}$ и ${{t}_{m}}$ – решение задачи (19), полученное на $m$-м шаге явным методом Эйлера,

${{D}_{{\max }}} = \left\{ {\begin{array}{*{20}{l}} {{{D}_{1}},\quad a + \gamma \geqslant - \frac{{2(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}{{{))}}^{{3/2}}}}}{{{{h}_{\mu }}\exp (\gamma {{t}_{m}})}},} \\ {{{D}_{2}},\quad a + \gamma < - \frac{{2(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}{{{))}}^{{3/2}}}}}{{{{h}_{\mu }}\exp (\gamma {{t}_{m}})}},} \end{array}} \right.$

$\begin{array}{*{20}{l}} {{{D}_{1}} = a + \gamma + \sqrt {{{{(a - \gamma )}}^{2}} - 4{{a}^{3}}\gamma y_{m}^{2}\exp (2\gamma {{t}_{m}})} ,} \\ {{{D}_{2}} = a + \gamma - \sqrt {{{{(a - \gamma )}}^{2}} - 4{{a}^{3}}\gamma y_{m}^{2}\exp (2\gamma {{t}_{m}})} .} \end{array}$

Доказательство. Линеаризуем правые части уравнений задачи (19), используя разложение по формуле Тэйлора в окрестности точки с координатами $y = {{y}_{m}}$ и $t = {{t}_{m}}$ до членов первой степени:

$\begin{gathered} {{f}_{1}}(y,t) = \frac{{ay\exp \gamma t}}{{\sqrt {1 + {{a}^{2}}{{y}^{2}}\exp 2\gamma t} }} = \frac{{a{{y}_{m}}\exp (\gamma {{t}_{m}})}}{{\sqrt {1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}})} }} + \frac{{a\exp (\gamma {{t}_{m}})}}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}}(y - {{y}_{m}}) + \\ \, + \frac{{a{{y}_{m}}\gamma \exp (\gamma {{t}_{m}})}}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}}(t - {{t}_{m}}) + O((y - {{y}_{m}}{{)}^{2}} + {{(t - {{t}_{m}})}^{2}}), \\ \end{gathered} $

$\begin{gathered} {{f}_{2}}(y,t) = \frac{{\exp \gamma t}}{{\sqrt {1 + {{a}^{2}}{{y}^{2}}\exp 2\gamma t} }} = \frac{{\exp (\gamma {{t}_{m}})}}{{\sqrt {1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}})} }} - \frac{{{{y}_{m}}{{a}^{2}}\exp (3\gamma {{t}_{m}})}}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}}(y - {{y}_{m}}) + \\ \, + \frac{{\gamma \exp (\gamma {{t}_{m}})}}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}}(t - {{t}_{m}}) + O((y - {{y}_{m}}{{)}^{2}} + {{(t - {{t}_{m}})}^{2}}). \\ \end{gathered} $

Запишем линеаризованную систему уравнений для задачи (19):

(22)

$\begin{gathered} \frac{{dy}}{{d\mu }} = \frac{1}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}}\left[ {a\exp (\gamma {{t}_{m}})y + a{{y}_{m}}\gamma \exp (\gamma {{t}_{m}})t + {{a}^{3}}y_{m}^{3}\exp (3\gamma {{t}_{m}}) - {{a}^{2}}{{y}_{m}}{{t}_{m}}\gamma \exp (\gamma {{t}_{m}})} \right], \\ \frac{{dt}}{{d\mu }} = \frac{1}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}}\left[ { - {{y}_{m}}{{a}^{2}}\exp (3\gamma {{t}_{m}})y + \gamma \exp (\gamma {{t}_{m}})t + } \right. \\ \left. {\, + \exp (\gamma {{t}_{m}}) + 2{{a}^{2}}y_{m}^{2}\exp (3\gamma {{t}_{m}}) - \gamma {{t}_{m}}\exp (\gamma {{t}_{m}})} \right]. \\ \end{gathered} $

В линеаризованной системе проведем замену

$\bar {y} = y - y{\kern 1pt} * = y - {{y}_{m}},$

$\bar {t} = t - t{\kern 1pt} * = t - {{t}_{m}} - \frac{{{{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}) + 1}}{\gamma },$

где $(t_{*}^{{}},y_{*}^{{}})$ – точка покоя линеаризованной системы, координаты которой имеют значения

$y{\kern 1pt} * = {{y}_{m}},\quad t{\kern 1pt} * = \frac{{\gamma {{t}_{m}} - {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}) - 1}}{\gamma }.$

Получим новую линеаризованную систему уравнений

(23)

$\begin{gathered} \frac{{d\bar {y}}}{{d\mu }} = \frac{{a\exp (\gamma {{t}_{m}})}}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}}\bar {y} + \frac{{a\gamma {{y}_{m}}\exp (\gamma {{t}_{m}})}}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}}\bar {t}, \\ \frac{{d\bar {t}}}{{d\mu }} = \frac{{ - {{a}^{2}}{{y}_{m}}\exp (3\gamma {{t}_{m}})}}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}}\bar {y} + \frac{{\gamma \exp (\gamma {{t}_{m}})}}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}}\bar {t}. \\ \end{gathered} $

Поскольку линеаризованная система (23) получена использованием линейной замены, и правая часть системы уравнений задачи (19) ограничена, то условие абсолютной устойчивости разностной схемы (22) будет эквивалентно условию абсолютной устойчивости для линеаризованной системы (23).

Введем обозначения

${{\alpha }_{1}} = \frac{{a\exp (\gamma {{t}_{m}})}}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}},\quad {{\alpha }_{2}} = \frac{{a\gamma {{y}_{m}}\exp (\gamma {{t}_{m}})}}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}},$

${{\beta }_{1}} = \frac{{ - {{a}^{2}}{{y}_{m}}\exp (3\gamma {{t}_{m}})}}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}},\quad {{\beta }_{2}} = \frac{{\gamma \exp (\gamma {{t}_{m}})}}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}}.$

Тогда явная схема метода Эйлера для преобразованной системы (23) имеет вид

(24)

${{\bar {y}}_{{m + 1}}} = {{\bar {y}}_{m}} + {{h}_{\mu }}\left( {{{\alpha }_{1}}{{{\bar {y}}}_{m}} + {{\alpha }_{2}}{{{\bar {t}}}_{m}}} \right),\quad {{\bar {t}}_{{m + 1}}} = {{\bar {t}}_{m}} + {{h}_{\mu }}\left( {{{\beta }_{1}}{{{\bar {y}}}_{m}} + {{\beta }_{2}}{{{\bar {t}}}_{m}}} \right).$

Найдем собственные значения схемы (24) , рассмотрев характеристическое уравнение

$\left| {\begin{array}{*{20}{c}} {1 + {{\alpha }_{1}}{{h}_{\mu }} - \xi }&{{{\alpha }_{2}}{{h}_{\mu }}} \\ {{{\beta }_{1}}{{h}_{\mu }}}&{1 + {{\beta }_{2}}{{h}_{\mu }} - \xi } \end{array}} \right| = 0.$

Получаем следующие собственные значения:

${{\xi }_{{1,2}}} = 1 + \frac{1}{2}\frac{{\exp (\gamma {{t}_{m}})}}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}}{{h}_{\mu }}\left( {a + \gamma \pm \sqrt {{{{(a - \gamma )}}^{2}} - 4{{a}^{3}}\gamma y_{m}^{2}\exp (2\gamma {{t}_{m}})} } \right).$

Выражение для собственных значений зависит от значений двух параметров $a$ и $\gamma $. Проведение анализа на большее собственное значение может оказаться трудоемким, поскольку оно зависит от знаков входящих параметров.

Ограничимся случаем, когда параметры $a$ и $\gamma $ удовлетворяют условию (20), гарантирующему неотрицательность подкоренных выражений в ${{\xi }_{{1,2}}}$. Тогда наибольшее собственное значение целиком зависит от последнего множителя. Тогда обозначим

(25)

${{D}_{1}} = a + \gamma + \sqrt {{{{(a - \gamma )}}^{2}} - 4{{a}^{3}}\gamma y_{m}^{2}\exp (2\gamma {{t}_{m}})} ,\quad {{D}_{2}} = a + \gamma - \sqrt {{{{(a - \gamma )}}^{2}} - 4{{a}^{3}}\gamma y_{m}^{2}\exp (2\gamma {{t}_{m}})} .$

Вычислим максимальное значение

(26)

Таким образом, получаем, что абсолютная устойчивость разностной схемы явного метода Эйлера для преобразованной модифицированным наилучшим аргументом $\mu $ вида (5) системы (19) возможна только при выполнении неравенства

(27)

$\left| {1 + \frac{1}{2}\frac{{\exp (\gamma {{t}_{m}})}}{{{{{(1 + {{a}^{2}}y_{m}^{2}\exp (2\gamma {{t}_{m}}))}}^{{3/2}}}}}{{h}_{\mu }}{{D}_{{\max }}}} \right| \leqslant 1.$

В стоящем под модулем выражении знак может менять только ${{D}_{{\max }}}$. При значениях шага ${{h}_{\mu }} > 0$ абсолютная устойчивость схемы (24) возможна только при значениях параметра ${{D}_{{\max }}} < 0$. Разрешая неравенство (27) относительно шага ${{h}_{\mu }}$, полагая ${{D}_{{\max }}} = - \left| {{{D}_{{\max }}}} \right|$, получим

$0 < {{h}_{\mu }} \leqslant \frac{{4{{{\left( {1 + {{a}^{2}}y_{m}^{2}\exp 2\gamma {{t}_{m}}} \right)}}^{{3/2}}}}}{{\left| {{{D}_{{\max }}}} \right|\exp \gamma {{t}_{m}}}}.$

При отрицательных значениях шага ${{h}_{\mu }}$ абсолютная устойчивость схемы (24) возможна только при значениях параметра ${{D}_{{\max }}} > 0$. Разрешая неравенство (27) относительно шага ${{h}_{\mu }}$, полагая ${{D}_{{\max }}} = \left| {{{D}_{{\max }}}} \right|$, получим

$ - \frac{{4{{{\left( {1 + {{a}^{2}}y_{m}^{2}\exp 2\gamma {{t}_{m}}} \right)}}^{{3/2}}}}}{{\left| {{{D}_{{\max }}}} \right|\exp \gamma {{t}_{m}}}} \leqslant {{h}_{\mu }} < 0.$

Объединяя два этих неравенства, получим условие (21).

В тривиальных случаях, при ${{h}_{\mu }} = 0$ или ${{D}_{{\max }}} = 0$ (т.е. $a = \gamma = 0$) абсолютная устойчивость также имеет место. В обоих случаях максимальное собственное значение ${{\xi }_{{\max }}} = 1$. Однако при практических расчетах случай ${{h}_{\mu }} = 0$ исключается.

Сформулируем некоторые следствия из теоремы 3.

Следствие 1. При значении параметра $\gamma = 0$ экспоненциальный аргумент продолжения решения (5) становится наилучшим аргументом, принимая вид (3), а условия устойчивости разностных схем явного метода Эйлера для задач (19) и (10) будут совпадать.

Доказательство. Очевидно, что при $\gamma = 0$ задачи (19) и (10) совпадают. Значит, должны совпадать и условия абсолютной устойчивости разностной схемы явного метода Эйлера для этих задач. Докажем это.

При $\gamma = 0$ условие (21) примет вид

$\left| {{{h}_{\mu }}} \right| \leqslant \frac{{4{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}}{{\left| {{{D}_{{\max }}}} \right|}}.$

Учитывая выражения для ${{D}_{{\max }}}$ (26) и ${{D}_{1}},\;{{D}_{2}}$ (25), при $\gamma = 0$ получим

$\left| {{{D}_{{\max }}}} \right| = 2\left| a \right|,$

что доказывает следствие 1.

Замечание 5. Следствие 1 устанавливает связь между наилучшим аргументом и его модификацией. Экспоненциальный наилучший аргумент (5) является обобщением наилучшего аргумента (3). Экспоненциальный наилучший аргумент нацелен на повышение эффективности численного решения жестких начальных задач. Варьируя значение параметра $\gamma $, можно добиться увеличения области абсолютной устойчивости разностной схемы явного метода Эйлера для преобразованной задачи.

Следствие 2. При значении параметра $\gamma = - a$ условие устойчивости разностной схемы явного метода Эйлера для задачи $(19)$ принимает вид

(28)

$\left| {{{h}_{\mu }}} \right| \leqslant \frac{{2\left[ {\exp (a{{t}_{m}}) + {{a}^{2}}y_{m}^{2}\exp ( - a{{t}_{m}})} \right]}}{{\left| a \right|}}.$

Доказательство. При $\gamma = - a$ условие (21) примет вид

(29)

$\left| {{{h}_{\mu }}} \right| \leqslant \frac{{4{{{\left( {1 + {{a}^{2}}y_{m}^{2}\exp ( - 2a{{t}_{m}})} \right)}}^{{3/2}}}}}{{\left| {{{D}_{{\max }}}} \right|\exp ( - a{{t}_{m}})}}.$

Вычислим значение ${\text{|}}{{D}_{{\max }}}{\text{|}}$. При $\gamma = - a$ оно примет вид

$\left| {{{D}_{{\max }}}} \right| = {{D}_{1}} = \sqrt {{{{(a - \gamma )}}^{2}} - 4{{a}^{3}}\gamma y_{m}^{2}\exp (2\gamma {{t}_{m}})} = 2\left| a \right|\sqrt {1 + {{a}^{2}}y_{m}^{2}\exp ( - 2a{{t}_{m}})} .$

Подставляя найденное выражение для $\left| {{{D}_{{\max }}}} \right|$ в неравенство (29), перепишем его в виде

$\left| {{{h}_{\mu }}} \right| \leqslant \frac{{4{{{\left( {1 + {{a}^{2}}y_{m}^{2}\exp ( - 2a{{t}_{m}})} \right)}}^{{3/2}}}}}{{2\left| a \right|\sqrt {1 + {{a}^{2}}y_{m}^{2}\exp ( - 2a{{t}_{m}})} \exp ( - a{{t}_{m}})}},$

упростив правую часть неравенства, получим условие (28).

Замечание 6. Условие абсолютной устойчивости (28) зависит от значения параметра $a$. Рассмотрим функцию

$F(a) = \frac{{\left[ {\exp (a{{t}_{m}}) + {{a}^{2}}y_{m}^{2}\exp ( - a{{t}_{m}})} \right]}}{{{{{\left( {1 + {{a}^{2}}y_{m}^{2}} \right)}}^{{3/2}}}}},$

которая получена как частное правых частей неравенств (28) и (11). Функция $F(a)$ ограничена снизу, положительная и удовлетворяет условию

$\mathop {\lim }\limits_{a \to \pm \infty } F(a) = + \infty .$

Это означает, что существуют такие значения $a_{ - }^{*} < 0$ и $a_{ + }^{*} > 0$, для которых справедливо неравенство

$F(a) > 1\quad \forall a \in [a_{ + }^{*}, + \infty )\quad {\text{при}}\quad {\kern 1pt} a \geqslant 0\quad {\text{и}}\quad \forall a \in ( - \infty ;a_{ - }^{*}]\quad {\text{при}}\quad {\kern 1pt} a \leqslant 0,$

т.е. условие абсолютной устойчивости (28) при указанных условиях определяет большую область абсолютной устойчивости, чем при условии (11). Проанализируем этот факт.

При значениях $a_{ - }^{*}$ и $a_{ + }^{*}$, равных нулю, согласно следствию 1, $F(a) = 0$. Однако в окрестности нуля функция $F(a)$ может вести себя следующим образом: сначала возрастать до некоторого значения больше единицы, затем убывать до значения меньше единицы, после чего монотонно возрастать. Поэтому, как правило, значения $a_{ - }^{*}$ и $a_{ + }^{*}$ полагаются отличными от нуля. В зависимости от рассматриваемой точки $M({{t}_{m}},{{y}_{m}})$ значения $a_{ - }^{*}$ и $a_{ + }^{*}$ могут быть как умеренными, так и большими по модулю.

Следовательно, выбор $\gamma = - a$ хорошо подходит для решения жестких и сверхжестких начальных задач и менее пригоден для начальных задач с малыми значениями параметра $a$. Однако при больших значениях параметра $a$ возникают вычислительные трудности с проведениями расчетов, поскольку значения $\exp (a{{t}_{m}})$ могут достигать очень больших значений.

4. ЭКСПОНЕНЦИАЛЬНЫЙ ТЕСТ

В статье [37] на примере экспоненциального теста, предложенного Н.Н. Калиткиным и А.А. Беловым в [41], показаны преимущества использования наилучшего и модифицированного наилучшего аргументов. Рассмотрим экспоненциальную тестовую задачу

(30)

$\frac{{du}}{{dt}} = - \xi (t){\kern 1pt} u{\kern 1pt} \left( {{{u}^{2}} - {{a}^{2}}} \right),\quad u(0) = {{u}_{0}},\quad t \in \left[ {0,2\pi } \right],$

имеющую точное решение

(31)

$u(t) = \frac{{a{{u}_{0}}}}{{\sqrt {u_{0}^{2} + \left( {{{a}^{2}} - u_{0}^{2}} \right){\kern 1pt} \exp \left( { - 2{{a}^{2}}\Xi (t)} \right)} }},$

где $\xi (t) = {{\xi }_{0}}\cos (t)$, $\Xi (t) = {{\xi }_{0}}\sin (t)$, ${{u}_{0}} = 0.5$, $a = \pi $. Жесткость задачи регулируется параметром ${{\xi }_{0}}$: чем больше ${{\xi }_{0}}$, тем уже зоны пограничных и внутренних слоев (см. фиг. 1).

Фиг. 1.

Точное решение (31) экспоненциального теста (30) для ${{\xi }_{0}} = 1$, ${{\xi }_{0}} = 50$.

Проведем анализ результатов численного решения исходной задачи (30), а также задачи, преобразованной к наилучшему аргументу $\lambda $ вида

$d{{\lambda }^{2}} = d{{u}^{2}} + d{{t}^{2}},$

имеющей вид

(32)

$\begin{gathered} \frac{{du}}{{d\lambda }} = \frac{{ - \xi (t){\kern 1pt} u{\kern 1pt} \left( {{{u}^{2}} - {{a}^{2}}} \right)}}{{\sqrt {1 + {{\xi }^{2}}(t){\kern 1pt} {{u}^{2}}{\kern 1pt} {{{\left( {{{u}^{2}} - {{a}^{2}}} \right)}}^{2}}} }},\quad u(0) = {{u}_{0}}, \\ \frac{{dt}}{{d\lambda }} = \frac{1}{{\sqrt {1 + {{\xi }^{2}}(t){\kern 1pt} {{u}^{2}}{\kern 1pt} {{{\left( {{{u}^{2}} - {{a}^{2}}} \right)}}^{2}}} }}\quad t(0) = 0, \\ \end{gathered} $

и задачи

(33)

$\begin{gathered} \frac{{du}}{{d\mu }} = \frac{{ - \xi (t){\kern 1pt} u{\kern 1pt} \left( {{{u}^{2}} - {{a}^{2}}} \right){\kern 1pt} \exp \left( {\gamma t} \right)}}{{\sqrt {1 + {{\xi }^{2}}(t){\kern 1pt} {{u}^{2}}{\kern 1pt} {{{\left( {{{u}^{2}} - {{a}^{2}}} \right)}}^{2}}{\kern 1pt} \exp \left( {2\gamma t} \right)} }},\quad u(0) = {{u}_{0}}, \\ \frac{{dt}}{{d\mu }} = \frac{{\exp (\gamma t)}}{{\sqrt {1 + {{\xi }^{2}}(t){\kern 1pt} {{u}^{2}}{\kern 1pt} {{{\left( {{{u}^{2}} - {{a}^{2}}} \right)}}^{2}}{\kern 1pt} \exp \left( {2\gamma t} \right)} }},\quad t(0) = 0, \\ \end{gathered} $

преобразованной к модифицированному наилучшему аргументу $\mu $:

$d{{\mu }^{2}} = d{{u}^{2}} + \exp ( - 2\gamma t){\kern 1pt} d{{t}^{2}}.$

В [37] были представлены численные расчеты задач (30), (32) и (33) явным методом Эйлера с переменным шагом интегрирования, вычисляемым по правилу Рунге–Ромберга–Ричардсона. Эти результаты приведены в табл. 1, где $\theta $ – значение локальной погрешности. Начальный шаг имеет значение ${{h}_{0}}{{ = 10}^{{ - 5}}}$, а параметр модифицированного наилучшего аргумента $\gamma = - 100$.

Таблица 1.

Средняя ошибка и время счета явного метода Эйлера решения задач (30), (32) и (33) с переменным шагом интегрирования

${{\xi }_{0}}$	$\theta $	Исходная задача		Наилучший аргумент		Модифицированный наилучший аргумент
${{\xi }_{0}}$	$\theta $	${{\varepsilon }_{{{\text{a}}vg}}}$	$t,c$	${{\varepsilon }_{{{\text{a}}vg}}}$	$t,c$	${{\varepsilon }_{{{\text{a}}vg}}}$	$t,c$
1	10^–4	0.0072	0.01	0.02	0.02	$1.1 \times {{10}^{{ - 4}}}$	0.33
	10^–5	0.02	0.05	0.02	0.03	$2 \times {{10}^{{ - 5}}}$	0.8
	10^–6	0.006	0.11	0.01	0.18	$3.6 \times {{10}^{{ - 6}}}$	2.5
	10^–7	$8.5 \times {{10}^{{ - 4}}}$	0.3	0.002	0.4	$6.3 \times {{10}^{{ - 7}}}$	8.2
	10^–8	$1.4 \times {{10}^{{ - 4}}}$	0.8	$2.3 \times {{10}^{{ - 4}}}$	1.25	$1.1 \times {{10}^{{ - 7}}}$	25
10	10^–6	–	–	–	–	$7.9 \times {{10}^{{ - 4}}}$	2.3
	10^–7	–	–	–	–	$8.5 \times {{10}^{{ - 4}}}$	7.2
	10^–8	–	–	–	–	$2.1 \times {{10}^{{ - 4}}}$	22
50	10^–6	–	–	–	–	0.006	2
	10^–7	–	–	–	–	0.003	6.6
	10^–8	–	–	–	–	$1.6 \times {{10}^{{ - 4}}}$	22

Проведем анализ изменения шага интегрирования в процессе численного решения задач (30), (32) и (33) явным методом Эйлера.

На фиг. 2б изображены значения шага интегрирования в каждой точке решения задач (30), (32) и (33). Для исходной задачи (30) и задачи, преобразованной к наилучшему аргументу (32), наблюдается уменьшение шага интегрирования при прохождении пограничных и внутренних слоев. При решении задачи (33), преобразованной к модифицированному наилучшему аргументу, шаг практически не уменьшается на протяжении всего процесса вычисления, только приостанавливая свое увеличение при прохождении участков быстрого изменения интегральной кривой. Отметим, что параметр модифицированного наилучшего аргумента здесь имеет значение $\gamma = - 10$, в отличие от расчетов, приведенных в табл. 1. Это сделано, чтобы более наглядно представить полученные результаты. Согласно полученным авторами результатам, при $\gamma = - 100$ шаг интегрирования задачи (33) увеличивается еще значительнее, что делает графики изменения шага интегрирования для задач (30), (32) менее информативными.

Фиг. 2.

(а) –Численное решение задач (30), (32) и (33) для ${{\xi }_{0}} = 1$ явным методом Эйлера с переменным шагом, рассчитанным по правилу Рунге–Ромберга–Ричардсона с допустимой погрешностью $\theta {{ = 10}^{{ - 8}}}$, начальным шагом ${{h}_{0}}{{ = 10}^{{ - 5}}}$ и параметром модифицированного наилучшего аргумента $\gamma = - 10$. (б) – Соответствующие значения шага интегрирования в каждой точке решения задач (30), (32) и (33).

Также стоит обратить внимание: на фиг. 2а видно, что решение задачи, преобразованной к модифицированному наилучшему аргументу, наиболее близко к точному решению (31). И хотя численное решение заняло больше времени – 22 с, чем при решении задач (30) – 6 с и (32) – 4 с, его точность на несколько порядков выше. Увеличение времени счета для задачи (33) связано с необходимостью вычисления экспоненты с большими значениями аргументов.

Таким образом, преобразование исходной задачи к наилучшему и модифицированному наилучшему аргументам позволяет увеличить шаг интегрирования при использовании явного метода Эйлера, что соответствует полученным теоретическим результатам. При этом применение модифицированного наилучшего аргумента в некоторых случаях практически освобождает от необходимости подбора начального шага при использовании явного метода Эйлера с переменным шагом.

5. ЗАКЛЮЧЕНИЕ

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

В настоящей работе исследовалась область абсолютной устойчивости явного численного метода Эйлера при решении задач, преобразованных к наилучшему аргументу и его модификации. Исследование проводилось на тестовой начальной задаче Далквиста, так как она моделирует локальное поведение интегральной кривой в окрестности каждой точки линеаризованной решаемой задачи.

Теорема об увеличении области абсолютной устойчивости явного метода Эйлера для задачи, преобразованной к наилучшему аргументу, была сформулирована и доказана ранее в работах Е.Б. Кузнецова и В.И. Шалашилина. Однако в их доказательстве были сделаны следующие допущения, которые были устранены в текущей работе: шаг интегрирования принимался исключительно больше нуля и последнее уравнение преобразованной системы не рассматривалось, так как считалось, что оно не оказывает существенного влияния на размер области абсолютной устойчивости.

Для экспоненциального модифицированного наилучшего аргумента подобного исследования ранее не проводилось. Полученные теоретические результаты показывают, что преобразование к такому аргументу позволяет расширить область абсолютной устойчивости явного метода Эйлера при корректном выборе параметра $\gamma $. Однако полученные теоретические результаты трудно применять при выборе значения параметра $\gamma $ на практике. Разработка более удобного способа оценки значения параметра $\gamma $, входящего в экспоненциальный аргумент, является направлением дальнейших исследований.

Список литературы

Чанг К., Хауэс Ф. Нелинейные сингулярно возмущенные краевые задачи: теория и приложения. М.: Мир, 1988. 247 с.
Прандтль Л. Теория несущего крыла. Ч. 1. Движение жидкости с очень малым трением. М.; Л.: ГНТИ, 1931. С. 5–11.
Тихонов А.Н. Системы дифференциальных уравнений, содержащие малые параметры при производных // Матем. сб. 1952. Т. 31 (73). № 3. С. 575–586.
Васильева А.Б. О дифференцировании решений систем дифференциальных уравнений, содержащих малый параметр // Матем. сб. 1951. Т. 28 (70). № 1. С. 131–146.
Бутузов В.Ф., Васильева А.Б., Нефедов Н.Н. Асимптотическая теория контрастных структур (обзор) // Автомат. и телемех. 1997. № 7. С. 4–32.
Бутузов В.Ф., Левашова Н.Т., Мельникова А.А. Контрастная структура типа ступеньки в сингулярно возмущенной системе эллиптических уравнений // Ж. вычисл. матем. и матем. физ. 2013. Т. 53. № 9. С. 1427–1447.
Бутузов В.Ф. О контрастных структурах с многозонным внутренним слоем // Моделирование и анализ информационных систем. 2017. Т. 24. № 3. С. 288–308.
Curtiss C.F., Hirschfelder J.O. Integration of stiff equations // Proc. Nat. Acad. Sci. USA. 1952. V. 38. P. 235–243.
Bjurel G., Dahlquist G., Lindberg B., Linde S., Oden L. Survey of stiff ordinary differential equations // Royal Institute of Technology, Stockholm, Depart. Informat. Proces., Comput. Sci. Rep. NA 70.11. 1970.
Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М.: Мир, 1999.
Деккер К., Вервер Я. Устойчивость методов Рунге–Кутты для жестких нелинейных дифференциальных уравнений. М.: Мир, 1988.
Gear C.W. Numerical initial value problems in ordinary differential equations. Prentice-Hall, Upper Saddle River, 1971.
Белов А.А., Калиткин Н.Н. Проблема нелинейности при численном решении сверхжестких задач Коши // Матем. моделирование. 2016. Т. 28. № 4. С. 16–32. (Перевод: Belov A.A., Kalitkin N.N. Nonlinearity problem in the numerical solution of super stiff Cauchy problems // Math. Model. Comput. Simulat. 2013. V. 8. Iss. 6. P. 638–650).
Калиткин Н.Н., Пошивайло И.П. Вычисления с использованием обратных схем Рунге–Кутты // Матем. моделирование. 2013. Т. 25. № 10. С. 79–96 (Перевод: Kalitkin N.N., Poshivaylo I.P. Computations with inverse Runge–Kutta schemes // Math. Model. Comput. Simulat. 2013. V. 6. Iss. 3. P. 272–285).
Куликов Г.Ю. Вложенные симметричные неявные гнездовые методы Рунге–Кутты типов Гаусса и Лобатто для решения жестких обыкновенных дифференциальных уравнений и гамильтоновых систем // Ж. вычисл. матем. и матем. физ. 2015. Т. 55. № 6. С. 986–1007 (Перевод: 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. Iss. 6. P. 983–1003).
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. No. 3. P. A1593–A1613.
Kulikov G.Yu., Shindin S.K. Adaptive nested implicit Runge–Kutta formulas of Gauss type // Appl. Numeric. Math. 2009. V. 59. № 3–4. P. 707–722.
Куликов Г.Ю. Об устойчивости симметричных формул Рунге–Кутты // Докл. АН. 2003. Т. 389. № 2. С. 164–168 (Перевод: Kulikov G.Yu. On the stability of symmetric Runge–Kutta formulas // Doklady. Mathematics. 2003. V. 67. Iss. 2. P. 184–188).
Kulikov G.Yu. Symmetric Runge–Kutta methods and their stability // Russian J. Numer. Anal. Math. Model. 2003. V. 18. Iss. 1. P. 13–41.
Куликов Г.Ю. Неявные гнездовые методы Рунге–Кутты типа Гаусса и Лобатто с локальным и глобальным контролем точности для жестких обыкновенных дифференциальных уравнений // Ж. вычисл. матем. и матем. физ. 2020. Т. 60. № 7. С. 1170–1192 (Перевод: Kulikov G.Yu. Nested implicit Runge–Kutta pairs of Gauss and Lobatto types with local and global error controls for stiff ordinary differential equations // Comput. Math. Math. Phys. 2020. V. 60. Iss. 7. P. 1134–1154).
Скворцов Л.М. О неявных методах Рунге–Кутты, полученных в результате обращения явных методов // Матем. моделирование. 2017. Т. 29. № 1. С. 3–19 (Перевод: Skvortsov L.M. On implicit Runge–Kutta methods received as a result of inversion of explicit methods // Math. Model. Comput. Simulat. 2017 V. 9. Iss. 4. P. 498–510).
Скворцов Л.М., Козлов О.С. Эффективная реализация диагонально-неявных методов Рунге–Кутты // Матем. моделирование. 2014. Т. 26. № 1. С. 96–108 (Перевод: Skvortsov L.M., Kozlov O.S. Efficient implementation of diagonally implicit Runge–Kutta methods // Math. Model. Comput. Simulat. 2014 V. 6. Iss. 4. P. 415–424).
Скворцов Л.М. Однократно неявные диагонально расширенные методы Рунге–Кутты четвертого порядка // Ж. вычисл. матем. и матем. физ. 2014. Т. 54. № 5. С. 755–765 (Перевод: Skvortsov L.M. Singly implicit diagonally extended Runge–Kutta methods of fourth order // Comput. Math. Math. Phys. 2014. V. 54. Iss. 5. P. 775–784).
Скворцов Л.М. Эффективная реализация неявных методов Рунге–Кутты второго порядка // Матем. моделирование. 2013. Т. 25. № 5. С. 15–28 (Перевод: Skvortsov L.M. Efficient implementation of second order implicit Runge–Kutta methods // Math. Model. Comput. Simulat. 2013. V. 5. Iss. 6. P. 565–574).
Новиков Е.А. Явные методы для жестких систем. Новосибирск: Наука. Сибирское предприятие РАН, 1997.
Новиков Е.А., Шорников Ю.В. Моделирование жестких гибридных систем. СПб: Изд-во “Лань”, 2019.
Novikov E.A., Rybkov M.V. Application of explicit methods with extended stability regions for solving stiff problems // J. Siberian Federal University. Math. Phys. 2016. V. 9. No. 2. P. 209–219.
Новиков Е.А. Алгоритм интегрирования жестких задач с помощью явных и неявных методов // Изв. Сарат. ун-та. Новая серия. Сер. “Математика. Механика. Информатика”. 2012. Т. 12. № 4. С. 19–27.
Новиков А.Е., Новиков Е.А. Численное решение жестких задач с небольшой точностью // Матем. моделирование. 2010. Т. 22. № 1. С. 46–56 (Перевод: Novikov A.E., Novikov E.A. Numerical integration of stiff systems with low accuracy // Math. Model. Comput. Simulat. 2010. V. 2. Iss. 4. P. 443–452).
Скворцов Л.М. Явные адаптивные методы Рунге–Кутты для жестких и колебательных задач // Ж. вычисл. матем. и матем. физ. 2011. Т. 51. № 8. С. 1434–1448 (Перевод: Skvortsov L.M. Explicit adaptive Runge–Kutta methods for stiff and oscillation problems // Comput. Math. Math. Phys. 2011. V. 51. Iss. 8. P. 1339–1352).
Скворцов Л.М. Явные адаптивные методы Рунге–Кутты // Матем. моделирование. 2011. Т. 23. № 7. С. 73–87 (Перевод: Skvortsov L.M. Explicit adaptive Runge–Kutta methods // Math. Model. Comput. Simulat. 2011. V. 4. Iss. 1. P. 82–91).
Лебедев В.И. Как решать явными методами жесткие системы дифференциальных уравнений // Вычисл. процессы и системы. Вып. 8. М: Наука, 1991.
Лебедев В.И. Явные разностные схемы для решения жестких задач с комплексным или разделимым спектром // Ж. вычисл. матем и матем. физ. 2000. Т. 40. № 12. С. 1801–1812.
Шалашилин В.И., Кузнецов Е.Б. Метод продолжения решения по параметру и наилучшая параметризация в прикладной математике и механике. М.: Эдиториал УРСС, 1999. 224 с. (Shalashilin V.I., Kuznetsov E.B. Parametric continuation and optimal parametrization in applied mathematics and mechanics. Kluwer Acad. Publ. Dordrecht, Boston, London, 2003. 236 p.)
Kulikov G.Yu. On quasi-consistent integration by Nordsieck methods // J. Comput. Appl. Math. 2009. V. 225. Iss. 1. P. 268–287.
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. Iss. 9. P. 2351–2364.
Kuznetsov E.B., Leonov S.S., Tsapko E.D. A new numerical approach for solving initial value problems with exponential growth integral curves // IOP Conf. Ser.: Mater. Sci. Eng. 2020. V. 927.
Кузнецов Е.Б., Шалашилин В.И. Задача Коши как задача продолжения решения по параметру // Ж. вычисл. матем. и матем. физ. 1993. Т. 33. № 12. С. 1792–1805.
Dahlquist G. A special stability problem for linear multistep methods // BIT. 1963. № 3. P. 27–43.
Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989.
Белов А.А., Калиткин Н.Н. Особенности расчета контрастных структур в задачах Коши // Матем. моделирование. 2016. Т. 28. № 10. С. 97–109.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал