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

2.1. Модельная задача 1

В качестве первой модельной задачи выбрана задача об одномерном движении материальной точки единичной массы в поле кубического потенциала ${{U}_{C}}(q) = {{q}^{3}}{\text{/}}3 - {{q}^{2}}{\text{/}}2$ [9], [11]–[14]. В данном случае число степеней свободы равно 1, импульс материальной точки $p(t)$ и координата $q(t)$ являются скалярами. Гамильтониан $H(p,q) = {{p}^{2}}{\text{/}}2 + {{U}_{C}}(q)$.

В настоящей работе используются только периодические решения модельной задачи 1, которые существуют при $0 < {{q}_{0}} < 1.5$, ${\text{|}}{{p}_{0}}{\text{|}} < {{q}_{0}}{{(1 - 2{{q}_{0}}{\text{/}}3)}^{{0.5}}}$. Период решения зависит от начального значения гамильтониана ${{H}_{0}}$ и вычисляется по формуле

где ${{q}_{{{\text{min}}}}}({{H}_{0}})$, ${{q}_{{{\text{max}}}}}({{H}_{0}})$ – координаты точек возврата, которые являются корнями уравнения ${{H}_{0}} - U(q) = 0$ и определяются следующим образом: В этих точках $p = 0$. В малой окрестности минимума потенциала ($q = 1$) движения приближенно описываются гармоническими колебаниями с периодом ${{T}_{{\text{Г}}}} = 2\pi $ [14].

У периодического решения задачи для каждого $q \in [{{q}_{{{\text{min}}}}},{{q}_{{{\text{max}}}}}]$ значение импульса вычисляется по формуле $p = \pm {{\left( {2{{H}_{0}} - 2{{q}^{3}}{\text{/}}3 + {{q}^{2}}} \right)}^{{0.5}}}$.

2.2. Модельная задача 2

В качестве второй модельной задачи выбрана задача о плоском движении двух тел единичной массы (задача Кеплера) с гамильтонианом [3, с. 9]

Здесь ${\mathbf{p}} = {{({{p}_{1}},{{p}_{2}})}^{{\text{T}}}}$ – импульс второго тела, ${\mathbf{q}} = {{({{q}_{1}},{{q}_{2}})}^{{\text{T}}}}$ – координаты второго тела относительно первого. Эта задача имеет также другой первый интеграл (помимо $H({\mathbf{p}},{\mathbf{q}}))$: угловой момент $L({\mathbf{p}},{\mathbf{q}}) = {{q}_{1}}{{p}_{2}} - {{q}_{2}}{{p}_{1}}$.

Для рассматриваемой задачи известно аналитическое решение [15, с. 56], [3], [5, с. 49]. Возможны три случая: эллиптическая, параболическая и гиперболическая орбиты. Для эллиптической орбиты формулы аналитического решения имеют следующий вид:

где $a$ – большая полуось орбиты, $e \in [0,1)$ – эксцентриситет орбиты, $E$ – эксцентрическая аномалия. $E$ находится как решение уравнения Кеплера где ${{t}_{{\text{П}}}}$ – момент времени, в который второе тело проходит через перицентр.

3.1. Двухпараметрическое семейство трехстадийных симметрично-симплектических методов Рунге–Кутты

Рассмотрим двухпараметрическое семейство трехстадийных симметрично-симплектических методов Рунге–Кутты с таблицей Бутчера следующего вида [11], [16], [17]:

Здесь ${{b}_{1}}$ и ${{s}_{{12}}}$ – свободные параметры семейства, ${{b}_{1}} > 1{\text{/}}6$, ${{\tilde {c}}_{1}}({{b}_{1}}) = 0.5{{(6{{b}_{1}})}^{{ - 0.5}}}$. Методы являются неявными. Значениям ${{b}_{1}} = 5{\text{/}}18$, ${{s}_{{12}}} = s_{{12}}^{*} = 0.75 \times {{0.6}^{{0.5}}} \approx 0.580948$ соответствует единственный метод шестого порядка аппроксимации – метод Кунцмана–Бутчера [4, с. 209]. Остальным значениям параметров соответствуют методы четвертого порядка. При ${{b}_{1}} = 0.5$ независимо от значения ${{s}_{{12}}}$ методы вырождаются в двухстадийный симметрично-симплектический метод Рунге–Кутты – метод Хаммера–Холлингсуорта четвертого порядка [4, с. 207], [17].

Применительно к модельным задачам формулы методов имеют вид

где ${{{\mathbf{p}}}_{n}}$, ${{{\mathbf{q}}}_{n}}$ – векторы импульсов и координат на $n$-м шаге метода; $\tau > 0$ – величина шага метода; ${{\overline b }_{i}} = {{b}_{1}}{{a}_{{1i}}} + (1 - 2{{b}_{1}}){{a}_{{2i}}} + {{b}_{1}}{{a}_{{3i}}}$ (${{a}_{{ij}}}$ – коэффициенты таблицы Бутчера); ${{{\mathbf{Z}}}_{{in}}}$ являются решениями редуцированной системы разрешающих уравнений Здесь ${{\bar {a}}_{{ij}}} = \sum\nolimits_{k = 1}^3 \,{{a}_{{ik}}}{{a}_{{kj}}}$; ${{c}_{i}}$ – параметры из левого столбца таблицы Бутчера.

Для решения (3.3) далее в работе использована модификация метода простой итерации [18, с. 191], [19, с. 381] – метод покоординатных итераций [20, с. 321], [21, с. 285]. На $(k + 1)$-й итерации приближение ${\mathbf{Z}}_{{in}}^{{(k + 1)}}$ решения системы (3.3) вычисляется по следующим формулам:

${\mathbf{Z}}_{{in}}^{{(0)}}$ – начальные приближения.

Решение существует и единственно в некоторой окрестности $Q$ точки ${{{\mathbf{q}}}_{n}}$ при выполнении следующего условия [4, с. 206]:

Согласно проведенным исследованиям, выполняется неравенство $\mathop {max}\limits_i \sum\nolimits_{k = 1}^3 \,{\text{|}}{{\bar {a}}_{{ik}}}{\text{|}} < 1.2$ для ${{b}_{1}} = 5{\text{/}}18$, ${{s}_{{12}}} \in [ - 0.5,1.5]$ и ${{b}_{1}} = 0.5$, ${{s}_{{12}}} = 0$.

В качестве начальных приближений при расчетах, описанных далее в настоящей работе, использованы ${\mathbf{Z}}_{{in}}^{{(0)}} = \tau {{c}_{i}}{{{\mathbf{p}}}_{n}}$ [3]. Как показали исследования для модельной задачи 1, в подавляющем большинстве случаев эти начальные приближения обеспечили меньшее число итераций, чем приближения ${\mathbf{Z}}_{{in}}^{{(0)}} = 0$. Итерационный процесс останавливался [22], если

где ${\mathbf{Z}}_{n}^{{(k)}} = \mathop {\left( {Z_{{1n,1}}^{{(k)}}, \ldots ,Z_{{1n,d}}^{{(k)}},\;Z_{{2n,1}}^{{(k)}}, \ldots ,Z_{{2n,d}}^{{(k)}},\;Z_{{3n,1}}^{{(k)}}, \ldots ,Z_{{3n,d}}^{{(k)}}} \right)}\nolimits^{\text{T}} $ есть полученное на k-й итерации приближенное решение (3.3), ${\mathbf{Z}}_{n}^{{(0)}}$ – начальное приближение, ${{\varepsilon }_{{{\text{abs}}}}} \geqslant 0$ и ${{\varepsilon }_{{{\text{rel}}}}} \geqslant 0$ – параметры метода решения системы разрешающих уравнений.

Кубический потенциал является простейшим потенциалом, при котором любой метод рассматриваемого семейства не сохраняет полную энергию при фиксированных значениях свободных параметров ${{b}_{1}}$ и ${{s}_{{12}}}$. В случаях линейного и квадратичного потенциалов методы рассматриваемого семейства сохраняют полную энергию, как и все симплектические методы Рунге–Кутты [3, с. 101]. Далее будем искать связь между этими параметрами, при которой полная энергия сохраняется.

3.2. Метод нулевого дисбаланса полной энергии

В работе [11] предложено семейство консервативных вычислительных методов решения задачи Коши для гамильтоновых систем. Методы основаны на двухпараметрическом семействе методов Рунге–Кутты (3.1)–(3.3). Система разрешающих уравнений (3.3) методов Рунге–Кутты дополняется уравнением нулевого дисбаланса полной энергии

где при этом параметры семейства ${{b}_{1}}$ и ${{s}_{{12}}}$ рассматриваются как дополнительные неизвестные. Предложенные методы работоспособны, если полученная расширенная система разрешающих уравнений имеет решение. Как показали исследования, для модельной задачи 1 расширенная система имеет однопараметрическое множество решений $\{ {{b}_{1}}(r),{{s}_{{12}}}(r)\} $ для любых $({{p}_{n}},{{q}_{n}})$, лежащих на фазовых траекториях, соответствующих $H \in (0, - 1{\text{/}}6)$, и $\tau \leqslant 0.05{{T}_{{\text{Г}}}}$. Для выбора конкретного метода из предложенного семейства можно зафиксировать значение одного из неизвестных, например ${{b}_{1}} = 5{\text{/}}18$, и искать ${{s}_{{12}}}$.

Расширенная система разрешающих уравнений (3.3), (3.5) в описанных ниже расчетах решалась с помощью двух вложенных друг в друга итерационных процессов. На внешних итерациях решалось уравнение (3.5), т.е. для фиксированного допустимого значения ${{b}_{1}} = 5{\text{/}}18$ находилось значение ${{s}_{{12}}}$, обеспечивающее сохранение полной энергии. Для определения ${{{\mathbf{Z}}}_{{in}}}$ и, следовательно, ${{{\mathbf{p}}}_{{n + 1}}}$ и ${{{\mathbf{q}}}_{{n + 1}}}$ по известным ${{b}_{1}}$ и ${{s}_{{12}}}$ использован внутренний итерационный процесс, описанный в предыдущем подразделе 3.1. Таким образом, на каждом временном шаге система (3.3) при фиксированных ${{b}_{1}}$, ${{s}_{{12}}}$ могла решаться несколько раз, каждый раз – для другого ${{s}_{{12}}}$. При втором и следующих на текущем временном шаге решениях (3.3) в качестве начальных приближений использовались решения (3.3), полученные при предыдущем решении этой системы (т.е. для другого значения ${{s}_{{12}}}$).

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

или $\left| {s_{{12}}^{{(j + 1)}} - s_{{12}}^{{(j)}}} \right| \leqslant {{\varepsilon }_{s}}$, где $s_{{12}}^{{(j)}}$, $s_{{12}}^{{(j + 1)}}$ – полученные на предыдущей и текущей итерациях приближенные решения (3.5); ${{\varepsilon }_{{\Delta H}}} \geqslant 0$, ${{\varepsilon }_{s}} \geqslant 0$ – заданные параметры [24, с. 43]; $s_{{12}}^{{(j)}}$ также может являться начальным приближением. Если для одного из начальных приближений условие (3.6) выполнялось, это значение принималось в качестве решения.

3.3. Исследование эффективности методов решения уравнения нулевого дисбаланса полной энергии

Для ответа на вопрос, какой вычислительный метод решения уравнения (3.5) эффективнее (метод секущих или метод Мюллера), проведены две серии вычислительных экспериментов с модельной задачей 1. В первой серии использовались числа с плавающей запятой двойной точности (длина мантиссы 53 бит), во второй серии – числа с плавающей запятой повышенной точности (длина мантиссы 448 бита). Для работы с числами повышенной точности применялась программная библиотека GNU MPFR [25]. Расчеты проводились для множества начальных значений ${{p}_{0}} = 0$, ${{q}_{0}} \in {{Q}_{0}}$, величин шага $\tau \in \Theta $ и отрезков интегрирования длиной $10T$. Период точного решения $T$ разный для разных начальных значений гамильтониана ${{H}_{0}}$. Здесь и далее

Для метода секущих использовались следующие начальные приближения: $s_{{12}}^{{(1)}} = s_{{12}}^{*}$, $s_{{12}}^{{(2)}} = s_{{12}}^{*} + 4 \times {{10}^{{ - 4}}}$, для метода Мюллера: $s_{{12}}^{{(1)}} = s_{{12}}^{*}$, $s_{{12}}^{{(2)}} = s_{{12}}^{*} + 4 \times {{10}^{{ - 4}}}$, $s_{{12}}^{{(3)}} = 0.5\left( {s_{{12}}^{{(1)}} + s_{{12}}^{{(2)}}} \right)$.

Найдены $ma{{x}_{{{{p}_{0}} = 0,{{q}_{0}} \in {{Q}_{0}},\tau \in \Theta }}}{{N}_{{{\text{внешн}},{\text{ср}}}}}$ и $ma{{x}_{{{{p}_{0}} = 0,{{q}_{0}} \in {{Q}_{0}},\tau \in \Theta }}}{{N}_{{{\text{внутр}},{\text{ср}}}}}$. Здесь ${{N}_{{{\text{внешн}},{\text{ср}}}}}$ – среднее на один временной шаг число внешних итераций, ${{N}_{{{\text{внутр}},{\text{ср}}}}} = \sum\nolimits_{n = 1}^N \,{{N}_{{{\text{внутр}},n}}}{\text{/}}(N)$ – среднее на один временной шаг  число  внутренних  итераций, ${{N}_{{{\text{внутр}},n}}}$ –  суммарное число внутренних итераций на $n$-м шаге по времени (учитывались все внешние итерации), N – число шагов по времени. Результаты расчетов представлены в табл. 1. Равенство нулю ${{\varepsilon }_{{\Delta H}}}$ означает, что условие остановки внешних итераций (3.6) не использовалось (за исключением случаев точного равенства $\Delta H\left( {s_{{12}}^{{(j)}}} \right) = 0$).

Согласно представленным результатам, величина $ma{{x}_{{{{p}_{0}} = 0,{{q}_{0}} \in {{Q}_{0}},\tau \in \Theta }}}{{N}_{{{\text{внутр}},{\text{ср}}}}}$ немного меньше для метода секущих в случае двойной точности и меньше для метода Мюллера в случае повышенной точности. Поэтому для решения уравнения (3.5) при вычислениях, рассматриваемых в разд. 6, 7, применялся метод Мюллера с указанными выше начальными приближениями.

6.1. Симплектичность

Для исследования симплектичности рассматриваемых методов введем матрицы глобального дефекта симплектичности

и локального дефекта симплектичности Матрицы $\Delta {{{\mathbf{S}}}_{n}}$ и $\delta {{{\mathbf{S}}}_{n}}$ являются кососимметричными матрицами. В случае модельной задачи 1 достаточно проанализировать только зависимости элементов от времени на приближенных решениях. Близость к нулю этих элементов характеризует реальную симплектичность метода.

Для вычисления частных производных использовались правые разностные производные [18, с. 34].

Приведем результаты расчетов повышенной точности. Приращение аргумента в формуле правой разностной производной в этом случае равно ${{10}^{{ - 30}}}$. Исследования проводились на малых промежутках времени, сравнимых с периодом, и на больших промежутках времени длиной $1000T$.

Метод нулевого дисбаланса полной энергии. Рассмотрим сначала $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ для этого метода на малых промежутках времени. На фиг. 1 изображены графики зависимостей $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ от времени при различных $({{p}_{0}},{{q}_{0}})$ и $\tau = 0.05{{T}_{{\text{Г}}}}$. Согласно графикам, эти зависимости $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ являются почти периодическими.

Теперь рассмотрим $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ на больших отрезках времени. В табл. 3 представлены максимальные абсолютные значения $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ на отрезках $t \in [0,1000T]$ для разных начальных значений $({{p}_{0}},{{q}_{0}})$ и величин шага $\tau $. Согласно полученным результатам для ${{p}_{0}} = 0$, ${{q}_{0}} \in {{Q}_{0}}$ и $\tau \in \Theta $, максимальные абсолютные значения $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ убывают при уменьшении $\tau $.

Трехстадийные симметрично-симплектические методы Рунге–Кутты. Для этих методов Рунге–Кутты (постоянные ${{b}_{1}}$, ${{s}_{{12}}}$) $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ в случае точного решения системы разрешающих уравнений (3.3), точного вычисления производных и точной арифметики равны нулю. Однако значения этих величин, вычисленные с помощью правых разностей, в общем случае не равны нулю из-за погрешности машинной арифметики и погрешностей решения системы разрешающих уравнений и вычисления производных. Расчеты показывают, что значения дефектов симплектичности удовлетворяют неравенствам ${\text{|}}\Delta {{S}_{{n,12}}}{\text{|}} < {{10}^{{ - 20}}}$, ${\text{|}}\delta {{S}_{{n,12}}}{\text{|}} < {{10}^{{ - 31}}}$ на промежутках времени $t \in [0,1000T].$

Метод дискретного градиента. На фиг. 2 показаны графики зависимостей $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ от времени на малом промежутке времени. В табл. 3 представлены максимальные абсолютные значения $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ на отрезках $t \in [0,1000T]$ для разных начальных значений $({{p}_{0}},{{q}_{0}})$ и величин шага $\tau $. Как показали исследования для ${{p}_{0}} = 0$ и ${{q}_{0}} \in {{Q}_{0}}$, максимальные абсолютные значения $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ убывают при уменьшении ${{H}_{0}}$ и $\tau $.

Неявные гнездовые методы Рунге–Кутты. На фиг. 3–5 изображены графики зависимостей $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ от времени на малом промежутке времени для НГРК-методов. Согласно проведенным исследованиям для ${{p}_{0}} = 0$ и ${{q}_{0}} \in {{Q}_{0}}$, эти зависимости для всех трех рассматриваемых НГРК-методов являются почти периодическими. Форма графиков, соответствующих одинаковым ${{q}_{0}}$, на фиг. 3 и 5 совпадает. В табл. 4 представлены максимальные абсолютные значения $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ на отрезках $t \in [0,1000T]$ для НРГК-методов типа Гаусса при разных начальных значениях $({{p}_{0}},{{q}_{0}})$ и величинах шага $\tau $. Максимальные абсолютные значения $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ для метода Лобатто IIIA примерно в 1.7 раза меньше соответствующих значений для двухуровневого НГРК-метода типа Гаусса. При ${{p}_{0}} = 0$ и ${{q}_{0}} \in {{Q}_{0}}$ максимальные абсолютные значения $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ для всех рассматриваемых НГРК-методов убывают при уменьшении $\tau $.

Из анализа графиков на фиг. 1–5, табл. 3, 4 и не приведенных в настоящей работе данных для остальных начальных значений ${{p}_{0}} = 0$, ${{q}_{0}} \in {{Q}_{0}}$ следует, что абсолютные значения компонент дефектов симплектичности всех рассматриваемых методов не имеют тенденции к возрастанию на больших промежутках времени.

Сравнение фиг. 1 и 2, данные табл. 3 и данные для остальных $({{p}_{0}},{{q}_{0}})$ указывают на то, что $\mathop {max}\limits_n {\text{|}}\Delta {{S}_{{n,12}}}{\text{|}}$ и $\mathop {max}\limits_n {\text{|}}\delta {{S}_{{n,12}}}{\text{|}}$ для нового консервативного метода нулевого дисбаланса полной энергии при одинаковых значениях шага по крайней мере на 3 порядка меньше, чем соответствующие значения для метода дискретного градиента. Компоненты дефектов симплектичности обоих методов имеют сравнимые величины, если шаг в методе дискретного градиента на 2 порядка меньше.

Максимальные абсолютные значения $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ для двухуровневого НРГК-метода типа Гаусса и метода Лобатто IIIA на порядок или несколько порядков больше, чем соответствующие величины для метода нулевого дисбаланса (при фиксированных ${{q}_{0}}$, $\tau $).

Максимальные абсолютные значения $\Delta {{S}_{{n,12}}}$ и $\delta {{S}_{{n,12}}}$ для трехуровневого НГРК-метода типа Гаусса имеют сопоставимые величины (отличаются не более чем в 3 раза) с соответствующими значениями для метода нулевого дисбаланса, за исключением начальных значений ${{p}_{0}} = 0$, ${{q}_{0}} = 0.99$. В последнем случае трехуровневый НГРК-метод типа Гаусса дает меньшие $\mathop {max}\limits_n {\text{|}}\Delta {{S}_{{n,12}}}{\text{|}}$ и $\mathop {max}\limits_n {\text{|}}\delta {{S}_{{n,12}}}{\text{|}}$. При ${{p}_{0}} = 0$, ${{q}_{0}} \in \{ 0.7,0.8,0.9\} $ метод нулевого дисбаланса дает меньшие $\mathop {max}\limits_n {\text{|}}\Delta {{S}_{{n,12}}}{\text{|}}$, но большие $\mathop {max}\limits_n {\text{|}}\delta {{S}_{{n,12}}}{\text{|}}$, чем трехуровневый НГРК-метод типа Гаусса. При ${{p}_{0}} = 0$, ${{q}_{0}} \in \{ 0.05,0.1,0.2,0.3,0.4,0.5,0.6\} $ метод нулевого дисбаланса дает меньшие $\mathop {max}\limits_n {\text{|}}\Delta {{S}_{{n,12}}}{\text{|}}$ и $\mathop {max}\limits_n {\text{|}}\delta {{S}_{{n,12}}}{\text{|}}$.

Вычисления на двойной точности. При численном вычислении частных производных приращение аргумента должно быть не меньше величины, линейно зависящей от квадратного корня из величины погрешности вычисления функции [18, с. 186], [20], [21]. Из-за этого двойной точности может быть недостаточно для вычисления $\Delta {{S}_{{n,12}}}$, $\delta {{S}_{{n,12}}}$. Как показали расчеты, зависимости $\Delta {{S}_{{n,12}}}$, $\delta {{S}_{{n,12}}}$ в таком случае могут значительно отличаться от зависимостей, полученных в случае повышенной точности, например, не быть периодическими и нарастать со временем. Отметим, что в случае метода нулевого дисбаланса для установления почти периодической зависимости $\Delta {{S}_{{n,12}}}$ от времени при ${{p}_{0}} = 0$, ${{q}_{0}} = 0.05$, $\tau = 0.05{{T}_{{\text{Г}}}}$, $t \in [0,1000T]$ нужно использовать приращение аргумента в правой разностной производной менее ${{10}^{{ - 18}}}$, которое слишком мало для вычислений на двойной точности.

Выводы. Итак, строго говоря, метод нулевого дисбаланса полной энергии, метод дискретного градиента и НГРК-методы не являются симплектическими. Однако максимальные абсолютные значения компонент дефектов симплектичности $\mathop {max}\limits_n {\text{|}}\Delta {{S}_{{n,12}}}{\text{|}}$, $\mathop {max}\limits_n {\text{|}}\delta {{S}_{{n,12}}}{\text{|}}$ для метода нулевого дисбаланса на порядки меньше соответствующих значений для метода дискретного градиента. Двухуровневые НГРК-методы (типа Гаусса и Лобатто) также дают бóльшие значения $\mathop {max}\limits_n {\text{|}}\Delta {{S}_{{n,12}}}{\text{|}}$, $\mathop {max}\limits_n {\text{|}}\delta {{S}_{{n,12}}}{\text{|}}$, чем метод нулевого дисбаланса. Только трехуровневый НГРК-метод типа Гаусса при ${{p}_{0}} = 0$, ${{q}_{0}} \in \{ 0.7,0.8,0.9,0.99\} $ (из рассмотренных ${{p}_{0}} = 0$, ${{q}_{0}} \in {{Q}_{0}}$) дает меньшие значения $\mathop {max}\limits_n {\text{|}}\delta {{S}_{{n,12}}}{\text{|}}$ и при ${{p}_{0}} = 0$, ${{q}_{0}} = 0.99$ – меньшие значения $\mathop {max}\limits_n {\text{|}}\Delta {{S}_{{n,12}}}{\text{|}}$, чем метод нулевого дисбаланса.

6.2. Обратимость

Для исследования обратимости численного решения задачи Коши введем векторы $\Delta {{\rho }_{n}}$, $\delta {{\rho }_{n}}$ – глобального и локального дефектов обратимости на $n$-м шаге. В случае модельной задачи 1 $\Delta {{\rho }_{n}} \in {{R}^{2}}$, $\delta {{\rho }_{n}} \in {{R}^{2}}$. Глобальный дефект обратимости $\Delta {{\rho }_{n}}$ вычисляется следующим образом. Пусть имеется значение численного решения $({{p}_{n}},{{q}_{n}})$, найденное исследуемым методом для начальных значений $({{p}_{0}},{{q}_{0}})$. Для начальных значений $( - {{p}_{n}},{{q}_{n}})$ тем же методом находится значение решения $(\mathop {\underline p }\nolimits_{2n} ,\mathop {\underline q }\nolimits_{2n} )$ с помощью $n$ шагов. Тогда имеем

Локальный дефект обратимости $\delta {{\rho }_{n}}$ вычисляется следующим способом. Возьмем приближенное значение решения $({{p}_{n}},{{q}_{n}})$, найденное исследуемым методом. Сделаем один шаг метода для начальных значений $( - {{p}_{n}},{{q}_{n}})$, получим $(\mathop {\underline p }\nolimits_{n + 1} ,\mathop {\underline q }\nolimits_{n + 1} )$. Значение $\delta {{\rho }_{n}}$ вычисляется по формуле

Максимальные по ${{q}_{0}}$, $\tau $ и $n$ значения компонент глобального и локального дефектов обратимости ($ma{{x}_{{{{p}_{0}} = 0,\,{{q}_{0}} \in {{Q}_{0}},\,\tau \in \Theta ,n}}}{{\left\| {\Delta {{\rho }_{n}}} \right\|}_{c}}$ и $ma{{x}_{{{{p}_{0}} = 0,{{q}_{0}} \in {{Q}_{0}},\tau \in \Theta ,n}}}{{\left\| {\delta {{\rho }_{n}}} \right\|}_{c}}$) указаны в табл. 5. Глобальный дефект обратимости вычислялся на временны́х отрезках $t \in [0,200T]$ в случае двойной точности и $t \in [0,50T]$ в случае повышенной; локальный дефект обратимости – на отрезках $t \in [0,1000T]$. Вычисление $\Delta {{\rho }_{n}}$ на более коротких временны́х отрезках обусловлено длительностью расчетов.

Для метода дискретного градиента при $\tau = 0.0002{{T}_{{\text{Г}}}}$ получены равенства: $ma{{x}_{{{{p}_{0}} = 0,{{q}_{0}} \in {{Q}_{0}},n}}}{{\left\| {\Delta {{\rho }_{n}}} \right\|}_{c}} = $ $ = 4.87092 \times {{10}^{{ - 10}}}$, $ma{{x}_{{{{p}_{0}} = 0,{{q}_{0}} \in {{Q}_{0}},n}}}{{\left\| {\delta {{\rho }_{n}}} \right\|}_{c}} = 9.86548 \times {{10}^{{ - 8}}}$ для двойной точности; $ma{{x}_{{{{p}_{0}} = 0,{{q}_{0}} \in {{Q}_{0}},n}}}{{\left\| {\Delta {{\rho }_{n}}} \right\|}_{c}} = $ $ = 1.85364 \times {{10}^{{ - 72}}}$, $ma{{x}_{{{{p}_{0}} = 0,{{q}_{0}} \in {{Q}_{0}},n}}}{{\left\| {\delta {{\rho }_{n}}} \right\|}_{c}} = 1.12959 \times {{10}^{{ - 76}}}$ для повышенной точности. В этих случаях расчеты $\Delta {{\rho }_{n}}$ проводились для $t \in [0,10T]$.

Сравним глобальный и локальный дефекты обратимости используя максимальные по ${{q}_{0}}$, $\tau $ и $n$ абсолютные значения компонент. Глобальный и локальный дефекты обратимости значительно меньше в случае повышенной точности. Методы нулевого дисбаланса и трехстадийные симметрично-симплектические методы Рунге–Кутты имеют примерно одинаковые глобальные дефекты обратимости в случае двойной точности. В случае повышенной точности глобальный дефект обратимости для метода нулевого дисбаланса на три порядка меньше, чем глобальные дефекты обратимости для трехстадийных симметрично-симплектических методов Рунге–Кутты. Локальные дефекты обратимости трехстадийных симметрично-симплектических методов Рунге–Кутты меньше на шесть (в случае двойной точности) или два (в случае повышенной точности) порядка, чем соответствующие дефекты метода нулевого дисбаланса. Метод дискретного градиента дает глобальный дефект обратимости, который на несколько порядков больше соответствующего дефекта метода нулевого дисбаланса. Локальный дефект обратимости метода дискретного градиента имеет примерно такой же порядок, как и локальный дефект обратимости метода нулевого дисбаланса. Дефекты обратимости НГРК-методов отличаются от соответствующих дефектов обратимости трехстадийных симметрично-симплектических методов Рунге–Кутты не более чем на один порядок. Для всех восьми вычислительных методов компоненты глобального дефекта обратимости могут нарастать со временем, компоненты локального дефекта обратимости колеблются в окрестности 0.

Таким образом, все восемь вычислительных методов с большой точностью являются обратимыми.

6.3. Дисбаланс полной энергии

Рассмотрим значения дисбаланса полной энергии $\Delta {{H}_{n}} = H({{p}_{n}},{{q}_{n}}) - {{H}_{0}}$, полученные при использовании исследуемых методов решения задачи Коши. В табл. 6 приведены максимальные абсолютные значения дисбаланса полной энергии на отрезках $t \in [0,1000T]$ для различных $({{p}_{0}},{{q}_{0}})$ и $\tau $ при двойной и повышенной точности.

Метод нулевого дисбаланса полной энергии. Как показали вычисления для ${{p}_{0}} = 0$, ${{q}_{0}} \in {{Q}_{0}}$ и $\tau \in \Theta $, абсолютные значения дисбаланса полной энергии на приближенном решении модельной задачи, полученном методом нулевого дисбаланса полной энергии при повышенной точности, не превышают величины порядка ${{10}^{{ - 123}}}$. Для двойной точности абсолютные значения дисбаланса полной энергии не превышают величины порядка ${{10}^{{ - 14}}}$ и практически не зависят от рассматриваемых величин шага.

Метод дискретного градиента. Абсолютные значения дисбаланса полной энергии для метода дискретного градиента не превышают величины порядка ${{10}^{{ - 69}}}$ в случае повышенной точности. Для двойной точности абсолютные значения дисбаланса полной энергии не превышают величины порядка ${{10}^{{ - 9}}}$.

Трехстадийные симметрично-симплектические методы Рунге–Кутты. Для трехстадийных симметрично-симплектических методов Рунге–Кутты при повышенной точности зависимости дисбаланса полной энергии $\Delta {{H}_{n}}$ от времени являются почти периодическими. В случае двойной точности при ${{b}_{1}} = 5{\text{/}}18$, ${{s}_{{12}}} = 0$ и ${{b}_{1}} = 0.5$, ${{s}_{{12}}} = 0$ зависимости $\Delta {{H}_{n}}$ также являются почти периодическими. Однако для двойной точности зависимости $\Delta {{H}_{n}}$ при ${{b}_{1}} = 5{\text{/}}18$, ${{s}_{{12}}} = s_{{12}}^{*}$ и некоторых комбинациях ${{p}_{0}}$, ${{q}_{0}}$ и $\tau $ (например, при ${{p}_{0}} = 0$, ${{q}_{0}} = 0.99$, $\tau = 0.01{{T}_{{\text{Г}}}}$) не являются периодическими.

Максимальные абсолютные значения дисбаланса полной энергии $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$ для повышенной точности на отрезках $t \in [0,1000T]$ при различных начальных значениях $({{p}_{0}},{{q}_{0}})$, значениях параметров ${{b}_{1}}$, ${{s}_{{12}}}$ и величинах шага $\tau $ приведены в табл. 7. Для двойной точности величины $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$ при указанных в табл. 7 значениях $({{p}_{0}},{{q}_{0}})$ имеют тот же порядок и отличаются не более чем на один знак после запятой, за исключением значения $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}} = 6.21725 \times {{10}^{{ - 15}}}$ для ${{p}_{0}} = 0$, ${{q}_{0}} = 0.9$, $\tau = 0.01{{T}_{{\text{Г}}}}$, ${{b}_{1}} = 5{\text{/}}18$, ${{s}_{{12}}} = s_{{12}}^{*}$, а также значений $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}} < 3 \times {{10}^{{15}}}$ для ${{p}_{0}} = 0$, ${{q}_{0}} = 0.99$, $\tau \in \{ 0.01{{T}_{Г}},0.02{{T}_{Г}},0.03{{T}_{Г}}\} $, b₁ = 5/18, ${{s}_{{12}}} = s_{{12}}^{*}$. Согласно данным расчетов для ${{p}_{0}} = 0$, ${{q}_{0}} \in {{Q}_{0}}$ и $\tau \in \Theta $, максимальные абсолютные значения дисбаланса полной энергии $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$ убывают при уменьшении $\tau $.

Итак, выбор значений ${{s}_{{12}}}$ в соответствии с уравнением нулевого дисбаланса (3.5), как правило, уменьшает абсолютное значение дисбаланса и делает его практически не зависимым от изменений величины шага в разумных пределах. Только в случае двойной точности при некоторых из рассмотренных комбинаций ${{p}_{0}}$, ${{q}_{0}}$, $\tau $ (при ${{p}_{0}} = 0$, ${{q}_{0}} = 0.9$, $\tau = 0.01{{T}_{{\text{Г}}}}$; ${{p}_{0}} = 0$, ${{q}_{0}} = 0.99$, $\tau \in \{ 0.03{{T}_{Г}},0.04{{T}_{Г}}\} $) значения $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$ при постоянных ${{b}_{1}} = 5{\text{/}}18$, ${{s}_{{12}}} = s_{{12}}^{*}$ меньше, чем для метода нулевого дисбаланса. Отметим, что при ${{p}_{0}} = 0$, ${{q}_{0}} = 0.99$, $\tau \in \{ 0.01{{T}_{Г}},0.02{{T}_{Г}}\} $ в случае двойной точности значения $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$ для постоянных ${{b}_{1}} = 5{\text{/}}18$, ${{s}_{{12}}} = s_{{12}}^{*}$ и для метода нулевого дисбаланса совпадают точно. Это вызвано тем, что условие остановки внешних итераций (3.6) метода нулевого дисбаланса удовлетворяется уже для начального приближения $s_{{12}}^{{(0)}} = s_{{12}}^{*}$ на всех шагах метода нулевого дисбаланса, т.е. внешние итерации при расчетах фактически не выполнялись. При ${{b}_{1}} = 5{\text{/}}18$, ${{s}_{{12}}} = 0$ и ${{b}_{1}} = 0.5$, ${{s}_{{12}}} = 0$ для двойной точности, а также при всех рассмотренных комбинациях постоянных ${{b}_{1}}$, ${{s}_{{12}}}$ для повышенной точности получаются бóльшие $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$, чем при использовании метода нулевого дисбаланса (по данным расчетов для ${{p}_{0}} = 0$, ${{q}_{0}} \in {{Q}_{0}}$, $\tau \in \Theta $).

На фиг. 6 изображены графики зависимостей значения ${{s}_{{12}}}$ от времени для метода нулевого дисбаланса. Зависимости являются почти периодическими в случае повышенной точности вычислений при ${{p}_{0}} = 0$, ${{q}_{0}} \in {{Q}_{0}}$ и $\tau \in \Theta $. Величины ${{s}_{{12}}}$ отличаются от $s_{{12}}^{*}$ не более чем в третьем знаке после запятой. Такие поправки позволяют уменьшить дисбаланс полной энергии, особенно в случае повышенной точности.

Метод дискретного градиента в случае повышенной точности также обеспечивает при ${{p}_{0}} = 0$, ${{q}_{0}} \in {{Q}_{0}}$ и $\tau \in \Theta $ гораздо меньшие значения $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$, чем трехстадийные симметрично-симплектические методы Рунге–Кутты. Однако для двойной точности метод Кунцмана–Бутчера обеспечивает для большинства комбинаций ${{q}_{0}} \in {{Q}_{0}}$, $\tau \in \Theta $ меньшие $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$. При ${{b}_{1}} = 5{\text{/}}18,$ ${{b}_{1}} = 0.{\text{5}}$ и ${{s}_{{12}}} = 0$ трехстадийные симметрично-симплектические методы Рунге–Кутты дают бóльшие значения $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$, чем метод дискретного градиента, для подавляющего большинства комбинаций ${{q}_{0}} \in {{Q}_{0}}$, $\tau \in \Theta $ в случае двойной точности.

Неявные гнездовые методы Рунге–Кутты. Характеристики зависимостей дисбаланса полной энергии $\Delta {{H}_{n}}$ от времени для НГРК-методов аналогичны характеристикам соответствующих зависимостей для трехстадийных симметрично-симплектических методов Рунге–Кутты. При повышенной точности зависимости $\Delta {{H}_{n}}$ от времени являются почти периодическими. Для двойной точности – зависимости почти периодические в случае двухуровневого метода типа Гаусса и метода Лобатто IIIA. Для трехуровневого метода типа Гаусса зависимости $\Delta {{H}_{n}}$ могут не являться периодическими при двойной точности.

Максимальные абсолютные величины дисбаланса полной энергии $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$ в случае повышенной точности на отрезках $t \in [0,1000T]$ для НГРК-методов при различных начальных значениях $({{p}_{0}},{{q}_{0}})$ представлены в табл. 8. Для двойной точности величины $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$ при указанных в этой таблице $({{p}_{0}},{{q}_{0}})$ имеют тот же порядок и отличаются не более чем на один знак после запятой, за исключением значения $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}} = 4.16334 \times {{10}^{{ - 16}}}$ для трехуровневого НГРК-метода типа Гаусса при ${{p}_{0}} = 0$, ${{q}_{0}} = 0.99$, $\tau = 0.01{{T}_{{\text{Г}}}}$. По крайней мере для ${{p}_{0}} = 0$, ${{q}_{0}} \in {{Q}_{0}}$ и $\tau \in \Theta $ максимальные абсолютные значения дисбаланса полной энергии $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$ убывают при уменьшении $\tau $.

Как показали расчеты для ${{p}_{0}} = 0$, ${{q}_{0}} \in {{Q}_{0}}$, $\tau \in \Theta $, в случае повышенной точности значения $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$ для НГРК-методов значительно больше соответствующих значений для метода нулевого дисбаланса полной энергии и метода дискретного градиента. В случае двойной точности при ${{p}_{0}} = 0$, ${{q}_{0}} \in {{Q}_{0}}$, $\tau \in \Theta $ значения $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$ для НГРК-методов больше соответствующих значений для метода нулевого дисбаланса. Двухуровневый НГРК-метод типа Гаусса дает бóльшие значения $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$, чем метод дискретного градиента при всех ${{q}_{0}} \in {{Q}_{0}}$, $\tau \in \Theta $, кроме ${{q}_{0}} = 0.99$, $\tau \in \{ 0.01{{T}_{Г}},0.02{{T}_{Г}}\} $. Метод Лобатто IIIA также дает бóльшие значения $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$, чем метод дискретного градиента при всех ${{q}_{0}} \in {{Q}_{0}}$, $\tau \in \Theta $, кроме ${{q}_{0}} = 0.99$, $\tau \in \{ 0.01{{T}_{Г}},0.02{{T}_{Г}},0.04{{T}_{Г}}\} $. Только трехуровневый НГРК-метод типа Гаусса дает для большиства комбинаций ${{q}_{0}} \in {{Q}_{0}}$, $\tau \in \Theta $ (для 28 комбинаций из 55) меньшие значения $\mathop {max}\limits_n {\text{|}}\Delta {{H}_{n}}{\text{|}}$, чем метод дискретного градиента.