A fast normal splitting preconditioner for attractive coupled nonlinear Schroedinger equations with fractional Laplacian

全文:

Введение. Одним из самых популярных методов расчета в квантовой механике [1, 2] является интеграл Фейнмана по броуновским путям, который приводит к стандартным уравнениям Шредингера. Если вместо броуновского движения рассматривать полет Леви, то возникает система дробных уравнений Шредингера (fractional Schroedinger equations, FSE) [3]. Вместо лапласиана в стандартных уравнениях Шредингера [4] FSE включает дробную производную порядка α (1 < α < 2). При α = 2 FSE сводится к стандартному случаю. Кроме того, FSE имеет важные приложения в квантовой механике, полупроводниках и других областях [5]. Из-за нелокальной природы оператора дробного дифференциала точное решение FSE получить трудно. Для изучения природы FSE было разработано большое количество численных методов, таких, как методы конечных элементов [6, 7], спектральные методы [8, 9], методы коллокации [10, 11], методы конечных разностей [12–16] и пр.

В статье рассматриваются следующие одномерные (1D) пространственные дробно-связанные нелинейные уравнения Шредингера (coupled nonlinear Schroedinger, CNLS):

$\left\{\begin{array}{l}r{u}_t-\gamma {\left(-\Delta \right)}^{\frac{\alpha }{2}}u+\rho \left({\left|u\right|}^2+\beta {\left|v\right|}^2\right)u=\mathrm{0,}\\ r{u}_t-\gamma {\left(-\Delta \right)}^{\frac{\alpha }{2}}v+\rho \left({\left|v\right|}^2+\beta {\left|u\right|}^2\right)u=\mathrm{0,}\end{array}\right.x\in ℝ,0<t\le T,$  (0.1)

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

$u\left(x\mathrm{,0}\right)=u_0\left(x\right),u\left(x\mathrm{,0}\right)=v_0\left(x\right),x\in ℝ,$(0.2)

где 1 = $\sqrt{-1}$, параметры γ > 0, β ≥ 0, ρ – вещественные константы, дробный лапласиан

${\left(-\Delta \right)}^{\frac{\alpha }{2}}=-\frac{{\partial }^{\alpha }}{\partial {\left|x\right|}^{\alpha }}$  (0.3)

определяется дробной производной одномерного пространства Рисса [17, 18] при 1 < α < 2. При параметре β = 0 можно получить несвязанные нелинейные уравнения Шредингера (decoupled nonlinear Schroedinger, DNLS) [19]. Что касается параметра ρ, то он определяет три различных варианта нелинейного члена FSE. При ρ = 0 нелинейные члены исчезают и система (0.1) описывает свободные частицы [17, 20]. При ρ < 0 нелинейные члены в системе (0.1) представляют собой отталкивающее взаимодействие частиц [21–23]. При ρ > 0 нелинейные члены в системе (0.1) представляют собой притягивающее взаимодействие частиц [21, 24, 25]. В данной статье рассматривается только ρ > 0, в этом случае матрица коэффициентов дискретизированной линейной системы является комплексной симметричной неопределенной. Случай ρ < 0, приводящий к сложным симметричным отрицательно определенным линейным системам, обсуждается в [26–30].

Теоретически решение FSE сохраняет массу и энергию, например система (0.1). В частности, решения u(x,t),v(x,t) задачи (0.1) удовлетворяют сохранению массы [13]:

$\begin{array}{l}{||u\left(\cdot ,t\right)||}_{L_2}^2=\underset{ℝ}{\int }{\left|u\left(x,t\right)\right|}^{2}dx={||u\left(\cdot ,0\right)||}_{L_2}^2,\\ {||v\left(\cdot ,t\right)||}_{L_2}^2=\underset{ℝ}{\int }{\left|v\left(x,t\right)\right|}^{2}dx={||u\left(\cdot ,0\right)||}_{L_2}^2,\end{array}$(0.4)

и сохранению энергии:

$E\left(t\right)=\frac{\gamma }{2}\underset{ℝ}{\int }\left(\overline{u}{\left(-\Delta \right)}^{\frac{\alpha }{2}}u+\overline{v}{\left(-\Delta \right)}^{\frac{\alpha }{2}}v\right)dx-\frac{\rho }{4}\underset{ℝ}{\int }\left({\left|u\right|}^4+{\left|v\right|}^4+2\beta {\left|u\right|}^2{\left|v\right|}^2\right)dx$, (0.5)

где функционал энергии удовлетворяет E(t) = E(0). Вышеуказанные законы сохранения остаются в силе и при дискретизации системы (0.1).

В этой статье вместо всей вещественной оси ℝ рассмотрение ведется на ограниченном интервале a ≤ x ≤ b и приняты граничные условия Дирихле:

$u\left(a,t\right)=u\left(b,t\right)=\mathrm{0,}v\left(a,t\right)=v\left(b,t\right)=\mathrm{0,}0\le t\le T.$

Для дискретизации дробных пространственных уравнений CNLS (0.1) используется схема линейно-неявной консервативной разности (linearly implicit conservative difference, LICD), предложенная в [13]. На каждом временном уровне необходимо решать сложные симметричные линейные системы с теплиц-плюс-диагональной структурой (D − T + гI) u = b, где D – диагональная, а T – теплицева симметричная положительно определенная матрица. При ρ > 0 матрица D положительно полуопределена, а D − T неопределена. Тогда D − T + гI является комплексно-симметричной неопределенной.

Существует множество численных методов для работы со сложными симметричными линейными системами. Первый класс методов – это методы неявной итерации с чередующимся направлением [31–35]. В [31, 32] были предложены итерационные методы эрмитова и косоэрмитова разложения (skew-hermitian splitting, HSS) для работы с общими неэрмитовыми положительно определенными линейными системами. Кроме того в этих публикациях были разработаны модифицированные итерационные методы HSS (modified HSS, MHSS) [33] и предварительно обусловленные итерационные методы MHSS (preconditioned MHSS, PMHSS) [34, 35] с учетом симметричной структуры матрицы системы. Второй класс методов – итерационные C-to-R [36], которые требуют высококачественных предобусловливателей для дополнения Шура для обеспечения устойчивости. Построение дополнения Шура затруднено, поскольку матрица дополнения Шура плотна и имеет вид S +${\tilde{S}}^{-1}$, где S, $\tilde{S}$ ∈ ℝ^M×M – плотные матрицы с теплиц-плюс-диагональной структурой. К сожалению, методы типа HSS и итерации C-to-R не могут быть применены для прямой обработки комплексной симметричной неопределенной матрицы D − T + гI, поскольку они предназначены для неэрмитовой положительно определенной матрицы.

Для линейных систем с теплиц-плюс-диагональной структурой быстрые прямые решатели не разработаны. Следовательно, лучшим вариантом является использование предварительно обусловленных методов итерации подпространства Крылова [37, 38]. Реализация умножения матрицы на вектор и решение линейной системы с обобщенной невязкой (generalized residual, GR) всегда необходимы на каждой итерации предобусловливаемых методов итерации подпространства Крылова. Если существуют быстрые реализации умножения матрицы на вектор и можно эффективно сконструировать и реализовать высококачественные предобусловливатели, итерации подпространства Крылова будут работать очень хорошо. Многие эффективные предобусловливатели работают с эрмитовыми положительно определенными теплиц-плюс-диагональными системами. Это, например, класс полосовых (banded) предобусловливателей для эрмитовых систем «теплиц-плюс-полоса» [39]; класс приближенных предобусловливателей обратного циркулянта-плюс-диагонали (approximate inverse circulant-plus-diagonal, AICD) [40]; предобусловливатели диагонального и теплицевого разложения (diagonal and toeplitz splitting, DTS), а также предобусловливатели диагонального и циркулянтного разложения (diagonal and circulant splitting, DCS), предложенные в [41] (1D-случай) и [42] (многомерный случай). Тем не менее, никакие предобусловливатели, описанные выше, не могут быть напрямую применены для решения сложных линейных систем (D − T + гI) u = b, поскольку они предназначены для эрмитовых положительно определенных систем, но не для неопределенных симметричных теплиц-плюс-диагональных.

Для разработки эффективных итерационных методов рассматривается вещественная блочная (два на два) форма линейной системы (D − T + гI) u = b. Матрица блочной системы раскладывается на нормальную блочную матрицу и антисимметричную блочную матрицу и строится класс итерационных методов нормального и антисимметричного разложения (normal and anti-symmetric splitting, NASS), основанный на подходе неявных итерационных методов чередующегося направления (alternating direction implicit, ADI) [43, 44]. На каждом шаге итерационного метода NASS используются две линейные подсистемы с матрицами коэффициентов:

$\left[\begin{array}{cc}\left(\omega +1\right)I& T\\ -T& \left(\omega +1\right)I\end{array}\right],\left[\begin{array}{cc}\omega I& -D\\ D& \omega I\end{array}\right],$ (0.6)

которые необходимо решить. Первое можно решить итеративно с помощью метода предобусловленного обобщенного метода минимальной невязки (generalized minimal residual method, GMRES) с блочным предобусловливателем:

$\left[\begin{array}{cc}\left(\omega +1\right)I& C\\ -C& \left(\omega +1\right)I\end{array}\right]$ (0.7)

с циркулянтной аппроксимацией C. Каждая итерация предобусловленного GMRES может быть выполнена за 𝒪(M log M) операций с плавающей точкой (флоп, flops), если применяется быстрый алгоритм, например быстрое преобразование Фурье (БПФ) [37]. Второе можно решить непосредственно с помощью разреженного исключения Гаусса (Gaussian elimination, GE) за 𝒪(M) флоп. Теоретически итерационный метод NASS безусловно сходится к единственному решению. Получена верхняя оценка коэффициента сжатия итерационного метода NASS, зависящая только от спектров симметричной положительно определенной матрицы Теплица T ∈ ℝ^M×M. Оптимальное значение параметра итерации, минимизирующее верхнюю границу коэффициента сжатия итерационного метода NASS, определяется нижней и верхней границами собственных значений T.

Класс предобусловливателей разложения матрицы, называемый предобусловливателями NASS, может быть получен естественным путем из итерационного метода NASS. Чтобы уменьшить вычислительную сложность предобусловливателей NASS, для замены блоков Теплица используются циркулянтные аппроксимации, что приводит к созданию класса более экономичных циркулянтных улучшенных нормальных и антисимметричных предобусловливателей (circulant improved normal and anti-symmetric, CNAS). Реализация предобусловливателя CNAS требует разрешения двух линейных подсистем. Матрицы подсистем имеют вид

$\left[\begin{array}{cc}\left(\omega +1\right)I& C\\ -C& \left(\omega +1\right)I\end{array}\right],\left[\begin{array}{cc}\omega I& -D\\ D& \omega I\end{array}\right].$ (0.8)

Первое можно решить за 𝒪(M log M) флоп на основе БПФ, а второе можно решить непосредственно с помощью разреженного GE за 𝒪(M) флоп. Теоретический анализ показывает, что собственные значения матрицы предварительно обусловленной системы NASS группируются около 1. Кроме того, можно доказать, что предварительное обусловливание системы CNAS представляет собой небольшое возмущение нормы и коррекцию низкого ранга матрицы предварительно обусловленной системы NASS. Это значит, что большинство собственных значений матрицы предобусловленной системы CNAS сгруппированы вокруг аналогов для NASS, т.е. могут быть расположены около 1. Численные эксперименты показывают, что предобусловливатель CNAS значительно повышает вычислительную эффективность методов итерации подпространства Крылова (таких, как GMRES), а соответствующий предобусловленный метод (PGMRES) ведет себя независимо от размера пространственной сетки и нечувствителен к дробному порядку. Кроме того, замечено, что дискретная масса и энергия сохраняются по схеме LICD в сочетании с предварительно обусловленным методом CNAS GMRES с точки зрения заданной точности.

Изложение построено следующим образом. В разд. 1 схема LICD применяется для дискретизации дробно-пространственных уравнений CNLS, что приводит к неопределенной комплексной симметричной линейной системе (D − T + гI) u = b. В разд. 2 представлен итерационный метод NASS и установлена соответствующая асимптотическая теория сходимости. В разд. 3 создаются и анализируются новые предобусловливатели. В разд. 4 рассматриваются и обсуждаются реализации и вычислительные сложности новых предобусловливателей. В разд. 5 приведены числовые результаты и даются заключительные замечания.

1. Дискретизация и линейная система. Дробные по пространству уравнения CNLS (0.1) дискретизируются по схеме LICD [13]. Пусть N и M – целые положительные числа, через τ = T∕N и h = (b − a)∕(M + 1) заданы размеры временных шагов и пространственных сеток соответственно. Тогда t_n = nτ для n = 0,N и x_j = a + jh для j = 0,M + 1 представляют временные уровни и узлы пространственных ячеек.

Пусть u_jⁿ ≈ u(x_j,t_n) и v_jⁿ ≈ v(x_j,tⁿ) – аппроксимации сеточной функции. Одномерный дробный лапласиан (−Δ)можно дискретизировать дробной центрированной разностью [45, 46] в ограниченном пространственном интервале [a, b] как

${\left(-\Delta \right)}^{\frac{\alpha }{2}}u\left(x_j,t\right)=\frac{1}{h^{\alpha }}{\Delta }_h^{\alpha }{u}_j+\mathcal{O}\left(h^2\right).$

Здесь 

${\Delta }_h^{\alpha }{u}_j=\sum _{k=1}^M{c}_{j-k}{u}_k$ (1.1)

с коэффициентами $c_k={\left(-1\right)}^k\Gamma \left(\alpha +1\right)/\left[\Gamma \left(\alpha /2-k+1\right)\Gamma \left(\alpha /2+k+1\right)\right],\forall k\in \mathbb{Z},$

где Γ(⋅) – гамма-функция. Кроме того, коэффициенты c_k удовлетворяют следующим свойствам [45]:

$c_0\ge \mathrm{0,}{c}_k=c_{-k}\le 0\left(\forall k\ge 1\right),\sum _{k=-\infty ,k\ne 0}^{+\infty }\left|c_k\right|=c_0.$

Схема LICD [13] для системы (0.1) выглядит следующим образом:

$\left\{\begin{array}{l}l\frac{u_j^{n+1}-u_j^{n-1}}{2\tau }-\frac{\gamma }{h^{\alpha }}{\Delta }_h^{\alpha }{\widehat{u}}_j^n+\rho \left({\left|u_j^n\right|}^2+\beta {\left|v_j^n\right|}^2\right){\widehat{u}}_j^n=\mathrm{0,}\\ l\frac{v_j^{n+1}-v_j^{n-1}}{2\tau }-\frac{\gamma }{h^{\alpha }}{\Delta }_h^{\alpha }{\widehat{v}}_j^n+\rho \left({\left|v_j^n\right|}^2+\beta {\left|u_j^n\right|}^2\right){\widehat{v}}_j^n=\mathrm{0,}\end{array}\right.$ (1.2)

где û_jⁿ = (u_jⁿ⁺¹ + u_jⁿ⁻¹)∕2, û_jⁿ = (v_jⁿ⁺¹ + v_jⁿ⁻¹)∕2, для j = 1, 2,…,M,n = 1, 2,…,N − 1. Начальные и граничные условия: u_j⁰ = u₀(x_j), v_j⁰ = v₀(x_j), u₀ⁿ = u_M+1ⁿ = 0, v₀ⁿ = v_M+1ⁿ = 0. Эквивалентно первая строка (1.2) приводит к сложной линейной системе:

$A^{n+1}{u}^{n+1}=b^{n+1},\forall n\ge \mathrm{1,}$ (1.3)

где Aⁿ⁺¹ = Dⁿ⁺¹ − T + гI ∈ ℂ^M×M, Dⁿ⁺¹ = diag{d₁ⁿ⁺¹,d₂ⁿ⁺¹,…,d_Mⁿ⁺¹} ∈ ℝ^M×M диагонально с d_jⁿ⁺¹= = ρτ(|u_jⁿ|² + β|v_jⁿ|²), I ∈ ℝ^M×M – единичная матрица, а T ∈ ℝ^M×M является теплицевой:

$T=\mu \left[\begin{array}{ccccc}{c}_0& c_{-1}& \cdots & c_{2-M}& c_{1-M}\\ c_1& c_0& \ddots & \ddots & c_{2-M}\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ c_{M-2}& \ddots & \ddots & c_0& c_{-1}\\ c_{M-1}& c_{M-2}& \cdots & c_1& c_0\end{array}\right]$ (1.4)

с μ = γτ∕h^α. Очевидно, что вторая строка (1.2) также допускает связанную линейную систему, такую же, как (1.3).

Ввиду того, что ρ > 0, γ > 0, β ≥ 0 и свойств коэффициентов c_k, Dⁿ⁺¹ – неотрицательная диагональная матрица, а T строго диагонально доминантна и симметрична. Следовательно, Dⁿ⁺¹ − T симметричная неопределенная, а Aⁿ⁺¹ комплексно симметричная и неопределенная.

В [13] доказано, что дискретизированные связанные линейные системы (1.2) сохраняют дискретную массу и энергию, т.е.

$Q_1^n=Q_1^0,Q_2^n=Q_2^0,E^n=E^0,1\le n\le N,$ (1.5)

где Q₁ⁿ = (∥uⁿ⁺¹∥² + ∥uⁿ∥²)∕2 и Q₂ⁿ = (∥vⁿ⁺¹∥² + ∥vⁿ∥²)∕2 с ∥u∥² = ⟨u,u⟩ и

$〈u,v〉=h\sum _{k=1}^{M-1}{u}_j{\overline{v}}_j,$

$\begin{array}{c}{E}^n=\frac{\gamma n}{4h^{\alpha }}\sum _{j=1}^{M-1}\left({\overline{u}}_j^{n+1}{\Delta }_h^{\alpha }{u}_j^{n+1}+{\overline{u}}_j^n{\Delta }_h^{\alpha }{u}_j^n+{\overline{v}}_j^{n+1}{\Delta }_h^{\alpha }{v}_j^{n+1}+{\overline{v}}_j^n{\Delta }_h^{\alpha }{v}_j^n\right)-\\ -\frac{\rho h}{4}\sum _{j=1}^{M-1}\left[\left({\left|u_j^n\right|}^2{\left|u_j^{n+1}\right|}^2+{\left|v_j^n\right|}^2{\left|v_j^{n+1}\right|}^2\right)+\beta {\left|u_j^n\right|}^2{\left|v_j^{n+1}\right|}^2+{\left|v_j^n\right|}^2{\left|u_j^{n+1}\right|}^2\right].\end{array}$ (1.6)

2. Итеративный метод NASS. Рассматривается построение итерационного метода для сложных линейных систем вида (1.3). Если игнорировать верхние индексы, он записывается как

$Au=b,$(2.1)

где A = D − T + гI ∈ ℂ^M×M комплексно-симметричная и неопределенная, причем D ∈ ℝ^M×M – положительная полуопределенная диагональная матрица и T ∈ ℝ^M×M – симметричная положительно определенная теплицевая матрица. Комплексный неизвестный вектор имеет вид u = y + гz ∈ ℂ^M, а комплексный вектор b = p + гq ∈ ℂ^M – правая часть, где y, z, p, q ∈ ℝ^M – вещественные векторы. Затем сложную линейную систему (2.1) можно преобразовать в следующую вещественную блочную линейную систему:

$\widehat{\mathcal{R}}\widehat{x}\equiv \left[\begin{array}{cc}D-T& -I\\ I& D-T\end{array}\right]\left[\begin{array}{c}y\\ z\end{array}\right]=\left[\begin{array}{c}p\\ q\end{array}\right]\equiv \widehat{f}$ (2.2)

где $\widehat{\mathcal{R}}$ ∈ ℝ^2M×2M несимметрично и неопределенно.

Определим перестановки блоков

$\mathcal{P}=\left[\begin{array}{cc}-I& 0\\ 0& I\end{array}\right],\mathcal{Q}=\left[\begin{array}{cc}0& I\\ I& 0\end{array}\right]$ (2.3)

и пусть $\mathcal{R}=\mathcal{P}\widehat{\mathcal{R}}\mathcal{Q},x=\mathcal{Q}\widehat{x},f=\mathcal{P}\widehat{\widehat{f}}.$

Тогда вещественная блочная линейная система (2.2) эквивалентна вещественной несимметричной положительно определенной блочной линейной системе следующим образом:

$\mathcal{R}x\equiv \left[\begin{array}{cc}I& T-D\\ D-T& I\end{array}\right]\left[\begin{array}{c}z\\ y\end{array}\right]=\left[\begin{array}{c}-p\\ q\end{array}\right]\equiv f.$ (2.4)

Матрица ℛ ∈ ℝ^2M×2M допускает следующее нормальное и антисимметричное (NAS) разложение:

$\mathcal{R}=\mathcal{T}+\mathcal{D},$ (2.5)

$\mathcal{T}=\left[\begin{array}{cc}I& T\\ -T& I\end{array}\right],\mathcal{D}=\left[\begin{array}{cc}0& -D\\ D& 0\end{array}\right].$ (2.6)

Основываясь на разложении NAS (2.5) и духе неявной итерации с чередующимся направлением (ADI) [43, 44], можно построить метод итерации NASS следующим образом.

Итерационный метод NASS. Пусть x⁽⁰⁾ ∈ ℝ^2M – произвольное начальное предположение. Для k = 0, 1, 2,…, пока последовательность итераций {x^(k)} k≥0 не сойдется, вычислить следующую итерацию x^(k+1), согласно следующей процедуре:

$\left\{\begin{array}{l}\left(\omega I+\mathcal{T}\right)x^{\left(k+\frac{1}{2}\right)}=\left(\omega I-\mathcal{D}\right)x^{\left(k\right)}+f,\\ \left(\omega I+\mathcal{D}\right)x^{\left(k+1\right)}=\left(\omega I-\mathcal{T}\right)x^{\left(k+\frac{1}{2}\right)}+f,\end{array}\right.$(2.7)

где ω > 0 – произвольный параметр итерации.

Итерационная схема (2.7) может быть переформулирована для итерации с фиксированной точкой как x^(k+1) = ℒ_ωx^(k) + ℱ_ω⁻¹f, где матрица итераций это

$L_{\omega }={\mathcal{F}}_{\omega }^{-1}{\mathcal{G}}_{\omega }.$ (2.8)

Здесь

${\mathcal{F}}_{\omega }=\frac{1}{2\omega }\left(\omega I+\mathcal{T}\right)\left(\omega I+\mathcal{D}\right)$ (2.9)

и

${\mathcal{G}}_{\omega }=\frac{1}{2\omega }\left(\omega I-\mathcal{T}\right)\left(\omega I-\mathcal{D}\right).$

Действительно, матрицы ℱ_ω и образуют разложение ℛ, т.е.

${\mathcal{R}}_{\omega }={\mathcal{F}}_{\omega }-{\mathcal{G}}_{\omega }.$ (2.10)

В следующей теореме обеспечивается свойство сходимости итерационного метода NASS.

Теорема 1. Пусть ℛ ∈ ℝ^2M×2M – матрица системы (2.4), 𝒯 ∈ ℝ^2M×2M и 𝒟∈ ℝ^2M×2M – блочные матрицы, образующие разбиение ℛ в (2.5). Если ω – положительная константа, то итерационный метод NASS (2.7) корректно определен, а матрица итераций ℒ_ω задается (2.8). Спектральный радиус ρ(ℒ_ω) ограничен:

$\sigma \left(\omega \right)=\underset{{\lambda }_i\in \lambda \left(T\right)}{\max }\sqrt{\frac{{\left(\omega -1\right)}^2+{\lambda }_i^2}{{\left(\omega +1\right)}^2+{\lambda }_i^2}},$ (2.11)

где λ(T) – спектр T. Таким образом, можно записать

$\rho \left({\mathcal{L}}_{\omega }\right)\le \sigma \left(\omega \right)<\mathrm{1,}\forall \omega >\mathrm{0,}$ (2.12)

т. е. итерация NASS сходится к единственному решению блочной линейной системы (2.4).

Доказательство. Поскольку 𝒯 положительно определен, 𝒟 антисимметричен и ω > 0, то ωI + 𝒯 и ωI + 𝒟 обратимы и положительно определены. Таким образом, итерация NASS (2.7) четко определена. Заметим, что итерационная матрица ℒ_ω удовлетворяет условию ℒ_ω = = (ωI + 𝒟)⁻¹𝒰_ω𝒱_ω(ωI + 𝒟), где 𝒰_ω = (ωI + 𝒯 )⁻¹(ωI −𝒯) и 𝒱_omega = (ωI −𝒟)(ωI + 𝒟)⁻¹, можно установить следующее соотношение:

$\rho \left({\mathcal{L}}_{\omega }\right)=\rho \left({\mathcal{U}}_{\omega }{\mathcal{V}}_{\omega }\right)\le {||{\mathcal{U}}_{\omega }||}_2{||{\mathcal{V}}_{\omega }||}_2.$ (2.13)

Теперь сосредоточимся на оценке ∥𝒰_ω∥₂ и ∥𝒱_ω∥₂. Во-первых, легко проверить, что 𝒯^⊤𝒯 = 𝒯𝒯^⊤, тогда имеем

$\begin{array}{c}{{\mathcal{U}}_{\omega }}^{\top }{\mathcal{V}}_{\omega }={\left(\omega I-\mathcal{T}\right)}^{\top }{{\left(\omega I+\mathcal{T}\right)}^{-}}^{\top }{\left(\omega I+\mathcal{T}\right)}^{-1}\left(\omega I-\mathcal{T}\right)=\\ ={{\left(\omega I+\mathcal{T}\right)}^{-}}^{\top }{\left(\omega I+\mathcal{T}\right)}^{-1}{\left(\omega I-\mathcal{T}\right)}^{\top }\left(\omega I-\mathcal{T}\right)=\\ ={\left[\left(\omega I+\mathcal{T}\right){\left(\omega I+\mathcal{T}\right)}^{\top }\right]}^{-1}{\left(\omega I-\mathcal{T}\right)}^{\top }\left(\omega I-\mathcal{T}\right)=\\ ={\left[\begin{array}{cc}{\left(\omega +1\right)}^{2}I+T^2& 0\\ 0& {\left(\omega +1\right)}^{2}I+T^2\end{array}\right]}_{-1}\left[\begin{array}{cc}{\left(\omega -1\right)}^{2}I+T^2& 0\\ 0& {\left(\omega -1\right)}^{2}I+T^2\end{array}\right].\end{array}$ (2.14)

Ввиду того, что T симметрична положительно определена, справедливо равенство

${||{\mathcal{U}}_{\omega }||}_2=\underset{{\lambda }_i\in \lambda \left(T\right)}{\max }\sqrt{\frac{{\left(\omega -1\right)}^2+{\lambda }_i^2}{{\left(\omega +1\right)}^2+{\lambda }_i^2}}<1.$ (2.15)

Во-вторых, рассмотрим оценку ∥𝒱_ω∥₂. Поскольку 𝒟^⊤𝒟 = 𝒟𝒟^⊤, то следует, что

$\begin{array}{c}{\mathcal{V}}_{\omega }^{{}^{\top }}{\mathcal{V}}_{\omega }={\left(\omega I+{\mathcal{D}}^{\top }\right)}^{-1}\left(\omega I-{\mathcal{D}}^{\top }\right)\left(\omega I-\mathcal{D}\right){\left(\omega I+\mathcal{D}\right)}^{-1}=\\ =\left(\omega I-{\mathcal{D}}^{\top }\right)\left(\omega I-\mathcal{D}\right){\left[\left(\omega I+\mathcal{D}\right)\left(\omega I+{\mathcal{D}}^{\top }\right)\right]}^{-1}=\\ =\left[\begin{array}{cc}{\omega }^{2}I+D^2& 0\\ 0& {\omega }^{2}I+D^2\end{array}\right]{\left[\begin{array}{cc}{\omega }^{2}I+D^2& 0\\ 0& {\omega }^{2}I+D^2\end{array}\right]}^{-1}=I.\end{array}$

Тогда

${||{\mathcal{V}}_{\omega }||}_2=1.$ (2.16)

На основании (2.13)−(2.16) можно получить оценку (2.12).

Замечание 1. Для дальнейшего понимания теоретического результата, изложенного в приведенной выше теореме, сделаем следующие замечания.

1. Скорость сходимости итерации NASS ограничена величиной σ(ω) и зависит только от спектра теплицевой матрицы T.
2. Если ввести взвешенную векторную норму ∥x∥_𝒟 = ∥(ωI + 𝒟)x∥₂ для x ∈ ℝ^2M и соответствующую норму индуцированной взвешенной матрицы ∥𝒳∥_𝒟 = ∥(ωI + 𝒟)𝒳(ωI + 𝒟)⁻¹∥₂ для 𝒳 ∈ ℝ^2M×2M имеет место равенство ${||{\mathcal{L}}_{\omega }||}_{\mathcal{D}}={||{\mathcal{U}}_{\omega }{\mathcal{V}}_{\omega }||}_2<\sigma \left(\omega \right),\forall \omega >0.$

Более того, на основе формы итерации с фиксированной точкой итерации NASS можно проверить, что последовательность итераций {x^(k)} k ≥ 0 удовлетворяет условию ${||x^{\left(k+1\right)}-x^{\ast }||}_{\mathcal{D}}\le \sigma \left(\omega \right){||x^{\left(k\right)}-x^{\ast }||}_{\mathcal{D}},\forall k\ge 0.$

Следовательно, σ(ω) также служит верхней границей коэффициента сжатия итерации NASS в смысле ∥⋅∥_𝒟-нормы. Замечено, что при TD = DT выполняется 𝒰_ω𝒱_ω = 𝒱_ω𝒰_ω, что приводит к факту ρ(ℒ_ω) = ∥ℒ_ω∥_𝒟 = σ(ω). Тогда эти величины можно минимизировать при том же оптимальном значении ω.

3. Верхняя граница σ(ω) может быть дополнительно ограничена величиной (ω), т.е. $\sigma \left(\omega \right)\le \underset{{\lambda }_{\min }\le \lambda \le {\lambda }_{\max }}{\max }\sqrt{\frac{{\left(\omega -1\right)}^2+{\lambda }_i^2}{{\left(\omega +1\right)}^2+{\lambda }_i^2}}\equiv \widehat{\sigma }\left(\omega \right),$

где λ_min > 0 и λ_max > 0 – минимальное и максимальное собственные значения T. Поскольку функция g(ω; λ) = $\sqrt{\left[{\left(\omega -1\right)}^2+{\lambda }_{}^2\right]/\left[{\left(\omega +1\right)}^2+{\lambda }_{}^2\right]}$ монотонно возрастает при λ > 0, она имеет вид

$\widehat{\sigma }\left(\omega \right)=\sqrt{\frac{{\left(\omega -1\right)}^2+{\lambda }_{\max }^2}{{\left(\omega +1\right)}^2+{\lambda }_{\max }^2}}.$

Пусть

$\frac{\left(\widehat{\sigma }\left(\omega \right)\right)}{\omega }=\frac{2\left({\omega }^2-1-{\lambda }_{max}^2\right)}{{\left(\omega -1\right)}^2+{\lambda }_{max}^2 {}^{{\displaystyle \frac{1}{ 2 }} }{\left(\omega +1\right)}^{2 }+{\lambda }_{max}^2 {}^{\frac{3}{ 2 } }}=0$  (2.17)

тогда $\widehat{\sigma }\left(\omega \right)$ минимизируется в ω⋆ = $\sqrt{{\lambda }_{\max }^2+1}$. Если применяется параметр ω = ω⋆, скорость сходимости итерации NASS для решения блочной линейной системы (2.4) ограничивается следующим образом:

$\sigma \left({\omega }^{\ast }\right)\le \widehat{\sigma }\left({\omega }^{\ast }\right)=\frac{{\lambda }_{\max }}{1+\sqrt{{\lambda }_{\max }^2+1}}.$

3. Предварительное обусловливание. Итерационный метод NASS естественным образом вызывает матрицу предварительной обработки ℱ_ω, определенную в (2.9) для матрицы ℛ в (2.4). Соответствующая предобусловленная линейная система имеет вид

${ℱ}_{\omega }^{-1}ℛx={ℱ}_{\omega }^{-1}f$(3.1)

предобусловливатель ℱ_ω называется предобусловливателем NASS. Очевидно, ℱ_ω есть произведение блочных нормальных матриц ωI + 𝒯, ωI + 𝒟 и скалярного коэффициента 1∕(2ω).

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

${ℱ}_{\omega }^{}{z}^{\left(k\right)}=r^{\left(k\right)}, \forall k\ge$, (3.2)

где r^(k) – текущий вектор невязки, а z^(k) – вектор обобщенной невязки. Разрешение (3.2) состоит из решения линейных подсистем с матрицами коэффициентов ωI + 𝒯 и ωI + 𝒟. Чтобы снизить стоимость процесса предварительной обработки, рассматриваются циркулянтные аппроксимации для замены матрицы Теплица T в ℱ_ω [47–49]. Теоретическое рассматрение дается только для циркулянта Стрэнга C для T [47]. Для простоты анализа пусть M четное, в дальнейшем будем записывать как C = μ[ς_ij] ∈ ℝ^M×M с ς_ij = ς_i−j, ς_k = ς_−k для 0 ≤ k < M, ς= 0 и

${\varsigma }_k=\left\{\begin{array}{c}{c}_k,\text{\hspace{1em}if\hspace{0.33em}0}\le k<\frac{M}{2},\\ c_{M-k},\text{\hspace{1em}if\hspace{0.33em}}\frac{M}{2}<k<M,\end{array}\right.$

т.е.

$C=\mu \left[\begin{array}{cccccccc}{c}_0& c_{-1}& \cdots & c_{-\left(\frac{M}{2}-1\right)}& 0& c_{-\left(\frac{M}{2}-1\right)}& \cdots & c_{-1}\\ c_{-1}& c_0& \ddots & \ddots & c_{-\left(\frac{M}{2}-1\right)}& 0& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & c_{-\left(\frac{M}{2}-1\right)}\\ c_{-\left(\frac{M}{2}-1\right)}& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & 0\\ 0& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & c_{-\left(\frac{M}{2}-1\right)}\\ c_{-\left(\frac{M}{2}-1\right)}& 0& \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & c_{-1}\\ c_{-1}& c_{-2}& \cdots & 0& c_{-\left(\frac{M}{2}-1\right)}& \ddots & \cdots & c_0\end{array}\right],$  (3.3)

которая является вещественной симметричной. Затем строится улучшенный эффективный предобусловливатель на основе циркулянта ${\tilde{\mathcal{F}}}_{\omega }$:

${\tilde{\mathcal{F}}}_{\omega }=\frac{1}{2\omega }\left(\omega I+\mathcal{C}\right)\left(\omega I+\mathcal{D}\right),\mathcal{C}=\left[\begin{array}{cc}I& C\\ -C& I\end{array}\right].$

Очевидно, 𝒞 нормален, а 𝒟 антисимметричен. Таким образом, предобусловливатель ${\tilde{\mathcal{F}}}_{\omega }$ называется циркулянтным улучшенным нормальным и антисимметричным предобусловливателем (CNAS).

Теперь изучим свойство кластеризации собственных значений предварительно обусловленной матрицы системы ${\tilde{\mathcal{F}}}_{\omega }$⁻¹ℛ. Благодаря тому факту, что

${\tilde{\mathcal{F}}}_{\omega }^{-1}\mathcal{R}=\underset{}{\underbrace{{\tilde{\mathcal{F}}}_{\omega }^{-1}{\mathcal{F}}_{\omega }}}\underset{}{\underbrace{{\tilde{\mathcal{F}}}_{\omega }^{-1}\mathcal{R}}}$ (3.4)

это рассмотрение сводится к анализу свойств ℱ_ω⁻¹ℛ и ${\tilde{\mathcal{F}}}_{\omega }^{-1}$ ℱ_ω соответственно.

3.1. Cвойство кластеризации ℱ_ω⁻¹ℛ. Рассмотрим свойство кластеризации ℱ_ω⁻¹ℛ.

Теорема 2. Пусть ℛ ∈ ℝ^2M×2M – матрица системы в (2.4), 𝒯 ∈ ℝ^22M×2M и 𝒟 ∈ ℝ^22M×2M – блочные матрицы, образующие разбиение ℛ в (2.5), ω – положительная константа и σ(ω) определяется (2.11).

Когда предобусловливатель NASS ℱ_ω применяется к блочной линейной системе (2.4), собственные значения предобусловленной матрицы ℱ_ω⁻¹ℛ в (3.1) расположены в круге радиуса σ(ω) < 1 с центром в точке 1.

Доказательство. Согласно (2.8) и (2.10), оно записывается как ℒ_ω = I −ℱ_ω⁻¹ℛ. Таким образом, пусть η – собственное значение ℱ_ω⁻¹ℛ, тогда 1 − η – собственное значение ℒ_ω. Кроме того, теорема 1 показывает, что ρ(ℒ_ω) ≤ σ(ω), т. е. |1 −η| ≤ σ(ω) для всех η ∈ λ(ℱ_ω⁻¹ℛ), что и является результатом этой теоремы.

Теорема 3. Пусть ℛ ∈ ℝ^22M×2M – матрица системы в (2.4), 𝒯 ∈ ℝ^22M×2M и 𝒟∈ ℝ^22M×2M – блочные матрицы, образующие разбиение ℛ в (2.5), ω – положительная константа. Тогда собственные значения предвобусловленной матрицы ℱ_ω⁻¹ℛ группируются в точке 0₊ на правой полукомплексной плоскости, если ω стремится к положительной бесконечности.

Доказательство. Пусть λ – собственное значение ℒ_ω, а x – соответствующий собственный единичный вектор, тогда оно записывается как ℒ_ωx = λx, т.е. $\left(\omega I-\mathcal{T}\right)\left(\omega I-\mathcal{D}\right)x=\lambda \left(\omega I+\mathcal{T}\right)\left(\omega I+\mathcal{D}\right)x$.

Обозначим

$\widehat{\mathcal{T}}=\left[\begin{array}{cc}I& T\\ -T& I\end{array}\right].$ (3.5)

Если верно, что 𝒯 = I + $\widehat{\mathcal{T}}$, то имеем $\left[\left(\omega -1\right)I-\widehat{\mathcal{T}}\right]\left(\omega I-\mathcal{D}\right)x=\lambda \left[\left(\omega +1\right)I+\widehat{\mathcal{T}}\right]\left(\omega I+\mathcal{D}\right)x.$

Умножая x^∗ слева на обе части приведенного выше уравнения, получим, что

$\lambda =\frac{\omega \left(\omega -1\right)-\omega x^{\ast }\widehat{\mathcal{T}}x-\left(\omega -1\right)x^{\ast }\mathcal{D}x+x^{\ast }\widehat{\mathcal{T}}\mathcal{D}x}{\omega \left(\omega +1\right)+\omega x^{\ast }\widehat{\mathcal{T}}x+\left(\omega +1\right)x^{\ast }\mathcal{D}x+x^{\ast }\widehat{\mathcal{T}}\mathcal{D}x}.$

Поскольку x^∗$\widehat{\mathcal{T}}$x и x^∗𝒟 x являются чисто мнимыми числами в силу того, что $\widehat{\mathcal{T}}$ и 𝒟 вещественны и антисимметричны, запишем x^∗$\widehat{\mathcal{T}}$x = γ_Tı и x^∗𝒟x = γ_Dı с γ_T,γ_D ∈ ℝ. Кроме того, обозначим x^∗$\widehat{\mathcal{T}}$𝒟x = γ_TD^(R) + ıγ_TD^(I) с γ_TD^(R),γ_TD^(I)∈ ℝ. Из чего следует

$\lambda =\frac{\omega \left(\omega -1\right)+{\gamma }_{TD}^{\left(R\right)}+l\left[-\omega {\gamma }_T-\left(\omega -1\right){\gamma }_D+{\gamma }_{TD}^{\left(I\right)}\right]}{\omega \left(\omega -1\right)+{\gamma }_{TD}^{\left(R\right)}+l\left[\omega {\gamma }_T+\left(\omega -1\right){\gamma }_D+{\gamma }_{TD}^{\left(I\right)}\right]}.$

Поскольку ℱ_ω⁻¹ℛ = I −ℒ_ω, то 1-λ является собственным значением ℱ_ω⁻¹ℛ и получается

$1-\lambda =\frac{2}{\omega }\frac{1+l\left({\gamma }_T+{\gamma }_D\right)+\mathcal{O}\left(1/\omega \right)}{{\left[1+1/\omega +{\gamma }_{TD}^{\left(R\right)}/{\omega }^2\right]}^2+{\left[{\gamma }_T/\omega +\left(\omega +1\right){\gamma }_D/{\omega }^2+{\gamma }_{TD}^{\left(I\right)}/{\omega }^2\right]}^2}$

Очевидно, что вещественная часть 1-λ стремится к 0₊, а мнимая часть 1-λ – к 0, когда параметр ω стремится к положительной бесконечности.

3.2. Свойство $\tilde{\mathcal{F}}$_ω⁻¹ℱ_ω. Теперь обсудим свойство $\tilde{\mathcal{F}}$_ω⁻¹ℱ_ω. Для удобства определим следующие константы:

$\begin{array}{c}\theta =\frac{{\left(1-\frac{1+\alpha }{5+\alpha /2}\right)}^{5+\frac{\alpha }{2}}{e}^{1+\alpha }\Gamma \left(\alpha +1\right)\sin \left(\frac{\pi \alpha }{2}\right)}{\pi \alpha }\\ {\theta }_0=\frac{{\sqrt{2e}}^{13/12}\Gamma \left(\alpha +1\right)\sin \left(\frac{\pi \alpha }{2}\right)}{\pi \alpha }.\end{array}$ (3.6)

Далее вводятся две леммы следующего содержания.

Лемма 1. Эта и следующая леммы приведены в [30]. Пусть c_j = (−1)^jΓ(α + 1)∕[Γ(α∕2 − j + 1)Γ(α∕2 + j + 1)], k₀ ≥ 3, и 1 < α < 2, тогда

$\frac{\theta }{{\left(k_0+1/2\right)}^{\alpha }}<\sum _{j=k_0+1}\left|c_j\right|<\frac{{\theta }_0}{{\left(k_0-1\right)}^{\alpha }}$. (3.7)

На основе леммы 1 оценки собственных значений матрицы Теплица T и ее циркулянтной аппроксимации C представлены в следующей лемме.

Лемма 2. Пусть T – матрица Теплица в (1.4), C – циркулянт Стрэнга в (3.3) и M четно, тогда справедливо равенство

$\frac{2\gamma \tau \theta }{{\left(b-a\right)}^{\alpha }}<{\lambda }_T<\frac{2\gamma \tau }{h^{\alpha }}\left[\frac{\Gamma \left(\alpha +1\right)}{\Gamma \left(\alpha /2+1\right)}-\frac{\theta h^{\alpha }}{{\left(b-a\right)}^{\alpha }}\right],M\ge \mathrm{4,}$ (3.8)

и

$\frac{2^{\alpha +1}\gamma \tau \theta }{{\left(b-a\right)}^{\alpha }}<{\lambda }_C<\frac{2\gamma \tau }{h^{\alpha }}\left[\frac{\Gamma \left(\alpha +1\right)}{\Gamma {\left(\alpha /2+1\right)}^2}-\frac{2^{\alpha }\theta h^{\alpha }}{{\left(b-a\right)}^{\alpha }}\right],M\ge \mathrm{8,}$  (3.9)

где λ_T и λ_C – собственные значения T и C.

На основе лемм 1 и 2 свойство $\tilde{\mathcal{F}}$_ω⁻¹ℱ_ω равно суммированы в следующей теореме.

Теорема 4. Пусть 1 < α < 2 и M ≥ 8 четно, 𝜖 > 0 – небольшая константа и

$k_0=\left[{\left(\frac{\mu {\theta }_0}{ϵ}\right)}^{\frac{1}{\alpha }}\right]+\mathrm{1,}$ (3.10)

где ⌈.⌉ представляет собой округление вещественного числа до положительной бесконечности. Тогда существуют две матрицы $\tilde{\mathcal{E}}$ ∈ ℝ^2M×2M и $\widehat{\mathcal{F}}$ ∈ ℝ^2M×2M, удовлетворяющее rank $\left(\tilde{\mathcal{E}}\right)$= 4k₀:

${||\tilde{\mathcal{E}}||}_2<\frac{{\left({\omega }^2+v^2\right)}^{\frac{1}{2}}{M}^{\frac{1}{2}}\mu }{\omega \sqrt{{\left(\omega +1\right)}^2+{\left[\frac{2^{\alpha +1}\gamma \tau \theta }{{\left(b-a\right)}^{\alpha }}\right]}^2}}\left[\frac{c_0}{2}-\frac{\theta }{{\left(M-\frac{1}{2}\right)}^{\alpha }}\right],$ (3.11)

${||\tilde{\mathcal{F}}||}_2<\frac{{\left({\omega }^2+v^2\right)}^{\frac{1}{2}}{M}^{\frac{1}{2}}ϵ}{\omega \sqrt{{\left(\omega +1\right)}^2+{\left[\frac{2^{\alpha +1}\gamma \tau \theta }{{\left(b-a\right)}^{\alpha }}\right]}^2}},$ (3.12)

где ν = max_μi∈λ(D)|μ_i|, такой, что

${\tilde{\mathcal{F}}}_{\omega }^{-1}{\mathcal{F}}_{\omega }=I+\tilde{\mathcal{E}}+\tilde{\mathcal{F}}.$ (3.13)

Доказательство. Поскольку M четно, оно записывается как

$T-C=\mu \left[\begin{array}{ccc}0& {\widehat{F}}_{12}& {\widehat{E}}_{13}\\ {\widehat{\mathcal{F}}}_{12}^{\top }& 0& 0\\ {\widehat{\mathcal{F}}}_{13}^{\top }& 0& 0\end{array}\right],$

где $\widehat{F}$₁₂ ∈ ℝ^×(−k₀⁾ и $\widehat{E}$₁₃ ∈ ℝ^×k₀ имеют вид

$\begin{array}{l}{\widehat{F}}_{12}=\left[\begin{array}{cccc}{c}_{\frac{M}{2}}& c_{\frac{M}{2}+1}-c_{\frac{M}{2}-1}& \cdots & c_{M-\left(k_0+1\right)}-c_{k_0+1}\\ 0& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & c_{\frac{M}{2}+1}-c_{\frac{M}{2}-1}\\ ⋮& \ddots & \ddots & c_{\frac{M}{2}}\\ 0& \cdots & \cdots & 0\\ ⋮& & & ⋮\\ 0& \cdots & \cdots & 0\end{array}\right],\\ {\widehat{E}}_{13}=\left[\begin{array}{cccc}{c}_{M-k_0}-c_{k_0}& \cdots & \cdots & c_{M-1}-c_1\\ ⋮& \ddots & \ddots & ⋮\\ c_{\frac{M}{2}+1}-c_{\frac{M}{2}-1}& \ddots & \ddots & ⋮\\ c_{\frac{M}{2}}& \ddots & \ddots & c_{M-k_0}-c_{k_0}\\ 0& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & c_{\frac{M}{2}+1}-c_{\frac{M}{2}-1}\\ 0& \cdots & 0& c_{\frac{M}{2}}\end{array}\right].\end{array}$

Следовательно, выполняется

$T-C=\widehat{E}+\widehat{F},$  (3.14)

где

$\widehat{E}=\mu \left(\begin{array}{ccc}0& 0& {\widehat{E}}_{13}\\ 0& 0& 0\\ {\widehat{E}}_{13}^{\top }& 0& 0\end{array}\right),\widehat{F}=\mu \left(\begin{array}{ccc}0& {\widehat{F}}_{12}& 0\\ {\widehat{F}}_{12}^{\top }& 0& 0\\ 0& 0& 0\end{array}\right).$

Очевидно, оно записывается как rank $\widehat{E}$= 2k₀. Более того, ввиду лемм 1 и 2, а также структуры $\widehat{E}$ и $\widehat{F}$ можно получить следующие оценки:

$\begin{array}{c}{||\widehat{E}||}_{\infty }=\mu \max \left\{{||{\widehat{E}}_{13}||}_{\infty },{||{\widehat{E}}_{13}^{\top }||}_{\infty }\right\}=\mu {||{\widehat{E}}_{13}^{\top }||}_{\infty }\le \\ \le \mu \sum _{j=1}^{M-1}\left|c_j\right|=\mu \left(\frac{c_0}{2}-\sum _{j=M}^{\infty }\left|c_j\right|\right)<\\ <\mu \left[\frac{c_0}{2}-\frac{\theta }{{\left(M-\frac{1}{2}\right)}^{\alpha }}\right],\end{array}$ (3.15)

$\begin{array}{c}{||\widehat{F}||}_{\infty }=\mu \max \left\{{||{\widehat{F}}_{12}||}_{\infty },{||{\widehat{F}}_{12}^{\top }||}_{\infty }\right\}=\mu {||{\widehat{F}}_{12}||}_{\infty }\le \\ \le \mu \sum _{j=k_0+1}^{M-k_0-1}{c}_j<\mu \sum _{j=k_0+1}^{\infty }\left|c_j\right|<\\ <\frac{\mu {\theta }_0}{{\left(k_0-1\right)}^{\alpha }}<ϵ.\end{array}$  (3.16)

Можно проверить, что $\widehat{\mathcal{F}}$_ω⁻¹ℱ_ω − I и (ωI + 𝒞)⁻¹(𝒯 −𝒞) подобны, т.е.

$\begin{array}{c}{\tilde{\mathcal{F}}}_{\omega }^{-1}{\mathcal{F}}_{\omega }-I={\left(\omega I+\mathcal{D}\right)}^{-1}\left[{\left(\omega I+\mathcal{C}\right)}^{-1}\left(\omega I+\mathcal{T}\right)-I\right]\left(\omega I+\mathcal{D}\right)=\\ {\left(\omega I+\mathcal{D}\right)}^{-1}{\left(\omega I+\mathcal{C}\right)}^{-1}\left(\mathcal{T}-\mathcal{C}\right)\left(\omega I+\mathcal{D}\right),\end{array}$

причем

${\left(\omega I+\mathcal{C}\right)}^{-1}\left(\mathcal{T}-\mathcal{C}\right)={\left[\begin{array}{cc}\left(\omega +1\right)I& C\\ -C& \left(\omega +1\right)I\end{array}\right]}^{-1}\left[\begin{array}{cc}0& T-C\\ C-T& 0\end{array}\right].$  (3.17)

Учитывая (3.14) запишем ${\left(\omega I+\mathcal{C}\right)}^{-1}\left(\mathcal{T}-\mathcal{C}\right)=\widehat{\mathcal{E}}+\widehat{\mathcal{F}},$

где

$\begin{array}{c}\widehat{{\rm E}}={\left[\begin{array}{cc}\left(\omega +1\right)I& C\\ -C& \left(\omega +1\right)I\end{array}\right]}^{-1}\left[\begin{array}{cc}0& \widehat{E}\\ -\widehat{E}& 0\end{array}\right],\\ \widehat{\mathcal{F}}={\left[\begin{array}{cc}\left(\omega +1\right)I& C\\ -C& \left(\omega +1\right)I\end{array}\right]}^{-1}\left[\begin{array}{cc}0& \widehat{F}\\ -\widehat{F}& 0\end{array}\right].\end{array}$

Таким образом, выполняется ${\mathcal{F}}_{\omega }^{-1}{\mathcal{F}}_{\omega }-I=\tilde{\mathcal{E}}+\tilde{\mathcal{F}},$

где

$\tilde{\mathcal{E}}={\left(\omega I+\mathcal{D}\right)}^{-1}\tilde{\mathcal{E}}\left(\omega I+\mathcal{D}\right)$ (3.18)

и

$\tilde{\mathcal{F}}={\left(\omega I+\mathcal{D}\right)}^{-1}\tilde{\mathcal{F}}\left(\omega I+\mathcal{D}\right).$ (3.19)

Дополнительно имеем

$\begin{array}{c}{||\tilde{\mathcal{E}}||}_2\le {||{\left[\begin{array}{cc}\omega I& -D\\ D& \omega I\end{array}\right]}^{-1}||}_2{||\left[\begin{array}{cc}\omega I& -D\\ D& \omega I\end{array}\right]||}_2{||{\left[\begin{array}{cc}\left(\omega +1\right)I& C\\ -C& \left(\omega +1\right)I\end{array}\right]}^{-1}||}_2{||\widehat{E}||}_2=\\ =\frac{{\max }_{{\mu }_i\in \lambda \left(D\right)}\sqrt{{\omega }^2+{\mu }_i^2}}{{\min }_{{\mu }_i\in \lambda \left(D\right)}\sqrt{{\omega }^2+{\mu }_i^2}}\frac{1}{{\min }_{{\lambda }_i^{\left(c\right)}\in \lambda \left(C\right)}\sqrt{{\left(\omega +1\right)}^2+{\left({\lambda }_i^{\left(c\right)}\right)}^2}}{||\widehat{E}||}_2\le \\ \le \frac{{\max }_{{\mu }_i\in \lambda \left(D\right)}\sqrt{{\omega }^2+{\mu }_i^2}}{{\min }_{{\mu }_i\in \lambda \left(D\right)}\sqrt{{\omega }^2+{\mu }_i^2}}\frac{1}{{\min }_{{\lambda }_i^{\left(c\right)}\in \lambda \left(C\right)}\sqrt{{\left(\omega +1\right)}^2+{\left({\lambda }_i^{\left(c\right)}\right)}^2}}{M}^{\frac{1}{2}}{||\widehat{E}||}_{\infty }<\\ \end{array}$

$<\frac{{\left({\omega }^2+v^2\right)}^{\frac{1}{2}}{M}^{\frac{1}{2}}\mu }{\omega \sqrt{{\left(\omega +1\right)}^2+}{\left[\frac{2^{\alpha +1}\gamma \tau \theta }{{\left(b-a\right)}^{\alpha }}\right]}^2}\left[\frac{c_0}{2}-\frac{\theta }{{\left(M-\frac{1}{2}\right)}^{\alpha }}\right].$

Первое “=” связано с тем, что собственные значения

$\left[\begin{array}{cc}\omega I& -D\\ D& \omega I\end{array}\right]$ (3.20)

равны ω ± μ_iı с μ_i ∈ λ(D) и г = $\sqrt{-1}$, а собственные значения

$\left[\begin{array}{cc}\left(\omega +1\right)I& C\\ -C& \left(\omega +1\right)I\end{array}\right]$  (3.21)

представляют собой ω + 1 ± λ_i^(c)ı с λ_i^(c) ∈ λ(C). Второй знак “≤” обусловлен связью между ∥⋅∥₂ и ∥⋅∥_∞, а последнее строгое неравенство возникает вследствие леммы 2, оценки $||\widehat{E}||$_∞ и учета того, что ω > 0 и 0 ≤ μ_i ≤ ν. Аналогично имеем

${||\widehat{\mathcal{F}}||}_2<\frac{{\left({\omega }^2+v^2\right)}^{\frac{1}{2}}{M}^{\frac{1}{2}}ϵ}{\omega \sqrt{{\left(\omega +1\right)}^2+}{\left[\frac{2^{\alpha +1}\gamma \tau \theta }{{\left(b-a\right)}^{\alpha }}\right]}^2}.$

Замечание 2. Рассмотрим следующее соотношение, сформулированное в теореме 4:

${\widehat{\mathcal{F}}}_{\omega }^{-1}{\mathcal{F}}_{\omega }=\underset{}{\underbrace{I+\tilde{\mathcal{F}}}}+\tilde{\mathcal{E}}.$

С одной стороны, матрица I +$\widehat{\mathcal{F}}$ является малым возмущением единичной матрицы I. В частности, поскольку ∥$\widehat{\mathcal{F}}$∥₂ ≤ 𝜖, теорема Бауэра-Фике [50] приводит к тому, что 

$\left|\xi -1\right|\le ϵ,\forall \xi \in \lambda \left(I+\tilde{\mathcal{F}}\right),$

т.е. собственные значения I +$\tilde{\mathcal{F}}$ сгруппированы в небольшом диске с центром в точке 1 и радиусом 𝜖. С другой стороны, поскольку матрица $\tilde{\mathcal{E}}$ имеет ограниченную ℓ₂-норму и низкий ранг, матрицу $\tilde{\mathcal{F}}$_ω⁻¹ℱ_ω можно рассматривать как низкоранговую модификацию I +$\tilde{\mathcal{F}}$ . Тогда ожидаем, что собственные значения $\tilde{\mathcal{F}}$_ω⁻¹ℱ_ω также расположены в том же диске с центром в точке 1 и радиусом 𝜖, за исключением небольшого количества выбросов.

Наконец, следующая теорема устанавливает связь между $\tilde{\mathcal{F}}$_ω⁻¹ℛ и ℱ_ω⁻¹ℛ.

Теорема 5. Пусть 1 < α < 2, M ≥ 8 четно, 𝜖 – небольшая положительная константа, такая, что

$\frac{2^{\alpha }\mu {\theta }_0}{{\left(M-2\right)}^{\alpha }}<ϵ\le \mu {\theta }_0,k_0=\left[{\left(\frac{\mu {\theta }_0}{ϵ}\right)}^{\frac{1}{\alpha }}\right]+\mathrm{1,}v=\underset{{\mu }_i\in \lambda \left(D\right)}{\max }\left|{\mu }_i\right|.$ (3.22)

Тогда существуют две матрицы 𝒫_ω ∈ ℝ^2M×2M и 𝒬_ω ∈ ℝ^2M×2M, удовлетворяющие

$\text{rank}\left({\mathcal{P}}_{\omega }\right)=4k_0,$  (3.23)

${||{\mathcal{P}}_{\omega }||}_2<\frac{2\left({\omega }^2+v^2\right)M^{\frac{1}{2}}\mu }{\omega \sqrt{{\left(\omega +1\right)}^2+}{\left[\frac{2^{\alpha +1}\gamma \tau \theta }{{\left(b-a\right)}^{\alpha }}\right]}^2}\left[\frac{c_0}{2}-\frac{\theta }{{\left(M-\frac{1}{2}\right)}^{\alpha }}\right],$ (3.24)

${||Q_{\omega }||}_2<\frac{2\left({\omega }^2+v^2\right)M^{{\displaystyle \frac{1}{2}}}\epsilon }{\omega \sqrt{{\left(\omega +1\right)}^2+}{\left[{\displaystyle \frac{2^{\alpha +1}\gamma \tau \theta }{{\left(b-a\right)}^{\alpha }}}\right]}^2}$ (3.25)

при этом

${\tilde{\mathcal{F}}}_{\omega }^{-1}\mathcal{R}={\mathcal{F}}_{\omega }^{-1}\mathcal{R}+{\mathcal{P}}_{\omega }+{\mathcal{Q}}_{\omega }$  (3.26)

Доказательство. Отметим, что

$\begin{array}{c}{\mathcal{F}}_{\omega }^{-1}\mathcal{R}=I-{\tilde{\mathcal{F}}}_{\omega }^{-1}{\mathcal{Q}}_{\omega }=\\ ={\left[\begin{array}{cc}\omega I& -D\\ D& \omega I\end{array}\right]}^{-1}\left(I-{\mathcal{U}}_{\omega }{\mathcal{V}}_{\omega }\right)\left[\begin{array}{cc}\omega I& -D\\ D& \omega I\end{array}\right],\end{array}$

где 𝒰_ω и 𝒱_ω определены в доказательстве теоремы 1. Согласно (2.15) и (2.16), известно, что $||{\mathcal{U}}_{\omega }{\mathcal{V}}_{\omega }||$₂ < 1, и в этом случае $||I-{\mathcal{U}}_{\omega }{\mathcal{V}}_{\omega }||$_{2 ≤ ∥I∥2 +} $||{\mathcal{U}}_{\omega }{\mathcal{V}}_{\omega }||$₂ < 2.

Далее, согласно доказательству теоремы 4, имеем

${||{\left[\begin{array}{cc}\omega I& -D\\ D& \omega I\end{array}\right]}^{-1}||}_2{||\left[\begin{array}{cc}\omega I& -D\\ D& \omega I\end{array}\right]||}_2<\frac{\sqrt{{\omega }^2+v^2}}{\omega }.$

Таким образом $||{\mathcal{F}}_{\omega }^{-1}\mathcal{R}||$₂ < 2$\sqrt{{\omega }^2+v^2}$∕ω.

Ввиду (3.4) и (3.13) выполняется

${\tilde{\mathcal{F}}}_{\omega }^{-1}\mathcal{R}=\left(I-\widetilde{\mathcal{E}+}\tilde{\mathcal{F}} \right){\mathcal{F}}_{\omega }^{-1}\mathcal{R}={\mathcal{F}}_{\omega }^{-1}\mathcal{R}+{\mathcal{P}}_{\omega }+Q_{\omega }$(3.27)

где 𝒫_ω = $\tilde{\mathcal{E}}$ℱ_ω⁻¹ℛ и 𝒬_ω = $\tilde{\mathcal{F}}$ℱ_ω⁻¹ℛ. Из теоремы 4 следует rank $\left({\mathcal{P}}_{\omega }\right)$= 4k₀, 

Замечание 3. Теорема 5 показывает, что 𝒬_ω – матрица малой нормы, поэтому собственные значения ℱ_ω⁻¹ℛ + 𝒬_ω можно рассматривать как малые возмущения собственных значений ℱ_ω⁻¹ℛ. Кроме того, 𝒫_ω имеет ограниченную ℓ₂-норму и низкий ранг, поэтому можно ожидать, что большинство собственных значений $\tilde{\mathcal{F}}$_ω⁻¹ℛ распределены вблизи собственных значений ℱ_ω⁻¹ℛ + 𝒬 _ω. Таким образом, собственные значения $\tilde{\mathcal{F}}$_ω⁻¹ℛ группируются вокруг собственных значений ℱ_ω⁻¹ℛ, за исключением нескольких выбросов.

4. Реализация и сложность. В предыдущих разделах были предложены и проанализированы новые предобусловливатели для блочной линейной системы (2.4), а именно NASS ℱ_ω и CNAS $\tilde{\mathcal{F}}$_ω, реализация которых подробно обсуждается в этом разделе. Фактически можно игнорировать скалярный множитель 1∕(2ω) для ℱ_ω и $\tilde{\mathcal{F}}$_ω по той причине, что он не влияет на свойства предобусловленных матриц системы. Тогда упрощенные предобусловливатели записываются как

$\begin{array}{c}{\mathcal{F}}_{\text{NASS}}=\left[\begin{array}{cc}\widehat{\omega }I& T\\ -T& \widehat{\omega }I\end{array}\right]\left[\begin{array}{cc}\omega I& -D\\ D& \omega I\end{array}\right]=\\ =\left[\begin{array}{cc}I& 0\\ {\widehat{L}}_{21}& I\end{array}\right]\left[\begin{array}{cc}\widehat{\omega }I& T\\ 0& {\widehat{U}}_{22}\end{array}\right]\left[\begin{array}{cc}I& 0\\ L_{21}& I\end{array}\right]\left[\begin{array}{cc}\omega I& -D\\ 0& U_{22}\end{array}\right],\\ {\mathcal{F}}_{\text{CNAS}}=\left[\begin{array}{cc}\widehat{\omega }I& C\\ -C& \widehat{\omega }I\end{array}\right]\left[\begin{array}{cc}\omega I& -D\\ D& \omega I\end{array}\right]=\\ ={\left[\begin{array}{cc}F& 0\\ 0& F\end{array}\right]}^{-1}\left[\begin{array}{cc}I& 0\\ {\tilde{L}}_{21}& I\end{array}\right]\left[\begin{array}{cc}\widehat{\omega }I& \Lambda \\ 0& {\tilde{U}}_{22}\end{array}\right]\left[\begin{array}{cc}F& 0\\ 0& F\end{array}\right]\left[\begin{array}{cc}I& 0\\ L_{21}& I\end{array}\right]\left[\begin{array}{cc}\omega I& -D\\ 0& U_{22}\end{array}\right],\end{array}$

$\begin{array}{c}{||\tilde{\mathcal{E}}||}_2\le {||{\left[\begin{array}{cc}\omega I& -D\\ D& \omega I\end{array}\right]}^{-1}||}_2{||\left[\begin{array}{cc}\omega I& -D\\ D& \omega I\end{array}\right]||}_2{||{\left[\begin{array}{cc}\left(\omega +1\right)I& C\\ -C& \left(\omega +1\right)I\end{array}\right]}^{-1}||}_2{||\widehat{E}||}_2=\\ =\frac{{\max }_{{\mu }_i\in \lambda \left(D\right)}\sqrt{{\omega }^2+{\mu }_i^2}}{{\min }_{{\mu }_i\in \lambda \left(D\right)}\sqrt{{\omega }^2+{\mu }_i^2}}\frac{1}{{\min }_{{\lambda }_i^{\left(c\right)}\in \lambda \left(C\right)}\sqrt{{\left(\omega +1\right)}^2+{\left({\lambda }_i^{\left(c\right)}\right)}^2}}{||\widehat{E}||}_2\le \\ \le \frac{{\max }_{{\mu }_i\in \lambda \left(D\right)}\sqrt{{\omega }^2+{\mu }_i^2}}{{\min }_{{\mu }_i\in \lambda \left(D\right)}\sqrt{{\omega }^2+{\mu }_i^2}}\frac{1}{{\min }_{{\lambda }_i^{\left(c\right)}\in \lambda \left(C\right)}\sqrt{{\left(\omega +1\right)}^2+{\left({\lambda }_i^{\left(c\right)}\right)}^2}}{M}^{\frac{1}{2}}{||\widehat{E}||}_{\infty }<\\ <\frac{{\left({\omega }^2+v^2\right)}^{\frac{1}{2}}{M}^{\frac{1}{2}}\mu }{\omega \sqrt{{\left(\omega +1\right)}^2+}{\left[\frac{2^{\alpha +1}\gamma \tau \theta }{{\left(b-a\right)}^{\alpha }}\right]}^2}\left[\frac{c_0}{2}-\frac{\theta }{{\left(M-\frac{1}{2}\right)}^{\alpha }}\right].\end{array}$

где ω > 0 служит параметром ℱ_NASS и ℱ_CNAS, ℱ ∈ ℂ^M×M – дискретное преобразование Фурье (ДПФ), Λ = diag(Fc) ∈ ℂ^M×M, где c ∈ ℂ^M является первым столбцом C, $\widehat{\omega }$ = ω + 1, $\widehat{L}$₂₁ = −T∕$\widehat{\omega }$, ${\widehat{U}}_{22}$ = $\widehat{\omega }$I + T²∕$\widehat{\omega }$, L₂₁ = D∕omega, U₂₂ = ωI + D²∕ω, $\tilde{L}$₂₁ = −Λ ∕$\widehat{\omega }$, $\tilde{U}$₂₂ = $\widehat{\omega }$I + Λ²∕$\widehat{\omega }$. Диагональные матрицы L₂₁,U₂₂,$\tilde{L}$₂₁,$\tilde{U}$₂₂ можно вычислить заранее.

В методах предварительно обусловленных подпространств Крылова вектор обобщенной невязки вычисляется путем решения соответствующего уравнения, связанного с новыми предобусловливателями, на каждой итерации. Пусть x = (x₁^⊤,x₂^⊤)^⊤ ∈ ℝ^2M с x₁, x₂ ∈ ℝ^M и r = (r₁^⊤,r₂^⊤)^⊤ ∈ ℝ^2M с r₁, r₂ ∈ ℝ^M. Подробности для реализация предварительной подготовки NASS и CNAS приведены в описании алгоритмов 1 и 2.

Стоимость алгоритма 2 складывается из операций быстрого преобразования Фурье (БПФ) и его обратного (ОБПФ) размерности M и векторных операций размерности M (включая сложение векторов, диагональное умножение матрицы на вектор и решение диагональной линейной системы). Согласно [51], БПФ/ОБПФ могут быть выполнены за 𝒪(MlogM) флоп. Все операции с векторами размерности M могут быть произведены за 𝒪(M) флоп. Таким образом, в вычислениях метода ℱ_CNAS преобладают операции БПФ/ОБПФ. Это означает, что вектор обобщенной невязки может быть вычислен за 𝒪(MlogM) флоп на каждой итерации предобусловленных методов подпространств Крылова.

Стоимость алгоритма 1 складывается из векторных операций размерности M (около 𝒪(M) флоп), умножений теплицевой матрицы M × M на вектор (около 𝒪(MlogM) флопс, поскольку его можно реализовать на основе БПФ/ОБПФ) и однократного решения плотной линейной системы размерности M × M с матрицей коэффициентов (ω + 1)I + T²∕(ω + 1). Решение линейной системы может быть реализовано на основе прямого метода (около 𝒪(M³)) или метода предварительно обусловленного сопряженного градиента (PCG) с предобусловливателем (ω + 1)I + C²∕(ω + 1) (около 𝒪(kMlogM) флоп, где k – количество итераций PCG). Очевидно, что вычислительная нагрузка по реализации ℱ_NASS во многом зависит от стоимости решения плотной линейной системы, что намного дороже, чем реализация ℱ_CNAS. Поэтому настоятельно рекомендуется использовать новый предобусловливатель ℱ_CNAS.

Алгоритм 1. Решение уравнения с обобщенной невязкой ℱ_NASS x = r.

1: x₁ = r₁, x₂ = r₂ −$\widehat{L}$₂₁r₁;

2: $\widehat{U}$₂₂x₂ = x₂ и x₁ = (r₁ − Tx₂) ∕$\widehat{\omega }$;

3: x₂ = x₂ − L₂₁x₁;

4: U₂₂x₂ = x₂ и x₁ = (x₁ + Dx₂)∕ω.

На первом шаге этого алгоритма используется одно умножение теплицевой матрицы M ×M на вектор и две векторные операциии размерности M; на втором – одно решение плотной линейной системы M ×M , одно умножение теплицевой матрицы на вектор и две векторные операции; на третьем – две векторные операции; на четвертом – четыре векторные операции.

Алгоритм 2. Решение уравнение обобщенной невязки ℱ_CNAS x = r.

1: x_j = БПФ(r_j), j = 1, 2;

2: x₂ = x₂ −$\tilde{L}$₂₁x₁;

3: $\tilde{U}$₂₂x₂ = x₂ и x₁ = (x₁ − Λx₂)∕$\widehat{\omega }$;

4: x_j = ОБПФ(x_j), j = 1, 2;

5: x₂ = x₂ − L₂₁x₁;

6: U₂₂x₂ = x₂ и x₁ = (x₁ + Dx₂) ∕ω.

На первом шаге этого алгоритма осуществляются два БПФ размерности M; на втором и пятом шагах – по две векторные операции; на третьем и шестом шагах – по четыре векторные операции; на четвертом – два ОБПФ. размерности M.

5. Численные эксперименты. В этом разделе представлены доказательства свойств новых предобусловливателей ℱ_NASS и ℱ_CNAS. Метод GMRES применяется с предобусловливателем CNAS ℱ_CNAS для решения дробных уравнений CNLS в дискретизированном пространстве в случае притягивающего взаимодействия частиц. Дробные дискретизированные уравнения CNLS получены на основе схемы LICD. Для запуска схемы LICD требуется начальное значение на начальном временном уровне и приблизительное значение второго или более высокого порядка на первом временном уровне. Например, приближенное значение второго порядка можно получить с помощью неявной консервативной схемы второго порядка [12].

Во всех численных экспериментах блочная линейная система (2.4) на втором временном уровне дискретизированных дробных уравнений CNLS выбирается в качестве тестируемой линейной системы, а начальное предположение предварительно обусловленного метода GMRES устанавливается равным нулевому вектору. Кроме того, метод GMRES (предварительно обусловленный) выполняется без перезапуска и завершает либо относительную невязку ℓ₂-нормы тестируемой линейной системы, упавшую ниже 10⁻⁶, либо количество итераций, превышающее 3000.

5.1. Случай DNLS. Пусть β = 0, тогда система (0.1) расцеплена. Рассмотрим следующую усеченную систему:

$l{u}_t-\gamma {\left(-\Delta \right)}^{\frac{\alpha }{2}}u+\rho {\left|u\right|}^{2}u=\mathrm{0,}-20\le x\le \mathrm{20,}0<t\le T,$ (5.1)

при соблюдении начальных и граничных условий

$u\left(x\mathrm{,0}\right)=\mathrm{sech}\left(x\right)e^{2lx},u\left(-\mathrm{20,}t\right)=u\left(\mathrm{20,}t\right)=0.$ (5.2)

Здесь взяты параметры γ = 1, ρ = 2, 1 < α ≤ 2. Дискретизированные дробные уравнения DNLS получаются путем применения схемы LICD к (5.1) и (5.2). Требуется решить сложную симметричную линейную систему вида (2.1) на каждом временном уровне t_n для 1 < n ≤ N, что эквивалентно решению блочной линейной системы (2.4). В частности, блочная линейная система (2.4) решается методом предварительно обусловленной GMRES CNAS (CNAS-GMRES). В этом подразделе «IT» представляет количество итераций тестируемого метода.

На рис. 1 показаны кривые IT в зависимости от параметра ω ∈ (0, 8] CNAS-GMRES при α = 1.1 : 0.2 : 1.9, M = 6400, N = 200. Видно, что IT быстро увеличивается, когда ω стремится к нулю. Однако когда ω растет, IT быстро достигает минимума, а затем растет очень медленно, а оптимальное значение ω составляет примерно [2.2, 2.8] для всех случаев α. Таким образом, сходимость CNAS-GMRES менее чувствительна к параметру ω, пока ω не слишком близок к нулю, что делает предобусловливатель CNAS ℱ_CNAS прост в использовании. Кроме того, больший α приводит к большему IT, а это означает, что с увеличением α систему становится сложнее решить.

Рис. 1. Кривые IT от параметра ω ∈ (0, 8] CNAS-GMRES при α = 1.1 : 0.2 : 1.9, M = 6400, N = 200

На рис. 2 показаны кривые IT CNAS-GMRES в зависимости от размера пространственной сетки M схемы LICD, примененной к дробным пространственным уравнениям DNLS, когда α = 1.1 : 0.2 : 1.9, N = 200. Выбрано эмпирическое оптимальное значение ω. Показано, что IT CNAS-GMRES остается практически неизменным с увеличением количества пространственных дискретных точек, что указывает на независимость свойства сходимости CNAS-GMRES от размера пространственной сетки. Кроме того, больший дробный порядок α приводит к более высокому положению кривой на графике.

Рис. 2. Кривые IT CNAS-GMRES от размера пространственной сетки M при α = 1.1 : 0.2 : 1.9, N = 200

В табл. 1 перечислены IT CNAS-GMRES в сочетании с различными циркулянтными матрицами, а также эмпирическое оптимальное значение параметра ω CNAS-GMRES при α = 1.9, M = 6400, N = 200. Интервал поиска параметра ω равен (0, 4]. Перечислены только результаты нескольких репрезентативных циркулянтных матриц [51], включая циркулянтную матрицу Т.Чана, циркулянтную матрицу Стрэнга, циркулянтную матрицу Р.Чана, циркулянтные матрицы, построенные на основе некоторых известных ядер (например, модифицированного ядра Дирихле, ядра фон Ханна, ядра Хэмминга), и супероптимальная циркулянтная матрица. За исключением супероптимальной циркулянтной матрицы, IT CNAS-GMRES со всеми остальными циркулянтными матрицами меньше 10, а оптимальное эмпирическое значение ω практически одинаково, т.е. [0.21, 1.01] ([0.41, 0.81] для циркулянтной матрицы Т.Чана). Эти эксперименты показывают, что большинство циркулянтных матриц в литературе эффективны, когда они применяются к CNAS-GMRES, что позволяет легко выбрать циркулянтное приближение в CNAS-GMRES. Супероптимальная циркулянтная матрица плохо работает в наших экспериментах, и возможная причина заключается в том, что оптимальное значение ω CNAS-GMRES в этом случае остается за пределами интервала поиска (0, 4].

Таблица 1. IT CNAS-GMRES с различными циркулянтными матрицами и эмпирическое оптимальное значение параметра ω CNAS-GMRES при α = 1.9, M = 6400,N = 200

На рис. 3−5 показано распределение собственных значений матрицы ℛ, предобусловленной матрицы NASS ℱ_NASS⁻¹ℛ и предобусловленной матрицы CNAS ℱ_CNAS⁻¹ℛ, когда α = 1.1 : 0.4 : 1.9, M = 1600, 3200. Левые графики относятся к случаю M = 1600, а правые – к случаю M = 3200. На рис. 3 вещественные части собственных значений ℛ равны 1, а мнимые части распределены от −1.3 до 1.3 при M = 1600 и от −2.7 до 2.7 при M = 3200. На рис. 4 вещественные части собственных значений ℛ равны 1, а мнимые части распределены от −8 до 8 при M = 1600 и от −20 до 20 при M = 3200. На рис. 5 вещественные части собственных начений ℛ равны 1, а мнимые части распределены от −42 до 42 при M = 1600 и от −150 до 150 при M = 3200. При этом вещественные части собственных значений ℱ_NASS⁻¹ℛ и ℱ_CNAS⁻¹ℛ распределены от 1.3 до 2, а мнимые – от −0.5 до 0.5 на всех графиках. Очевидно, собственные значения ℱ_NASS⁻¹ℛ и ℱ_CNAS⁻¹ℛ более кластеризованы, чем ℛ (особенно для больших α). Кроме того, большинство собственных значений ℱ_CNAS⁻¹ℛ группируются вокруг значений ℱ_NASS⁻¹ℛи только несколько собственных значений ℱ_CNAS⁻¹ℛ отклоняются сильнее, что соответствует предсказанию замечаний 2, 3. Распределение собственных значений предварительно обусловленных случаев остается почти таким же, когда M увеличивается с 1600 до 3200, что указывает на свойство сходимости NASS-GMRES и CNAS-GMRES, не зависящее от размера пространственной сетки.

Рис. 3. Распределение собственных значений ℛ, ℱ_NASS⁻¹ℛ и ℱ_CNAS⁻¹ℛ, при α = 1.1, M = 1600, N = 200 (слева) и M = 3200, N = 200 (справа)

Рис. 4. Распределение собственных значений ℛ, ℱ_NASS⁻¹ℛ и ℱ_CNAS⁻¹ℛ, при α = 1.5, M = 1600,N = 200 (слева) и M = 3200, N = 200 (справа)

Рис. 5. Распределение собственных значений ℛ, ℱ_NASS⁻¹ℛ и ℱ_CNAS⁻¹ℛ, при α = 1.9, M = 1600, N = 200 (слева) и M = 3200, N = 200 (справа)

На рис. 6−9 слева показано численное решение u_CNAS, полученное с помощью CNAS-GMRES, а справа – его ошибка err_u = |u_CNAS − u_GE|, т.е. отклонение от решения, полученного в схеме LICD с использованием решения плотной линейной системы методом Гаусса, для дробно-пространственных уравнений DNLS (5.1) при M = 800, N = 200, α = 1.1 : 0.4 : 1.9 и α = 2. Замечено, что дробный порядок α будет влиять на форму волнового фронта как по высоте, так и по ширине. Чем меньше α, тем быстрее будет меняться форма волнового фронта. При стремлении α к 2 волновой фронт пространственного дробного уравнения DNLS сходится к волновому фронту недробного уравнения, и высота волнового фронта стабильна. Кроме того, ошибка между численным решением и точным решением схемы LICD остается всего лишь около 10⁻⁴ во всей задействованной пространственно-временной области, следовательно, численное решение, полученное с помощью CNAS-GMRES, является надежным.

Рис. 6. Численное решение (слева) и его ошибка (справа) дробного пространственного уравнения DNLS (5.1) при α = 1.1, M = 800, N = 200

Рис. 7. Численное решение (слева) и его ошибка (справа) дробного пространственного уравнения DNLS (5.1) при α = 1.5, M = 800, N = 200

Рис. 8. Численное решение (слева) и его ошибка (справа) дробного пространственного уравнения DNLS (5.1) при α = 1.9, M = 800, N = 200

Рис. 9. Численное решение (слева) и его ошибка (справа) дробного пространственного уравнения DNLS (5.1) при α = 2, M = 800, N = 200

В табл. 2 и на рис. 10 показаны относительные погрешности дискретной массы и энергии для различных значений α соответственно. Размер пространственной сетки и размер временного шага составляют h = 0.2 и τ = 0.05 для табл. 2 и h = 0.2 и τ = 0.001 для рис. 10. На каждом временном уровне блочная линейная система решается с помощью CNAS-GMRES и завершается, когда относительная невязка по норме ℓ₂ блочной линейной системы снижается ниже 10⁻¹⁵. Указанные ошибки остаются очень небольшими, а это означает, что линейный решатель CNAS-GMRES сохраняет устойчивость схемы LICD.

Рис. 10. Относительные погрешности дискретной энергии, т.е. |(Eⁿ −E⁰)∕E⁰|, при h = 0.2, τ = 0.001

Таблица 2. Относительные ошибки дискретной массы, т. е. |(Qⁿ − Q⁰)∕Q⁰| при h = 0.2, τ = 0.05

5.2. Случай CNLS. Эксперименты здесь проводятся в связанном случае, т. е. усеченной системе уравнений CNLS:

$\left\{\begin{array}{c}l{u}_t-\gamma {\left(-\Delta \right)}^{\frac{\alpha }{2}}u+\rho \left({\left|u\right|}^2+\beta {\left|v\right|}^2\right)u=\mathrm{0,}\\ l{u}_t-\gamma {\left(-\Delta \right)}^{\frac{\alpha }{2}}v+\rho \left({\left|v\right|}^2+\beta {\left|u\right|}^2\right)v=\mathrm{0,}\end{array}\right.-20\le x\le \mathrm{20,}0<t\le T,$ (5.3)