Журнал вычислительной математики и математической физики, 2021, T. 61, № 5, стр. 706-722

2.1. Методы подпространства Крылова и НВ перезапуск

Предположим, что по данным $A \in {{\mathbb{R}}^{{n \times n}}}$, ${v} \in {{\mathbb{R}}^{n}}$ и $t > 0$ нужно вычислить действие матричной экспоненты матрицы $ - tA$ на вектор ${v}$, т.е.

Эта задача эквивалентна решению задачи Коши где, слегка пренебрегая точностью обозначений, мы используем $t$ и как независимую переменную в (3), и для обозначения длины интервала в (2). Метод подпространства Крылова для вычисления действия матричной экспоненты можно рассматривать как галеркинскую проекцию задачи Коши (3) на подпространство Крылова Сначала ортонормальный базис подпространства ${{\mathcal{K}}_{k}}(A,{v})$ вычисляется обычным процессом Арнольди (или, если $A = {{A}^{{\text{т}}}}$, процессом Ланцоша) (см. [24], [28]–[30]) и сохраняется в виде столбцов ${{{v}}_{1}}$, …, ${{{v}}_{k}}$ матрицы ${{V}_{k}} = [{{{v}}_{1}}\; \ldots \;{{{v}}_{k}}] \in {{\mathbb{R}}^{{n \times k}}}$, так что выполняется так называемое разложение Арнольди: где ${{e}_{k}} = {{(0,\; \ldots ,\;0,\;1)}^{{\text{т}}}} \in {{\mathbb{R}}^{k}}$, ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{H} }_{k}} \in {{\mathbb{R}}^{{(k + 1) \times k}}}$ – верхняя хессенбергова матрица, а матрица ${{H}_{k}} \in {{R}^{{k \times k}}}$ состоит из первых $k$ строк матрицы ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{H} }_{k}}$. После этого крыловская аппроксимация ${{y}_{k}}(t) \approx exp( - tA){v}$ определяется как (см. [31]–[33]) где $u(t):\mathbb{R} \to {{\mathbb{R}}^{k}}$ решает задачу Коши с проецированной матрицей ${{H}_{k}} = V_{k}^{{\text{т}}}A{{V}_{k}}$: Здесь $\beta = \left\| {v} \right\|$ и ${{e}_{1}} = {{(1,\;0,\; \ldots ,\;0)}^{{\text{т}}}} \in {{\mathbb{R}}^{k}}$. Заметим, что задача Коши (7) – это галеркинская проекция задачи Коши (3) на подпространство Крылова и что $u(t)$ может быть вычислена как Если $k$ не слишком велико, то (8) – предпочтительный способ вычисления $u$, который может быть реализован многими стандартными методами линейной алгебры (см. [24]–[26]). Вычисление $exp( - t{{H}_{k}})$ обычно является более простой операцией, чем решение системы дифференциальных уравнений в (7), что требует выбора подходящего решателя (в частности, жесткого или нежесткого) и его параметров. Метод подпространства Крылова (5)–(8) иногда называют полиномиальным методом подпространства Крылова, чтобы подчеркнуть тот факт, что вектора подпространства (4) – многочлены от $A$.

Естественным способом контроля (неизвестной) ошибки крыловского приближения (6) является отслеживание невязки ${{r}_{k}}(t)$ этого приближения ${{y}_{k}}(t)$ по отношению к системе дифференциальных уравнений $y{\text{'}} = - Ay$, а именно (см. [12], [27], [34]),

Невязка ${{r}_{k}}(t)$ доступна в ходе процесса Арнольди и может быть вычислена по формуле (см. [12], [27], [34]) Как видим, ${{r}_{k}}(t)$ – скалярная функция, помноженная на ${{{v}}_{{k + 1}}}$. Следовательно, $V_{k}^{{\text{т}}}{{r}_{k}}(t) = 0$ для всех $t > 0$, а (6) – действительно галеркинская проекция на ${{\mathcal{K}}_{k}}(A,{v})$. Некоторые результаты по сходимости невязки и ее связи с ошибкой могут быть найдены, например, в [22], [27].

Метод подпространства Крылова СО (“сдвиг–обращение”) (см. [7], [8]) для вычисления (2) отличается от стандартного полиномиального метода подпространства Крылова, описанного выше, тем, что подпространство Крылова строится для сдвинутой и обращенной матрицы ${{(I + \gamma A)}^{{ - 1}}}$, а не для $A$, где параметр $\gamma > 0$ фиксирован. Это делается для ускорения сходимости: процесс Арнольди имеет тенденцию лучше приближать наибольшие по модулю собственные числа, а для матрицы ${{(I + \gamma A)}^{{ - 1}}}$ они соответствуют малым собственным значениям матрицы $A$. Именно эти собственные числа важны для матричной экспоненты (компоненты, соответствующие большим собственным числам, не столь важны, они угасают экспоненциально быстро). Цена этого ускорения в сходимости – необходимость решать линейные системы с матрицей $I + \gamma A$ на каждой крыловской итерации. Разложение Арнольди (5) для сдвинутой и обращенной (СО) матрицы ${{(I + \gamma A)}^{{ - 1}}}$ принимает вид

Это соотношение удобнее использовать после преобразования где обозначение $\tilde { \cdot }$ указывает, что проекция построена для СО матрицы ${{(I + \gamma A)}^{{ - 1}}}$, а матрица ${{H}_{k}}$ определяется как преобразование, обратное к сдвигу и обращению: Подчеркнем, что матрицы ${{V}_{{k + 1}}}$ и ${{H}_{k}}$ здесь отличаются от матриц ${{V}_{{k + 1}}}$ и ${{H}_{k}}$ из соотношения (5). Подробный анализ метода подпространства Крылова СО и родственных методов можно найти в (см. [7], [8], [10]).

В методе подпространства Крылова СО невязка легко может быть вычислена следующим образом (см. [27]):

Здесь $u(t)$ – решение спроецированной задачи Коши (7), где ${{H}_{k}}$ – матрица, получающаяся обратным преобразованием (12).

Величина сдвига $\gamma $ обычно выбирается в соответствии с длиной $t$ временного интервала $[0,t]$ (см. [8]), при этом часто используется величина $\gamma = t{\text{/}}10$. Таким образом, изменение параметра $\gamma $ означает изменение $t$. Стандартный полиномиальный метод подпространства Крылова (5)–(8) имеет привлекательное свойство инвариантности относительно $t$: коль скоро ${{V}_{{k + 1}}}$ и $\mathop {\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{H} }\nolimits_k $ вычислены, они могут быть использованы для любых времен $t$ (хотя качество аппроксимации ${{y}_{k}}(t) \approx y(t)$, вообще говоря, ухудшается с ростом $t$). К сожалению, это свойство полностью не распространяется на метод подпространства Крылова СО: в этом методе матрицы ${{V}_{{k + 1}}}$ и $\mathop {\tilde {\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{H} }}\nolimits_k $ зависят от величины $\gamma $, которая, в свою очередь, зависит от $t$. Тем не менее на практике можно использовать вычисленные матрицы Арнольди ${{V}_{{k + 1}}}$ и $\mathop {\tilde {\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{H} }}\nolimits_k $ для определенного диапазона $t$, не вычисляя их заново.

Недавно предложенный НВ перезапуск основан на том факте, что невязка как функция от $t$ обычно является для обычного метода подпространства Крылова неубывающей функцией. Следовательно, когда выполнено максимально допустимое число ${{k}_{{max}}}$ крыловских итераций (так что хранение большего числа базисных векторов Крылова и работа с ними слишком затратны), мы можем найти подынтервал $[0,\delta ]$, $\delta < t$, на котором невязка уже достаточно мала по норме. Тогда мы можем перезапустить метод, полагая ${v}: = {{y}_{{{{k}_{{max}}}}}}(\delta )$, уменьшая временной интервал $t: = t - \delta $ и выполняя следующие ${{k}_{{max}}}$ крыловских итераций для задачи (2) с обновленными ${v}$ и $t$. Схема НВ перезапуска представлена на фиг. 1.

2.2. ТНВ перезапуск: идеи и алгоритм

В методе подпространства Крылова СО (5)–(8) невязка как функция от $t$ проявляет гораздо более нерегулярное поведение, чем в обычном полиномиальном методе подпространства Крылова (фиг. 2). Если для метода подпространства Крылова СО применяется НВ перезапуск, то может случиться, что такое $\delta $, что $\left\| {{{r}_{k}}(s)} \right\|$ не превышает заданной точности для $s \in [0,\delta ]$, не может быть найдено или слишком мало для практически эффективного перезапуска. Разумеется, можно устроить перезапуск, положив $\delta $ равной любой точке $s$, где $\left\| {{{r}_{k}}(s)} \right\|$ достаточно мала (фиг. 2). Однако гарантии, что $mi{{n}_{{s \in [0,t]}}}\left\| {{{r}_{k}}(s)} \right\|$ не превосходит заданной точности, нет, и в таком случае перезапуск с любым $\delta \in [0,t]$ неизбежно приводит к потере точности.

В данной работе предлагается подход, позволяющий устранить этот недостаток НВ перезапуска для метода подпространства Крылова СО. Подход наш основан на следующих двух наблюдениях.

1. Поскольку $\gamma $ обычно выбирается пропорционально $t$, меньшая величина сдвига $\gamma $ означает более короткий временной интервал $[0,t]$. Для несимметричных матриц $A$ невязка ${{r}_{k}}(s)$ в методе подпространства Крылова СО становится меньше по норме с уменьшением $\gamma $ на некотором подынтервале $s \in [0,\delta ]$, $0 < \delta < t$ (фиг. 2, 3).

2. Как уже обсуждалось выше, чтобы изменить $\gamma $ в подпространстве Крылова СО, необходимо заново выполнить весь процесс Арнольди. Однако если линейная система с матрицей $I + \gamma A$ решена для определенного значения сдвига $\gamma $, то часть уже проделанной вычислительной работы может быть использована для решения линейных систем с меньшим значением сдвига $\tilde {\gamma } \leqslant \gamma $. В частности, если для некоторого сдвига $\gamma $ вычислено (разреженное) LU разложение, то оно может быть успешно использовано как предобусловливатель для решения линейных систем с новым значением сдвига, с матрицами $I + \tilde {\gamma }A$, $\tilde {\gamma } \leqslant \gamma $ (см. утверждение 2).

На основе этих наблюдений для метода подпространства Крылова СО предлагается организовать НВ перезапуск без потери точности следующим образом. Предположим, что может быть выполнено не более ${{k}_{{max}}}$ шагов процесса Арнольди или Ланцоша, так как хранение и использование более ${{k}_{{max}}}$ крыловских векторов слишком затратны. Тогда выполняются итерации $k = 1,\;2,\; \ldots ,\;{{k}_{{max}}}$ с контролем на каждой итерации нормы невязки $\left\| {{{r}_{k}}(t)} \right\|$ (см. (13)). Если норма невязки меньше заданной точности, процесс успешно оканчивается. Иначе после выполнения шага $k = {{k}_{{max}}}$, значения функции $\left\| {{{r}_{k}}(s)} \right\|$ анализируются на отрезке $s \in [0,t]$. Если не может быть найдено точки $s$ такой, что норма $\left\| {{{r}_{k}}(s)} \right\|$ достаточно мала, то мы уменьшаем $\gamma $ в два раза и повторяем шаги метода $k = 1,\;2,\; \ldots ,\;{{k}_{{max}}}$. Затем процедура перезапуска (подробно описанная на фиг. 4) повторяется до тех пор, пока норма невязки не будет достаточно малой в конечной точке заданного временного интервала.

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

Отметим, что для симметричных матриц $A$ описанное выше поведение невязки в методе подпространства Крылова СО с уменьшением $\gamma $ не наблюдается. Это подтверждается довольно точной оценкой леммы 3.1 из [8] (где полагаем $\mu = 0$ и выбираем $\gamma $ пропорционально $t = \tau $).

Следующие лемма и утверждение показывают, что норма невязки ${{r}_{k}}(s)$ метода подпространства Крылова СО является функцией, ограниченной по времени. Оценка, дающая ограничение, зависит от $\gamma $.

Лемма 1. Пусть для $A \in {{\mathbb{R}}^{{n \times n}}}$ выполняется соотношение (1) и пусть ${{H}_{k}}$ – матрица, полученная в методе подпространства Крылова СО (см. (11), (12)). Тогда существует такая константа ${{\omega }_{k}} \geqslant 0$, что

Доказательство. Пусть $\omega = mi{{n}_{{x \in {{\mathbb{C}}^{n}},\,\left\| x \right\| = 1}}}{\text{Re}}(x{\text{*}}Ax)$. Известно (см., например, [35, Теорема 2.4]), что

В силу (1) эти два эквивалентные соотношения выполняются для некоторой константы $\omega \geqslant 0$. Кроме того, эти соотношения эквивалентны неравенству (см. [35, теорема 2.13]) которое выполняется для всех $\gamma > 0$ и всех $\omega \in \mathbb{R}$ при условии, что $1 + \gamma \omega > 0$. Пусть $\gamma > 0$ так же, как и в методе подпространства Крылова СО. Определим так, чтобы $\left\| {{{{\tilde {H}}}_{k}}} \right\| = 1{\text{/}}(1 + \gamma {{\omega }_{k}})$. Получаем откуда следует, что $0 \leqslant \omega \leqslant {{\omega }_{k}}$. Поскольку ${{\tilde {H}}_{k}} = {{(I + \gamma {{H}_{k}})}^{{ - 1}}}$ (см. (12)), видим, снова используя (см. [35, теорема 2.13]), что Это равносильно искомому неравенству (14). Доказательство завершено.

Определим функцию $\varphi (z)$ (см., например, [1]):

Утверждение 1. Пусть для $A \in {{\mathbb{R}}^{{n \times n}}}$ выполняется соотношение (1) и пусть ${{r}_{k}}(t)$ – невязка метода подпространства Крылова СО (13) при решении задачи (2). Тогда для всех $t \geqslant 0$ справедливо

где функция $u(t)$ определена в (8), (12), ${{\omega }_{k}} \geqslant 0$ – константа из (14), а ${{\tilde {h}}_{{k + 1,k}}}$ и ${{h}_{{k,1}}}$ – элементы матриц ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\tilde {H}} }_{k}} \in {{\mathbb{R}}^{{(k + 1) \times k}}}$ и ${{H}_{k}} \in {{\mathbb{R}}^{{k \times k}}}$ соответственно (см. (11), (12)). Минимум в соотношении (18) берется по двум элементам множества, обозначенного ${\text{\{ }} \ldots {\text{\} }}$. Подчеркнем, что ${{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\tilde {H}} }_{k}}$ в оценке выше зависит от $\gamma $ (а следовательно, также и ${{H}_{k}}$, $u(t)$, ${{\omega }_{k}}$).

Доказательство. Соотношение (17) идентично (13) (см. (12)), а доказательство (13) можно найти в [27]. В силу (7), (12) имеем

Далее нетрудно проверить, что (см. (16)) Поэтому, учитывая $u(0) = \beta {{e}_{1}}$ и $\tilde {H}_{k}^{{ - 1}} = I + \gamma {{H}_{k}}$, можно оценить Здесь используется неравенство справедливое в силу (14), (15) (см., например, [1, доказательство леммы 2.4]). Оценка (19) особенно полезна для малых $t$. Получим теперь альтернативную оценку, которая может быть точнее для больших $t$: где используется неравенство (14). Из (19), (20) следует, что Это позволяет оценить что и дает (18). Доказательство завершено.

Подчеркнем, что оценка (18), к сожалению, не настолько точна, чтобы отразить зависимость $\left\| {{{r}_{k}}(t)} \right\|$ от $\gamma $ (см. фиг. 2 и 3). Однако следует сделать следующее

Замечание 1. Численные эксперименты показывают, что величина $\left| {{{h}_{{k,1}}}} \right|$ (вспомним, что $\left| {{{\beta }_{k}}(0)} \right| = \beta {{\tilde {h}}_{{k + 1,k}}}\left| {{{h}_{{k,1}}}} \right|$), появляющаяся в (18), обычно мала, много порядков меньше, чем другой член $\tfrac{1}{\gamma }min{\text{\{ }} \ldots {\text{\} }}$, участвующий в правой части (18). Если $\left| {{{h}_{{k,1}}}} \right| = 0$, то ${{\beta }_{k}}(0) = 0$. Таким образом, оценка (18) показывает, что для любого $k$ и для любой точности $\varepsilon > 0$ можно найти такой временной интервал $[0,\;\delta ]$, что $\left\| {{{r}_{k}}(s)} \right\| \leqslant \varepsilon $ для $s \in [0,\delta ]$. В этом случае временной интервал можно сократить (фиг. 1), поэтому как НВ, так и ТНВ перезапуски гарантируют сходимость метода подпространства Крылова СО для любой длины перезапуска. Разумеется, указанное $\delta $ может быть слишком мало, чтобы использовать его на практике. Поэтому подстройка параметра $\gamma $, как это делается в ТНВ перезапуске, может быть необходима для успешной работы.

2.3. Решение сдвинутых линейных систем

Покажем теперь, что LU разложение, вычисленное для матрицы $I + \gamma A$, может быть успешно использовано как предобусловливатель для сдвинутой матрицы $I + \tilde {\gamma }A$ с уменьшенным значением сдвига $\tilde {\gamma }$, $0 < \tilde {\gamma } \leqslant \gamma $. Говоря точнее, для решений сдвинутой линейной системы

будем применять предобусловливатель Тогда нетрудно показать (см. утверждение 2 ниже), что даже метод простых итераций с таким предобусловливанием, точнее сходится безусловно, т.е. для спектрального радиуса $\rho (\tilde {G})$ матрицы перехода $\tilde {G}$ выполняется $\rho (\tilde {G}) < 1$. Следовательно, собственные числа предобусловленной матрицы ${{\mathcal{M}}^{{ - 1}}}\mathcal{A}$ расположены на комплексной плоскости внутри круга единичного радиуса с центром в точке $1 + 0i$, ${{i}^{2}} = - 1$. Это означает, что современные итерационные методы подпространства Крылова, такие как GMRES, BiCGSTAB, QMR или аналогичные им (см. [29], [30], [36]), будут успешно решать предобусловленную линейную систему (21).

Однако чем меньше значение сдвига $\tilde {\gamma }$, тем лучше обусловлена сдвинутая матрица $I + \tilde {\gamma }A$. Следовательно, для малого $\tilde {\gamma }$ может оказаться, что непредобусловленный итерационный метод сходится достаточно быстро. Поэтому в утверждении 2 нами дается условие, достаточное для того, чтобы предобусловленный метод простых итераций сходился быстрее, чем непредобусловленный метод.

Утверждение 2. Пусть для $A \in {{\mathbb{R}}^{{n \times n}}}$ выполняется соотношение (1), $0 < \tilde {\gamma } \leqslant \gamma $ и пусть линейная система с матрицей $I + \tilde {\gamma }A$ решается итерационно. Тогда предобусловленный метод простых итераций (22) с матрицей предобусловливателя $\mathcal{M} = I + \gamma A$ сходится.

Кроме того, пусть метод простых итераций без предобусловливания также сходится. Тогда предобусловленный метод простых итераций (22) с матрицей предобусловливателя $\mathcal{M} = I + \gamma A$ сходится быстрее, чем метод простых итераций без предобусловливания при условии, что

где $\rho (A)$ – спектральный радиус матрицы $A$.

Доказательство. Пусть $\lambda $ – некоторое собственное число матрицы $A$. Тогда собственные числа предобусловленной матрицы ${{(I + \gamma A)}^{{ - 1}}}(I + \tilde {\gamma }A)$ имеют вид

Предобусловленный метод простых итераций сходится тогда и только тогда, когда собственные числа матрицы перехода $\tilde {G} = I - {{(I + \gamma A)}^{{ - 1}}}(I + \tilde {\gamma }A)$ по модулю меньше единицы, т.е. если Левую часть этого неравенства можно оценить сверху: где второе неравенство выполняется потому, что собственные числа матрицы $A$ имеют неотрицательную вещественную часть (см. (1)). Следовательно, предобусловленный метод простых итераций сходится.

Далее заметим, что матрица перехода $G$ метода простых итераций без предобусловливания – это  матрица  $G = I - (I + \tilde {\gamma }A) = - \tilde {\gamma }A$.  Предобусловленный метод сходится быстрее, чем непредобусловленный, при условии $\rho (\tilde {G}) < \rho (G)$, т.е. если

Левую часть здесь можно оценить величиной $1 - \tfrac{{\tilde {\gamma }}}{\gamma }$, так что неравенство выполняется, если Легко проверить, что это неравенство равносильно (23). Доказательство завершено.

3.1. Подробности экспериментов

Мы реализовали ТНВ перезапуск так, как показано на фиг. 4, со следующими модификациями.

1. Находимое алгоритмом приближенное наименьшее значение нормы невязки зависит от числа точек, в которых вычисляется норма невязки. Поэтому это число точек, обозначенное в описании алгоритма (фиг. 4) ${{n}_{{{\text{steps}}}}}$, задается в соответствии с требуемой точностью $tol$:

2. Если величина $\delta $ (т.е. такая величина, что $\left\| {{{r}_{k}}(\delta )} \right\| \leqslant {\text{tol}}$, см. фиг. 4) остается равной нулю после двух последовательных уменьшений $\gamma $, то $\gamma $ получает свое изначальное значение $4\gamma $, помноженное на 0.8, т.е. $\gamma : = 0.8 \times 4 \times \gamma $, а число тестовых точек ${{n}_{{{\text{steps}}}}}$ удваивается. После этого обычная работа алгоритма продолжается. Это означает, что если изначально $\gamma = 1$, $\gamma $ в алгоритме последовательно принимает значения 1, 0.5, 0.25, 0.8, 0.4, 0.2, 0.64, 0.32, … . Как только найдено положительное $\delta $, т.е. перезапуск успешен, число тестовых точек ${{n}_{{{\text{steps}}}}}$ уменьшается до 500.

3. Как только в ходе выполнения алгоритма $\gamma $ меняется, пропорционально меняется и длина временного интервала, на котором по ${{n}_{{{\text{steps}}}}}$ тестовым точкам ищется приближенное минимальное значение нормы невязки. Например, если значение $\gamma $ уменьшается вдвое, то временной интервал поиска меняется от $[0,\;t]$ к $[0,t{\text{/}}2]$. Это делается потому, что $\left\| {{{r}_{k}}(s)} \right\|$ навряд ли будет мала для $s > t{\text{/}}2$ (фиг. 3). Как только перезапуск оказался успешным, т.е. $\delta > 0$ и временной интервал уменьшен (строка алгоритма $t: = t - \delta $), мы увеличиваем интервал поиска до $[0,\;t]$.

4. Последняя модификация состоит в том, что в нашей реализации итерационный метод GMRES(10) (см. [37]) может быть использован вместо LU разложения не только тогда, когда сдвиг $\gamma $ уменьшен. В качестве предобусловливателя с GMRES(10) может быть использовано неполное LU разложение ILUT($\varepsilon $) (см. [38, гл. 10]. Оно вычисляется один раз и используется для всех значений $\gamma $. Мы используем реализацию GMRES, доступную на сайте www.netlib.org/templates/, из [36]. Предобусловливатель применяется справа, а в качестве остановочного критерия итераций GMRES(10) при решении линейных систем “сдвиг–обращение” $(I + \gamma A)x = b$ используется такое условие на невязку res_i линейной системы:

где ${\text{tol}}$ – требуемая заданная точность вычисления действия матричной экспоненты. При малых значениях $\gamma $ может быть разумным отключить предобусловливатель (см. утверждение 2).

Начальное значение $\gamma $ может быть задано как необязательный параметр нашей процедуры ТНВ перезапуска. По умолчанию, если начальное значение $\gamma $ не задано пользователем, оно устанавливается равным $t{\text{/}}20$. Заметим, что значение $\gamma = t{\text{/}}10$ рекомендуется в [8] для умеренных величин допустимой точности ${\text{tol}} \approx {{10}^{{ - 6}}}$ для симметричных матриц. Выбор начального значения $\gamma : = t{\text{/}}20$ в ТНВ перезапуске представляется разумным выбором, поскольку, согласно нашему опыту, оптимальное значение $\gamma $ для несимметричных матриц обычно меньше, чем $t{\text{/}}10$.

В экспериментах, представленных ниже, в рамках метода подпространства Крылова СО мы сравниваем ТНВ перезапуск с НВ перезапуском. В работе [22] НВ перезапуск был всесторонне протестирован и сравнен с другими тремя методами перезапуска: перезапуском пакета EXPOKIT (см. [39]), перезапуском Нихоффа–Хохбрук (см. [23, гл. 3]) и невязочным перезапуском (см. [12], [27]). Значения ошибок, представленные в этом разделе, получены для численного решения ${{y}_{k}}(t)$, как

где ${{y}_{{{\text{ref}}}}}(t)$ – референтное решение, вычисленное с высокой точностью функцией ${\text{phiv}}$ пакета EXPOKIT (см. [39]). Численные тесты были проведены в Матлабе на линукс-компьютере с 8 процессорами Intel Xeon E5504 2.00GHz.

3.2. Задача конвекции-диффузии

В этой задаче матрица $A$ получается стандартной пятиточечной конечно-разностной аппроксимацией оператора конвекции-диффузии, определенного для функций $u(x,y)$ с $(x,y) \in \Omega = [0,1] \times [0,1]$ и ${{\left. u \right|}_{{\partial \Omega }}} = 0$. Оператор имеет вид

где Pe – число Пекле. Здесь конвективные члены (первые производные) записаны в специальном виде так, чтобы их вклады в матрицу A представляли собой кососимметричную матрицу (см. [40]). В экспериментах использовалась равномерная сетка $802 \times 802$ и значения числа Пекле Pe = 200 и Pe = 1000. Размер задачи для этой сетки – $n = {{800}^{2}} = 640\,000$. Для обоих значений числа Пекле ${{\left\| {\tfrac{1}{2}\left( {A + {{A}^{{\text{т}}}}} \right)} \right\|}_{2}} \approx 6000$, в то время как ${{\left\| {\tfrac{1}{2}\left( {A - {{A}^{{\text{т}}}}} \right)} \right\|}_{2}} \approx 0.5$ для ${\text{Pe}} = 200$ и ${{\left\| {\tfrac{1}{2}\left( {A - {{A}^{{\text{т}}}}} \right)} \right\|}_{2}} \approx 2.5$ для ${\text{Pe}} = 1000$. Следовательно, в обоих случаях матрицы можно считать слабо несимметричными. Значения функции $sin(\pi x)sin(\pi y)$ на конечно-разностной сетке присваивались начальному вектору ${v}$, который затем нормализовывался ${v}: = {v}{\text{/}}\left\| {v} \right\|$. Задавалось конечное время $t = 1$.

В экспериментах начальное значение сдвига $\gamma $ в ТНВ перезапуске не задавалось и было по умолчанию $t{\text{/}}20$, а в НВ перезапуске использовалось обычное значение $\gamma = t{\text{/}}10$. Это не обязательно дает преимущество ТНВ перезапуску, потому что оптимальные значения $\gamma $, находимые в ТНВ перезапуске, все равно были меньше, чем $t{\text{/}}20$.

Результаты для этой тестовой задачи представлены в табл. 1. Как можно увидеть в первых двух строках таблицы, НВ перезапуск не в состоянии дать меньшую ошибку при уменьшенном значении допустимой точности, несмотря на возросший объем вычислений (27 вместо 20 итераций). В то же время ТНВ перезапуск, где используется тот же самый линейный решатель – разреженное LU разложение – справляется с поставленной задачей, хотя процессорное время и увеличивается в 10 раз (см. строку 3 табл. 1). Отметим, что процессорное время, измеренное средствами Матлаба, не всегда правильно отражает вычислительную работу. В данном случае оно не соответствует числу крыловских шагов метода (73 шага с ТНВ перезапуском вместо 27 шагов c НВ перезапуском). Причина здесь в том, что прямые методы решения линейных систем (LU разложение и оператор действия обратной матрицы “обратная дробная черта” \) реализованы в среде Матлаб весьма эффективно, чего нельзя сказать об итерационных решателях. Из-за этого процессорное время вычисления действий предобусловливателя внутри GMRES(10) оказывается значительным. Однако, как только алгоритм определил подходящее значение $\gamma $, это значение успешно может быть использовано для повторных вычислений действия матричной экспоненты: в этом случае, как видим в строке 4 табл. 1, мы получаем требуемую высокую точность при небольшом увеличении расчетного времени.

Мы также тестируем наш подход на этой задаче с итерационным линейным решателем – итерационным методом GMRES(10) с предобусловливателем ILUT($\varepsilon = {{10}^{{ - 3}}}$). В строке 5 табл. 1 показано, что  НВ перезапуск в комбинации с предобусловленным GMRES требует 36 шагов Арнольди (вместо 27 шагов для НВ перезапуска в комбинации с прямым методом решения, см. табл. 1, строка 2), потому что невязки в этих двух реализациях метода (с LU разложением и с методом GMRES(10)) слегка отличаются. ТНВ перезапуск с тем же самым предобусловленным итерационным методом требует в два раза больше процессорного времени, чем НВ перезапуск, но дает меньшую ошибку (см. табл. 1, строка 6). Наконец, в последней строке таблицы мы видим, что коль скоро правильное значение сдвига определено ТНВ алгоритмом, ТНВ перезапуск позволяет получать лучшую точность при сравнимых вычислительных затратах.

В нижней части табл. 1 представлены результаты для большего числа Пекле. При требуемой точности $tol = 1e - 06$ перезапуск НВ дает результат с точностью $3.36e - 07$, что вполне удовлетворительно. Однако при требуемой точности $tol = 1e - 08$ метод оказывается не в состоянии получить меньшую ошибку, хотя вычислительные затраты выросли с 14 до 27 шагов. В следующих двух строках таблицы показаны результаты для ТНВ перезапуска. ТНВ перезапуск позволяет получить более точный результат при возросшем процессорном времени. Из последней строки табл. 1 следует, что как только оптимальное значение $\gamma $ найдено, аналогичная точность может быть получена за примерно то же процессорное время.

3.3. Уравнения Максвелла в среде без потерь

Рассмотрим уравнения Максвелла в трехмерной области, заполненной непроводящей средой без источников электромагнитного поля:

Здесь $\varepsilon $ и $\mu $ – относительные диэлектрическая и магнитная проницаемости соответственно, являющиеся скалярными функциями переменных $(x,y,z)$, а магнитное поле ${\mathbf{H}}$ и электрическое поле ${\mathbf{E}}$ – неизвестные вектор-функции переменных $(x,y,z,t)$. Краевые условия задают на границе области нулевые тангенциальные компоненты электрического поля, что физически означает идеально проводящую границу области или так называемые краевые условия “большого бака” (см. [41], [42]). Данная модельная задача взята из [42]: в пространственной области $[ - 6.05,6.05] \times [ - 6.05,6.05] \times [ - 6.05,6.05]$, наполненной воздухом (относительная диэлектрическая проницаемость ${{\varepsilon }_{r}} = 1$), помещен образец из материала с относительной диэлектрической проницаемостью ${{\varepsilon }_{r}} = 5.0$, занимающий подобласть $[ - 4.55,4.55] \times [ - 4.55,4.55] \times [ - 4.55,4.55]$. В образце имеются 27 сферических отверстий (${{\varepsilon }_{r}} = 1$) радиуса $1.4$, центры отверстий расположены в точках $({{x}_{i}},{{y}_{j}},{{z}_{k}}) = (3.03i,3.03j,3.03k)$, $i,j,k = - 1,\;0,\;1$. Задаются нулевые начальные значения для всех компонент обоих полей ${\mathbf{H}}$ и ${\mathbf{E}}$, кроме компонент $x$ и $y$ поля ${\mathbf{E}}$. Последние не равны нулю в центре области и представляют собой световой импульс.

Дискретизация по пространству центральными, разнесенными по сетке, конечными разностями (схема Йи (Yee)) приводит в системе дифференциальных уравнений вида (3). Пространственная сетка в этом тесте состоит из $40 \times 40 \times 40$ или $80 \times 80 \times 80$ ячеек Йи, так что размер задачи $n = 413\,526$ или $n = 3\,188\,646$ соответственно. После дискретизации вектор начальных значений ${v} \in {{\mathbb{R}}^{n}}$ нормализуется ${v}: = {v}{\text{/}}\left\| {v} \right\|$. Сравнение результатов, полученных на этих двух пространственных сетках, показывает, что разрешение сетки достаточно для этого теста. Конечное время $t = 1$.

Этот тест является сложным для метода подпространства Крылова СО, потому что матрица $A$ сильно несимметрична (можно выбрать такую диагональную матрицу $D$, что ${{D}^{{ - 1}}}AD$ – кососимметричная). Для сильно несимметричных задач, таких как дискретизированные уравнения Максвелла в непроводящей среде, методы подпространства Крылова СО зачастую не являются эффективными (см. [43]). Действительно, другие методы подпространства Крылова с перезапуском оказываются более эффективными в этой тестовой задаче (см. [22]). Кроме того, это – трехмерная векторная задача, где 6 переменных (компоненты $x$, $y$ и $z$ обоих полей) соответствуют каждой ячейкe вычислительной сетки. Поэтому в зависимости от конкретных значений параметров задачи решать линейные системы со сдвинутой матрицей $I + \gamma A$ может быть весьма непросто. Тем не менее для данной конкретной задачи оказывается, что величина $\gamma \left\| A \right\|$ достаточно мала, так что условие (23) не выполняется и даже непредобусловленный метод простых итераций может быть успешно использован для решения сдвинутых линейных систем. В представленных здесь экспериментах для этой цели был использован итерационный метод GMRES(10). Таким образом, мы рассматриваем эту тестовую задачу, чтобы показать возможности предложенного ТНВ перезапуска.

Опыт показывает, что для успешной работы с сильно несимметричными матрицами $A$ методы подпространства Крылова СО должны использовать гораздо меньшие величины сдвига $\gamma $, чем обычно используемое значение $t{\text{/}}10$ (см. [44]). Поэтому мы задаем для $\gamma $ значение $t{\text{/}}80 = 1{\text{/}}80$ в обоих методах перезапуска (при этом ТНВ перезапуск при необходимости может уменьшить это значение). Результаты представлены в табл. 2 и на фиг. 5 и 6. ТНВ перезапуск очевидно превосходит НВ перезапуск как по вычислительным затратам, так и по полученной точности. Потеря точности в НВ перезапуске происходит на первых перезапусках, когда норма невязки оказывается больше требуемой точности по всему временному интервалу. Во избежание потери точности ТНВ перезапуск уменьшает $\gamma $, что не только восстанавливает точность, но и приводит к более быстрому решению сдвинутых линейных систем (напомним, что при меньших значениях $\gamma $ непредобусловленный GMRES сходится быстрее). Более того, уменьшенные значения $\gamma $ при этом приводят к дальнейшему выигрышу эффективности в ТНВ перезапуске. Как видно на фиг. 6, этот выигрыш получается из-за больших временных интервалов $[0,\;\delta ]$, на которых норма невязки оказывается меньше заданной точности (напомним, что временной интервал на каждом перезапуске сокращается с $[0,\;t]$ до $[0,\;t - \delta ]$).