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

1. ВВЕДЕНИЕ

Граничные методы решения задач математической физики давно и успешно используются для решения, в основном, эллиптических задач. Главное преимущество граничной постановки – снижение размерности исходной задачи на единицу путем перехода к уравнениям для неизвестных величин, заданных лишь на границе области, в которой требуется найти решение [1]. Непосредственное обобщение такого подхода на нестационарные задачи наталкивается на очевидную трудность: граница для эволюционной задачи представляет собой пространственно-временнóе многообразие – прямое произведение пространственной границы области и отрезка времени, на котором требуется найти решение, – причем с ростом времени интегрирования эта граница разрастается, заставляя хранить и решать все большее количество переменных и уравнений [2]. На практике подобный подход оказывается неприменим уже при малых временах интегрирования.

В данной работе представлен граничный метод решения некоторых гиперболических уравнений, свободный от указанного недостатка. Мы рассматриваем уравнения, для которых выполняется принцип Гюйгенса, т.е., иными словами, решения которых обладают конечной памятью “назад” по времени. Последнее эквивалентно наличию лакун в решении данных уравнений, т.е. областей пространства-времени, где решения обращаются в ноль по окончании работы источника (правой части уравнения), порождающего данное решение. Практически важными примерами подобных уравнений являются трехмерное волновое уравнение с постоянными коэффициентами и система Максвелла в вакууме. Ограниченная во времени зависимость решений этих уравнений позволяет разбить интервал интегрирования на небольшие части и продвигаться по времени шагами фиксированной длительности, последовательно обновляя решение некоторых граничных уравнений. Тем самым удается совместить преимущества граничной формулировки (снижение размерности по пространственным переменным на единицу) с алгоритмом обновления решения в духе обычной разностной схемы для эволюционного уравнения. За счет снижения размерности построенный алгоритм оказывается более эффективным с точки зрения вычислительной сложности, чем традиционное решение эволюционного уравнения с помощью явной разностной схемы (это свойство выполняется начиная с определенной длительности расчета).

Другой особенностью предлагаемого подхода является способ сведения и особенно решения граничной задачи. Традиционный подход базируется на дискретизации интегралов, взятых по границе, и решении возникающей алгебраической системы. Для областей со сложной геометрией указанный прием может вызвать значительные вычислительные сложности. Дополнительная трудность связана с сингулярностью ядер в фигурирующих свертках. В нашем изложении мы опираемся на метод разностных потенциалов [3], позволяющий вычислять возникающие свертки путем решения вспомогательной задачи, поставленной в геометрически простой области (скажем, кубе), с помощью любой подходящей разностной схемы на простейшей прямоугольной сетке. Последнее обстоятельство обеспечивает легкость реализации алгоритма, возможность контролировать и увеличивать точность расчета путем использования схем повышенной точности и другие преимущества, присущие конечно-разностному подходу. При этом никакой потери точности, связанной с аппроксимацией криволинейных границ/граничных условий на прямоугольной сетке, не происходит.

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

В разд. 2 представлена общая идея метода в непрерывной формулировке. Основным моментом здесь является схема продвижения во времени небольшими шагами. Разд. 3 содержит детали дискретизации алгоритма. В разд. 4 приводятся результаты тестового расчета.

2. ИДЕЯ МЕТОДА

Рассмотрим внешнюю задачу для трехмерного волнового уравнения с постоянными коэффициентами

где ${{\square }_{c}}\; \equiv \tfrac{1}{{{{c}^{2}}}}\tfrac{{{{\partial }^{2}}}}{{\partial {{t}^{2}}}} - \Delta $ есть оператор Даламбера, а $c$ – скорость распространения взаимодействия. Решение $u({\mathbf{x}},t)$ требуется найти снаружи ограниченной области $\Omega \subset {{\mathbb{R}}^{3}}$ в течение времени $t \in (0,T]$. Подобная постановка задачи характерна для (нестационарного) процесса рассеяния волн на препятствии $\Omega $, причем функция $u({\mathbf{x}},t)$ играет роль отраженного поля. В такой трактовке нулевые начальные условия (2.3) отвечают ситуации, когда падающее излучение только достигло (в момент $t = 0$) препятствия $\Omega $. Оператор ${{{\mathbf{l}}}_{\Gamma }}$ в (2.2) задает граничные условия на поверхности рассеивателя и определяет тип рассеяния. Простейшими примерами граничных условий на $\partial \Omega $ могут служить условия Дирихле ${{{\mathbf{l}}}_{\Gamma }}u \equiv u = - {{u}^{{{\text{inc}}}}}$ или Неймана ${{{\mathbf{l}}}_{\Gamma }}u \equiv \tfrac{{\partial u}}{{\partial {\mathbf{n}}}} = - \tfrac{{\partial {{u}^{{{\text{inc}}}}}}}{{\partial {\mathbf{n}}}}$ (${\mathbf{n}}$ – внешняя нормаль к границе $\Gamma $ области $\Omega $), записанные в терминах падающего излучения ${{u}^{{{\text{inc}}}}}$. В общем случае граничные условия могут быть гораздо более сложными, в т.ч. нелокальными, обеспечивающими, однако, корректность поставленной задачи (2.1)–(2.3).

Рассмотрим несколько бóльшую вспомогательную область $\Omega {\text{'}} \supset \Omega $ такую, что разность (черта над $\Omega $ означает замыкание) содержит в себе, в частности, расчетную область, где требуется найти решение поставленной задачи. Без потери общности будем предполагать, что $\Omega {\text{'}}$ представляет собой параллелепипед. Выберем также некоторую запаздывающую функцию Грина (фундаментальное решение) $G({\mathbf{x}},t)$ уравнения (2.1). Применим формулу Грина к четырехмерной области , $t \leqslant T$ (см. [5]):

выражающую решение $u\left( {{\mathbf{x}},t} \right)$ внутри данной области через значения решения на границе. Данная (трехмерная) граница состоит из кусков, отвечающих начальному моменту времени, , и ненулевым временам ${{\Gamma }_{t}} \equiv \partial \Omega \times (0,t]$ и $\Gamma _{t}^{'} \equiv \partial \Omega {\text{'}} \times (0,t]$, причем $\Gamma _{t}^{'}$ отвечает пространственно-временнóй границе вспомогательной области $\Omega {\text{'}}$.

Первый член в (2.4) равняется нулю из-за однородных граничных условий (2.3). Последний интеграл в формуле (2.4) берется по пространственно-временнóй границе вспомогательной области $\Gamma _{t}^{'} \equiv \partial \Omega {\text{'}} \times (0,t]$ и равняется нулю при условии, что выбранная функция Грина есть фундаментальное решение трехмерного волнового уравнения, распространяющееся в пустом пространстве, т.е.

В этом случае граничные значения решения и его производной на $\Gamma _{t}^{'}$ не дают вклада в решение $u({\mathbf{x}},t)$ во внутренних точках $\Omega {\text{'}}$, поскольку все, что достигло пространственной границы $\partial \Omega {\text{'}}$ есть уходящие вовне $\Omega {\text{'}}$ волны. Таким образом, при указанном выборе функции Грина в выражении (2.4) остается лишь второе слагаемое, которое позволяет выразить значение решения в расчетной области через его значения (и значения его нормальной производной) на границе рассеивателя $\Omega $ во все предыдущие моменты времени.

По аналогии с формулой (2.4) вводится потенциал Кальдерона с плотностью ${{\xi }_{{{{\Gamma }_{t}}}}} = ({{\xi }_{0}},{{\xi }_{1}})$, определенной на ${{\Gamma }_{t}}$:

где ${\mathbf{x}} \in \tilde {\Omega } = {{\mathbb{R}}^{3}}{\backslash }\bar {\Omega }$. Отметим, что потенциал Кальдерона (2.6) удовлетворяет волновому уравнению на $\tilde {\Omega }$, т.е. ${{\square }_{c}}{{{\mathbf{P}}}_{{\tilde {\Omega }}}}{{\xi }_{{{{\Gamma }_{t}}}}}({\mathbf{x}},t) = 0$, для произвольной плотности ${{{\mathbf{\xi }}}_{{{{\Gamma }_{t}}}}}$.

Сравнивая формулы (2.4) и (2.6), нетрудно заметить, что функции ${{\xi }_{0}}$ и ${{\xi }_{1}}$ можно интерпретировать как след на $\partial \Omega $ уходящего из $\Omega {\text{'}}$ решения и его нормальной производной. Доказывается, см. [3], что пара ${{\xi }_{{{{\Gamma }_{t}}}}} = ({{\xi }_{0}},{{\xi }_{1}})$ действительно является следом уходящего решения тогда и только тогда, когда ${{\xi }_{{{{\Gamma }_{t}}}}}$ удовлетворяет граничному уравнению с проектором (ГУП):

где ${{{\mathbf{P}}}_{{{{\Gamma }_{t}}}}} \equiv {\mathbf{T}}{{{\mathbf{r}}}_{{{{\Gamma }_{t}}}}}{{{\mathbf{P}}}_{{\tilde {\Omega }}}}{{\xi }_{{{{\Gamma }_{t}}}}}$, а оператор взятия следа определяется как ${\mathbf{T}}{{{\mathbf{r}}}_{{{{\Gamma }_{t}}}}}f({\mathbf{x}},t) = \mathop {\left. {\left( {f,\tfrac{{\partial f}}{{\partial n}}} \right)} \right|}\nolimits_{{{\Gamma }_{t}}} $. Очевидно, что существует (бесконечное) множество решений ГУП (2.7), поскольку существует множество распространяющихся вовне решений исходного уравнения (2.1). Чтобы выделить то из них, которое отвечает типу рассеяния, заданному граничным условием (2.2), нужно решить ГУП (2.7) совместно с граничным условием, связывающим решение и его нормальную производную на границе рассеивателя: Подчеркнем, что граничное условие (2.2) может определять одну из компонент плотности, ${{\xi }_{0}}$ (задача Дирихле с заданной на границе функцией) или ${{\xi }_{1}}$ (задача Неймана с заданной нормальной производной), или связь между ними (задача третьего рода и более сложные нелокальные условия). Однако никогда обе компоненты плотности (и функция, и нормальная производная на границе) не фиксируются одновременно, поскольку это привело бы к переопределенной задаче, не имеющей однозначного решения.

Таким образом, ГУП (2.7) совместно с граничным условием (2.8) являются эквивалентной граничной формулировкой исходной системы (2.2) с искомой величиной ${{\xi }_{{{{\Gamma }_{t}}}}}$ (возможно, только одной неизвестной компонентой в случае условий Дирихле и Неймана). Размерность задачи при этом снижается на единицу по сравнению с исходной системой. По найденной плотности ${{\xi }_{{{{\Gamma }_{t}}}}}$ решение $u({\mathbf{x}},t)$ в трехмерной области $\tilde {\Omega } = {{\mathbb{R}}^{3}}{\backslash }\Omega $ восстанавливается по формуле (2.6).

Подчеркнем, что интегрирование в формуле (2.6) проводится по (2 + 1)-мерному многообразию и в случае геометрически сложного рассеивателя $\Omega $ может представлять вычислительно нетривиальную задачу. В этой связи удобно использовать другое, эквивалентное представление потенциала Кальдерона. Рассмотрим функцию $w({\mathbf{x}},t)$, определенную на ${{\mathbb{R}}^{3}}{\backslash }\Omega \times (0,t]$ и такую, что ${\mathbf{T}}{{{\mathbf{r}}}_{{{{\Gamma }_{t}}}}}w = {{\xi }_{{{{\Gamma }_{t}}}}}$, $w({\mathbf{x}},0) = \tfrac{{\partial w({\mathbf{x}},0)}}{{\partial t}} = 0$. Иными словами, от функции $w({\mathbf{x}},t)$ требуется лишь то, чтобы ее след и след ее нормальной производной на ${{\Gamma }_{t}}$ совпадал с ${{\xi }_{{{{\Gamma }_{t}}}}}$, а также выполнялись указанные начальные условия. В остальном эта функция произвольна, в частности, она не обязана быть решением волнового уравнения (2.1), т.е. в общем случае ${{\square }_{c}}w \ne 0$ на ${{\mathbb{R}}^{3}}{\backslash }\Omega \times (0,t]$. На практике функцию $w({\mathbf{x}},t)$ можно построить продолжением с границы ${{\Gamma }_{t}}$ с помощью ряда Тейлора, где необходимые производные второго и более высоких порядков вычисляются с помощью исходного уравнения (2.1), см. разд. 2. Применяя формулу Грина к данной функции на и разрешая получающееся соотношение относительно интересующей нас величины (2.6), нетрудно убедиться, что потенциал Кальдерона можно представить в следующем виде:

где функция $v({\mathbf{x}},t) = G * \,\square \mathop {\left. {_{c}w} \right|}\nolimits_{\tilde {\Omega }} $ определяется через свертку функции Грина (2.5) c величиной $\square \mathop {\left. {_{c}w} \right|}\nolimits_{\tilde {\Omega }} $, доопределенной нулем внутри рассеивателя, т.е. на $\Omega $. Отметим, что эта свертка берется по всему трехмерному пространству ${{\mathbb{R}}^{3}}$ и по временнóму интервалу $(0,t]$. Эту свертку удобно вычислять как решение вспомогательной задачи (ВЗ): поставленной формально в неограниченной области ${{\mathbb{R}}^{3}}$ с однородными начальными условиями. На практике эта задача решается, конечно, в ограниченной области $\Omega {\text{'}}$, но на границе $\Omega {\text{'}}$ ставятся условия излучения. Последнее эквивалентно выбору требуемой функции Грина (2.5). Переход от расчета потенциала Кальдерона через интеграл (2.6) к расчету в виде решения ВЗ (2.10) играет ключевую роль в удобстве, гибкости и простоте реализации излагаемого подхода.

Подчеркнем еще раз, что ВЗ (2.10) решается в простой (прямоугольной) вспомогательной области $\Omega {\text{'}}$ без “дырки”, отвечающей рассеивателю $\Omega $; вся информация о геометрии этого объекта (возможно, сложной и не совпадающей с линиями декартовой сетки), а также граничных условиях на нем содержится в функции $w(x,t)$ и переходит в правую часть (2.10) и далее в формулу (2.9). Данное обстоятельство позволяет использовать конечно-разностный метод на простейших прямоугольных сетках для расчета объектов сложной геометрической формы без какой-либо потери точности, связанной с несовпадением границ объекта с линиями сетки.

Представим теперь, что нам требуется найти решение задачи указанным методом в конечный момент времени $t = T$. С ростом $T$ граница ${{\Gamma }_{T}} = \partial \Omega \times (0,T]$, где определена плотность ${{\xi }_{{{{\Gamma }_{T}}}}}$, увеличивается линейно, и уже при умеренных $T$ метод становится практически неприменимым из-за необходимости хранения и решения уравнения относительно разрастающейся со временем величины ${{\xi }_{{{{\Gamma }_{T}}}}}$.

Для того, чтобы снять это ограничение, разобъем интервал $(0,T]$ на $K$ частей длины ${{T}_{0}}$, $K{{T}_{0}} = T$, и на каждой парциальной пространственно-временной границе

отвечающей указанному разбиению, введем соответствующие плотности ${{\xi }_{{{{\Gamma }_{k}}}}}$, см. фиг. 2а. На первый взгляд, подобное формальное разбиение не решает проблемы разрастания границы ${{\Gamma }_{T}}$ со временем, поскольку в силу причинности плотность ${{\xi }_{{{{\Gamma }_{k}}}}}$ на любом $k$-м элементе разбиения должна зависеть от плотностей на всех предыдущих участках ${{\xi }_{{{{\Gamma }_{s}}}}}$, $s = 1,...,k - 1$. Казалось бы, это не позволяет решать систему (2.7), (2.8) для ${{\xi }_{{{{\Gamma }_{k}}}}}$ изолированно, не задействуя все предыдущие величины ${{\xi }_{{{{\Gamma }_{s}}}}}$.

Однако легко заметить, что для выбранной функции Грина (2.5) интегрирование в потенциале Кальдерона (2.6) осуществляется по пересечению светового конуса с вершиной в точке границы ${{\Gamma }_{k}}$ и поверхности рассеивателя $\Omega $ в предыдущие моменты времени, см. фиг. 2а, где для демонстрации положено $k = K$. Рассматривая все точки из ${{\Gamma }_{K}}$, где мы интересуемся плотностью ${{\xi }_{{{{\Gamma }_{K}}}}}$, данное пересечение не пусто максимум в течение времени $2{{T}_{0}}$, отсчитывая назад от последнего (наиболее позднего) момента времени множества ${{\Gamma }_{K}}$ (т.е. для случая, изображенного на рисунке, $t = T$), при условии, что размер разбиения принят равным

Иными словами, при условии (2.11) плотности ${{\xi }_{{{{\Gamma }_{s}}}}}$, отвечающие более ранним элементам разбиения ${{\Gamma }_{s}}$, $s = 1,...,K - 2$, никак не влияют на плотность, определенную на элементе ${{\Gamma }_{K}}$. То же самое можно объяснить на языке лакун решения волнового уравнения в свободном пространстве, когда решение, порожденное источником (правой частью), локализованным в пространстве-времени на границе предыдущего элемента разбиения (${{\Gamma }_{1}}$ на фиг. 2б), заведомо уходит из области $\Omega $ за время (2.11). Таким образом, глубина зависимости решения уравнения (2.7) ограничена назад одним элементом разбиения, что позволяет построить рекуррентную схему обновления плотности: где и $({\mathbf{x}},t) \in {{\Gamma }_{k}}$ в (2.13). Отметим, что интеграл (2.13) берется по предыдущему элементу разбиения ${{\Gamma }_{{k - 1}}}$, но его значение рассчитывается в точке из следующего элемента ${{\Gamma }_{k}}$.

Парциальные плотности ${{\xi }_{{{{\Gamma }_{k}}}}}$ и ${{\xi }_{{{{\Gamma }_{{k - 1}}}}}}$ удобно разложить в ряд на ${{\Gamma }_{k}}$ и ${{\Gamma }_{{k - 1}}}$ по системе каких-либо базисных функций

где двухкомпонентные базисные функции ${{\psi }_{{0,s}}} = ({{\psi }_{s}},0)$ и ${{\psi }_{{1,s}}} = (0,{{\psi }_{s}})$ используются для представления решения и его нормальной производной на границе. Естественно также использовать одни и те же базисные функции для любого элемента разбиения $k$. Например, в случае сферического рассеивателя $\Omega $ в качестве ${{\psi }_{s}}({\mathbf{x}},t)$ можно взять прямое произведение сферических гармоник (по угловым переменным) и полиномов Чебышева (по времени). В зависимости от наложенного граничного условия (2.2) часть коэффициентов в указанных разложениях может быть известна; например, в случае условия Дирихле заданная на границе функция однозначно определяет величины $c_{{0,s}}^{{({\text{I}})}}$, $c_{{0,s}}^{{({\text{II}})}}$; условие Неймана фиксирует $c_{{1,s}}^{{({\text{I}})}}$, $c_{{1,s}}^{{({\text{II}})}}$, а условие III рода устанавливает связь между этими системами коэффициентов.

Укажем также, что ядро интегралов (2.6) и (2.13) с учетом явного вида (2.5) инвариантно относительно сдвигов во времени. Поэтому с учетом использования одного и того же базиса для всех элементов разбиения результат действия операторов (2.6) и (2.13) на базисные элементы на множестве ${{\Gamma }_{k}}$ (левая часть уравнения (2.12)) можно рассчитать на любых двух смежных элементах разбиения длительностью ${{T}_{0}}$, например, ${{\Gamma }_{1}}$ и ${{\Gamma }_{2}}$, и обозначить результат через ${{{\mathbf{P}}}_{{{{\Gamma }_{{{{T}_{0}}}}}}}}{{\psi }_{s}}$ и ${{{\mathbf{R}}}_{{{{\Gamma }_{{{{T}_{0}}}}}}}}{{\psi }_{s}}$. При этом уравнение (2.7) примет вид

Подчеркнем, что действие этих операторов вычисляется через вспомогательную задачу (2.10), причем ВЗ для вычисления ${{{\mathbf{R}}}_{{{{\Gamma }_{{{{T}_{0}}}}}}}}{{\psi }_{s}}$ решается в течение времени $t \in [0,2{{T}_{0}}]$, а для ${{{\mathbf{P}}}_{{{{\Gamma }_{{{{T}_{0}}}}}}}}{{\psi }_{s}}$ – в течение $t \in [{{T}_{0}},2{{T}_{0}}]$. Число таких ВЗ равно числу базисных элементов в разложении (2.14). После указанного расчета уравнение (2.15) превращается в уравнение для коэффициентов разложения плотности на двух смежных элементах разбиения ${{\Gamma }_{k}}$ и ${{\Gamma }_{{k - 1}}}$, где через I обозначается тождественный оператор, причем $c_{{0,s}}^{{({\text{I}})}}$ и $c_{{1,s}}^{{({\text{I}})}}$ известны с предыдущего шага. Для выбора нужного решения к уравнению (2.16) следует добавить граничное условие (2.8), записанное в терминах коэффициентов разложения:

3. ДИСКРЕТИЗАЦИЯ

Введем на вспомогательной области $\Omega {\text{'}}$ (параллелепипиде) прямоугольную сетку с шагом $h$ по каждому пространственному измерению. Еще раз подчеркнем, что хотя разностное решение задачи определено и требуется лишь вне препятствия $\Omega $, мы тем не менее рассматриваем дискретизацию всей области $\Omega $. Заметим также, что граница объекта $\partial \Omega $ может не совпадать с линиями введенной прямоугольной сетки. Решением задачи будем интересоваться в дискретные моменты времени с шагом $\tau $, т.е., иными словами, введем сетку на четырехмерной области $\Omega {\text{'}} \times [0,T]$ с указанными шагами. Для дискретизации уравнения (2.2) воспользуемся компактной разностной схемой четвертого порядка точности с шаблоном, состоящим из $3 \times 3 \times 3 \times 3 = 81$ точек:

с параметром $\theta = \tfrac{1}{{12}}$. Подчеркнем, что уравнение (3.1) описывает лишь дискретизацию по времени; отметим, что эта схема является неявной. На временнóм слое $n + 1$ уравнение (3.1) далее аппроксимируется с четвертым порядком по пространственным переменным на шаблоне из $3 \times 3 \times 3$ точек, используя общую идею компактных схем, см. [6].

Обозначим через ${{\mathbb{N}}_{m}}$ точки сетки, принадлежащие данному 81-точечному шаблону с центром в узле сетки $m$. Введем, далее, несколько требуемых в дальнейшем подмножеств пространственно-временной прямоугольной сетки, покрывающей $\Omega {\text{'}} \times [0,T]$. Именно, обозначим через ${{\mathbb{M}}^{ - }}$ и ${{\mathbb{M}}^{ + }}$ подмножества узлов сетки, попадающих в области и $\Omega \times [0,T]$ соответственно. Таким образом, ${{\mathbb{M}}^{ - }}$ и ${{\mathbb{M}}^{ + }}$ представляют собой узлы сетки, лежащие вне и внутри объекта $\Omega $ во все моменты времени $t \in [0,T]$; узлы, попадающие на границу препятствия $\partial \Omega $, включаются в ${{\mathbb{M}}^{ + }}$. При этом центр шаблона ${{\mathbb{N}}_{m}}$ находится внутри областей и $\Omega \times [0,T]$, если $m \in {{\mathbb{M}}^{ - }}$ и $m \in {{\mathbb{M}}^{ + }}$ соответственно.

Пусть ${{\gamma }^{ + }}$ обозначает подмножество узлов, удовлетворяющих одновременно двум условиям: эти узлы принадлежат шаблону ${{\mathbb{N}}_{m}}$ с центром $m \in {{\mathbb{M}}^{ - }}$ и лежат внутри $\Omega $. Аналогично, подмножество ${{\gamma }^{ - }}$ обозначает узлы, принадлежащие шаблону ${{\mathbb{N}}_{m}}$ с центром $m \in {{\mathbb{M}}^{ + }}$ и лежащие вне $\Omega $. Объединение ${{\gamma }^{ + }} \cup {{\gamma }^{ - }} \equiv \gamma $ образует тонкий слой узлов, охватывающий пространственно-временную границу $\Gamma $ рассеивателя $\Omega $ и заменяющий $\Gamma $ в дискретной постановке, см. фиг. 3. Рассмотрим также множества ${{\mathbb{N}}^{ + }} = \bigcup\nolimits_{m \in {{\mathbb{M}}^{ + }}}^{} {{{\mathbb{N}}_{m}}} $ и ${{\mathbb{N}}^{ - }} = \bigcup\nolimits_{m \in {{\mathbb{M}}^{ - }}}^{} {{{\mathbb{N}}_{m}}} $, включающие все точки шаблона ${{\mathbb{N}}_{m}}$, когда центр последнего пробегает все точки из ${{\mathbb{M}}^{ + }}$ или ${{\mathbb{M}}^{ - }}$ соответственно. В частности, ${{\mathbb{N}}^{ - }}$ содержит точки на границе $\partial \Omega {\text{'}}$, первый и последний слои по времени $t = \{ 0,T\} $ внутри области , а также подмножество ${{\gamma }^{ + }}$. Объединение ${{\mathbb{N}}^{0}} = {{\mathbb{N}}^{ + }} \cup {{\mathbb{N}}^{ - }}$ есть набор всех точек сетки, которые попадают в шаблон ${{\mathbb{N}}_{m}}$, когда $m \in {{\mathbb{M}}^{ + }} \cup {{\mathbb{M}}^{ - }}$.

Все введенные выше (под)множества узлов отвечали сетке, покрывающей всю пространственно-временную область $\Omega {\text{'}} \times [0,T]$, отвечающую конечному моменту времени $T$. Очевидно, что все эти конструкции можно определить и на любой части $\Omega {\text{'}} \times [0,t]$, $t < T$, снабдив их подстрочным индексом $t$. В частности, множество ${{\gamma }_{t}}$ является дискретным аналогом границы ${{\Gamma }_{t}}$ из разд. 2.

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

Для демонстрации работы описанного метода рассмотрим задачу рассеяния плоской волны на сфере радиуса ${{R}_{0}} = 1$ см. Эта классическая задача имеет аналитическое решение в виде ряда, см. [5], с которым мы и будем сравнивать численное решение. В качестве падающего излучения возьмем не одну, а сумму двух плоских волн с несоизмеримыми угловыми частотами, ${{\omega }_{1}} = 3$ c^–1 и ${{\omega }_{2}} = {{\omega }_{1}}{\text{/}}\sqrt 2 $, получая тем самым апериодическое во времени решение. Отметим, что для избранных частот длина волны падающего излучения является величиной порядка диаметра рассеивателя $2{{R}_{0}}$. Решением будем интересоваться в сферическом слое вокруг рассеивателя ${{R}_{0}} < r < {{R}_{1}}$, где ${{R}_{1}} = 1.5$ см.

В качестве вспомогательной области $\Omega {\text{'}}$ возьмем куб ${{[ - {{R}_{2}},{{R}_{2}}]}^{3}}$, ${{R}_{2}} > {{R}_{1}}$, такой, что за максимальное время расчета вспомогательной задачи $2{{T}_{0}}$ волны, отраженные от внешней границы $\partial \Omega {\text{'}}$, не успеют достичь расчетной области ${{R}_{0}} < r < {{R}_{1}}$, см. фиг. 2б. Нетрудно проверить, что для этого величина ${{R}_{2}}$ должна удовлетворять условию $2{{R}_{2}} - {{R}_{0}} - {{R}_{1}} > 2c{{T}_{0}}$. Такой выбор размера вспомогательной области позволяет поставить произвольные граничные условия на $\partial \Omega {\text{'}}$ вместо каких-либо приближенных условий излучения и получить решение ВЗ (3.3) в интервале времени $[0,2{{T}_{0}}]$ на ${{R}_{0}} < r < {{R}_{1}}$ как если бы ВЗ решалась в неограниченной области. Отметим, что подобный прием позволяет точно, без какого-либо приближения, смоделировать процесс ухода отраженных волн в бесконечность и является существенным преимуществом излагаемого метода. Возможность этого обусловлена коротким временем расчета ВЗ (максимум $2{{T}_{0}}$ секунд), необходимым для формирования матриц (3.11), (3.12).

Построенная таким образом область $\Omega {\text{'}}$ покрывается прямоугольной сеткой, на которой алгоритм реализуется в рамках упомянутой разностной схемы четвертого порядка [6]. Уменьшая шаг сетки, можно проследить за сходимостью алгоритма, см. фиг. 4. Из данного рисунка видно, что двоичный логарифм ошибки уменьшается на четыре при измельчении сетки вдвое, что подтверждает правильную работу алгоритма. Заметим также, что на данном рисунке представлен результат достаточно длительного расчета: конечный момент ${{T}_{{{\text{sim}}}}}$ отвечает времени, за которое сигнал пересечет шар $r \leqslant {{R}_{0}}$ около 6000 раз. При этом никаких признаков возрастания ошибки со временем не наблюдается.

Отметим также, что для построения базиса ${{\psi }_{{0,s}}}$, ${{\psi }_{{1,s}}}$ используется прямое произведение полиномов Чебышева ${{T}_{n}}(t)$ (по времени) и сферических гармоник ${{Y}_{{lm}}}(\vartheta ,\varphi )$ (для разложения по границе $\Omega $, т.е. по сфере). Как показал эксперимент, максимального порядка $n = 10$ и $l = 12$ хватает с большим запасом для результатов с уровнем точности, представленным на фиг. 4.

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

В работе предложен алгоритм расчета нестационарных процессов распространения волн, позволяющий: (i) понизить размерность задачи на единицу за счет сведения к граничной формулировке (ii) производить расчеты для областей сложной геометрической формы, пользуясь обычными разностными схемами на простейших прямоугольных сетках без какой-либо потери точности, связанной с аппроксимацией границ (iii) точно учитывать условия излучения для внешних задач (iv) эффективно использовать возможности параллелизации вычислений.

Обсудим следствия (i) и (iv) более подробно. Очевидно, что обновление решения на элементе разбиения длины ${{T}_{0}}$ согласно (3.10) занимает меньше времени, чем продвижение решения на это же время с помощью какой-либо явной разностной схемы, в силу того, что размерность задачи в первом случае на единицу меньше. Однако чтобы начать решать рекуррентное соотношение (3.10), требуется рассчитать матрицы, фигурирующие в (3.10). Это наиболее вычислительно длительная часть алгоритма, т.к. на этом этапе приходится решать $4S$ задач на протяжении времени максимум $2{{T}_{0}}$. Из-за этого обстоятельства при небольших временах расчета обычная явная схема справится с заданием быстрее, чем предложенный алгоритм. При более длительных временах затраты времени на вычисление матриц ${\mathbf{Q}}_{{{{\gamma }_{{{{T}_{0}}}}}}}^{{(0)}}$, ${\mathbf{Q}}_{{{{\gamma }_{{{{T}_{0}}}}}}}^{{(1)}}$, ${\mathbf{R}}_{{{{\gamma }_{{{{T}_{0}}}}}}}^{{(0)}}$, ${\mathbf{R}}_{{{{\gamma }_{{{{T}_{0}}}}}}}^{{(1)}}$ будут скомпенсированы более быстрым, чем в случае явной схемы, продвижением согласно (3.10). Таким образом, существует время расчета, начиная с которого предложенный алгоритм становится выгоднее самого очевидного, простого и быстрого решения – продвижения с помощью явной разностной схемы. Также важно отметить, что матрицы ${\mathbf{Q}}_{{{{\gamma }_{{{{T}_{0}}}}}}}^{{(0)}}$, ${\mathbf{Q}}_{{{{\gamma }_{{{{T}_{0}}}}}}}^{{(1)}}$, ${\mathbf{R}}_{{{{\gamma }_{{{{T}_{0}}}}}}}^{{(0)}}$, ${\mathbf{R}}_{{{{\gamma }_{{{{T}_{0}}}}}}}^{{(1)}}$ рассчитываются параллельно, т.к. их столбцы независимы друг от друга. Чем больше ядер, узлов и т.д. задействовано в этом расчете, тем меньшее время потребуется на данный расчет, и тем скорее (при меньших временах расчета) описанный алгоритм окажется эффективнее явной схемы. Укажем также, что подробное исследование описанного метода применительно к внутренним задачам представлено в работе [7].

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