Журнал вычислительной математики и математической физики, 2022, T. 62, № 5, стр. 809-822

ВВЕДЕНИЕ

Многие приложения в физике связаны с решением трехмерного уравнения Гельмгольца в неограниченной области. Применение стандартных технологий решения этой проблемы таких, как методы конечных элементов или конечных объемов, сопряжены со сложностями, связанными с неограниченностью расчeтной области и необходимостью удовлетворять условию излучения на бесконечности. Существует большое количество численных методов решения внешних краевых задач в неограниченной области [1]. Одним из наиболее известных является представление искомого решения в виде комбинации потенциалов простого и двойного слоя, в которых плотности потенциалов находятся решением интегральных уравнений Фредгольма на границе области [2]–[4]. Другой подход для решения внешних краевых задач состоит в комбинированном применении метода конечных элементов и интегрального метода путeм декомпозиции области на составляющие, допускающие применение в них этих методов. В 80-х годах прошлого века K. Feng и D. Yu предложили разбить всю область решения внешней краевой задачи на две подобласти: без пересечения или с пересечением (см. [5]). При декомпозиции без пересечения (см. фиг. 1a) вся область вне поверхности тела $\partial {\kern 1pt} D$ разбивается на ограниченную область S₊ между сферой $\partial {\kern 1pt} S$ поверхностью $\partial {\kern 1pt} D$, и на область S_–, ограниченную только снизу сферой $\partial {\kern 1pt} S$. При декомпозиции с пересечением (фиг. 1б) ограниченная область ${{\Omega }_{ + }}$ находится между поверхностью $\partial \,\Omega $ и поверхностью $\partial {\kern 1pt} D$, а неограниченная область ${{\Omega }_{ - }}$, так же как и в декомпозиции без пересечения, ограничена только снизу сферой $\partial {\kern 1pt} S$.

Для решения исходной задачи был предложен альтернирующий метод Шварца, в котором внутренняя задача решалась методом конечных элементов, а условия согласования ставились на внешней поверхности. В [5] и других работах [6]–[8] для численного решения задачи было отдано предпочтение декомпозиции области без пересечения (фиг. 1а). В предложенном методе нет необходимости искать решение внешней задачи в области вне сферической поверхности. Достаточно найти только значения нормальной производной на этой поверхности, если на ней заданы значения функции. Производные находятся вычислением значений оператора Пуанкаре–Стеклова, которые могут быть определены в виде ряда [6].

Необходимо отметить, что при решении внешней краевой задачи для искусственной поверхности $\partial {\kern 1pt} S$ необходимо вычислять интегралы с сингулярными ядрами, что вызывает необходимость использовать для корректного численного вычисления таких интегралов квадратуры специального вида. В настоящей работе для решения задачи предлагается использовать декомпозицию области с пересечением (фиг. 1б). В этом случае при решении внешней краевой задачи все вычисляемые интегралы не имеют сингулярностей, что значительно упрощает использование квадратур для их численного вычисления.

Предлагаемый в настоящей работе подход к решению внешней краевой задачи Дирихле для уравнения Гельмгольца состоит в нахождении приближeнных значений искомой функции и еe нормальной производной на поверхности вспомогательной сферы $\partial {\kern 1pt} S$, заключающей в себя исходную внутреннюю границу $\partial {\kern 1pt} D$, поскольку тогда можно найти значения функции в другой точке области, используя формулу Грина. При этом внешняя граница $\partial \,\Omega $ ограниченной подобласти ${{\Omega }_{ + }}$ не обязательно является сферой и выбирается из условия экономии вычислительных ресурсов и получения наиболее точного решения на сфере $\partial \,S$. Для решения задачи применяется альтернирующий метод Шварца, в котором последовательно пересчитываются граничные условия на сфере $\partial \,S$ и внешней вспомогательной поверхности $\partial \,\Omega $. На каждой итерации метода решаются внутренняя и внешняя краевые задачи. Решение внутренней краевой задачи производится в области, ограниченной внешней вспомогательной поверхностью $\partial \,\Omega $ и исходной границей $\partial {\kern 1pt} D$ (область ${{\Omega }_{ + }}$ на фиг. 1б), и находятся приближeнные значения искомой функции и еe нормальной производной на поверхности сферы $\partial \,S$. Для решения внешней задачи в области ${{\Omega }_{ - }}$, ограниченной снизу сферой $\partial \,S$, в отличие от рассмотренных выше методов, предлагается использовать формулу Грина. Такой подход имеет следующие преимущества. Во-первых, формула Грина является более универсальной, и позволяет находить решения уравнения Гельмгольца не только при положительном коэффициенте в уравнении, но и при произвольном, в том числе комплексном коэффициенте. Во-вторых, эта формула не представлена в виде ряда, как в рассмотренных выше методах, и поэтому нет необходимости аппроксимировать этот ряд частичной суммой, внося тем самым дополнительную погрешность в численный метод. Для применения формулы Грина необходимо знать значения функции и еe нормальной производной на сфере $\partial \,S$, а применение метода конечных объeмов при решении задачи в области ${{\Omega }_{ + }}$ позволяет получить и те, и другие значения на сфере $\partial \,S$. Отметим, что в предлагаемом итерационном алгоритме искомыми значениями при решении внутренней краевой задачи будут значения функции и еe нормальной производной только на сфере $\partial \,S$, а при решении внешней краевой задачи – только значения функции на поверхности $\partial \,\Omega $. В статье получены достаточные условия сходимости предложенного итерационного метода и оценка для производной решения уравнения Гельмгольца при отрицательном коэффициенте в этом уравнении.

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

1. ПОСТАНОВКА ЗАДАЧИ И МЕТОД РЕШЕНИЯ

Рассмотрим открытую область D в пространстве ${{\mathbb{R}}^{3}}$, ограниченную поверхностью $\partial {\kern 1pt} D$. Внешняя краевая задача для уравнения Гельмгольца состоит в нахождении функции $u\, \in \,{{C}^{1}}({{\mathbb{R}}^{3}}{{\backslash }}D)\, \cap \,{{C}^{2}}({{\mathbb{R}}^{3}}{{\backslash }}\bar {D})$, удовлетворяющей уравнению

краевому условию и условию излучения на бесконечности

Задача (1)–(3) предполагается решeнной, если удастся найти значения искомой функции ${{u}_{S}}$ и еe нормальной производной u_n на сфере $\partial {\kern 1pt} S$, являющейся границей шара $\bar {S} = S \cup \partial \,S$, $\bar {D} \subset S$. Действительно, в этом случае можно найти значения функции в произвольной точке пространства $r \in {{\mathbb{R}}^{3}}{{\backslash }}\bar {D}$ по формуле Грина:

где r₀ – радиус сферы $\partial {\kern 1pt} S$, ${{r}_{0}} \in \partial {\kern 1pt} S$, $R = \left| {\,r - {{r}_{0}}} \right| = \sqrt {{{r}^{2}} - 2{{r}_{{\,0}}}r\cos \gamma + r_{{\,0}}^{2}} $, γ – угол между векторами r и r₀, $\cos \gamma = \cos \theta \cos {{\theta }_{0}} + \sin \theta \sin {{\theta }_{0}}\cos \left( {\varphi - {{\varphi }_{{\,0}}}} \right)$, $\frac{\partial }{{\partial \,n}}$ – производная по нормали к поверхности в точке r₀.

Формула (4) является основой для предлагаемого итерационного метода решения задачи (1)–(3), который заключается в последовательном решении вспомогательных краевых задач в ограниченной и неограниченной области.

Введeм дополнительную ограниченную область ${{\Omega }_{0}}$ с границей $\partial \,{{\Omega }_{0}}$ такую, что $\bar {S} \subset {{\Omega }_{0}}$. Будем решать методом конечных объeмов внутреннюю краевую задачу в области ${{\Omega }_{0}}{{\backslash }}\bar {D}$ при заданных граничных условиях на $\partial {\kern 1pt} D$ и пересчитываемых на каждой итерации краевых условиях на $\partial \,{{\Omega }_{0}}$. На первой итерации метода задаeм, например, нулевые граничные условия на поверхности $\partial \,{{\Omega }_{0}}$ и решаем внутреннюю краевую задачу в области ${{\Omega }_{0}}{{\backslash }}\bar {D}$. При этом нас будут интересовать полученные значения для функции и еe нормальной производной только на поверхности $\partial \,S$. По этим значениям определим новые граничные условия на внешней границе расчeтной области $\partial \,{{\Omega }_{0}}$ по формуле (4), и проведeм аналогично вторую и последующие итерации. Отметим, что на разных шагах этого процесса расчeтные области для внутренних краевых задач можно выбирать динамическим образом, т.е. на j-й итерации вместо ${{\Omega }_{0}}$ определить область ${{\Omega }_{{j + 1}}}$ с границей $\partial \,{{\Omega }_{{j + 1}}}$, $j = 0,1, \ldots ,J$. Количество итераций J определяется по условию сходимости метода с заданной точностью $\varepsilon \ll 1$. Выбор области ${{\Omega }_{j}}$, заключающей в себя шар $\bar {S}$, на второй и всех последующих итерациях производится из соображений экономии вычислительных ресурсов и возможности получения наиболее точного решения на поверхности шара. При этом условие включения $S \subset {{\Omega }_{j}}$ остаeтся в силе. На всех последующих итерациях решается внутренняя краевая задача в области ${{\Omega }_{j}}{{\backslash }}\bar {D}$ и по формуле (4) восстанавливаются новые граничные условия Дирихле на поверхности $\partial \,{{\Omega }_{{\,j}}}$.

2. СХОДИМОСТЬ ИТЕРАЦИОННОГО МЕТОДА

Исследуем сходимость предложенного метода для решения уравнения

с граничными условиями (2), (3), где k – вещественное число.

2.1. Итерационные формулы для искомого приближeнного решения и его погрешности

Запишем формулу (4) в виде

где $G(r,{{r}_{0}}) = \frac{1}{{4\pi }}\frac{{{{e}^{{ - kR}}}}}{R}$. Тогда итерационный процесс определяется следующим образом: где

Определим погрешность метода как

Тогда получим

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

3. АНАЛИЗ СХОДИМОСТИ ЧАСТНОГО РЕШЕНИЯ С ПРОИЗВОЛЬНЫМ ВОЛНОВЫМ ЧИСЛОМ

В предыдущем разделе исследована сходимость метода для решения частного случая задачи (1)–(3) с отрицательным коэффициентом в уравнении Гельмгольца. В этом разделе будет исследована сходимость частного решения с произвольным волновым числом, что позволит сделать вывод о нецелесообразности использования предложенного метода декомпозиции с пересечением для решения задач с произвольными волновыми числами.

Проведeм анализ сходимости альтернирующего метода Шварца для решения задачи (1)–(3) для частного случая, когда поверхность $\partial \,D$ является сферой радиуса ${{r}_{D}}$, функция $f$ в формуле (2) – константа, а решением задачи (1)–(3) является функция $u(r) = \frac{{{{e}^{{\,i\,k\,r}}}}}{r}$.

Пусть на сфере $\partial \,S$ задано некоторое начальное значение искомой функции, являющейся константой: ${{u}^{0}}\left( {{{r}_{0}}} \right) = V_{R}^{0} + i\,V_{I}^{0} = const$. Выберем в качестве области ${{\Omega }_{0}}$ шар, ограниченный сферой радиуса ${{r}_{\Omega }}$, ${{r}_{D}} < {{r}_{0}} < {{r}_{\Omega }}$. Исследование сходимости сведeтся к следующим шагам.

Шаг 1. По заданному начальному значению искомой функции на сфере $\partial \,S$ решим внешнюю задачу и найдeм значение функции на сфере $\partial \,{{\Omega }_{0}}$.

Шаг 2. Решим внутреннюю задачу в области ${{\Omega }_{{\text{0}}}}{{\backslash }}\bar {D}$ с учeтом полученного значения на $\partial \,{{\Omega }_{0}}$ и заданного на $\partial \,D$, и найдeм новое значение для искомой функции на сфере $\partial \,S$.

Шаг 3. Сравним новые и начальные значения на сфере $\partial \,S$ с точными значениями и определим условия для убывания погрешности решения.

Детализируем предложенные шаги в данной постановке задачи.

Шаг 1. Решением внешней задачи для уравнения Гельмгольца вне сферы $\partial \,S$ является функция $\left( {{{C}_{R}} + i{{C}_{I}}} \right)\frac{{{{e}^{{\,i\,k\,r}}}}}{r}$, где

поэтому где $\phi = k\left( {{{r}_{\Omega }} - {{r}_{0}}} \right)$.

Шаг 2. Решением внутренней задачи для уравнения Гельмгольца в области ${{\Omega }_{0}}{{\backslash }}\bar {D}$ с граничными условиями ${{u}^{{\,1}}}\left( {{{r}_{\Omega }}} \right) = {{U}_{R}} + i{{U}_{I}}$ и ${{u}^{{\,1}}}\left( {{{r}_{D}}} \right) = {{u}_{R}} + i\,{{u}_{I}}$ является функция

где

Если значения U_R и U_I определены формулой (26), а значения функции на сфере $\partial \,D$ равны

то функция u¹(r) из формулы (27) примет на сфере $\partial \,S$ значение где

Шаг 3. Перейдeм к уравнениям для погрешностей решения. Введем обозначения:

Тогда $\varepsilon _{R}^{{\,1}} = {{a}_{R}}\,\varepsilon _{R}^{0} + {{b}_{R}}\,\varepsilon _{I}^{0}$, $\varepsilon _{I}^{{\,1}} = {{a}_{I}}\,\varepsilon _{R}^{0} + {{b}_{I}}\,\varepsilon _{I}^{0}$, или

где

Формула (30) задаeт преобразование вращения, помноженное на постоянный множитель χ, который и определяет сходимость метода для рассмотренного частного случая.

Условие сходимости $\left| \chi \right| < 1$ для рассмотренной простейшей задачи имеет непростую структуру ввиду наличия бесконечного количества интервалов расходимости. Для более общего случая условия сходимости, скорее всего, примут ещe более сложный вид. По этой причине авторы не рекомендуют решать уравнение Гельмгольца с произвольным волновым числом, отличное от вида (5), предложенным методом.

Рассмотрим условие сходимости модельной задачи с модификацией предложенного метода вида (22). Пусть ${{\eta }_{j}} = \eta = const$. Формулы (28) для данной модификации метода будут:

Перейдя к погрешностям решения $\varepsilon _{R}^{{\,i}}$ и $\varepsilon _{I}^{{\,i}}$ в формулах (29), получим

Выберем параметр η так, чтобы модуль погрешности ${\text{|}}{{\varepsilon }^{1}}{\text{|}} = \sqrt {{{{(\varepsilon _{R}^{1})}}^{2}} + {{{(\varepsilon _{I}^{1})}}^{2}}} $ был минимальным. Нетрудно показать, что оптимальное значение параметра ${{\eta }_{{opt}}}$ не зависит от |ε⁰| и равно

Преобразовав формулу (32) с учeтом явного вида параметров χ и $\phi $, получим

Рассмотрим убывание погрешности метода. Условие $\rho = {{{\text{|}}{{\varepsilon }^{1}}{\text{|}}} \mathord{\left/ {\vphantom {{{\text{|}}{{\varepsilon }^{1}}{\text{|}}} {{\text{|}}{{\varepsilon }^{0}}{\text{|}}}}} \right. \kern-0em} {{\text{|}}{{\varepsilon }^{0}}{\text{|}}}} < 1$ равносильно

Выбрав в качестве параметра η в формуле (33) значение ${{\eta }_{{opt}}}$ из (32), получим

Воспользовавшись явным видом параметра χ в формуле (31), получим окончательно

Таким образом, преимущество более общего подхода (22) для решения задачи состоит в том, что в нeм существует оптимальное значение параметра релаксации $\eta = \eta \left( {k,{{r}_{D}},{{r}_{0}},{{r}_{\Omega }}} \right)$ в (22), обеспечивающее сходимость метода при любых значениях входных параметров $k,\;{{r}_{D}},\;{{r}_{0}},\;{{r}_{\Omega }}$, за исключением случаев, когда $k\left( {{{r}_{0}} - {{r}_{D}}} \right) = {{\pi \,n} \mathord{\left/ {\vphantom {{\pi \,n} 2}} \right. \kern-0em} 2}$, где n – любое целое нечeтное число.

Полученные значения интервалов и скорости сходимости метода позволяют предположить, что в случае произвольной исходной поверхности решение внешней задачи предложенными методами имеет аналогичный характер или удовлетворяет ещe более сложным условиям. А именно: сходимость метода (6), (7) для решения задачи (1)–(3) имеет место только на интервалах, где выполняются определeнные соотношения между входными параметрами, характеризующими поверхности $\partial \,D$ и $\partial \,\Omega $, волновым числом k, и радиусом вспомогательной сферы r₀. Для модификации метода с релаксацией (22) сходимость будет иметь место при некоторых значениях параметров ${{\eta }_{j}}$ этого метода. Выбор соотношений для исходных параметров в первом варианте метода и оптимального значения параметра для второго варианта, необходимых для их сходимости в общем случае, является пока открытой проблемой. Таким образом, рекомендация о нецелесообразности использования также и модифицированного метода с релаксацией для решения задачи (1)–(3) с произвольным волновым числом остаeтся в силе.

4. АЛГОРИТМ НАХОЖДЕНИЯ ПРИБЛИЖЕННОГО РЕШЕНИЯ

Приведем принципиальную структуру программы для нахождения приближенного решения задачи (5), (2), (3).

Шаг 1. Задание начальных данных:

– задание трехмерной неравномерной структурированной тетраэдральной сетки, приближенно описывающей трехмерную замкнутую область ${{\Omega }_{j}}{{\backslash }}D$;

– задание равномерного разбиения сферы $\partial \,S$;

– задание критерия останова для итерационного решения внутренней краевой задачи ${{\partial }_{i}}$;

– задание критерия останова для внешнего итерационного процесса ${{\partial }_{0}}$;

Шаг 2. По заданной сетке проводится аппроксимация внутренней краевой задачи (6) хорошо известным методом барицентрических конечных объемов [12], [13], которая приводит к системе линейных алгебраических уравнений относительно приближeнных значений искомой функции u_h в открытой области ${{\Omega }_{j}}{{\backslash }}\bar {D}{{\backslash }}\partial {{\Omega }_{j}}$, значений ${{u}_{{\partial \,{{\Omega }_{j}}}}}$ в узлах сетки на поверхности $\partial \,{{\Omega }_{j}}$, и заданных значений ${{u}_{{\partial D}}} = {{\left( {f\left( r \right)} \right)}_{h}}$ в узлах сетки на поверхности $\partial \,D$:

где А – матрица баланса для внутренних узлов сетки, а ${{{\text{B}}}_{{\partial D}}}$ и ${{{\text{B}}}_{{\partial \,{{\Omega }_{j}}}}}$ – матрицы учета краевых условий типа Дирихле от значений в узлах, лежащих на границах $\partial \,D$ и $\partial \,{{\Omega }_{j}}$. Данное представление СЛАУ обусловлено тем, что позволяет отказаться от аппроксимации внутренней краевой задачи на каждом шаге внешнего итерационного процесса, ограничившись лишь пересчетом вклада значений функций в узлах, лежащих на границе $\partial \,{{\Omega }_{j}}$ в правой части СЛАУ (34).

Шаг 3. Внешний итерационный процесс:

– методом сопряженных градиентов в подпространствах Крылова [14] решается СЛАУ для внутренней краевой задачи

где $F_{h}^{j} = {{B}_{{\partial \,{{\Omega }_{j}}}}}\,u_{{\partial \,{{\Omega }_{j}}}}^{{j - 1}} + {{B}_{{\partial D}}}\,{{u}_{{\partial D}}}$, $u_{{\partial \,{{\Omega }_{{\,0}}}}}^{{ - 1}} \equiv 0$. По полученному приближенному решению $u_{h}^{j}$ проводится интерполяция приближeнных значений функции $u_{{S,h}}^{j}$ и еe нормальной производной $u_{{n,h}}^{j}$ на сфере $\partial \,S$ в точках, необходимых для решения внешней краевой задачи;

– по полученным значениям функции $u_{{S,h}}^{j}$ и ее нормальной производной $u_{{n,h}}^{j}$ численно находятся приближенные значения $u_{{\partial \,{{\Omega }_{{j + 1}}}}}^{j}$на границе $\partial \,{{\Omega }_{{j + 1}}}$, для чего интеграл (4), записанный в сферических координатах в виде двойного интеграла, аппроксимируется квадратурной формулой Симпсона по каждой из переменных. По этим значениям на границе $\partial \,{{\Omega }_{{j + 1}}}$ рассчитывается правая часть $F_{h}^{{j + 1}}$ для СЛАУ внутренней краевой задачи;

– выполняется проверка на критерий останова внешнего итерационного процесса:

если выполняется:

то внешний итерационный процесс завершается, иначе повторяем шаг 3.

5. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ И ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Будем искать численное решение внешней краевой задачи (5) с граничными условиями (2), (3), где k – вещественное число, имеющее следующее точное решение

${\text{или}}$где

$r = \sqrt {\,{{x}^{2}} + {{y}^{2}} + {{z}^{2}}} $,    ${{r}_{{\,1}}} = \sqrt {\,{{{\left( {x - {{x}_{0}}} \right)}}^{2}} + {{y}^{2}} + {{z}^{2}}} $,    ${{r}_{{\,2}}} = \sqrt {\,{{{\left( {x + {{x}_{0}}} \right)}}^{2}} + {{y}^{2}} + {{z}^{2}}} $,    ${{x}_{0}} = 0.1$.

Зададим на гранях куба с рeбрами, равными 0.4, граничные условия в соответствии с точным решением (36) или (37). Будем полагать, что куб, а также и все последующие в рассмотрении области, ограниченные выбранными поверхностями, имеют центр симметрии в начале координат. Выберем границу внешней расчeтной области $\partial \,{{\Omega }_{0}}$, которую оставим неизменной на каждой итерации, в виде поверхности куба с рeбрами, равными 2, на поверхности которого зададим нулевые граничные условия. В области ${{\Omega }_{0}}{{\backslash }}\bar {D}$, ограниченной поверхностями этих кубов, решается внутренняя задача Дирихле и находятся значения функции и еe нормальной производной на сфере $\partial \,S$, лежащей в области ${{\Omega }_{0}}{{\backslash }}\bar {D}$. По этим значениям определяются новые граничные условия на поверхности $\partial \,{{\Omega }_{0}}$ по формуле (4). Далее снова решается задача Дирихле в расчeтной области ${{\Omega }_{0}}{{\backslash }}\bar {D}$, и определяются новые приближeнные значения функции и еe нормальной производной на сфере. Поверхность сферы, на которой ищется приближeнное решение, является неизменной на всех итерациях в каждом численном эксперименте, однако может варьироваться в разных экспериментах.

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

редкая сетка L = 29 666, Ltet = 160 704;

средняя сетка L = 225 650, Ltet = 1 285 632;

густая сетка L = 1 759 794, Ltet = 1 028 5056.

Обозначим через ${{\varepsilon }_{{v}}}$, ${{\varepsilon }_{d}}$ и ${{\delta }_{{v}}}$, ${{\delta }_{d}}$ среднеквадратичные и максимальные погрешности отклонения от точного решения и от точной нормальной производной на сфере $\partial \,S$, вычисленные по формулам

где

являются точными и приближeнными значениями искомой функции и еe нормальной производной в узлах ${{\theta }_{j}},{{\varphi }_{{\,i}}}$ сетки в сферических координатах на сфере.

В табл. 1 представлены погрешности вычисления функции и еe нормальной производной на сфере радиуса ${{r}_{{\,0}}} = 0.5$ при ${{k}^{2}} = 1$, ${{N}_{\theta }} = {{N}_{\varphi }} = 17$. Число внешних итераций N_it для достижения критерия останова (35) итерационного процесса при ${{\partial }_{0}} = {{10}^{{ - 7}}}$ было ${{N}_{{it}}} = 3$. Точное решение определялось формулой (37). Вычисления проводились на редкой, средней и подробной сетках, в которых длина ребра элементарного тетраэдра в методе конечных объeмов уменьшалась приблизительно в 2 раза при переходе к более подробной сетке. Заметим, что увеличение количества точек на сфере не приводило к существенному уменьшению вычисляемых погрешностей, однако требовало дополнительных временных затрат. Это обстоятельство иллюстрирует то, что погрешность решения задачи обусловлена, прежде всего, погрешностью решения еe внутренней части.

Результаты, приведeнные в таблице, демонстрируют квадратичную сходимость полученного приближeнного решения и линейную сходимость для его производной при сгущении сетки и отсутствие зависимости числа внешних итераций от шага сетки.

В следующем численном эксперименте определялись погрешности метода при разных значениях радиуса вспомогательной сферы r₀ и разных коэффициентах k² в уравнении Гельмгольца. Вычисления проводились на подробной сетке, точное решение было задано формулой (36). Максимальные погрешности приближeнного значения функции ${{\delta }_{{v}}}$, ее нормальной производной на сфере δ_d и количество внешних итераций ${{N}_{{it}}}$, полученные при достижении условия (35), где ${{\partial }_{0}} = {{10}^{{ - 7}}}$ для разных значений радиуса вспомогательной сферы, приведены в табл. 2.

Для анализа полученных значений воспользуемся формулой (20) для коэффициента уменьшения погрешности ${\rm M}\left( {k,d,{{r}_{0}},{{x}_{0}}} \right)$ при переходе к следующей итерации метода. При заданных значениях параметров исходной задачи этот коэффициент будет меньше единицы при $0.487 \leqslant {{r}_{0}} \leqslant 0.696$ для k² = 1, и при $0.362 < {{r}_{0}} < 0.846$ для k² = 10. Полученные значения для погрешностей позволяют сделать вывод, что сходимость метода имеет место при больших значениях радиуса вспомогательной сферы, вплоть до ${{r}_{0}} < 0.95$. Это обстоятельство не противоречит условию сходимости метода (21), поскольку оно является достаточным.

Для обеспечения надeжности численных вычислений необходимо выбрать параметры метода, характеризующие расположение внешней поверхности и вспомогательной сферы таким образом, чтобы условие сходимости метода (21) было выполнено. Если при некотором выборе внешней границы расчeтной области ∂Ω_j условие ${\text{M}} < 1$ не выполняется для любых радиусов вспомогательных сфер, расположенных между поверхностями $\partial \,D$ и ∂Ω_j, то необходимо переместить внешнюю поверхность ∂Ω_j дальше от исходной поверхности ∂D. Действительно, множитель e^–kd/2d в формуле (20) быстро убывает с увеличением параметра d, где d – расстояние от начала координат до поверхности ∂Ω_j. Однако рост этого параметра влечeт увеличение расчeтной области для решения внутренней краевой задачи, что является нежелательным ввиду использования большего объeма вычислительных ресурсов для решения задачи. Другим параметром декомпозиции является радиус вспомогательной сферы r₀. Этот параметр связан со вторым множителем в формуле (20). При расположении сферы вблизи поверхности ∂D или ∂Ω_j параметры ${{\rho }_{0}}$ и ${{x}_{0}}$ уменьшаются, что приводит к увеличению коэффициента $\beta $, значит, и всего коэффициента M. По этой причине целесообразно размещать вспомогательную сферу $\partial \,S$ где-нибудь посередине между поверхностями ∂D и ∂Ω_j для уменьшения коэффициента $\beta $, который монотонно убывает при $0 < x < 2.08$.

Итак, подытожим сказанное о выборе параметров метода. Для уменьшения времени счeта и экономии вычислительных ресурсов выбираем внешнюю поверхность ∂Ω_j, расположенную относительно близко к исходной поверхности ∂D. После этого по формуле (20) производим вычисления коэффициента M при разных значениях радиуса вспомогательной сферы r₀. Если условие (21) при выбранном расположении вспомогательной поверхности ∂Ω_j не выполнено для всех сфер, расположенных между поверхностями ∂D и ∂Ω_j, то отодвигаем поверхность ∂Ω_j дальше от поверхности ∂D и снова производим вычисления коэффициента M для нового набора значений радиусов сфер r₀ между этими областями. Такую процедуру необходимо проводить до выполнения условия (21), после чего можно приступать к численным расчeтам по решению задачи (5) с условиями (2) и (3).