Журнал вычислительной математики и математической физики, 2022, T. 62, № 1, стр. 3-11

Моделирование динамических процессов в длинных джозефсоновских переходах. Проблема вычисления Вах. Численный метод оценки скорости роста ошибок округления

М. И. Зуев^1, *, С. И. Сердюкова^1, **

¹ Объединенный институт ядерных исследований
141980 М.о., Дубна, ул. Жолио-Кюри, 6, Россия

^* E-mail: zuevmax@jinr.ru
^** E-mail: sis@jinr.ru

Поступила в редакцию 23.03.2021
После доработки 23.03.2021
Принята к публикации 17.09.2021

EDN: LZVSQU
DOI: 10.31857/S0044466922010124

Аннотация

При численных расчетах вольт-амперных характеристик обычно используется схема Рунге–Кутты четвертого порядка точности. Расчеты проводятся на больших интервалах времени и на каждом шаге по времени проводится четырехкратный пересчет правых частей уравнений. Чтобы сократить время счета, предложено использовать вместо схемы Рунге–Кутты “явную” схему второго порядка точности. При $\tau = h$ доказаны оценки $\left\| {{{G}^{n}}} \right\|$ для всех $n$, гарантирующие ограниченность скорости роста ошибок округления, где $G$ – оператор перехода от слоя к слою; $\tau ,\;h$ – шаги сетки по $t,\;x$ соответственно. Разработан численно-аналитический алгоритм оценки ошибок округления для всех $\tau \leqslant h$. Установлена их ограниченность на всем интервале вычисления ВАХ длинных джозефсоновских переходов при использовании предлагаемой схемы. Расчеты проводились на суперкомпьютере GOVORUN с использованием системы REDUCE. Библ. 7. Фиг. 3.

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

ВВЕДЕНИЕ

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

Вычисление вольт-амперных характеристик систем $n$ длинных джозефсоновских переходов (ДДП) связано (см. [1, (7)]) с решением системы $n$ существенно нелинейных дифференциальных уравнений. Разработан алгоритм, позволяющий свести задачу к решению одного уравнения ${{u}_{{tt}}} = - \beta {{u}_{t}} + {{u}_{{xx}}} - sin(u) + I.$ Для решения последнего предлагается использовать разностную аппроксимацию второго порядка точности:

$\frac{{u_{\nu }^{{n + 1}} - 2u_{\nu }^{n} + u_{\nu }^{{n - 1}}}}{{{{\tau }^{2}}}} = - \beta \frac{{u_{\nu }^{{n + 1}} - u_{\nu }^{{n - 1}}}}{{2\tau }} + \frac{{u_{{\nu + 1}}^{n} - 2u_{\nu }^{n} + u_{{\nu - 1}}^{n}}}{{{{h}^{2}}}} - sin(u_{\nu }^{n}) + I(n),$

где $\tau ,\;h$ – шаги сетки по $t,\;x$ соответственно, ниже используются обозначения $\gamma = \tau {\text{/}}h,$ $\delta = \beta \tau {\text{/}}2$. Линейная часть без $I(n)$ (см. замечание 2 в конце работы) рассматриваемого разностного уравнения эквивалентна системе разностных уравнений

(1)

$\begin{gathered} u_{\nu }^{{n + 1}} - u_{\nu }^{n} = v_{\nu }^{{n + 1}}, \\ \frac{{v_{\nu }^{{n + 1}} - v_{\nu }^{n}}}{{{{\tau }^{2}}}} = - \beta \frac{{v_{\nu }^{{n + 1}} + v_{\nu }^{n}}}{{2\tau }} + \frac{{u_{{\nu + 1}}^{n} - 2u_{\nu }^{n} + u_{{\nu - 1}}^{n}}}{{{{h}^{2}}}}. \\ \end{gathered} $

Решается задача Коши, заданы начальные данные $u_{\nu }^{0},\;v_{\nu }^{0}$. После преобразования Фурье

${{U}^{n}} = \sum\limits_{\nu = - \infty }^\infty \,{{e}^{{i\nu \phi }}}[u_{\nu }^{n},v_{\nu }^{n}]{\kern 1pt} *$

(знак $ * $ означает преобразование вектора-строки в вектор-столбец) и некоторых алгебраических манипуляций получаем векторную систему

${{U}^{{n + 1}}} = A({{e}^{{i\phi }}}){{U}^{n}},\quad {{U}^{n}} = {{A}^{n}}({{e}^{{i\phi }}}){{U}^{0}}.$

Здесь $A({{e}^{{i\phi }}})$ – характеристическая матрица

$A({{e}^{{i\phi }}}) = \left[ {\begin{array}{*{20}{c}} {(1 - (4{{\gamma }^{2}}si{{n}^{2}}(\phi {\text{/}}2)){\text{/}}(1 + \delta )}&{(1 - \delta ){\text{/}}(1 + \delta )} \\ { - (4{{\gamma }^{2}}si{{n}^{2}}(\phi {\text{/}}2)){\text{/}}(1 + \delta )}&{(1 - \delta ){\text{/}}(1 + \delta )} \end{array}} \right].$

В случае $\tau = h$ получены (см. [2]) оценки норм $\left\| {{{A}^{n}}({{e}^{{i\phi }}})} \right\|$ в пространстве ${{L}_{2}}$ для всех $n$:

$\left\| {{{A}^{n}}({{e}^{{i\phi }}})} \right\| \leqslant 5\sqrt {\frac{2}{\pi }} {\kern 1pt} \sqrt {2 + \frac{{3.78}}{{\beta \tau }} + O(exp( - \beta \tau n))} ,\quad t = \tau n.$

Кроме того, было доказано, что при больших $\beta \tau n$ норма характеристической матрицы медленно убывает к некоторой константе:

$\left\| {{{A}^{n}}({{e}^{{i\phi }}})} \right\| \leqslant 5\sqrt {\frac{2}{\pi }} \sqrt {2 + \frac{{5.12}}{{\beta \tau \sqrt {\beta \tau n} }} + O(exp( - \beta \tau n))} .$

Полученные оценки гарантируют ограниченность роста ошибок округления в ${{l}_{2}}$-метрике для всех $n$. Особую трудность представляла оценка $\left\| {(\lambda _{1}^{n} - \lambda _{2}^{n}){\text{/}}({{\lambda }_{1}} - {{\lambda }_{2}})} \right\|,$ где ${{\lambda }_{1}},\;{{\lambda }_{2}}$ – собственные значения характеристической матрицы.

В настоящей работе разработан алгоритм прямого вычисления $\left\| {(\lambda _{1}^{n} - \lambda _{2}^{n}){\text{/}}({{\lambda }_{1}} - {{\lambda }_{2}})} \right\|$ при всех $\tau \leqslant h$ на всем интервале вычисления вольт-амперных характеристик ДДП. Это позволило выяснить качественное изменение $\left\| {{{A}^{n}}} \right\|$ в зависимости от $n$. Опишем этот алгоритм.

1. СВЕДЕНИЕ ОЦЕНКИ $\left\| {{{A}^{n}}({{e}^{{i\phi }}})} \right\|$ К ОЦЕНКЕ $\left\| {(\lambda _{1}^{n} - \lambda _{2}^{n}){\text{/}}({{\lambda }_{1}} - {{\lambda }_{2}})} \right\|$

Собственные значения характеристической матрицы $A$ удовлетворяют уравнению

${{\lambda }^{2}} - 2\frac{{1 - 2{{\gamma }^{2}}si{{n}^{2}}(\phi {\text{/}}2)}}{{1 + \delta }}\lambda + \frac{{1 - \delta }}{{1 + \delta }} = 0.$

Решая это уравнение, находим

(2)

${{\lambda }_{{1,2}}} = \frac{{1 - 2{{\gamma }^{2}}si{{n}^{2}}(\phi {\text{/}}2) \pm \sqrt {{{\delta }^{2}} - 4{{\gamma }^{2}}si{{n}^{2}}(\phi {\text{/}}2) + 4{{\gamma }^{4}}si{{n}^{4}}(\phi {\text{/}}2)} }}{{1 + \delta }}.$

Доказано (см. [2]), что $\left| {{{\lambda }_{1}}} \right| \leqslant 1$ и $\left| {{{\lambda }_{2}}} \right| \leqslant 1$ для всех $\gamma \leqslant 1$ и всех $0 \leqslant \phi \leqslant 2\pi $. Это необходимые условия устойчивости рассматриваемой задачи Коши в пространстве ${{L}_{2}}$, гарантирующие отсутствие неустойчивости экспоненциального типа. Но при наличии кратных собственных значений $\lambda ({{e}^{{i\phi }}})$ может возникать неустойчивость степенного типа (см. [3]) $\left\| {{{G}^{n}}} \right\| \asymp {{n}^{\theta }},$ $\theta > 0$.

Матрица преобразования подобия

$T = \left[ {\begin{array}{*{20}{c}} {{{\lambda }_{1}}}&1 \\ {{{\lambda }_{1}} - 1}&1 \end{array}} \right],\quad \det (T) = 1,$

приводит матрицу $A$ к треугольному виду:

(3)

$\begin{gathered} {{T}^{{ - 1}}}AT = B = \left[ {\begin{array}{*{20}{c}} {{{\lambda }_{1}}}&1 \\ 0&{{{\lambda }_{2}}} \end{array}} \right],\quad {{B}^{n}} = \left[ {\begin{array}{*{20}{c}} {\lambda _{1}^{n}}&{{{D}_{n}}} \\ 0&{\lambda _{2}^{n}} \end{array}} \right], \\ {{D}_{n}} = \frac{{\lambda _{1}^{n} - \lambda _{2}^{n}}}{{{{\lambda }_{1}} - {{\lambda }_{2}}}}. \\ \end{gathered} $

Сделав обратное преобразование, получаем

${{A}^{n}} = \left[ {\begin{array}{*{20}{c}} {\lambda _{1}^{{n + 1}} + (1 - {{\lambda }_{1}})\lambda _{2}^{n}}&{ - {{\lambda }_{1}}(\lambda _{1}^{n} - \lambda _{2}^{n})} \\ {({{\lambda }_{1}} - 1)(\lambda _{1}^{n} - \lambda _{2}^{n})}&{(1 - {{\lambda }_{1}})\lambda _{1}^{n} + {{\lambda }_{1}}\lambda _{2}^{n}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{{\lambda }_{1}}(1 - {{\lambda }_{1}})}&{\lambda _{1}^{2}} \\ { - {{{(1 - {{\lambda }_{1}})}}^{2}}}&{{{\lambda }_{1}}({{\lambda }_{1}} - 1)} \end{array}} \right]{{D}_{n}}.$

Элементы матрицы ${{A}^{n}}({{e}^{{i\phi }}})$ являются тригонометрическими полиномами $p({{e}^{{i\phi }}}) = \sum\nolimits_{k = - n}^n \,{{a}_{k}}{{e}^{{ik\phi }}}.$ Норма $p$ есть

$\left\| {p({{e}^{{i\phi }}})} \right\| = \sqrt {\sum\limits_{k = - n}^n \,{{{\left| {{{a}_{k}}} \right|}}^{2}}} = \sqrt {\frac{1}{{2\pi }}\int\limits_{ - \pi }^\pi \,{{{\left| {p({{e}^{{i\phi }}})} \right|}}^{2}}d\phi } .$

Соответствуюшая матричная норма (см. [4]) есть

${{\left\| {{{A}^{n}}({{e}^{{i\phi }}})} \right\|}^{2}} = \sqrt {\sum\limits_{k,j} \,{{{\left\| {a_{{k,j}}^{n}} \right\|}}^{2}}} = \sqrt {\sum\limits_{k,j} \,\frac{1}{{2\pi }}\int\limits_{ - \pi }^\pi \,{{{\left| {a_{{k,j}}^{n}} \right|}}^{2}}d\phi } .$

Так как $\left| {{{\lambda }_{i}}({{e}^{{i\phi }}})} \right| \leqslant 1$ для всех $0 \leqslant \phi \leqslant 2\pi $, справедлива следующая оценка:

$\left\| {{{A}^{n}}({{e}^{{i\phi }}})} \right\| \leqslant 7 + 5\left\| {{{D}_{n}}({{e}^{{i\phi }}})} \right\|.$

Oценка скорости роста ошибок округления свелась к оценке

$\left\| {{{D}_{n}}} \right\| = \sqrt {\frac{1}{{2\pi }}\int\limits_{ - \pi }^\pi \,\mathop {\left| {\frac{{\lambda _{1}^{n} - \lambda _{2}^{n}}}{{{{\lambda }_{1}} - {{\lambda }_{2}}}}} \right|}\nolimits^2 d\phi } = \sqrt {\frac{1}{\pi }\int\limits_0^\pi \,\mathop {\left| {\frac{{\lambda _{1}^{n} - \lambda _{2}^{n}}}{{{{\lambda }_{1}} - {{\lambda }_{2}}}}} \right|}\nolimits^2 d\phi } .$

Так как ${{\lambda }_{i}}$ являются функциями $co{{s}^{2}}(\phi {\text{/}}2)$, то

$\int\limits_{ - \pi }^\pi \,\left| {D_{n}^{2}} \right|d\phi = 2\int\limits_0^\pi \,\left| {D_{n}^{2}} \right|d\phi .$

2. СВЕДЕНИЕ ОЦЕНКИ $\left\| {{{D}_{n}}} \right\|$ К ВЫЧИСЛЕНИЮ КЛАССИЧЕСКОГО ИНТЕГРАЛА

В [2] рассматривалось только $\gamma = 1$. В этом случае собственные значения имеют более простую структуру:

${{\lambda }_{{1,2}}} = \frac{{cos(\phi ) \pm \sqrt {{{\delta }^{2}} - si{{n}^{2}}(\phi )} }}{{1 + \delta }}.$

Суть дела ясна в общем случае для $0 < \gamma \leqslant 1.$ Из оценок, полученных в [2], следует, что $\left\| {{{D}_{n}}} \right\| < 9.525$ для всех $n$.

В этой работе разработан численно-аналитический алгоритм для определения значений $\left\| {{{D}_{n}}} \right\|$ на всем интервале вычисления ВАХ: $t \leqslant {{t}_{{{\text{max}}}}}.$ По этому алгоритму (в случае $\tau = h = 0.1,$ ${{t}_{{{\text{max}}}}} = 1000$) были вычислены значения

$\left\| {{{D}_{n}}} \right\|,\;\;n = 3,7,14,21,35,....,5000,10\,000,$

по которым построен график (фиг. 1). Получена новая более точная оценка $\left\| {{{D}_{n}}} \right\| \leqslant 4.56$. Максимум достигается при $n = 63,$ $\left\| {{{D}_{{63}}}} \right\| = 4.5596...$ (фиг. 2, где представлен график $\left\| {{{D}_{n}}} \right\|$ в окрестности пика).

Фиг. 1.

График $\left\| {{{D}_{n}}} \right\|,\;n \leqslant 10\,000$, в случае $\tau = h = 0.1,$ ${{t}_{{{\text{max}}}}} = 1000$.

Фиг. 2.

В окрестности пика фиг. 1 звездочками отмечены значения $\left\| {{{D}_{n}}} \right\|$. Максимум $\left\| {{{D}_{n}}} \right\|$ достигается при n = 63.

Обе оценки гарантируют ограниченность ошибок округления при решении системы (1) и тем самым доказывают достoверность вычисленных вольт-амперных характеристик ДДП. На фиг. 3 представлены графики $\left\| {{{D}_{n}}} \right\|,\;n \leqslant 5000$, для различных τ, h. Сравнивая эти графики, видим, что быстрее всего $\left\| {{{D}_{n}}} \right\|$ растет в случае $\gamma = 1{\text{/}}5$. При $\gamma = 0$ собственные значения становятся равными единице, а $\left\| {{{D}_{n}}} \right\| = n.$

Фиг. 3.

Графики $\left\| {{{D}_{n}}} \right\|,\;n \leqslant 5000$: сплошная линия – для $\tau = h = 0.1$, штриховая линия – для $\tau = h = 0.05,$ пунктирная линия – для $h = 0.05,$ $\tau = h{\text{/}}5$.

Сделав в (2) замену переменных $z = 1 - {{\gamma }^{2}} + {{\gamma }^{2}}cos(\varphi )$, получаем

${{\lambda }_{{1,2}}} = \frac{{z \pm \sqrt {{{z}^{2}} + {{\delta }^{2}} - 1} }}{{1 + \delta }}.$

После соответствующей подстановки в (3) и ряда алгебраических манипуляций, находим, что ${{D}_{n}}$ есть многочлен степени $n - 1$ от $z$ с вещественными коэффициентами: ${{D}_{n}} = {{P}_{n}}(z){\text{/}}{{(1 + \delta )}^{{n - 1}}},$

(4)

${{P}_{n}} = \frac{{{{{(z + \sqrt {{{z}^{2}} + {{\delta }^{2}} - 1} )}}^{n}} - {{{(z - \sqrt {{{z}^{2}} + {{\delta }^{2}} - 1} )}}^{n}}}}{{2\sqrt {{{z}^{2}} + {{\delta }^{2}} - 1} }} = \sum\limits_{k = 1}^{nc} \,C_{n}^{{k - 1}}{{z}^{{n - (2k - 1)}}}{{({{z}^{2}} + {{\delta }^{2}} - 1)}^{{k - 1}}}.$

Здесь $nc$ есть целая часть $(n + 1){\text{/}}2,$ $nc = [(n + 1){\text{/}}2]$. Если $n$ – нечетное, то ${{P}_{n}}(z)$ и $P_{n}^{2}(z)$ – четные многочлены. Если $n$ – четное, $n - (2k - 1)$ – нечетное, ${{P}_{n}}(z)$ – нечетный, но $P_{n}^{2}(z)$ – четный многочлен.

Заметим, что $z(0) = 1,$ $z(\pi ) = 1 - 2{{\gamma }^{2}},$ $dz = - {{\gamma }^{2}}sin(\phi )d\phi ,$

${{\gamma }^{2}}cos(\phi ) = z - 1 + {{\gamma }^{2}},\quad {{\gamma }^{2}}sin(\phi ) = \sqrt {Y(z)} ,\quad d\phi = - \frac{{dz}}{{\sqrt {Y(z)} }},$

где

(5)

$Y(z) = {{\gamma }^{4}} - {{(z - 1 + {{\gamma }^{2}})}^{2}}.$

Задача свелась к вычислению интегралов

(6)

$\frac{1}{\pi }\int\limits_{1 - 2{{\gamma }^{2}}}^1 \,P_{n}^{2}(z)\frac{{dz}}{{\sqrt {Y(z)} }}.$

Заметим, что $Y(1) = Y(1 - 2{{\gamma }^{2}}) = 0.$ В расчетах был использован хорошо известный алгоритм (см., например, [5]):

$\int {\frac{{P_{n}^{2}(z)}}{{\sqrt {Y(z)} }}dz} = Q(z)\sqrt {Y(z)} + \lambda \int {\frac{{dz}}{{\sqrt {Y(z)} }}} ,$

где полином $Q(z)$ и константа $\lambda $ определяются соотношениeм

$P_{n}^{2}(z) = Q{\kern 1pt} '(z)Y(z) + \frac{1}{2}Q(z)Y{\kern 1pt} '(z) + \lambda .$

В рассматриваемом случае (см. (5)) соотношение принимает следующий вид:

(7)

$P_{n}^{2}(z) = Q{\kern 1pt} '(z)[ - {{z}^{2}} + 2z(1 - {{\gamma }^{2}}) - (1 - 2{{\gamma }^{2}})] + Q(z)[ - z + (1 - {{\gamma }^{2}})] + \lambda .$

В результате имеем

$ - \frac{1}{\pi }\int\limits_1^{1 - 2{{\gamma }^{2}}} \,P_{n}^{2}(z)\frac{{dz}}{{\sqrt {Y(z)} }} = \frac{1}{\pi }\left( {\mathop {\left. {Q(z)\sqrt {Y(z)} } \right|}\nolimits_{1 - 2{{\gamma }^{2}}}^1 + \lambda \int\limits_{1 - 2{{\gamma }^{2}}}^1 \frac{{dz}}{{\sqrt {Y(z)} }}} \right) = \lambda ,\quad \left\| {{{D}_{n}}} \right\| = \frac{{\sqrt \lambda }}{{{{{(1 + \delta )}}^{{n - 1}}}}}.$

В заключение заметим, что

$\frac{1}{\pi }\int\limits_{1 - 2{{\gamma }^{2}}}^1 \,\frac{{dz}}{{\sqrt {{{\gamma }^{4}} - {{{(z - 1 + {{\gamma }^{2}})}}^{2}}} }} = \frac{1}{\pi }\int\limits_{ - 1}^1 \,\frac{{du}}{{\sqrt {1 - {{u}^{2}}} }} = \frac{2}{\pi }\int\limits_0^1 \,\frac{{du}}{{\sqrt {1 - {{u}^{2}}} }} = \frac{2}{\pi }\arcsin (1) = 1.$

По ходу вычисления была сделана замена $z - 1 + {{\gamma }^{2}} = {{\gamma }^{2}}u.$ В случае больших $n$ трудность возникает уже при вычислении коэффициентов ${{P}_{n}}(z)$. Был разработан и реализован c использованием системы REDUCE (см. [6]) специальный алгоритм. При этом пришлось отдельно рассматривать случаи чeтных и нечетных $n$. Далее приводится алгоритм, реально реализованный на суперкомпьютере GOVORUN.

3. АЛГОРИТМ ВЫЧИСЛЕНИЯ $\left\| {{{D}_{n}}} \right\|$, СЛУЧАЙ ЧEТНЫХ $n$

В этом случае $n = 2nc$, (4) принимает вид

(8)

${{P}_{n}}(z) = z\sum\limits_{k = 1}^{nc} \,C_{n}^{{2k - 1}}{{z}^{{2(nc - k)}}}{{({{z}^{2}} + {{\delta }^{2}} - 1)}^{{k - 1}}}.$

Сделав подстановки ${{z}^{2}} = ({{\delta }^{2}} - 1){{w}^{2}}$, получим

${{P}_{n}}(z) = z{{({{\delta }^{2}} - 1)}^{{nc - 1}}}Pt(w),\quad {{P}_{n}}{{(z)}^{2}} = {{z}^{2}}{{({{\delta }^{2}} - 1)}^{{n - 2}}}Pt{{(w)}^{2}},$

где

$Pt(w) = \sum\limits_{k = 1}^{nc} \,C_{n}^{{2k - 1}}{{w}^{{2(nc - k)}}}{{({{w}^{2}} + 1)}^{{k - 1}}}.$

После замены $nc - k = \xi - 1$ имеем

$Pt(w) = \sum\limits_{\xi = 1}^{nc} \,C_{n}^{{2\xi - 1}}{{w}^{{2(\xi - 1)}}}{{({{w}^{2}} + 1)}^{{nc - \xi }}} = \sum\limits_{j = 1}^{nc} \,{{g}_{j}}{{w}^{{2(j - 1)}}},\quad {{g}_{j}} = \sum\limits_{k = 1}^j \,C_{n}^{{2k - 1}}C_{{nc - k}}^{{j - k}}.$

На первых шагах был реализован тривиальный алгоритм вычисления ${{g}_{j}}$:

${{g}_{j}} = \frac{{n{\kern 1pt} !}}{{(nc - j){\kern 1pt} !}}\sum\limits_{k = 1}^j \,\frac{{(nc - k){\kern 1pt} !}}{{(j - k){\kern 1pt} !{\kern 1pt} (2k - 1){\kern 1pt} !{\kern 1pt} (n - 2k + 1){\kern 1pt} !}},\quad 1 \leqslant j \leqslant nc.$

В процессе серии численных экспериментов был найден алгоритм, позволивший сократить расчетное время более, чем в 7 раз:

${{g}_{1}} = n;\quad {{g}_{j}} = \frac{{n{\kern 1pt} !}}{{(n - 2j + 1){\kern 1pt} !{\kern 1pt} (2j - 1){\kern 1pt} !}}\sum\limits_{k = 1}^j \,p(k),\quad 2 \leqslant j \leqslant nc - 1,$

$p(j) = 1,\quad p(k) = p(k + 1)\frac{{(2k + 1)k}}{{(j - k)(n - 2k + 1)}},\quad {{g}_{{nc}}} = {{2}^{{n - 1}}}.$

После того, как ${{g}_{j}}$ определены, коэффициенты $P_{n}^{2}(z)$ вычисляются следующим образом:

$P{{t}^{2}}(w) = \sum\limits_{k = 1}^{n - 1} \,{{b}_{k}}{{w}^{{2(k - 1)}}},\quad P_{n}^{2}(z) = \sum\limits_{k = 1}^{n - 1} \,{{b}_{k}}{{({{\delta }^{2}} - 1)}^{{n - k - 1}}}{{z}^{{2k}}}.$

Заметим, что ${{b}_{1}} = g_{1}^{2} = {{n}^{2}},$ ${{b}_{{n - 1}}} = g_{{nc}}^{2} = {{2}^{{2(n - 1)}}},$

${{b}_{k}} = \sum\limits_{l = 1}^k \,{{g}_{l}}{{g}_{{k - l + 1}}}\;\;{\text{для}}\;\;1 \leqslant k \leqslant nc,\quad {{b}_{k}} = \sum\limits_{l = k - nc + 1}^{nc} \,{{g}_{l}}{{g}_{{k - l + 1}}}\;\;{\text{для}}\;\;k > nc.$

Из определения (6) следует, что $Q(z)$ есть полином степени $2n - 3$:

$Q(z) = \sum\limits_{k = 1}^{n - 1} \,({{a}_{k}}{{z}^{{2k - 1}}} + {{c}_{k}}{{z}^{{2(k - 1)}}}),$

$Q{\kern 1pt} '(z) = \sum\limits_{k = 1}^{n - 1} \,(2k - 1){{a}_{k}}{{z}^{{2(k - 1)}}} + \sum\limits_{k = 2}^{n - 1} \,2(k - 1){{c}_{k}}{{z}^{{2k - 3}}}.$

Вычисляем шаг за шагом

${{a}_{{n - 1}}} = - \frac{{{{b}_{{n - 1}}}}}{{2(n - 1)}},\quad {{c}_{{n - 1}}} = \left( {2 + \frac{1}{{2n - 3}}} \right)(1 - {{\gamma }^{2}}){{a}_{{n - 1}}},$

${{a}_{{k - 1}}} = - {{a}_{k}}(1 - 2{{\gamma }^{2}}) - \frac{{{{b}_{k}}{{{({{\delta }^{2}} - 1)}}^{{n - k - 1}}} + {{a}_{k}}(1 - 2{{\gamma }^{2}})}}{{2(k - 1)}} + \left[ {2 + \frac{1}{{2(k - 1)}}} \right](1 - {{\gamma }^{2}}){{c}_{k}},$

${{c}_{{k - 1}}} = \left( {2 + \frac{1}{{2k - 3}}} \right)(1 - {{\gamma }^{2}}){{a}_{{k - 1}}} - \left( {1 + \frac{1}{{2k - 3}}} \right)(1 - 2{{\gamma }^{2}}){{c}_{k}},\quad k = n - 1,\;n - 2,...,2.$

После того, как величины ${{a}_{1}},\;{{c}_{1}}$ определены, вычисляем $\left\| {{{D}_{n}}} \right\|$. Учитывая, что

$P_{n}^{2}(0) = 0\;\;({\text{cм}}.\;(8)),\quad Q(0) = {{a}_{1}},\quad Q{\kern 1pt} '(0) = {{c}_{1}},\quad Y(0) = 1 - 2{{\gamma }^{2}},\quad Y{\kern 1pt} '(0) = 2(1 - {{\gamma }^{2}})\;\;({\text{см}}.\;(5)),$

получаем, что при $z = 0$ соотношение (7) принимает вид

$0 = - (1 - 2{{\gamma }^{2}}){{a}_{1}} + (1 - {{\gamma }^{2}}){{c}_{1}} + \lambda ,$

следовательно,

$\lambda = (1 - 2{{\gamma }^{2}}){{a}_{1}} - (1 - {{\gamma }^{2}}){{c}_{1}},\quad \left\| {{{D}_{n}}} \right\| = \sqrt \lambda {\text{/}}{{(1 + \delta )}^{{n - 1}}}.$

4. АЛГОРИТМ ВЫЧИСЛЕНИЯ $\left\| {{{D}_{n}}} \right\|$, CЛУЧАЙ НЕЧEТНЫХ $n$

В этом случае $n = 2nc - 1$, (4) принимает вид

(9)

${{P}_{n}}(z) = \sum\limits_{k = 1}^{nc} \,C_{n}^{{2k - 1}}{{z}^{{2(nc - k)}}}{{({{z}^{2}} + {{\delta }^{2}} - 1)}^{{k - 1}}}.$

Сделав подстановку ${{z}^{2}} = ({{\delta }^{2}} - 1){{w}^{2}}$, имеем

${{P}_{n}}(z) = {{({{\delta }^{2}} - 1)}^{{nc - 1}}}Pt(w),\quad P_{n}^{2}(z) = {{({{\delta }^{2}} - 1)}^{{n - 1}}}P{{t}^{2}}(w),$

$Pt(w) = \sum\limits_{k = 1}^{nc} \,C_{n}^{{2k - 1}}{{w}^{{2(nc - k)}}}{{(1 + {{w}^{2}})}^{{k - 1}}},\quad P{{t}^{2}}(w) = \sum\limits_{k = 1}^n \,{{b}_{j}}{{w}^{{2(j - 1)}}}.$

После замены $nc - k = \xi - 1$ получаем

$Pt(w) = \sum\limits_{\xi = 1}^{nc} \,C_{n}^{{n - 2(\xi - 1)}}{{w}^{{2(\xi - 1)}}}{{(1 + {{w}^{2}})}^{{nc - \xi }}} = \sum\limits_{j = 1}^{nc} \,{{g}_{j}}{{w}^{{2(j - 1)}}},\quad {{g}_{j}} = \sum\limits_{k = 1}^j \,C_{n}^{{n - 2(k - 1)}}C_{{nc - k}}^{\theta }.$

Здесь $\theta $ определяется соотношением $2(k - 1) + 2\theta = 2(j - 1),$ $\theta = j - k.$ Первоначально ${{g}_{j}}$ вычислялись напрямую:

${{g}_{j}} = \sum\limits_{k = 1}^j \,\frac{{n{\kern 1pt} !{\kern 1pt} (nc - k){\kern 1pt} {\kern 1pt} !}}{{(2(k - 1)){\kern 1pt} !{\kern 1pt} (n - 2(k - 1)){\kern 1pt} !{\kern 1pt} (j - k){\kern 1pt} !{\kern 1pt} (nc - j){\kern 1pt} !}}.$

${{g}_{1}} = 1,\quad {{g}_{j}} = \frac{{n{\kern 1pt} !}}{{(n - 2(j - 1)){\kern 1pt} !{\kern 1pt} (2(j - 1)){\kern 1pt} !}}\sum\limits_{k = 1}^j \,p(k),\quad 2 \leqslant j \leqslant nc - 1,\quad {{g}_{{nc}}} = {{2}^{{n - 1}}},$

$p(j) = 1,\quad p(k) = p(k + 1)\frac{{k(2k - 1)}}{{(j - k)[n - 2(k - 1)]}},\quad k = j - 1,\;j - 2,...,1.$

После того, как ${{g}_{j}}$ определены, коэффициенты $P_{n}^{2}(z)$ вычисляются следующим образом.

$P_{n}^{2}(z) = \sum\limits_{k = 1}^n \,{{b}_{k}}{{({{\delta }^{2}} - 1)}^{{n - k}}}{{z}^{{2(k - 1)}}},$

где

${{b}_{k}} = \sum\limits_{l = 1}^k \,{{g}_{l}}{{g}_{{k - l + 1}}}\;\;{\text{для}}\;\;k \leqslant nc,\quad {{b}_{k}} = \sum\limits_{l = k}^{nc} \,{{g}_{l}}{{g}_{{k - l + 1}}}\;\;{\text{для}}\;\;k > nc.$

В частности, имеем ${{b}_{1}} = 1,$ ${{b}_{n}} = {{2}^{{2(n - 1)}}}.$ Из определения (6) следует, что $Q(z)$ есть полином степени $2n - 3$:

$Q(z) = \sum\limits_{k = 1}^{n - 1} \,({{a}_{k}}{{z}^{{2k - 1}}} + {{c}_{k}}{{z}^{{2(k - 1)}}}),\quad Q{\kern 1pt} '(z) = \sum\limits_{k = 1}^{n - 1} \,(2k - 1){{a}_{k}}{{z}^{{2(k - 1)}}} + \sum\limits_{k = 2}^{n - 1} \,2(k - 1){{c}_{k}}{{z}^{{2k - 3}}}.$

И мы вычисляем последовательно

${{a}_{{n - 1}}} = - \frac{{{{b}_{n}}}}{{2(n - 1)}},\quad {{c}_{{n - 1}}} = {{a}_{{n - 1}}}(1 - {{\gamma }^{2}})\left( {2 + \frac{1}{{2n - 3}}} \right),$

${{a}_{{k - 1}}} = - {{a}_{k}}(1 - 2{{\gamma }^{2}}) - \frac{{{{b}_{k}}{{{({{\delta }^{2}} - 1)}}^{{n - k}}} + {{a}_{k}}(1 - 2{{\gamma }^{2}})}}{{2(k - 1)}} + {{c}_{k}}\left[ {2 + \frac{1}{{2(k - 1)}}} \right](1 - {{\gamma }^{2}}),$

${{c}_{{k - 1}}} = {{a}_{{k - 1}}}\left( {2 + \frac{1}{{2k - 3}}} \right)(1 - {{\gamma }^{2}}) - {{c}_{k}}\left( {1 + \frac{1}{{2k - 3}}} \right)(1 - 2{{\gamma }^{2}}),\quad k = n - 1,\;n - 2,...,2.$

После того, как ${{a}_{1}},\;{{c}_{1}}$ определены, вычисляем $\left\| {{{D}_{n}}} \right\|$. Принимая во внимание, что

$\begin{gathered} P_{n}^{2}(0) = {{({{\delta }^{2}} - 1)}^{{n - 1}}}\;\;({\text{см}}.\;(9)),\quad Q(0) = {{c}_{1}},\quad Q{\kern 1pt} '(0) = {{a}_{1}}, \\ Y(0) = - 1 + 2{{\gamma }^{2}},\quad Y{\kern 1pt} '(0) = 2(1 - {{\gamma }^{2}})\;\;({\text{см}}.\;(5)), \\ \end{gathered} $

получаем, что при $z = 0$ соотношение (7) принимает вид

$P_{n}^{2}(0) = ( - 1 + 2{{\gamma }^{2}}){{a}_{1}} + (1 - {{\gamma }^{2}}){{c}_{1}} + \lambda ,$

следовательно,

$\lambda = (1 - 2{{\gamma }^{2}}){{a}_{1}} - (1 - {{\gamma }^{2}}){{c}_{1}} + {{({{\delta }^{2}} - 1)}^{{n - 1}}},\quad \left\| {{{D}_{n}}} \right\| = \sqrt \lambda {\text{/}}{{(1 + \delta )}^{{n - 1}}}.$

Замечание 1. Большую часть времени занимает вычисление коэффициентов ${{g}_{j}}$. Использование многопроцессорных компьютеров упрощает проблему существенно. Мы можем одновремeнно вычислять ${{g}_{j}}$ для $j = 1 + m \times k,$ $2 + m \times k,...,m + m \times k \leqslant nc$ на $m$ процессорах.

Замечание 2. Ограниченность $\left\| {{{A}^{n}}({{e}^{{i\varphi }}})} \right\|$ гарантирует (в известной степени) корректность вычисления ВАХ длинных джозефсоновских переходов при использовании рассматриваемой разностной схемы с учетом отброшенных $ - sin(u_{\nu }^{n}) + I(n).$ Ошибки округления $\epsilon _{\nu }^{n},\;\eta _{\nu }^{n}$ (соответствующие $u_{\nu }^{n},\;{v}_{\nu }^{n}$) удовлетворяют системе уравнений

$\epsilon _{\nu }^{{n + 1}} = \epsilon _{\nu }^{n} + \frac{{1 - \delta }}{{1 + \delta }}\eta _{\nu }^{n} + \frac{{{{\gamma }^{2}}}}{{1 + \delta }}(\epsilon _{{\nu + 1}}^{n} - 2\epsilon _{\nu }^{n} + \epsilon _{{\nu - 1}}^{n}) + f_{\nu }^{n},$

$\eta _{\nu }^{{n + 1}} = \frac{{1 - \delta }}{{1 + \delta }}\eta _{\nu }^{n} + \frac{{{{\gamma }^{2}}}}{{1 + \delta }}(\epsilon _{{\nu + 1}}^{n} - 2\epsilon _{\nu }^{n} + \epsilon _{{\nu - 1}}^{n}) + f_{\nu }^{n}.$

Здесь $f_{\nu }^{n} = {{\tau }^{2}}[sin(u_{\nu }^{n} + \epsilon _{\nu }^{n}) - sin(u_{\nu }^{n})]{\text{/}}(1 + \delta ).$ Справедливы оценки

$\left| {f_{\nu }^{n}} \right| = \left| {{{\tau }^{2}}2sin(\epsilon _{\nu }^{n}{\text{/}}2)cos(u_{\nu }^{n} + \epsilon _{\nu }^{n}{\text{/}}2){\text{/}}(1 + \delta )} \right| \leqslant \left| {{{\tau }^{2}}\epsilon _{\nu }^{n}{\text{/}}(1 + \delta )} \right|,$

$\left\| {{{f}^{n}}} \right\| \leqslant {{\tau }^{2}}\left\| {{{\epsilon }^{n}}} \right\|{\text{/}}(1 + \delta ),\quad {{\left\| {{{f}^{n}}} \right\|}^{2}} = \sum\limits_\nu \,{{\left| {f_{\nu }^{n}} \right|}^{2}},\quad {{\left\| {{{\epsilon }^{n}}} \right\|}^{2}} = \sum\limits_\nu \,{{\left| {\varepsilon _{\nu }^{n}} \right|}^{2}}.$

Отсюда следует, что система (1) определяет (с точностью до $O({{\tau }^{2}})$) оператор перехода от слоя к слою ошибок округления, возникающих при использовании предлагаемой разностной схемы.

Благодарим Е.А. Колганову за бесценную помощь на первой стадии вычислений, А.М. Рапортиренко, установившему REDUCE на суперкомпьютер GOVORUN. Благодарим Г.М. Кобелькова за полезные критические замечания.

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

Башашин М.В., Земляная Е.В., Рахмонов И.Р., Шукринов Ю.М., Атанасова П.Х., Волохова А.В. Вычислительная схема и параллельная реализация для моделирования системы длинных джозефсоновских переходов // Компьют. исслед. и моделирование. 2016. Т. 6. № 4. С. 593–604.
Сердюкова С.И. Моделирование динамических процессов в длинных джозефсоновских переходах. Проблема вычисления вольт-амперных характеристик. Оценки скорости роста ошибок округления для разностной схемы второго порядка точности // Ж. вычисл. матем. и матем. физ. 2020. Т. 60. № 1. С. 159–166.
Урм В.Я. О необходимых и достаточных условиях устойчивости систем разностных уравнений // Докл. АН СССР. 1961. Т. 139. № 1. С. 40–43.
Ланкастер П. Теория матриц. Гл. 6. М.: Физматлит, 1982. C. 179.
Фихтенгольц Г.М. Курс дифференциального и интегрального исчисления. Т. II. М.: Физматлит, 1959. C. 67.
REDUCE User’s Guide for Unix Systems. Ver. 8 by Winfried Neun ZIB, 2004.
IT-ecosystem of the HybriLIT heterogeneous platform for high-performance computing and training of IT-specialists. Selected Papers of the 8th International Conference “Distributed Computing and Grid-technologies in Science and Education” (GRID 2018), Dubna, Russia, September 10–14, 2018, CEUR-WS.org/Vol-2267. by Gh. Adam, M. Bashashin, D. Belyakov et al.

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

Инструменты

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

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

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