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

1. ВВЕДЕНИЕ

Задачи рассеяния волн различной физической природы на неоднородных структурах, находящихся в трехмерной ограниченной области $Q$, имеют большое значение как с теоретической, так и с практической точек зрения. Указанные задачи могут быть описаны объемным интегральным уравнением следующего общего вида:

где $R = \left| {x - y} \right|$ – расстояние между точками $x = ({{x}_{1}},{{x}_{2}},{{x}_{3}})$ и $y = ({{y}_{1}},{{y}_{2}},{{y}_{3}})$, $\alpha ,\;\eta ,\;K,\;{{u}_{0}}$ – известные функции, причем $K(x - y)$ – дифференцируемая функция координат, $u$ – неизвестная функция. Приведем некоторые важные классы задач, сводящихся к уравнению (1).

• Рассеяние акустических волн на прозрачном неоднородном препятствии (см. [1]). В этом случае $d = 1$, $\alpha = 0$, а остальные функции, входящие в уравнение (1), являются скалярными.

• Рассеяние электромагнитных волн на неоднородном, в общем случае анизотропном, теле (см. [2]). В этом случае $d = 3$ и поэтому оператор в (1) будет сингулярным, функции $u$ и ${{u}_{0}}$ – векторными, а $\eta $ и $K$ – тензорными. Величина $\alpha $ является внеинтегральным членом сингулярного оператора и зависит от формы выделяемой особенности и ее центра. Например, если выделяемая особеность есть шар, то $\alpha = 1{\text{/}}3$.

Для решения уравнения (1), которое описывает реальные задачи, возможно использование только численных методов. В методе Галеркина или в методе коллокации уравнение (1) аппроксимируется системой линейных алгебраических уравнений. Для трехмерных задач размер $N$ получаемых систем достаточно велик.

Основными критериями эффективности численного алгоритма являются число арифметических операций $T$ и объем памяти $M$, необходимой для его реализации. При использовании прямых методов $T \sim {{N}^{3}}$ и $M \sim {{N}^{2}}$. Для итерационных методов

где ${{T}_{A}}$ – число арифметических операций, которое требуется для умножения матрицы на вектор, $L$ – количество итераций, необходимое для получения решения с заданной точностью, а ${{M}_{A}}$ – число параметров, определяющих матрицу. Используя метод коллокации на равномерной сетке и свойства ядер интегральных операторов, уравнение (1) можно аппроксимировать системой линейных уравнений с многоуровневой блочно тёплицевой матрицей и строить эффективные алгоритмы матрично-векторного умножения, в которых значения ${{T}_{A}}$ и ${{M}_{A}}$ практически пропорциональны размеру $N$ (см. [3]).

Однако во многих случаях целесообразно использование неравномерной сетки. Такая потребность возникает при аппроксимации сложной границы области $Q$ или при значительном изменении параметров среды в области неоднородности. К сожалению, при применении метода коллокации на неравномерной сетке матрица получаемой системы утрачивает свойство тёплицевости и поэтому непосредственное использование вышеуказанной техники невозможно.

В настоящей работе для аппроксимации интегрального уравнения (1) применяется метод коллокации на неравномерной сетке, а затем используется некоторая равномерная сетка, методы интерполяции функций и алгоритмы быстрого преобразования Фурье. Подчеркнем, что число узлов равномерной сетки оказывается сопоставимым с числом узлов исходной неравномерной сетки. Поэтому в итоге строятся эффективные алгоритмы численного решения уравнения (1) итерационными методами, в которых значения ${{T}_{A}}$ и ${{M}_{A}}$ практически пропорциональны размеру $N$, хотя коэффициент пропорциональности несколько больше по сравнению с алгоритмами на равномерной сетке

2. ПОСТАНОВКА ЗАДАЧИ

Для аппроксимации интегрального уравнения (1) будем использовать метод коллокации на неравномерной сетке. Представим область $Q$ в виде объединения ${{N}_{Q}}$ ячеек ${{\Omega }_{n}}$, $n = 1,2, \ldots ,{{N}_{Q}}$. Узловые точки в этих ячейках будем выбирать в их центрах, которые определяются формулами

где ${{x}^{c}} = (x_{1}^{c},x_{2}^{c},x_{3}^{c})$ – центр ячейки $\Omega $, а $\left| \Omega \right|$ – ее объем. Если в качестве ячеек выбираются тетраэдры произвольной формы, то можно достаточно точно описать многие сложные конфигурации области $Q$ и среды. Центр тетраэдра определяется простой формулой где $x_{1}^{{(k)}},\;x_{2}^{{(k)}},\;x_{3}^{{(k)}}$ – декартовы координаты $k$-й вершины тетраэдра.

Будем аппроксимировать интегральное уравнение (1) системой линейных алгебраических уравнений размера $ \sim {\kern 1pt} {{N}_{Q}}$ относительно значений неизвестного поля в узловых точках области $Q$, находящихся в центрах ячеек ${{\Omega }_{n}}$:

Для векторных задач тензор ${{\alpha }_{n}}$ определяется формой ячейки ${{\Omega }_{n}}$ и ее центром (см. [4]). Отметим, что, поскольку узловые точки находятся в центре ячеек, точность аппроксимации интегральных операторов $ \sim {\kern 1pt} {{h}^{2}}$, где $h$ – максимальный диаметр ячеек.

Для относительно небольших значений ${{N}_{Q}} \sim 10\,000$ можно решить систему уравнений (5) прямыми или итерационными методами. Ниже мы предлагаем эффективные алгоритмы, которые с использованием итерационных методов позволяют решать систему (5) со значительно большим размером ${{N}_{Q}}$.

3. ПОСТРОЕНИЕ АЛГОРИТМА РЕШЕНИЯ

Сначала рассмотрим следующий алгоритм интерполяции функций. Пусть известны значения дифференцируемой функции $f({{x}_{1}},{{x}_{2}},{{x}_{3}})$ в вершинах параллелепипеда. Тогда приближенное значение функции в любой точке параллелепипеда может быть представлено формулой

В формуле (6) ${{a}_{1}},\;{{b}_{1}},\;{{c}_{1}}$ – декартовы координаты левой нижней вершины параллелепипеда, а ${{h}_{1}},\;{{h}_{2}},\;{{h}_{3}}$ – длины ребер по осям ${{x}_{1}},\;{{x}_{2}},\;{{x}_{3}}$. Формула (6) дает точное значение в любой точке параллелепипеда, если функция имеет вид Точность интерполяции по формуле (6) имеет второй порядок по $h = max({{h}_{1}},{{h}_{2}},{{h}_{3}})$. Запишем (6) в виде где ${{f}^{t}}$ – значение функции $f(x)$ в $t$-й вершине параллелепипеда в нумерации, соответствующей (6), а значения $\beta (t,x)$ также определяются формулой (6).

В прямоугольной декартовой системе координат введем сетку так, чтобы область $Q$ целиком находилась в прямоугольном параллелепипеде $\Pi $, стороны которого имеют длины ${{N}_{i}}{{h}_{i}}$, $i = 1,2,3$, где ${{h}_{i}}$ – шаги сетки по декартовым координатам. Параллелепипед $\Pi $ разбивается данной сеткой на элементарные параллелепипеды $\Pi (p)$, $p = ({{p}_{1}},{{p}_{2}},{{p}_{3}})$, ${{p}_{i}} = 0,1, \ldots ,{{N}_{i}} - 1$. Количество элементарных параллелепипедов в $\Pi $ равно $N = {{N}_{1}}{{N}_{2}}{{N}_{3}}$. Центр параллелепипеда $\Pi (p)$ в единицах шагов сетки обозначим через

Определим узловые точки в параллелепипеде $\Pi $ формулами Узловую точку сетки будем обозначать $\tilde {p} = ({{\tilde {p}}_{1}},{{\tilde {p}}_{2}},{{\tilde {p}}_{3}})$. Понятно, что вершины элементарного параллелепипеда $\Pi (p)$ задаются восемью узловыми точками. Общее число узловых точек в параллелепипеде $\Pi $ равно $({{N}_{1}} + 1)({{N}_{2}} + 1)({{N}_{3}} + 1)$. Оно несколько больше, чем число элементарных параллелепипедов, но при больших значениях ${{N}_{1}},\;{{N}_{2}},\;{{N}_{3}}$ асимптотически такое же. Отметим, что узловая точка может быть вершиной от одного (вершины исходного параллелепипеда $\Pi $) до восьми (внутренние узловые точки) элементарных параллелепипедов. Запись $\tilde {p} \in \Pi (p)$ означает, что узловая точка $\tilde {p}$ является вершиной элементарного параллелепипеда $\Pi (p)$. Условие принадлежности можно сформулировать в виде равенства

Предположим, что на параллелепипеде $\Pi $ определена функция вида $F(x - y)$, которая зависит только от разности декартовых координат и является дифференцируемой функцией точек $x$ и $y$ при $x \in \Pi (p)$, $y \in \Pi (q)$, $p \ne q$. Рассмотрим сумму вида

где $k \geqslant 1$, а величины ${{a}_{l}}$ постоянны. Применим формулы интерполяции (6),(7) к функции $F(x - y)$ в точках $x$ и ${{y}_{l}}$ в параллелепипедах $\Pi (p)$ и $\Pi (q)$ соответственно. Получим Окончательно имеем В (12), (13) коэффициенты $\beta (\tilde {q},q,y)$, $y \in \Pi (q)$, и $\beta (\tilde {p},p,x)$, $x \in \Pi (p)$, определяются согласно формулам (6), (7). Формулы (12), (13) имеют второй порядок точности по $h$. Они будут использоваться для вычисления выражений вида (11), когда входящие в них точки ${{y}_{l}}$ находятся внутри нескольких элементарных параллелепипедов. В этом случае суммирование во второй сумме (12) проводится по всем вершинам элементарных параллелепипедов, содержащих точки ${{y}_{l}}$. При этом веса интерполяции суммируются в (13) в общих вершинах элементарных параллелепипедов.

Запишем уравнение (5) в виде, удобном для построения алгоритма. Будем полагать, что каждая ячейка ${{\Omega }_{n}} \subset {\text{Q}}$ находится внутри только одного из элементарных параллелепипедов, причем внутри него может быть несколько ячеек. Данное требование практически не накладывает ограничений на выбор ячеек: если ячейка принадлежит нескольким элементарным параллелепипедам, то ее можно разбить на несколько ячеек. Обозначим через $k(p)$ число ячеек в элементарном параллелепипеде $\Pi (p)$, через $x(l,p)$ центр $l$-й ячейки в $\Pi (p)$, а через $\eta (l,p)$ и $\alpha (l,p)$ значения ${{\eta }_{n}}$ и ${{\alpha }_{n}}$ для соответствующей ячейки. Обозначим через ${{\Pi }_{Q}}$ объединение всех элементарных параллелепипедов, внутри которых находятся ячейки, а через $\Omega (l,p)$ обозначим $l$-ю ячейку в параллелепипеде $\Pi (p)$. Тогда систему линейных уравнений (5) можно переписать в следующем виде:

Основные вычислительные затраты в итерационных алгоритмах решения системы (14) связаны с умножением матрицы на вектор, т.е. с вычислением сумм вида

Обозначим через ${{\Pi }_{0}}(p)$ объединение элементарных параллелепипедов в ${{\Pi }_{Q}}$, содержащих $\Pi (p)$ и расположенных в некоторой окрестности $\Pi (p)$. Опишем два возможных варианта определения ${{\Pi }_{0}}(p)$:

• область $\Pi _{0}^{1}(p)$ содержит параллелепипед $\Pi (p)$ и все элементарные параллелепипеды, соприкасающиеся с ним (максимальное число таких параллелепипедов равно $27$);

• область $\Pi _{0}^{2}(p)$ содержит параллелепипед $\Pi _{0}^{1}(p)$ и все элементарные параллелепипеды, соприкасающиеся с ним (максимальное количество таких параллелепипедов равно $125$).

Обозначим ${{\Pi }_{Q}}(p) = {{\Pi }_{Q}}{{\backslash }}{{\Pi }_{0}}(p)$. Тогда (15) можно представить в виде

Значения ${{W}_{1}}(l,p)$ будем вычислять непосредственно по вышеприведенной формуле. Число арифметических операций и память для этих вычислений пропорциональны ${{N}_{Q}}$.

Перейдем к построению алгоритма для вычисления значений ${{W}_{2}}(l,p)$. Определим функцию

Для вычисления интегралов (14) по ячейке в выражении для ${{W}_{2}}$ будем использовать приближенную формулу Формула (18) дает точное значение интеграла на линейных функциях вида $f(x) = {{a}_{0}} + {{a}_{1}}{{x}_{1}} + {{a}_{2}}{{x}_{2}} + {{a}_{3}}{{x}_{3}}$ $f(x) = {{a}_{0}} + {{a}_{1}}{{x}_{1}} + {{a}_{2}}{{x}_{2}} + {{a}_{3}}{{x}_{3}}$ и имеет второй порядок по $h$, где $h$ – диаметр ячейки $\Omega $. Из (14), (16)–(18) находим Отметим, что формулы вида (18) для вычисления интегралов, входящих в ${{W}_{1}}$, и применение интерполяции могут приводить к значительным погрешностям. Это обстоятельство необходимо учитывать при определении областей ${{\Pi }_{0}}(p)$ и шагов сетки.

В узловых точках параллелепипеда $\Pi $ определим функцию дискретного аргумента

Из свойств симметрии функции следует, что весь массив значений $B(\tilde {p} - \tilde {q})$, $\tilde {p},\tilde {q} \in \Pi $, определяется его частью, содержащей $({{N}_{1}} + 1)({{N}_{2}} + 1)({{N}_{3}} + 1)$ элементов. Функцию $B(\tilde {p} - \tilde {q})$ и формулы (12), (13) будем использовать для интерполяции функций $B(x(l,p) - y(m,q))$, входящих в (19). Положим где $\beta (\tilde {p},p,m) = \beta (\tilde {p},p,x(m,p))$. В (21) веса интерполяции значений $V$ в центрах ячеек параллелепипеда $\Pi (p)$ суммируются в узловых точках $\Pi (p)$. Тогда, учитывая (10), ${{W}_{2}}$ можно приближенно представить в виде

Далее, в узловых точках параллелепипеда $\Pi $ определим функцию

Рассмотрим сумму Сравнивая $\tilde {W}(\tilde {p})$ и ${{\tilde {W}}_{2}}(\tilde {p},p)$, при $\Pi (p) \subset {{\Pi }_{Q}}$ находим Суммы $\tilde {W}(\tilde {p})$ можно записать в следующем виде: Для вычисления этих сумм будем использовать технику умножения тёплицевой матрицы на вектор (см. [3]). Кратко опишем ее применительно к (28). Обозначим через ${{\Pi }_{2}}$ параллелепипед со сторонами $2{{N}_{1}}{{h}_{1}}$, $2{{N}_{2}}{{h}_{2}}$, $2{{N}_{3}}{{h}_{3}}$. Продолжим матричную функцию дискретного аргумента $B({{\tilde {p}}_{1}},{{\tilde {p}}_{2}},{{\tilde {p}}_{3}})$ на все целочисленные значения ${{\tilde {p}}_{1}},\;{{\tilde {p}}_{2}},\;{{\tilde {p}}_{3}}$, полагая ее периодической по каждой переменной с периодами соответственно $2({{N}_{1}} + 1)$, $2({{N}_{2}} + 1)$, $2({{N}_{3}} + 1)$. При этом доопределим функцию $B({{\tilde {p}}_{1}},{{\tilde {p}}_{2}},{{\tilde {p}}_{3}})$ нулем в тех точках, где она не определена. Доопределим вектор-функцию дискретного аргумента $V({{\tilde {p}}_{1}},{{\tilde {p}}_{2}},{{\tilde {p}}_{3}})$ нулем во всех узловых точках ${{\Pi }_{2}}$, не принадлежащих $\Pi $. Ясно, что при $\tilde {p} \in \Pi $ функция совпадает с $W({{\tilde {p}}_{1}},{{\tilde {p}}_{2}},{{\tilde {p}}_{3}})$ из (28). Выполнив дискретное преобразование Фурье по каждой переменной от обеих частей (29), получим равенство Таким образом, функцию $W(\tilde {p},\tilde {p}) \in \Pi $ можно вычислить, используя прямое и обратное быстрые дискретные преобразования Фурье.

Отметим, что точности интерполяции функций и аппроксимации интегралов имеют один и тот же порядок точности по $h$. Заметим, что точности интерполяции функции $B(x - y)$ и вычисления интегралов по формулам (3), (18) улучшаются с увеличением расстояния между точками $x$ и $y$. Это обстоятельство необходимо учитывать при выборе областей ${{\Pi }_{0}}(p)$.

Полученный результат можно сформулировать в матричном виде следующим образом.

Теорема. В случае неравномерной сетки матрица $A$ порядка ${{N}_{Q}}$ рассматриваемой системы линейных алгебраических уравнений допускает приближение

где $T$ – многоуровневая блочно-тёплицева матрица порядка $N = O({{N}_{Q}})$, соответствующая равномерной сетке, а ${{S}_{0}},\;{{S}_{1}},\;{{S}_{2}}$ – некоторые разреженные матрицы. Число ненулевых элементов в матрицах ${{S}_{1}}$ и ${{S}_{2}}$ имеет вид $O({{N}_{Q}})$. Число арифметических операций при умножении матрицы $\tilde {A}$ на вектор имеет вид $O({{N}_{Q}}log{{N}_{Q}})$.

4. ЗАКЛЮЧЕНИЕ

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