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

ВВЕДЕНИЕ

В работе представлены быстрые прямые алгоритмы реализации метода конечных элементов (МКЭ) на основе тензорных произведений 1D пространств МКЭ $n$-го порядка ($n \geqslant 2$) на прямоугольных параллелепипедах [1] для решения $N$-мерного уравнения типа Пуассона $ - \Delta u + \alpha u = f$ ($N \geqslant 2$) с краевыми условиями Дирихле. Алгоритмы основаны на хорошо известных Фурье-подходах (см., например, [2]–[5] и ссылки в них).

Ключевыми новыми элементами являются детальное описание собственных пар 1D задач на собственные значения для МКЭ высокого порядка и быстрые прямой и обратный алгоритмы разложения по соответствующим собственным векторам, использующие одновременно несколько версий БДПФ (быстрого дискретного преобразования Фурье), связанных не только с узлами сетки, но и с центрами элементов [6]. Это решает давно известную задачу (см., например, [5, c. 271]), и делает полные алгоритмы логарифмически оптимальными по отношению к числу элементов так же, как в случае полилинейных элементов ($n = 1$) или стандартных разностных схем. В краткой форме и в случае $N = 2$ алгоритмы представлены (без доказательств) в [7]. См. также методы другого типа для той же цели в [8], [9].

Отметим, что постоянная $\alpha $ может быть вещественной любого знака или даже комплексной, так что оператор $ - \Delta + \alpha $ может быть незнакоопределенным (в противоположность некоторым стандартным используемым итерационным методам).

Алгоритмы являются быстрыми с возрастанием количества элементов $K$ и на практике (причем более быстрыми, чем при стандартном теоретическом анализе, вследствие векторной обработки данных в обычных современных процессорах) и демонстрируют только слабое снижение эффективности с ростом $n$, начиная со стандартного случая $n = 1$. Для примера, в случае элементов 9-го порядка 2D МКЭ-система уравнений для $K = {{2}^{{20}}}$ элементов, содержащая почти $85 \times {{10}^{6}}$ неизвестных, и 3D МКЭ-система для $K = {{2}^{{18}}}$ элементов, содержащая более чем $190 \times {{10}^{6}}$ неизвестных, решаются за менее чем 2 и 15 мин соответственно на ординарном ноутбуке с использованием кода на Matlab R2016a, детали см. ниже.

Алгоритмы могут быть использованы для многочисленных других приложений, включая общие эллиптические уравнения 2-го порядка (для предобуславливания), для $N$-мерных уравнения теплопроводности, волнового уравнения, нестационарного уравнения Шрёдингера($\alpha $ может быть комплексным) и т.д., поскольку для их неявных дискретизаций по времени обычно получаются задачи рассматриваемого типа на верхнем слое по времени. В комбинации с другими методами они могут быть использованы для некоторых непрямоугольных областей и неравномерных сеток, в частности, с применением сеток, топологически эквивалентных равномерным прямоугольным [10]. Другие стандартные краевые условия также могут быть охвачены, см. краткое описание в [7]. Алгоритмы обладают высокой степенью параллелизуемости, так что они представляются полезными в научных вычислениях.

Подчеркнем также, что Фурье-структура алгоритмов является особенно ценной для решения некоторых задач волновой физики, в частности, использующих нелокальные граничные условия (см., например, [5], [11], [12]), откуда и возник наш собственный интерес.

Работа построена следующим образом. В разд. 1 даны формулировки стандартной 1D МКЭ-задачи на собственные значения и вспомогательных МКЭ-задач на собственные значения на эталонном элементе и внутри него. Основной разд. 2 посвящен описанию собственных значений и собственных векторов 1D МКЭ-задачи на собственные значения и быстрым прямому и обратному алгоритмам разложения по этим собственным векторам. Приложения к решению уравнения типа Пуассона в $N$-мерном прямоугольном параллелепипеде с краевым условием Дирихле описаны в разд. 3. Результаты численных экспериментов для $N = 2$ и 3 детально представлены в разд. 4; все они включают для сравнения известный случай $n = 1$.

1. ФОРМУЛИРОВКА 1D МКЭ ЗАДАЧИ НА СОБСТВЕННЫЕ ЗНАЧЕНИЯ

Сначала подробно рассмотрим МКЭ для простейшей 1D задачи на собственные значения для обыкновенного дифференциального уравнения (ОДУ)

Возьмем равномерную сетку $\mathop {\bar {\omega }}\nolimits_h $ с узлами ${{x}_{j}} = jh$, $j = \overline {0,K} $ (т.е. $0 \leqslant j \leqslant K$) и шагом $h = X{\text{/}}K$. Пусть $H_{h}^{{(n)}}[0,X]$ – МКЭ-пространство кусочно-полиномиальных функций $\varphi \in C[0,X]$ таких, что $\varphi (x) \in {{\mathcal{P}}_{n}}$ при $x \in [{{x}_{{j - 1}}},{{x}_{j}}]$, $j = \overline {1,K} $, с $\varphi (0) = \varphi (X) = 0$; здесь ${{\mathcal{P}}_{n}}$ – пространство многочленов степени не выше $n$, $n \geqslant 2$.

Пусть $S_{K}^{{(n)}}$ – пространство вектор-функций $w$ таких, что ${{w}_{j}} \in R$ при $j = \overline {0,K} $ с ${{w}_{0}} = {{w}_{K}} = 0$ и ${{w}_{{j - 1/2}}} \in {{\mathbb{R}}^{{n - 1}}}$, $j = \overline {1,K} $. Ясно, что $dimS_{K}^{{(n)}} = nK - 1$. Функция $\varphi \in H_{h}^{{(n)}}[0,X]$ однозначно определяется своими значениями в узлах сетки ${{\varphi }_{j}} = \varphi ({{x}_{j}})$, $j = \overline {0,K} $, с ${{\varphi }_{0}} = {{\varphi }_{K}} = 0$, и внутри элементов ${{\varphi }_{{j - 1/2}}} = \{ \varphi ({{x}_{{j - 1}}} + (l{\text{/}}n)h)\} _{{l = 1}}^{{n - 1}}$, $j = \overline {1,K} $, которые формируют элемент из $S_{K}^{{(n)}}$.

Воспользуемся следующей масштабированной операторной формой стандартной МКЭ-дискретизации задачи (1.1)

Здесь $\mathcal{A} = {{\mathcal{A}}^{{\text{T}}}} > 0$ и $\mathcal{C} = {{\mathcal{C}}^{{\text{T}}}} > 0$ – глобальные (масштабированные) операторы (матрицы) жесткости и масс, действующие в $S_{K}^{{(n)}}$ и вместе с $\lambda $ не зависящие от $h$; истинные приближенные собственные значения – это ${{\lambda }_{h}} = 4{{h}^{{ - 2}}}\lambda $.

Пусть $A = \{ {{A}_{{kl}}}\} _{{k,l = 0}}^{n}$ и $C = \{ {{C}_{{kl}}}\} _{{k,l = 0}}^{n}$ – локальные матрицы жесткости и масс, связанные с эталонным элементом ${{\sigma }_{0}} = [ - 1,1]$, со следующими элементами

где $\{ {{e}_{l}}\} _{{l = 0}}^{n}$ – лагранжев базис в ${{\mathcal{P}}_{n}}$ такой, что ${{e}_{l}}( - 1 + (2k){\text{/}}n) = {{\delta }_{{kl}}}$ при $k,l = \overline {0,n} $, а ${{\delta }_{{kl}}}$ – символ Кронекера. Матрицы $A$, $C$ и соответствующий матричный пучок $G(\lambda ): = A - \lambda C$ имеют следующую $3 \times 3$-блочную форму Здесь $\tilde {A}$, $\tilde {C}$ и $\tilde {G}(\lambda ) = \tilde {A} - \lambda \tilde {C}$ – квадратные матрицы порядка $n - 1$ с векторами-столбцами $a,c,g(\lambda ) = a - \lambda c \in {{\mathbb{R}}^{{n - 1}}}$, а $\mathop p\limits^ {{{\kern 1pt} }_{l}}\; \equiv {{(Pp)}_{l}} = {{p}_{{n - l}}}$, $l = \overline {1,n - 1} $ для $p \in {{\mathbb{R}}^{{n - 1}}}$. Матрицы $A$, $C$ и $G(\lambda )$ являются бисимметричными (т.е. симметричными по отношению к их основной и побочной диагоналям).

Отметим, что ${{P}_{{ij}}} = {{\delta }_{{i(n - j)}}}$ и ${{P}^{{\text{T}}}} = {{P}^{{ - 1}}} = P$. Пусть $\mathbb{R}_{e}^{{n - 1}}$ и $\mathbb{R}_{o}^{{n - 1}}$ – подпространства четных и нечетных векторов в ${{\mathbb{R}}^{{n - 1}}}$, т.е. таких, что соответственно $Pp = p$ и $Pp = - p$. Разложение ${{\mathbb{R}}^{{n - 1}}} = \mathbb{R}_{e}^{{n - 1}} \oplus \mathbb{R}_{o}^{{n - 1}}$ (при $n \geqslant 3$) реализуется формулами

Отметим, что $dim\mathbb{R}_{e}^{{n - 1}} = [n{\text{/}}2]$ и $\dim \mathbb{R}_{o}^{{n - 1}} = [(n - 1){\text{/}}2]$, с $\mathbb{R}_{o}^{{n - 1}} = \{ 0\} $ при $n = 2$. Ясно, что Здесь и ниже символ $ \cdot $ означает скалярное произведение векторов в ${{\mathbb{R}}^{{n - 1}}}$.

Теперь задачу (1.2) можно переписать в следующей явной форме

с ${{v}_{0}} = {{v}_{K}} = 0$, .

Рассмотрим также вспомогательные задачи на собственные значения на эталонном элементе ${{\sigma }_{0}}$ и внутри него

где ясно, что $A \geqslant 0$, $C > 0$ и $\tilde {A} = {{\tilde {A}}^{{\text{T}}}} > 0$, $\tilde {C} = {{\tilde {C}}^{{\text{T}}}} > 0$; см. некоторые их свойства в [11] (где была изучена задача, аналогичная (1.6), (1.7) на равномерной сетке на $[0,\infty )$ при $\lambda \in \mathbb{C}$). Обозначим через ${{S}_{n}}$ и ${{\tilde {S}}_{n}}$ их спектры. Пусть $\{ \lambda _{0}^{{(l)}},{{e}^{{(l)}}}\} _{{l = 1}}^{{n - 1}}$ – собственные пары задачи (1.9).

Лемма 1.1. 1. Подпространства $\mathbb{R}_{e}^{{n - 1}}$ и $\mathbb{R}_{e}^{{n - 1}}$ инвариантны по отношению к $\tilde {A}$ и $\tilde {C}$. Поэтому каждый собственный вектор в $\{ {{e}^{{(l)}}}\} _{{l = 1}}^{{n - 1}}$ может быть выбран четным или нечетным. Также $\lambda _{0}^{{(l)}} > 0$, $l = \overline {1,n - 1} $.

2. Аналогичные свойства верны для задачи (1.8) за исключением наличия одного простого нулевого собственного значения.

Доказательство. Любая бисимметричная матрица $B$ порядка $n - 1$ коммутирует с $P$, т.е.

откуда следует основной результат п. 1. Свойство $\lambda _{0}^{{(l)}} > 0$, $l = \overline {1,n - 1} $, хорошо известно.

Для п. 2 рассуждение аналогично с учетом того, что $A{{(1 \ldots 1)}^{{\text{T}}}} = 0$ (относительно простоты $\lambda = 0$ см. [11, утверждение 5]).

Прямым вычислением можно проверить, что все собственные значения в ${{S}_{n}}$ и ${{\tilde {S}}_{n}}$ являются простыми и ${{S}_{n}} \cap {{\tilde {S}}_{n}} = \not {0}$ по крайней мере при $1 \leqslant n \leqslant 9$, см. [11].

Для малых $n$ можно найти ${{S}_{n}}$ и ${{\tilde {S}}_{n}}$ аналитически (с помощью Mathematica), в частности, ${{\tilde {S}}_{2}} = \{ 2.5\} $, ${{\tilde {S}}_{3}} = \{ 2.5,10.5\} $, ${{\tilde {S}}_{4}} = \{ 14 \pm \sqrt {133} ,10.5\} $ и ${{\tilde {S}}_{5}} = \{ 14 \pm \sqrt {133} ,30 \pm 9\sqrt 5 \} $ (повторяемость собственных значений не случайна, см. [11]).

Выберем $\{ {{e}^{{(l)}}}\} _{{l = 1}}^{{n - 1}}$ как указано в лемме 1.1, с использованием масштабирования $\tilde {C}{{e}^{{(l)}}} \cdot {{e}^{{(l)}}} = 1$.

Лемма 1.2. Пусть $\tilde {G}(\lambda )p = - g(\lambda )$, где . Пусть векторы $a$ и $c$ разложены как

См. $\tilde {G}(\lambda )$, $g(\lambda )$, $a$ и $c$ в (1.3). Тогда верны следующие формулы:

Доказательство. Для разложений $p = \sum\nolimits_{l = 1}^{n - 1} \,{{p}_{l}}{{e}^{{(l)}}}$ и (1.11) имеем

откуда легко следует результат.

2. РЕШЕНИЕ 1D МКЭ-ЗАДАЧ НА СОБСТВЕННЫЕ ЗНАЧЕНИЯ И СООТВЕТСТВУЮЩИЕ АЛГОРИТМЫ, ОСНОВАННЫЕ НА БДПФ

Сделаем следующее

Предположение.

$(A)\;\;{\text{собственные}}\;{\text{значения}}\;{\text{в}}\;{{S}_{n}}\quad {\text{и}}\quad {{\tilde {S}}_{n}}\;{\text{простые}}\quad {\text{и}}\quad {{S}_{n}} \cap {{\tilde {S}}_{n}} = \not {0}.$

Напомним, что это верно по крайней мере при $2 \leqslant n \leqslant 9$ (ниже в разд. 4 это проверяется вплоть до $n = 21$).

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

где , с параметром $\theta $. Его решение эквивалентно нахождению корней многочлена степени не выше $n$, см. [11]. В частности, в силу леммы 1.2 это уравнение может быть переписано в виде Здесь $\mathop a\limits^ {{{\kern 1pt} }^{{(l)}}}\; = \mathop a\limits^ {\kern 1pt} \; \cdot {{e}^{{(l)}}}$ и $\mathop c\limits^ {{{\kern 1pt} }^{{(l)}}}\; = \mathop c\limits^ {\kern 1pt} \; \cdot {{e}^{{(l)}}}$. Более того, при $2 \leqslant n \leqslant 9$ вычисления помогают подтвердить, что векторы ${{e}^{{(l)}}}$ являются четными/нечетными соответственно для нечетных/четных $l$ при условии, что $\lambda _{0}^{{(1)}} < \ldots < \lambda _{0}^{{(n - 1)}}$; следовательно, тогда $\mathop a\limits^ {{{\kern 1pt} }^{{(l)}}}\; = {{( - 1)}^{l}}{{a}^{{(l)}}}$ и $\mathop c\limits^ {{{\kern 1pt} }^{{(l)}}}\; = {{( - 1)}^{l}}{{c}^{{(l)}}}$, $l = \overline {1,n - 1} $.

Введем простейшее скалярное произведение в $S_{K}^{{(n)}}$ и соответствующий квадрат $\mathcal{C}$-нормы

В следующей теореме описываются собственные значения и собственные векторы задачи (1.2).

Теорема 2.1. 1. Спектр задачи (1.2) состоит из ${{\tilde {S}}_{n}}$ и чисел $\{ \lambda _{k}^{{(l)}}\} _{{l = 1}}^{n}$, которые являются всеми $n$ (и  всеми положительными вещественными) решениями уравнения (2.2) с $\theta = {{\theta }_{k}}: = cos\tfrac{{\pi k}}{K}$ при $k = \overline {1,K - 1} $. Числа $\{ \lambda _{k}^{{(l)}}\} _{{l = 1}}^{n}$ отличаются от $\{ \lambda _{0}^{{(l)}}\} _{{l = 1}}^{{n - 1}}$ и различны при фиксированном $k$.

2. Собственному значению $\lambda _{0}^{{(l)}}$ отвечает собственный вектор

при $l = \overline {1,n - 1} $. Здесь ${{( - P)}^{{j - 1}}}e = {{( - 1)}^{{j - 1}}}e$ для четного $e$, ${{( - P)}^{{j - 1}}}e = e$ для нечетного $e$.

3. Собственному значению $\lambda _{k}^{{(l)}}$ отвечает собственный вектор

где ${{s}_{{k,j}}}: = sin\tfrac{{\pi kj}}{K}$, а $p_{k}^{{(l)}} \in {{\mathbb{R}}^{{n - 1}}}$ служит решением невырожденной алгебраической системы уравнений $\tilde {G}(\lambda _{k}^{{(l)}})p_{k}^{{(l)}} = - g(\lambda _{k}^{{(l)}})$ при $k = \overline {1,K - 1} $, $l = \overline {1,n} $.

4. Введенные собственные векторы являются $\mathcal{C}$-ортогональными, т.е.

для всех $k,\tilde {k} \in \overline {0,K - 1} $, $l \in \overline {1,n - {{\delta }_{{k0}}}} $, $\tilde {l} \in \overline {1,n - {{\delta }_{{\tilde {k}0}}}} $ таких, что $k \ne \tilde {k}$ и/или $l \ne \tilde {l}$.

Как следствие, они образуют базис в $S_{K}^{{(n)}}$, т.е. любой $w \in S_{K}^{{(n)}}$ может быть однозначно разложен как

Доказательство. Будем различать два случая.

1. Пусть сначала $\lambda \in {{\tilde {S}}_{n}}$ и $e$ удовлетворяет (1.9). Тогда для любого $j = \overline {1,K} $ с использованием уравнения (1.7) получим

Ясно, что Поскольку по предположению (A), то имеем $g(\lambda ) \cdot e \ne 0$ или $G(\lambda )\left( {\begin{array}{*{20}{c}} 0 \\ e \\ 0 \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {g(\lambda ) \cdot e} \\ 0 \\ {\mathop g\limits^ (\lambda ) \cdot e} \end{array}} \right).$. В силу предположения (A) и леммы 1.1 $\lambda $ является простым в ${{\tilde {S}}_{n}}$ и $e$ четен или нечетен. Соответственно или $\mathop g\limits^ (\lambda ) \cdot e = g(\lambda ) \cdot e \ne 0$ и ${{{v}}_{j}} = - {{{v}}_{{j - 1}}}$, или $\mathop g\limits^ (\lambda ) \cdot e = - g(\lambda ) \cdot e \ne 0$ и ${{{v}}_{j}} = {{{v}}_{{j - 1}}}$. Поскольку ${{{v}}_{0}} = 0$, то в обоих случаях получим Поэтому уравнение (1.7) сводится к $\tilde {G}{{{v}}_{{j - 1/2}}} = 0$, откуда следует, что ${{v}_{{j - 1/2}}} = {{c}_{{j - 1/2}}}e$.

Теперь уравнение (1.6) сводится к

поэтому ${{c}_{{j + 1/2}}} = - {{c}_{{j - 1/2}}}$ при четном $e$ или ${{c}_{{j + 1/2}}} = {{c}_{{j - 1/2}}}$ при нечетном $e$. Следовательно, искомый собственный вектор $v$ удовлетворяет (2.7) и ($v$ определен с точностью до ненулевого постоянного множителя). Таким образом приходим к собственным векторам (2.3).

2. Пусть теперь . Тогда из уравнения (1.7) получим

Подставляя эту формулу в уравнение (1.6), выведем трехточечное уравнение где

Непосредственно проверяются следующие равенства

Если ${{\hat {g}}_{n}}(\lambda ) = {{\hat {g}}_{0}}(\lambda ) = 0$, то эти равенства означают, что $\lambda $ является по крайней мере двукратным собственным значением задачи (1.8), что противоречит предположению (A).

Если $\mathop {\hat {g}}\nolimits_n (\lambda ) = 0$ и $\mathop {\hat {g}}\nolimits_0 (\lambda ) \ne 0$, то уравнение (2.9) вместе с (2.8) приводит к ${v} = 0$, поэтому такое $\lambda $ не удовлетворяет (1.2).

Значит, $\mathop {\hat {g}}\nolimits_n (\lambda ) \ne 0$ и уравнение (2.9) упрощается до

с функцией $\hat {\gamma }(\lambda ) = \mathop {\hat {g}}\nolimits_0 (\lambda ){\text{/}}\mathop {\hat {g}}\nolimits_n (\lambda )$ (см. ее также в (2.1)). Поскольку ${{{v}}_{0}} = {{{v}}_{n}} = 0$, то можно использовать разложение и ввести вектор его коэффициентов ${\mathbf{\tilde {v}}}: = ({{\tilde {v}}_{1}}, \ldots ,{{\tilde {v}}_{{K - 1}}})$. Применение этого разложения в (2.10) дает Ясно, что это равенство выполнено для некоторого ${\mathbf{\tilde {v}}} \ne 0$ тогда и только тогда, когда

Отметим, что ${\mathbf{\tilde {v}}} = 0$ эквивалентно ${v} = 0$ в $S_{K}^{{(n)}}$ (с учетом формулы (2.8)). Поэтому $\lambda $ удовлетворяет (2.12); более того, ${{\tilde {v}}_{j}} = {{\delta }_{{kj}}}$ и, как следствие, ${{{v}}_{j}} = {{s}_{{k,j}}}$, $j = \overline {0,K} $, а также

см. (2.8) (все последние три равенства выполнены с точностью до одинакового ненулевого множителя). Таким образом приходим к собственным векторам (2.4).

Общее количество собственных значений (с учетом возможной кратности) равно $dimS_{K}^{{(n)}} - (n - 1) = n(K - 1)$. Максимальное количество корней алгебраических уравнений (2.12) при всех $k$ то же самое, поэтому каждое уравнение (2.12) должно иметь ровно $n$ различных корней (при фиксированном $k$ записанный вектор $v$ определяется значением $\lambda $ однозначно).

3. Свойство (2.5) заведомо выполняется для собственных векторов $s_{k}^{{(l)}}$ и $s_{{\tilde {k}}}^{{(\tilde {l})}}$, отвечающих различным собственным значениям задачи (1.2), в частности, при $k = 0$ и $\tilde {k} \ne 0$, или $k = \tilde {k}$ и $l \ne \tilde {l}$.

Оставшийся случай будет рассмотрен ниже в следствии 2.1 связанной с этим леммы 2.1.

Отметим, что: 1) векторы $s_{0}^{{(l)}}$ используются только для описания алгоритма, и только векторы ${{e}^{{(l)}}}$ применяются в его реализации; 2) $s_{{k,j}}^{{(l)}}$ не зависят от $l$; 3) векторы $p_{k}^{{(l)}}$ не зависят от $j$ и могут быть также вычислены в силу леммы 1.2.

Лемма 2.1. Пусть $w \in S_{K}^{{(n)}}$ и ${{w}_{{j - 1/2}}} = q{{w}_{{j - 1}}} + \mathop q\limits^ {{w}_{j}}$, $j = \overline {1,K} $, для некоторого $q \in {{\mathbb{R}}^{{n - 1}}}$. Тогда

при $k = \overline {1,K - 1} $, $l = \overline {1,n} $, где ${{({{s}_{k}},w)}_{{{{\omega }_{h}}}}}: = \sum\nolimits_{j = 1}^{K - 1} \,{{s}_{{k,j}}}{{w}_{j}}$.

Доказательство. 1. Вспомнив обозначение (1.3), для любых $v,w \in S_{K}^{{(n)}}$, имеем

2. Согласно формулам (2.15) и (2.3) при четном ${{e}^{{(l)}}}$ получим

поскольку $(c - \mathop c\limits^ ) \cdot {{e}^{{(l)}}} = 0$ и $\tilde {C}{{e}^{{(l)}}} \cdot (q - \mathop q\limits^ ) = 0$ для всех $c,q \in {{\mathbb{R}}^{{n - 1}}}$, а также ${{w}_{0}} = {{w}_{K}} = 0$.

При нечетном ${{e}^{{(l)}}}$ аналогично получим

поскольку $(c + \mathop c\limits^ ) \cdot {{e}^{{(l)}}} = 0$ и $\tilde {C}{{e}^{{(l)}}} \cdot (q + \mathop q\limits^ ) = 0$ для всех $c,q \in {{\mathbb{R}}^{{n - 1}}}$, а также ${{w}_{0}} = {{w}_{K}} = 0$. Равенство (2.13) доказано.

3. Формула (2.15) вместе с (2.4) и (2.5) влекут равенство

В силу формулы

сначала выводим Затем с использованием формул (2.5) и (2.10) получим В силу свойств ${{s}_{{k,0}}} = {{s}_{{k,K}}} = 0$, ${{w}_{0}} = {{w}_{K}} = 0$ и формул (2.17) далее выводим Сложив (2.18) и (2.19), докажем (2.14).

Следствие 2.1. Верно свойство ортогональности (2.5) из теоремы 2.1, п. 4.

Доказательство. Остается рассмотреть случай $k,\tilde {k} \in \overline {1,K - 1} $ и $k \ne \tilde {k}$. Поскольку тогда ${{({{s}_{k}},{{s}_{{\tilde {k}}}})}_{{{{\omega }_{h}}}}} = 0$, то результат следует из (2.14).

Вычисление $w \in S_{K}^{{(n)}}$ по коэффициентам ${{w}_{{kl}}}$ разложения (2.6) назовем обратным ${{F}_{n}}$-преобразованием, а вычисление коэффициентов ${{w}_{{kl}}}$ по $w \in S_{K}^{{(n)}}$ – прямым ${{F}_{n}}$-преобразованием. Рассмотрим также связанное с этим разложение $y \in S_{K}^{{(n)}}$ вида

и вычисление коэффициентов ${{\tilde {y}}_{{kl}}}$ по $y \in S_{K}^{{(n)}}$, которое назовем прямым $F{{C}_{n}}$-преобразованием.

Опишем их быструю БДПФ-реализацию.

Теорема 2.2. 1. Обратное ${{F}_{n}}$-преобразование может быть реализовано в соответствии со следующими формулами

где ${{d}_{{k,e}}}$ и ${{d}_{{k,o}}}$ – соответственно четные и нечетные компоненты вектора ${{d}_{k}}: = \sum\nolimits_{l = 1}^n \,{{w}_{{kl}}}p_{k}^{{(l)}}$. Заметим, что ${{( - P)}^{{j - 1}}}e = e$ при нечетном $j$ и при четном $j$ для любого $e \in {{\mathbb{R}}^{{n - 1}}}$.

Набор $\{ {{w}_{j}}\} _{{j = 1}}^{{K - 1}}$ можно вычислить с помощью стандартного обратного БДПФ по синусам. Набор $\{ {{w}_{{j - 1/2}}}\} _{{j = 1}}^{K}$ можно вычислить с помощью $n - 1$ модифицированного обратного БДПФ, связанного с центрами элементов, в количестве $[n{\text{/}}2]$ по синусам и $[(n - 1){\text{/}}2]$ по косинусам с использованием продолжений ${{d}_{{K,e}}}: = 0$ и ${{d}_{{0,o}}}: = 0$, см. алгоритмы DST-I, DST-III и DCT-III в [6].

2. Прямое $F{{C}_{n}}$-преобразование можно реализовать на основе стандартной формулы

Здесь, во-первых, при $k = 0$, $l = \overline {1,n - 1} $ верны формулы

Во-вторых, при $k = \overline {1,K - 1} $, $l = \overline {1,n} $ верны следующие формулы:

Отметим, что суммы в формуле (2.25) не зависят от $l$. Набор всех этих коэффициентов можно вычислить с помощью $n$ стандартных прямых БДПФ по синусам.

3. Аналогично п. 2, прямое ${{F}_{n}}$-преобразование можно реализовать на основе стандартной формулы

Здесь при $k = 0$, $l = \overline {1,n - 1} $ верна следующая формула:

При $k = \overline {1,K - 1} $, $l = \overline {1,n} $ применима формула (2.25) с $y: = \mathcal{C}w$. С другой стороны, верна также следующая формула

где $q_{{k,e}}^{{(l)}}$ и $q_{{k,o}}^{{(l)}}$ – соответственно четная и нечетная компоненты вектора $q_{k}^{{(l)}}: = \tilde {C}p_{k}^{{(l)}} + c$. Снова все эти коэффициенты можно вычислить с помощью $n$ стандартных прямых БДПФ по синусам.

Доказательство. 1. Пусть коэффициенты ${{w}_{{kl}}}$ разложения (2.6) известны. В соответствии с первыми из формул (2.3) и (2.4) значения $w$ для целых индексов в (2.6) сводятся к (2.21).

Чтобы вычислить $w$ для полуцелых индексов, преобразуем вторую сумму в (2.6). В силу разложения (1.4) перепишем вторую формулу (2.4) в виде

Затем, применив также вторую формулу (2.3), получим формулу (2.22).

2. Теперь рассмотрим вычисление коэффициентов в разложении (2.20) при заданном $y \in {{S}_{K}}$. В силу свойства ортогональности (2.5), они сначала могут быть выражены в виде (2.23) при $k = 0$, $l = \overline {1,n - 1} $ и $k = \overline {1,K - 1} $, $l = \overline {1,n} $.

Формулы (2.15) и (2.3) влекут равенства

Лемма 2.1 немедленно влечет формулу (2.26), поскольку ${{({{s}_{k}},{{s}_{k}})}_{{{{\omega }_{h}}}}} = K{\text{/}}2$.

В силу формул (2.3) для числителя формулы (2.23) при $k = 0$ можно записать

В силу формул (2.4) для того же числителя при $k = \overline {1,K - 1} $ получим

Поэтому сдвинув на 1 индекс во второй из этих сумм, применив тождество ${{a}_{1}}{{b}_{1}} + {{a}_{2}}{{b}_{2}} = $ $ = 0.5({{a}_{1}} + {{a}_{2}})({{b}_{1}} + {{b}_{2}}) + $ $0.5({{a}_{1}} - {{a}_{2}})({{b}_{1}} - {{b}_{2}})$ и вспомнив разложение (1.4), выведем Поскольку также то в итоге установим формулу (2.25).

3. Благодаря свойству ортогональности (2.5) справедлива формула (2.27).

В силу формул (2.3), (2.15), а также ${{\tilde {C}}^{{\text{T}}}} = \tilde {C}$ и $P = {{P}^{{\text{T}}}}$ для ее числителя при $k = 0$ можно записать

В силу формул (2.16), (2.18) для того же числителя при $k = \overline {1,K - 1} $ получим

где $q_{k}^{{(l)}} = \tilde {C}p_{k}^{{(l)}} + c$. Преобразовав последнюю сумму таким же образом, как второе и третье слагаемые в (2.30), установим (2.29).

3. ПРИЛОЖЕНИЯ К ОБОБЩЕННОМУ УРАВНЕНИЮ ПУАССОНА

Для начала обратимся к простой 1D краевой задаче для ОДУ

где $\alpha \ne - {{(\pi i{\text{/}}X)}^{2}}$, $i \in \mathbb{N}$. Ее МКЭ-дискретизация имеет операторную форму где ${{f}^{h}} \in S_{K}^{{(n)}}$ – это МКЭ-усреднение $f$. Ее решение можно записать в виде разложения типа (2.6), где $\tilde {f}_{{kl}}^{h}$ – коэффициенты разложения типа (2.20) вектора ${{f}^{h}}$; напомним, что ${{\delta }_{{k0}}}$ – символ Кронекера. Формула корректна для невырожденных знаменателей, что имеет место при всех $h$, если $\alpha \in \mathbb{R}$ и $\alpha = {\text{const}} > - {{(\pi {\text{/}}X)}^{2}}$, или по крайней мере для достаточно малых $h$ в противном случае, или если $\alpha \in \mathbb{C}{\backslash }\mathbb{R}$. В последнем случае коэффициенты разложения становятся комплексными, а функцию $f$, а поэтому и вектор ${{f}^{h}}$ и его коэффициенты разложения $\tilde {f}_{{kl}}^{h}$ также можно считать комплексными. Подобное относится и к многомерному случаю ниже.

Теперь рассмотрим детально решение $N$-мерной ($N \geqslant 2$) краевой задачи

где $\Delta $ – оператор Лапласа, а $\alpha = {\text{const}}$ не принадлежит его спектру (при тех же краевых условиях).

Введем пространство $H_{{{{h}_{1}}}}^{{({{n}_{1}})}}[0,{{X}_{1}}] \otimes \ldots \otimes H_{{{{h}_{N}}}}^{{({{n}_{N}})}}[0,{{X}_{N}}]$ кусочно-полиномиальных функций в $\bar {\Omega }$, где ${{h}_{i}} = {{X}_{i}}{\text{/}}{{K}_{i}}$ и ${{n}_{i}} \geqslant 2$, $i = \overline {1,N} $. Пусть ${\mathbf{K}} = ({{K}_{1}}, \ldots ,{{K}_{N}})$ и ${\mathbf{n}} = ({{n}_{1}}, \ldots ,{{n}_{N}})$.

Введем также пространство вектор-функций $S_{{\mathbf{K}}}^{{({\mathbf{n}})}} = S_{{{{K}_{1}}}}^{{({{n}_{1}})}} \otimes \ldots \otimes S_{{{{K}_{N}}}}^{{({{n}_{N}})}}$. Например, при $N = 2$, эти функции являются числами для индексов $({{j}_{1}},{{j}_{2}})$, ${{j}_{1}} = \overline {0,{{K}_{1}}} $, ${{j}_{2}} = \overline {0,{{K}_{2}}} $ и векторами из ${{\mathbb{R}}^{{{{n}_{1}} - 1}}}$, ${{\mathbb{R}}^{{{{n}_{2}} - 1}}}$ и ${{\mathbb{R}}^{{({{n}_{1}} - 1) \times ({{n}_{2}} - 1)}}}$ соответственно для индексов

а также нулевыми векторами для ${{j}_{1}} = 0,{{K}_{1}}$ и ${{j}_{2}} = 0,{{K}_{2}}$. Аналогично 1D случаю существует естественный изоморфизм между функциями из $H_{{{{h}_{1}}}}^{{({{n}_{1}})}}[0,{{X}_{1}}] \otimes \ldots \otimes H_{{{{h}_{N}}}}^{{({{n}_{N}})}}[0,{{X}_{N}}]$ и векторами из $S_{{\mathbf{K}}}^{{({\mathbf{n}})}}$.

МКЭ-дискретизация задачи (3.4) может быть записана в следующей операторной форме

где ${{\mathcal{A}}_{i}}$ и ${{\mathcal{C}}_{i}}$ – это версии введенных выше операторов $\mathcal{A}$ и $\mathcal{C}$, действующие по переменной ${{x}_{i}}$ (зависящие от ${{K}_{i}}$ и ${{n}_{i}}$), $i = \overline {1,N} $, а ${{f}^{h}} \in S_{{\mathbf{K}}}^{{({\mathbf{n}})}}$ – МКЭ-усреднение $f$.

Напомним, что случай неоднородных краевых условий Дирихле $u(x) = b(x)$ на $\partial \Omega $ в (3.4) мог бы быть легко охвачен сведением к (3.5) с модифицированным вектором ${{f}^{h}}$, зависящим от аппроксимации ${{b}^{h}}$ для $b$.

Чтобы вычислить ее решение, ${{F}_{n}}$-преобразования из теоремы 2.2 могут быть применены двояко.

(a) Рассмотрим кратное разложение ${{f}^{h}} \in S_{{\mathbf{K}}}^{{({\mathbf{n}})}}$ типа (2.20)

Тогда разложение решения имеет следующий вид: Здесь $\{ \lambda _{{i,{{k}_{i}}}}^{{({{l}_{i}})}},s_{{i,{{k}_{i}}}}^{{({{l}_{i}})}}\} $ являются версиями введенных выше собственных пар $\{ \lambda _{k}^{{(l)}},s_{k}^{{(l)}}\} $ по отношению к ${{x}_{i}}$. Формула корректно определена при всех ${\mathbf{h}}: = ({{h}_{1}}, \ldots ,{{h}_{N}})$, если $\alpha \in \mathbb{R}$ и $\alpha > - {{\pi }^{2}}(X_{1}^{{ - 2}} + \ldots + X_{N}^{{ - 2}})$, или по крайней мере при достаточно малых ${\mathbf{h}}$ в противном случае, или если $\alpha \in \mathbb{C}{\backslash }\mathbb{R}$.

Алгоритм (a) состоит из двух достаточно стандартных шагов:

(1) нахождение коэффициентов разложения (3.6) для ${{f}^{h}}$ (посредством прямых $F{{C}_{n}}$-преобразований по ${{x}_{1}}$,…, ${{x}_{N}}$);

(2) нахождение $v$ по коэффициентам его разложения (3.7) (посредством обратных ${{F}_{n}}$-преобразований по ${{x}_{1}}$,…, ${{x}_{N}}$).

(б) Рассмотрим разложение ${{f}^{h}}$ типа (2.20) по ${{x}_{2}}$,…, ${{x}_{N}}$, т.е.

теперь с коэффициентами $\tilde {f}_{{{{k}_{2}}{{l}_{2}}, \ldots ,{{k}_{N}}{{l}_{N}}}}^{h} \in S_{{{{K}_{1}}}}^{{({{n}_{1}})}}$. Тогда коэффициенты ${{v}_{{kl}}} \in S_{{{{K}_{1}}}}^{{({{n}_{1}})}}$ аналогичного разложения решения $v \in S_{{\mathbf{K}}}^{{({\mathbf{n}})}}$, а именно служат решениями 1D задач по ${{x}_{1}}$Их матрицы симметричны и положительно определены, если $\alpha \in \mathbb{R}$ и $\alpha > - {{\pi }^{2}}(X_{1}^{{ - 2}} + \ldots + X_{N}^{{ - 2}})$, либо невырожденны по крайней мере при достаточно малых ${\mathbf{h}}$ в противном случае, или если $\alpha \in C{\backslash }\mathbb{R}$.

Алгоритм (б) состоит из трех достаточно стандартных шагов:

(1) нахождение коэффициентов разложения (3.8) для ${{f}^{h}}$ (посредством прямых $F{{C}_{n}}$-преобразований по ${{x}_{2}}$, …, ${{x}_{N}}$);

(2) решение набора независимых 1D задач (3.10) для коэффициентов разложения $v$;

(3) нахождение $v$ по коэффициентам его разложения (3.9) (посредством обратных ${{F}_{n}}$-преобразований по ${{x}_{2}}$, …, ${{x}_{N}}$).

Реализация алгоритмов (a) и (б) требует соответственно $O({{K}_{1}} \ldots {{K}_{N}}lo{{g}_{2}}({{K}_{1}} \ldots {{K}_{N}}))$ и $O({{K}_{1}} \ldots {{K}_{N}}lo{{g}_{2}}({{K}_{2}} \ldots {{K}_{N}}))$ арифметических действий. Здесь не анализируется зависимость от ${\mathbf{n}}$, но она умеренная согласно численным экспериментам в следующем разделе.

Важно, что они могут быть применены для решения разнообразных эволюционных уравнений в частных производных (УрЧП) таких, как уравнение теплопроводности, волновое уравнение или уравнение Шрёдингера, поскольку для их неявных по времени дискретизаций обычно получаются задачи типа (3.5) на верхнем слое по времени.

Кроме того, алгоритм (б) непосредственно обобщается на случаи более общих уравнений, чем в (3.4), с коэффициентами зависящими от ${{x}_{1}}$ (что существенно, в частности, в полярных или цилиндрических координатах), различных краевых условий при ${{x}_{1}} = 0,{{X}_{1}}$ и неравномерной сетки по ${{x}_{1}}$ [3]. Он может быть также применен для сведения 3D задач в цилиндрической области к набору независимых 2D задач в основании цилиндра.

4. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ

1. Сначала проверим, что собственные значения каждой из задач (1.8) и (1.9) хорошо разделены. Определим их спектральные щели как

где ${{S}_{n}} = :\{ {{\lambda }_{{0(l)}}}\} _{{l = 1}}^{{n + 1}}$, и представим ${{\delta }_{n}}$ и ${{\tilde {\delta }}_{n}}$ на фиг. 1 (слева). Слагаемые ${{\pi }^{2}}{\text{/}}4$ и $3{{\pi }^{2}}{\text{/}}4$ являются спектральными щелями (фактически щелями между двумя минимальными собственными значениями) соответствующих задач для ОДУ, см. [11]. Для обеих величин ${{\delta }_{n}}$ и $\mathop {\tilde {\delta }}\nolimits_n $ наблюдается убывание и быстрое стремление к 0 с возрастанием $n$. Проверено также, что ${{S}_{n}} \cap {{\tilde {S}}_{n}} = \not {0}$ при всех $2 \leqslant n \leqslant 21$.

Представим также спектральные числа обусловленности $\tilde {A}$ и $\tilde {C}$ на фиг. 1 (справа) и обратим внимание на их быстрый рост с возрастанием $n$ (к сожалению).

Отметим, что все вычисления были выполнены на ординарном ноутбуке ASUS-U36S с Intel Core i3-2350M CPU 2.3 GHz, 8 Gb, Win 10 x64. Для реализации алгоритмов были разработаны коды на Matlab R2016a, и подчеркнем, что при этом были применены несколько базовых и продвинутых методов векторизации кодов [13], чтобы заметно ускорить их выполнение.

В случае 1D задачи для ОДУ (3.2) при $\alpha = 1$ дадим числа обусловленности глобальной матрицы МКЭ $4{{h}^{{ - 2}}}\mathcal{A} + \alpha \mathcal{C}$ при $\alpha = 1$ и ее локальной версии в зависимости от $K = 2,4,...,1024$ при $n = \overline {1,9} $ на фиг. 2 для последующих ссылок. Отметим их быстрый рост с возрастанием $K$.

Далее проанализируем время выполнения ${{F}_{n}}$-преобразований в зависимости от выбора $K$ с его возрастанием до весьма высоких значений ${{2}^{{20}}}$. Рассмотрим три способа выбора $K$: степени двойки ($K = {{2}^{p}}$), простые числа и составные числа (не степени двойки). Простые и составные значения $K$ выбираются случайным образом между двумя последовательными степенями двойки, и сравнение выполняется в стиле [14]. Используются внешние функции comp_dst для DST-I и DST-III и comp_dct для DCT-III из Large Time-Frequency Analysis Toolbox [15], основанном на библиотеке FFTW [16], которая также используется в Matlab-функции fft для БДПФ.

Времена выполнения для DST-I и DST-III, а также наших прямого и обратного ${{F}_{n}}$-преобразований в случае $n = 5$ (для определенности) даны соответственно на фиг. 3 и 4. Время выполнения для нахождения пар $\{ \lambda _{k}^{{(l)}},p_{k}^{{(l)}}\} $ не принимается во внимание, поскольку оно становится пренебрежимо малым при многократном применении преобразований при фиксированном $K$, как это требуется ниже. Для наших преобразований времена выполнения соответственно выше (чем для DST-I и DST-III) из-за кратного использования БДПФ и некоторых дополнительных вычислений. Обратное ${{F}_{n}}$-преобразование требует меньше времени, чем прямое. Кроме того, наилучшие результаты получаются в основном для $K = {{2}^{p}}$, хотя это не так для DST-I (к счастью, оно применяется в оптимальном случае $K = {{2}^{p}} - 1$). Для двух других выборов $K$ времена выполнения похуже, но близки, и разница между всеми ними меньше, чем в случае как DST-I, так и DST-III. Эти результаты положительного плана. Отметим, что данные на фиг. 4 могут быть аппроксимированы линейными функциями при ${{2}^{5}} \leqslant K \leqslant {{2}^{{10}}}$ и ${{2}^{{10}}} \leqslant K \leqslant {{2}^{{20}}}$, но с довольно плоским наклоном в первом случае и заметно более крутым во втором.

Замечание 4.1. Последнее явление объясняется продвинутой архитектурой современных процессоров, включая кэш-память, потоковые SIMD-расширения, продвинутые векторные расширения набора команд и т.д., а также применением указанных выше реализаций БДПФ высокого качества.

2. Наши основные вычислительные результаты относятся к решению задачи (3.4) как для $N = 2$, так и $N = 3$, при $\alpha = 1$ и ${{X}_{i}} = 1$. Возьмем для простоты ${{K}_{i}} = K$ и ${{n}_{i}} = n$. Для вычисления ${{f}^{h}}$ применим кратные квадратурные формулы Гаусса с $n + 1$ узлами по ${{x}_{i}}$. Собственные пары ${\text{\{ }}\lambda _{k}^{{(l)}},p_{k}^{{(l)}}{\text{\} }}$ 1D задач вычисляются с четверной точностью (с использованием Mathematica) для повышения устойчивости по отношению к ошибкам округления. Отметим, что система (3.5) содержит ${{(Kn - 1)}^{N}}$ неизвестных. Для сравнения включим в анализ простейший и известный случай $n = 1$, реализованный в коде единым образом.

Рассмотрим сначала 2D случай ($N = 2$) и возьмем гладкое точное решение $u(x): = sin(2\pi {{x}_{1}})sin(3\pi {{x}_{2}})\cosh (\sqrt 2 {{x}_{1}} - {{x}_{2}})$. Вычислим МКЭ-решение для различных значений $n$ и $K$ для того, чтобы изучить поведение его (абсолютной) погрешности, см. фиг. 5 и табл. 1 (где значения ${{R}_{C}}$ меньше 3 опущены). Поведение погрешности в равномерной сеточной норме для алгоритма (a) стандартно: скорость убывания ${{R}_{C}}$ в основном пропорциональна ${{(n + 1)}^{2}}$, за исключением случая $n = 2$, когда она выше и аналогична $n = 3$ (такое исключение ранее уже отмечалось в [17]). Любопытно, что при $n = 9$ и минимальном $K = 2$ погрешность уже меньше, чем при $n = 1$ и максимальном $K = 1024$. Разумеется, погрешность не может быть лучше, чем уровень ошибок округления, который достигается тем быстрее, чем выше $5 \leqslant n \leqslant 9$. Наблюдается почти полное отсутствие влияния ошибок округления.

Для алгоритма (б) ситуация аналогична, но только до уровня погрешности $ \sim $${{10}^{{ - 11}}}$. Как только это значение достигается при фиксированном $3 \leqslant n \leqslant 9$, далее наблюдается рост погрешности с возрастанием $K$, что означает ощутимое влияние ошибок округления. Оно происходит из-за соответствующего роста чисел обусловленности в системе (3.10) с возрастанием $K$ или $n$, см. выше фиг. 2. Как следствие, алгоритм (a) предпочтительнее (б), если требуется очень высокая точность. Отметим, что при использовании только двойной точности вычислений уровень наилучшей погрешности становится $ \sim {\kern 1pt} {{10}^{{ - 12}}}{\kern 1pt} - {\kern 1pt} {{10}^{{ - 13}}}$ и результаты остаются устойчивыми для алгоритма (a), но они практически не меняются для алгоритма (б). Поэтому двойная точность всех вычислений возможна, если указанная точность достаточна (что имеет место в большинстве приложений). Отметим, что устойчивость метода должен повысить переход в разделе 1 от лагранжева базиса к ортогональным многочленам.

Проанализируем также времена выполнения обоих алгоритмов (с использованием кратного запуска кода и медианного времени их выполнения), см. фиг. 6 и табл. 2 и 3 для одних и тех же $K$ и $n$. Очевидно, что они не зависят от указанного выше специального выбора $u$. Время выполнения для вычисления пар $\{ \lambda _{k}^{{(l)}},p_{k}^{{(l)}}\} $ во внимание не принимается, поскольку рассматривается случай, когда они вычисляются заранее и хранятся (напомним, что они не зависят от данных задачи для УрЧП (3.4)). Обратим внимание на слабо сверхлинейную зависимость времен от $K$ и их слабый монотонный рост по $n$. Отметим, что все отношения последовательных времен исполнения в обеих таблицах даже меньше, чем нижняя граница 4 для соответствующих теоретических отношений (отношения меньше 1 опущены); см. в этой связи замечание 4.1. Вопреки теоретическим ожиданиям почти все отношения для алгоритма (a) меньше, чем для (б). Отношения растут с возрастанием $K$ и $n$, и для алгоритма (б) их наибольшее значение уже близко к 4. При максимальных $K = 1024$ и $n = 9$, система (3.5) содержит почти $85 \times {{10}^{6}}$ неизвестных, но требуется только менее 2 мин для ее решения, что является отличным результатом.

В заключение рассмотрим наиболее интересный 3D случай ($N = 3$) и возьмем точное решение $u(x): = sin(2\pi {{x}_{1}})sin(3\pi {{x}_{2}})sin(4\pi {{x}_{3}})cosh(\sqrt 2 {{x}_{1}} - {{x}_{2}} + {{x}_{3}}{\text{/}}\sqrt 3 )$. Опять вычислим МКЭ-решение для различных значений $K = 2,4, \ldots ,64$ и $n = \overline {1,9} $ и изучим поведение погрешности, см. фиг. 7 и табл. 4 (где значения ${{R}_{C}}$ меньше 3 опущены). Здесь можно сделать в основном те же самые выводы, что и в 2D случае. Отметим, что более плохие свойства устойчивости алгоритма (б) проявляются только при высоких $n = 8,9$, поскольку берется намного меньшее максимальное значение $K$.

Времена выполнения в 3D случае представлены на фиг. 8 и в табл. 5 и 6. Опять естественны выводы, аналогичные сделанным в 2D случае. Все отношения в обеих таблицах заметно ниже, чем нижняя граница 8 для соответствующих теоретических отношений. Важно, что для максимальных $K = 64$ и $n = 9$ система (3.5) содержит более, чем $190 \times {{10}^{6}}$ неизвестных, но требуется только менее 15 мин для ее решения, что является хорошим результатом (результат менее быстрый, чем в 2D случае, из-за неэффективной Matlab-реализации вложенных циклов).