Журнал вычислительной математики и математической физики, 2023, T. 63, № 2, стр. 230-244
Улучшенная квадратурная формула для потенциала простого слоя
П. А. Крутицкий 1, *, И. О. Резниченко 2, **
1 ИПМ им. М.В. Келдыша РАН
125047 Москва, Миусская пл. 4, Россия
2 МГУ им. М.В. Ломоносова, 1, стр. 2, физический факультет
119991 Москва, Ленинские горы, Россия
* E-mail: biem@mail.ru
** E-mail: io.reznichenko@physics.msu.ru
Поступила в редакцию 05.05.2022
После доработки 19.05.2022
Принята к публикации 04.08.2022
- EDN: BMYLDN
- DOI: 10.31857/S0044466923020114
Аннотация
Выводится улучшенная квадратурная формула для потенциала простого слоя с гладкой плотностью, заданной на замкнутой либо разомкнутой поверхности. Формула дает равномерную аппроксимацию потенциала вблизи поверхности и сохраняет свойство непрерывности потенциала при стремлении точки наблюдения из области к поверхности, что подтверждается численными тестами. Предложенная в работе квадратурная формула дает более высокую точность при вычислении потенциала вблизи поверхности, чем известные квадратурные формулы, что также подтверждается численными тестами. Кроме того, выводится квадратурная формула для прямого значения потенциала простого слоя на поверхности. Для этой формулы проведены численные тесты, подтверждающие ее эффективность и точность. Библ. 20. Табл. 5.
1. ВВЕДЕНИЕ
Метод потенциалов позволяет решать краевые задачи для уравнений Лапласа и Гельмгольца путем сведения их к граничным интегральным уравнениям относительно плотности потенциала. Находя плотность потенциала из интегрального уравнения, можно получить решение краевой задачи. Метод потенциалов описан в классических монографиях [1–3]. При численной реализации метода потенциалов необходимо иметь квадратурные формулы для потенциалов, позволяющие вычислять потенциалы с хорошей точностью во всей области, где решается краевая задача.
Стандартные квадратурные формулы для потенциала простого слоя для уравнений Лапласа и Гельмгольца, используемые в инженерных расчетах, не дают равномерной аппроксимации потенциала вблизи поверхности $\Gamma $, на которой задана плотность потенциала, и даже стремятся к бесконечности, когда точка, в которой вычисляется квадратурная формула, стремится к определенным точкам на поверхности $\Gamma $ (см. [4], гл. 2), тогда как сам потенциал непрерывен во всем пространстве, в том числе во всех точках на поверхности $\Gamma $. Следовательно, стандартные квадратурные формулы не сохраняют важнейшее свойство потенциала, а именно, его ограниченность и непрерывность на поверхности $\Gamma $. В [5], [6] предложена квадратурная формула, которая сохраняет указанное свойство потенциала простого слоя. В настоящей работе проведено дальнейшее улучшение квадратурной формулы из [5], [6] с учетом гладкости коэффициента, появляющегося при параметризации в потенциале. Благодаря сделанным преобразованиям, канонический интеграл, возникающий в квадратурной формуле, вычисляется с более высокой точностью, чем в [5], [6]. Это позволяет усовершенствовать квадратурную формулу из [5], [6] и предложить улучшенную формулу, обеспечивающую повышенную точность вычислений.
В двумерном случае улучшенная квадратурная формула для гармонического потенциала простого слоя с плотностью, заданной на разомкнутых кривых, предложена в [7], [8]. Эта формула может применяться при нахождении численных решений краевых задач для уравнений Лапласа и Гельмгольца вне разрезов и разомкнутых кривых на плоскости. Такие задачи изучались в [9–16].
2. ПОСТАНОВКА ЗАДАЧИ
Введем в пространстве декартову систему координат $x = ({{x}_{1}},{{x}_{2}},{{x}_{3}}) \in {{R}^{3}}$. Пусть $\Gamma $ – простая гладкая замкнутая либо ограниченная разомкнутая поверхность класса ${{C}^{2}}$, содержащая свои предельные точки (см. [17], гл. 14, § 1). Если поверхность $\Gamma $ замкнутая, то она должна ограничивать объемно-односвязную внутреннюю область (см. [18], c. 201). Предположим, что поверхность $\Gamma $ параметризована так, что на нее отображается прямоугольник:
(1)
$\begin{gathered} y = ({{y}_{1}},{{y}_{2}},{{y}_{3}}) \in \Gamma ,\quad {{y}_{1}} = {{y}_{1}}(u,v),\quad {{y}_{2}} = {{y}_{2}}(u,v),\quad {{y}_{3}} = {{y}_{3}}(u,v); \\ u \in [0,A],\quad v \in [0,B];\quad {{y}_{j}}(u,v) \in {{C}^{2}}([0,A] \times [0,B]),\quad j = 1,2,3. \\ \end{gathered} $Введем $N$ точек ${{u}_{n}}$ с шагом $h$ на отрезке $[0,A]$ и $M$ точек ${{v}_{m}}$ на отрезке $[0,B]$ и рассмотрим разбиение прямоугольника $[0,A] \times [0,B]$, который отображается на поверхность $\Gamma $:
Известно (см. [17], гл. 14, § 1), что компоненты вектора нормали (не единичного) $\eta (y) = ({{\eta }_{1}}(y),{{\eta }_{2}}(y),{{\eta }_{3}}(y))$ в точке поверхности $y = ({{y}_{1}},{{y}_{2}},{{y}_{3}}) \in \Gamma $ определяются через определители второго порядка формулами
(2)
${{\eta }_{1}} = \left| {\begin{array}{*{20}{c}} {{{{({{y}_{2}})}}_{u}}}&{{{{({{y}_{3}})}}_{u}}} \\ {{{{({{y}_{2}})}}_{v}}}&{{{{({{y}_{3}})}}_{v}}} \end{array}} \right|,\quad {{\eta }_{2}} = \left| {\begin{array}{*{20}{c}} {{{{({{y}_{3}})}}_{u}}}&{{{{({{y}_{1}})}}_{u}}} \\ {{{{({{y}_{3}})}}_{v}}}&{{{{({{y}_{1}})}}_{v}}} \end{array}} \right|,\quad {{\eta }_{3}} = \left| {\begin{array}{*{20}{c}} {{{{({{y}_{1}})}}_{u}}}&{{{{({{y}_{2}})}}_{u}}} \\ {{{{({{y}_{1}})}}_{v}}}&{{{{({{y}_{2}})}}_{v}}} \end{array}} \right|.$Заметим, что если ${\text{|}}\eta (y(u,v)){\kern 1pt} {\text{|}} = 0$ в некоторой точке, то функция ${\text{|}}\eta (y(u,{v})){\kern 1pt} {\text{|}}$ может быть недифференцируемой в этой точке. Поэтому дополнительно потребуем, чтобы
Кроме того, потребуем, чтобы
Из условия (4) следует, что ${\text{|}}\eta (y(u,v)){\kern 1pt} {\text{|}} \in {{C}^{1}}((0,A) \times (0,B))$, но условие (3) не следует.Потенциал простого слоя используется при решении краевых задач методом интегральных уравнений (см. [19]). Рассмотрим потенциал простого слоя для уравнения Гельмгольца с заданной на поверхности $\Gamma $ плотностью $\mu (y) \in {{C}^{0}}(\Gamma )$:
(5)
$\begin{gathered} {{\mathcal{V}}_{k}}[\mu ](x) = \frac{1}{{4\pi }}\int\limits_\Gamma \frac{{\mu (y){{e}^{{ik|x - y|}}}}}{{{\text{|}}x - y{\kern 1pt} {\text{|}}}}d{{s}_{y}} = \frac{1}{{4\pi }}\int\limits_0^A \,du\int\limits_0^B \,dv\frac{{\mu (y(u,v))\exp (ik{\kern 1pt} {\text{|}}x - y(u,v){\kern 1pt} {\text{|}}{\kern 1pt} )}}{{{\text{|}}x - y(u,v){\kern 1pt} {\text{|}}}}{\text{|}}\eta (y(u,v)){\kern 1pt} {\text{|}} = \\ = \frac{1}{{4\pi }}\sum\limits_{n = 0}^{N - 1} {\kern 1pt} {\kern 1pt} \sum\limits_{m = 0}^{M - 1} \,\int\limits_{{{u}_{n}} - h/2}^{{{u}_{n}} + h/2} \,\int\limits_{{{v}_{m}} - H/2}^{{{v}_{m}} + H/2} \frac{{\mu (y(u,v))\exp (ik{\text{|}}x - y(u,v){\text{|}})}}{{{\text{|}}x - y(u,v){\text{|}}}}{\text{|}}\eta (y(u,v)){\kern 1pt} {\text{|}}{\kern 1pt} dudv, \\ \end{gathered} $Пусть ${{\mu }_{{nm}}} = \mu (y({{u}_{n}},{{v}_{m}}))$, тогда
для $u \in [{{u}_{n}} - h{\text{/}}2,{{u}_{n}} + h{\text{/}}2]$ и $v \in [{{v}_{m}} - H{\text{/}}2,{{v}_{m}} + H{\text{/}}2]$.В [5], [6] показано, что при $u \in [{{u}_{n}} - h{\text{/}}2,{{u}_{n}} + h{\text{/}}2]$ и $v \in [{{v}_{m}} - H{\text{/}}2,{{v}_{m}} + H{\text{/}}2]$
Следовательно, используя (5), имеем
(7)
${{\mathcal{V}}_{k}}[\mu ](x) \approx \frac{1}{{4\pi }}\sum\limits_{n = 0}^{N - 1} {\kern 1pt} {\kern 1pt} \sum\limits_{m = 0}^{M - 1} \,{{\mu }_{{nm}}}\exp (ik{\kern 1pt} {\text{|}}x - y({{u}_{n}},{{{v}}_{m}}){\kern 1pt} {\text{|}}{\kern 1pt} )\int\limits_{{{u}_{n}} - h/2}^{{{u}_{n}} + h/2} du\int\limits_{{{{v}}_{m}} - H/2}^{{{{v}}_{m}} + H/2} d{v}\frac{{{\text{|}}\eta (y(u,{v})){\kern 1pt} {\text{|}}}}{{{\text{|}}x - y(u,{v}){\kern 1pt} {\text{|}}}},$Как видно из соотношения (7), для получения квадратурной формулы необходимо вычислить интеграл
который будем называть каноническим интегралом.3. КАНОНИЧЕСКИЙ ИНТЕГРАЛ
Пусть точка $x$ не принадлежит кусочку поверхности, вдоль которого изменяется точка $y = y(u,v)$ в интеграле (8). Центром этого кусочка является точка $y({{u}_{n}},{{v}_{m}})$. В интеграле (8) $(u - {{u}_{n}}) \in [ - h{\text{/}}2,h{\text{/}}2]$, $(v - {{v}_{m}}) \in [ - H{\text{/}}2,H{\text{/}}2]$. Разложим ${{y}_{j}}(u,v)$ по формуле Тейлора с центром в точке $({{u}_{n}},{{v}_{m}})$, тогда для $j = 1,2,3$ получим
где Все производные берутся в точке $({{u}_{n}},{{{v}}_{m}})$. Положим(9)
${{\alpha }^{2}}{{\beta }^{2}} - {{\delta }^{2}} = {\text{|}}\eta (y({{u}_{n}},{{v}_{m}})){\kern 1pt} {{{\text{|}}}^{2}}.$В силу (3) для всех возможных $n,\;m$ при $u \in [{{u}_{n}} - h{\text{/}}2,{{u}_{n}} + h{\text{/}}2]$ и $v \in [{{v}_{m}} - H{\text{/}}2,{{v}_{m}} + H{\text{/}}2]$ функция ${\text{|}}\eta (y(u,v)){\text{|}}$ может быть разложена в точке $({{u}_{n}},{{v}_{m}})$ по формуле Тейлора с остаточным членом в форме Пеано (см. [17], гл. 10, § 5.3):
(11)
$\begin{gathered} {\text{|}}\eta (y(u,v)){\kern 1pt} {\text{|}} = {\text{|}}\eta (y({{u}_{n}},{{v}_{m}})){\kern 1pt} {\text{|}}\; + \;{\text{|}}\eta {\kern 1pt} {\text{|}}_{u}^{'}(u - {{u}_{n}})\; + \;{\text{|}}\eta {\kern 1pt} {\text{|}}_{{v}}^{'}(v - {{v}_{m}}) + o\left( {\sqrt {{{{(u - {{u}_{n}})}}^{2}} + {{{(v - {{v}_{m}})}}^{2}}} } \right) = \\ \, = {\text{|}}\eta (y({{u}_{n}},{{v}_{m}})){\kern 1pt} {\text{|}}\, + \;{\text{|}}\eta {\kern 1pt} {\text{|}}_{u}^{'}U\; + \;{\text{|}}\eta {\kern 1pt} {\text{|}}_{{v}}^{'}V + o\left( {\sqrt {{{{(u - {{u}_{n}})}}^{2}} + {{{(v - {{v}_{m}})}}^{2}}} } \right). \\ \end{gathered} $Выражения для ${\text{|}}\eta {\kern 1pt} {\text{|}}_{u}^{'}$ и ${\text{|}}\eta {\kern 1pt} {\text{|}}_{{v}}^{'}$ можно найти с помощью формулы (9), эти выражения имеют вид
Целью этого раздела является вычисление интеграла (8), который можно записать в виде
(12)
$ = \;\int\limits_{ - h/2}^{h/2} dU\int\limits_{ - H/2}^{H/2} dV\frac{{{\text{|}}\eta (y({{u}_{n}},{{{v}}_{m}})){\kern 1pt} {\text{|}}\; - \;{\text{|}}\eta {\kern 1pt} {\text{|}}_{{v}}^{'}Q{\text{/}}{{\beta }^{2}} + U({\kern 1pt} {\text{|}}\eta {\kern 1pt} {\text{|}}_{u}^{'}\; - \;{\text{|}}\eta {\kern 1pt} {\text{|}}_{{v}}^{'}\delta {\text{/}}{{\beta }^{2}})\; + \;{\text{|}}\eta {\kern 1pt} {\text{|}}_{{v}}^{'}(V + \delta U{\text{/}}{{\beta }^{2}} + Q{\text{/}}{{\beta }^{2}})}}{{\beta \sqrt {{{{(V + \delta U{\text{/}}{{\beta }^{2}} + Q{\text{/}}{{\beta }^{2}})}}^{2}} + ( - {{{(\delta U + Q)}}^{2}}{\text{/}}{{\beta }^{2}} + {{\alpha }^{2}}{{U}^{2}} + 2PU + {{r}^{2}}){\text{/}}{{\beta }^{2}}} }} = $(13)
$\begin{gathered} {{J}_{1}}(H) = \int\limits_{ - h/2}^{h/2} U\ln \left| {\frac{H}{2} + \frac{{\delta U + Q}}{{{{\beta }^{2}}}} + \sqrt {{{{\left( {\frac{H}{2} + \frac{{\delta U + Q}}{{{{\beta }^{2}}}}} \right)}}^{2}} - \frac{{{{{(\delta U + Q)}}^{2}}}}{{{{\beta }^{4}}}} + \frac{{{{\alpha }^{2}}{{U}^{2}} + 2PU + {{r}^{2}}}}{{{{\beta }^{2}}}}} } \right|dU = \\ = \int\limits_{ - h/2}^{h/2} U\ln \left| {\varepsilon + {{\delta }_{0}}U + \sqrt {\alpha _{0}^{2}{{U}^{2}} + 2{{\beta }_{0}}U + {{\chi }_{0}}} } \right|dU, \\ \end{gathered} $(14)
${{J}_{2}}(H) = \int\limits_{ - h/2}^{h/2} \sqrt {{{{\left( {\frac{H}{2} + \frac{{\delta U + Q}}{{{{\beta }^{2}}}}} \right)}}^{2}} - \frac{{{{{(\delta U + Q)}}^{2}}}}{{{{\beta }^{4}}}} + \frac{{{{\alpha }^{2}}{{U}^{2}} + 2PU + {{r}^{2}}}}{{{{\beta }^{2}}}}} \,dU = \int\limits_{ - h/2}^{h/2} \sqrt {\alpha _{0}^{2}{{U}^{2}} + 2{{\beta }_{0}}U + {{\chi }_{0}}} \,dU.$3.1. Вычисление интеграла ${{J}_{1}}(H)$ при ${{\chi }_{1}} > 0$
Рассмотрим случай ${{\chi }_{1}} > 0$. Сделаем гиперболическую замену переменной в интеграле ${{J}_{1}}(H)$ по формулам
(15)
$\begin{gathered} {{J}_{{11}}} = \int\limits_{{{z}_{ - }}}^{{{z}_{ + }}} \frac{{({{z}^{2}} - 1)}}{{2z}}\,{\kern 1pt} \frac{{({{z}^{2}} - 1)({{\delta }_{1}}({{z}^{2}} + 1) + {{\chi }_{1}}({{z}^{2}} - 1))}}{{4{{\varepsilon }_{1}}{{z}^{2}} + 2z({{z}^{2}} - 1){{\delta }_{1}} + 2z({{z}^{2}} + 1){{\chi }_{1}}}}\,{\kern 1pt} \frac{{d{\kern 1pt} z}}{z} = \\ \, = \int\limits_{{{z}_{ - }}}^{{{z}_{ + }}} \frac{{({{z}^{2}} - 1)}}{{2z}}{\kern 1pt} \,\frac{{({{\delta }_{1}} + {{\chi }_{1}}){{z}^{4}} - 2{{\chi }_{1}}{{z}^{2}} + {{\chi }_{1}} - {{\delta }_{1}}}}{{2{{z}^{2}}\left( {({{\delta }_{1}} + {{\chi }_{1}}){{z}^{2}} + 2{{\varepsilon }_{1}}z + {{\chi }_{1}} - {{\delta }_{1}}} \right)}}d{\kern 1pt} z = \\ \end{gathered} $I. Пусть ${{\delta }_{ + }} \ne 0$.
Используя [20], п. 1.2.8.13, вычислим интеграл ${{J}_{0}}$. Если $\varepsilon _{1}^{2} - {{\delta }_{ + }}{{\delta }_{ - }} > 0$, то
Интеграл ${{J}_{{12}}}$ вычисляется, согласно [20], п. 1.2.8.21, и имеет вид
Интеграл ${{J}_{{13}}}$ вычисляется, согласно [20], пп. 1.2.8.19 и 1.2.4.17, и имеет вид
Для интеграла ${{J}_{{14}}}$ надо рассмотреть несколько случаев. Если $\varepsilon _{1}^{2} - {{\delta }_{ + }}{{\delta }_{ - }} \ne 0$ и ${{\delta }_{ - }} \ne 0$, то в соответствии с [20], п. 1.2.8.25 получаем
Интеграл ${{J}_{{15}}}$ находится в соответствии с [20], п. 1.2.8.27 и дается выражением
Замечание. Если $\varepsilon _{1}^{2} - {{\delta }_{ + }}{{\delta }_{ - }} = 0$ и ${{\delta }_{ - }} \ne 0$, то при вычислении интеграла ${{J}_{{11}}}$ линейная комбинация интегралов ${{J}_{0}}$ равна нулю.
II. Пусть ${{\delta }_{ + }} = 0$.
Как указано выше, в этом случае в формуле (15) отсутствует слагаемое с интегралом ${{J}_{{12}}}$, поэтому интеграл ${{J}_{{12}}}$ вычислять не нужно. Кроме того, ${{\delta }_{ - }} \ne 0$ если ${{\delta }_{ + }} = 0$, поэтому надо рассмотреть два случая: ${{\varepsilon }_{1}} \ne 0$ и ${{\varepsilon }_{1}} = 0$.
Если ${{\delta }_{ + }} = 0$ и ${{\varepsilon }_{1}} \ne 0$, то
Если ${{\delta }_{ + }} = 0$ и ${{\varepsilon }_{1}} = 0$, то
3.2. Вычисление интеграла ${{J}_{1}}(H)$ при ${{\chi }_{1}} = 0$
Поскольку ${{\alpha }_{0}} > 0$, подкоренное выражение в ${{J}_{1}}(H)$ принимает вид
1. Если $ - {{\beta }_{0}}{\text{/}}\alpha _{0}^{2} \in \left( { - h{\text{/}}2,h{\text{/}}2} \right)$, то
2. Если $ - {{\beta }_{0}}{\text{/}}\alpha _{0}^{2} \in \left[ {h{\text{/}}2, + \infty } \right)$, то
3. Если $ - {{\beta }_{0}}{\text{/}}\alpha _{0}^{2} \in \left( { - \infty , - h{\text{/}}2} \right]$, то
1. При ${{\kappa }_{1}} \ne 0$
2. При ${{\kappa }_{1}} = 0$
3.3. Вычисление интеграла ${{J}_{2}}(H)$
I. Пусть ${{\chi }_{1}} > 0.$
Интеграл ${{J}_{2}}$ вычисляется, согласно [20], п. 1.2.41.8, и имеет вид
II. Пусть ${{\chi }_{1}} = 0.$
Введем обозначения $t = {{\alpha }_{0}}U + {{\beta }_{0}}{\text{/}}{{\alpha }_{0}},$ ${{t}_{ \pm }} = {{\alpha }_{0}}( \pm h{\text{/}}2) + {{\beta }_{0}}{\text{/}}{{\alpha }_{0}}$. Поскольку ${{\alpha }_{0}} > 0$ и, очевидно, ${{t}_{ + }} > {{t}_{ - }},$ получаем
(16)
${{J}_{2}}(H) = \int\limits_{ - h/2}^{h/2} \sqrt {{{{({{\alpha }_{0}}U + {{\beta }_{0}}{\text{/}}{{\alpha }_{0}})}}^{2}}} \,dU = \frac{1}{{{{\alpha }_{0}}}}\int\limits_{{{t}_{ - }}}^{{{t}_{ + }}} \left| t \right|dt,\quad {{J}_{2}}(H) = \left\{ \begin{gathered} (t_{ + }^{2} - t_{ - }^{2}){\text{/}}(2{{\alpha }_{0}}),\quad 0 \leqslant {{t}_{ - }} < {{t}_{ + }}, \hfill \\ (t_{ + }^{2} + t_{ - }^{2}){\text{/}}(2{{\alpha }_{0}}),\quad {{t}_{ - }} < 0 < {{t}_{ + }}, \hfill \\ - (t_{ + }^{2} - t_{ - }^{2}){\text{/}}(2{{\alpha }_{0}}),\quad {{t}_{ - }} < {{t}_{ + }} \leqslant 0. \hfill \\ \end{gathered} \right.$4. ОСНОВНОЙ РЕЗУЛЬТАТ
Сформулируем основной результат этой работы в виде теоремы.
Теорема. Пусть $\Gamma $ – простая гладкая замкнутая поверхность класса C2, ограничивающая объёмно-односвязную внутреннюю область, либо простая гладкая ограниченная разомкнутая поверхность класса C2, содержащая свои предельные точки. Пусть Γ допускает параметризацию (1) со свойствами (3), (4), и $\mu (y) \in {{C}^{0}}(\Gamma )$. Тогда для потенциала простого слоя $(5)$ при $x \notin \Gamma $ имеет место квадратурная формула
(17)
${{\mathcal{V}}_{k}}[\mu ](x) \approx \frac{1}{{4\pi }}\sum\limits_{n = 0}^{N - 1} {\kern 1pt} {\kern 1pt} \sum\limits_{m = 0}^{M - 1} \,{{\mu }_{{nm}}}\exp (ik{\kern 1pt} {\text{|}}x - y({{u}_{n}},{{v}_{m}}){\kern 1pt} {\text{|}}{\kern 1pt} ){{\Theta }_{{nm}}}(x),\quad k \geqslant 0,$Как следует из приведенных ниже численных тестов, остаточный член в (17) можно оценить как $O({{H}^{2}})$ равномерно по $x \notin \Gamma $, т.е. формула (17) дает равномерную аппроксимацию и обеспечивает равномерную сходимость к потенциалу простого слоя для точек, расположенных вне $\Gamma $.
Если $k = 0$, то потенциал простого слоя для уравнения Гельмгольца переходит в потенциал простого слоя для уравнения Лапласа, соответственно, квадратурная формула (17) при $k = 0$ принимает вид квадратурной формулы для гармонического потенциала простого слоя.
5. ЧИСЛЕННЫЕ ТЕСТЫ
Тестирование квадратурных формул для потенциалов простого слоя для уравнений Лапласа и Гельмгольца проведено в случае, когда поверхность $\Gamma $ явлется сферой единичного радиуса, которая задана параметрически уравнениями
(18)
${{y}_{1}}(u,v) = \sin v\cos u,\quad {{y}_{2}}(u,v) = \sin v\sin u,\quad {{y}_{3}}(u,v) = \cos v,$В рассматриваемых тестовых примерах для потенциала простого слоя с заданной на единичной сфере плотностью известно явное выражение во всем пространстве, поэтому точные значения потенциала можно сравнить с приближенными, вычисленными по квадратурным формулам. Во всех тестах приближенное значение потенциала простого слоя вычислялось по квадратурной формуле (17) в некоторых точках на вспомогательных сферах, имеющих центры в начале координат и радиусы, равные $1 \pm \Delta R$. Тем самым вспомогательные сферы находятся либо внутри, либо снаружи сферы единичного радиуса, на которой задана плотность потенциала, на расстоянии $\Delta R$ от нее. Затем были рассчитаны значения абсолютных погрешностей в этих точках либо относительных погрешностей (когда абсолютная погрешность делится на модуль точного значения потенциала в данной точке), и для каждой вспомогательной сферы определялись максимумы значений этих погрешностей. Найденные погрешности сравниваются с погрешностями, полученными для стандартной и улучшенной квадратурных формул из [5], [6].
Координаты точек, которые использовались для оценки максимальной абсолютной либо относительной погрешности
(19)
$\begin{gathered} x_{j}^{{ql}} = R{{y}_{j}}({{u}_{q}},{{v}_{l}}),\quad j = 1,2,3, \\ {{u}_{q}} = \frac{{2\pi }}{{2N}}q,\quad q = 0,1,2;\quad {{v}_{l}} = \frac{\pi }{{2M}}l,\quad l = 0, \ldots ,2M, \\ \end{gathered} $Вычисления проводились для различных значений $M$ и $N$. Значения шагов определяются формулами $h = 2\pi {\text{/}}N,$ $H = \pi {\text{/}}M$. Если $N = M = 25$, то $h \approx 0.25$, $H \approx 0.13$; если $N = M = 50$, то $h \approx 0.126$, $H \approx 0.063$; если $N = M = 100$, то $h \approx 0.063$, $H \approx 0.031$.
В таблицах приведены рассчитанные максимальные значения погрешностей. В левом столбце указано отличие радиуса вспомогательной сферы от единицы: для внутренних сфер радиус равен $1 - \Delta R$, для внешних сфер он равен $1 + \Delta R$. В верхней строке указаны значения $M,\;N$. Первое число в ячейках таблиц – максимальная погрешность для стандартной квадратурной формулы из [5], [6] на данной вспомогательной сфере, второе число после точки с запятой – максимальная погрешность на данной сфере для улучшенной квадратурной формулы из [5], [6], а третье число – максимальная погрешность для квадратурной формулы (17).
Тест 1 для квадратурной формулы в случае уравнения Лапласа. В данном тесте использовалась плотность потенциала $\mu (y(u,v)) = 4\pi ,$ тогда гармонический потенциал простого слоя имеет вид
Таблица 1.
$\Delta R$ | $M = N = 25$ | $M = N = 50$ | $M = N = 100$ |
---|---|---|---|
Внутренние сферы | |||
0.1 | 4.0E–3; 4.0E–3; 0.0019 | 9.0E–4; 9.7E–4; 0.00044 | 2.2E–4; 2.4E–4; 0.00011 |
0.01 | 0.21; 0.040; 0.0035 | 0.042; 0.014; 0.00076 | 7.1E–3; 4.2E–3; 0.00015 |
0.001 | 2.5; 0.065; 0.0042 | 0.60; 0.038; 0.0011 | 0.14; 0.017; 0.00026 |
0.0001 | 25; 0.069; 0.0043 | 6.3; 0.045; 0.0012 | 1.56; 0.027; 0.0003 |
Внешние сферы | |||
0.1 | 4.9E–3; 5.8E–3; 0.0015 | 9.9E–4; 1.4E–3; 0.00035 | 2.4E–4; 3.4E–4; 8.6E–5 |
0.01 | 0.21; 0.048; 0.003 | 0.043; 0.016; 0.00073 | 7.2E–3; 4.3E–3; 0.00014 |
0.001 | 2.5; 0.075; 0.0035 | 0.60; 0.046; 0.00045 | 0.15; 0.019; 0.00026 |
0.0001 | 25; 0.070; 0.0043 | 6.3; 0.047; 0.0011 | 1.6; 0.031; 0.00021 |
Тест 2 для квадратурной формулы в случае уравнения Лапласа. В данном тесте использовалась плотность потенциала $\mu (y(u,v)) = \cos u\sin v$, при этом гармонический потенциал простого слоя имеет вид
Таблица 2.
$\Delta R$ | $M = N = 25$ | $M = N = 50$ | $M = N = 100$ |
---|---|---|---|
Внутренние сферы | |||
0.1 | 3.4E–3; 2.7E–3; 0.0024 | 1.1E–4; 6.4E–4; 5.7E–4 | 2.8E–7; 1.6E–4; 1.4E–4 |
0.01 | 0.21; 4.7E–3; 0.0044 | 0.042; 1.0E–3; 9.7E–4 | 7.1E–3; 2.2E–4; 2.0E–4 |
0.001 | 2.5; 5.5E–3; 0.0051 | 0.60; 1.4E–3; 0.0013 | 0.14; 3.4E–4; 3.1E–4 |
0.0001 | 25; 5.6E–3; 0.0052 | 6.2; 1.5E–3; 0.0014 | 1.6; 3.7E–4; 3.5E–4 |
Внешние сферы | |||
0.1 | 4.1E–3; 3.9E–4; 0.00066 | 1.6E–4; 8.5E–5; 1.4E–4 | 6.9E–7; 2.1E–5; 3.4E–5 |
0.01 | 0.21; 1.7E–3; 0.002 | 0.042; 4.2E–4; 5.1E–4 | 7.1E–3; 6.8E–5; 8.9E–5 |
0.001 | 2.5; 4.7E–3; 0.0044 | 0.60; 5.5E–4; 4.7E–4 | 0.14; 1.8E–4; 2.1E–4 |
0.0001 | 25; 5.5E–3; 0.0052 | 6.2; 1.4E–3; 0.0013 | 1.6; 2.8E–4; 2.6E–4 |
Тест 3 для квадратурной формулы в случае уравнения Гельмгольца. В данном тесте использовалась плотность потенциала $\mu (y(u,v)) = k.$ Тогда потенциал простого слоя для уравнения Гельмгольца имеет вид
Таблица 3.
$\Delta R$ | $M = N = 25$ | $M = N = 50$ | $M = N = 100$ |
---|---|---|---|
Внутренние сферы | |||
0.1 | 4.3E–3; 4.8E–3; 0.0029 | 9.7E–4; 1.2E–3; 6.7E–4 | 2.4E–4; 2.9E–4; 1.7E–4 |
0.01 | 0.24; 0.047; 0.0069 | 0.050; 0.017; 0.0017 | 8.4E–3; 4.9E–3; 3.8E–4 |
0.001 | 2.9; 0.077; 0.008 | 0.72; 0.045; 0.0022 | 0.17; 0.021; 5.7E–4 |
0.0001 | 30; 0.082; 0.0081 | 7.4; 0.053; 0.0023 | 1.9; 0.032; 6.1E–4 |
Внешние сферы | |||
0.1 | 5.4E–3; 6.9E–3; 0.0027 | 1.1E–3; 1.7E–3; 5.6E–4 | 2.7E–4; 4.1E–4; 1.4E–4 |
0.01 | 0.25; 0.057; 0.0063 | 0.051; 0.018; 0.0015 | 8.6E–3; 5.1E–3; 3.4E–4 |
0.001 | 2.9; 0.089; 0.0075 | 0.72; 0.054; 0.0017 | 0.17; 0.023; 5.2E–4 |
0.0001 | 30; 0.083; 0.0081 | 7.4; 0.056; 0.0022 | 1.9; 0.037; 5.4E–4 |
Тест 4 для квадратурной формулы в случае уравнения Гельмгольца. В данном тесте использовалась плотность потенциала $\mu (y(u,v)) = {{k}^{3}}\cos v.$ При этом потенциал простого слоя имеет вид
Таблица 4.
$\Delta R$ | $M = N = 25$ | $M = N = 50$ | $M = N = 100$ |
---|---|---|---|
Внутренние сферы | |||
0.1 | 3.9E–3; 4.4E–3; 0.0013 | 8.6E–4; 1.1E–3; 3.1E–4 | 2.1E–4; 2.7E–4; 7.8E–5 |
0.01 | 0.097; 0.040; 0.0023 | 0.019; 0.014; 5.2E–4 | 2.8E–3; 4.2E–3; 1.1E–4 |
0.001 | 1.2; 0.065; 0.0032 | 0.30; 0.038; 8.0E–4 | 0.071; 0.017; 1.9E–4 |
0.0001 | 13; 0.069; 0.0033 | 3.1; 0.045; 8.9E–4 | 0.78; 0.027; 2.3E–4 |
Внешние сферы | |||
0.1 | 3.9E–3; 5.4E–3; 7.0E–4 | 8.8E–4; 1.3E–3; 1.5E–4 | 2.1E–4; 3.2E–4; 3.8E–5 |
0.01 | 0.098; 0.048; 0.0018 | 0.019; 0.015; 4.8E–4 | 2.9E–3; 4.3E–3; 1.0E–4 |
0.001 | 1.2; 0.074; 0.0028 | 0.30; 0.046; 5.6E–4 | 0.071; 0.019; 1.7E–4 |
0.0001 | 13; 0.070; 0.0033 | 3.1; 0.047; 8.5E–4 | 0.78; 0.031; 1.9E–4 |
Выводы. Как указано выше, первое число в ячейках таблиц – погрешность стандартной квадратурной формулы, второе число – погрешность улучшенной квадратурной формулы из [5], [6], третье число – погрешность квадратурной формулы (17). Из таблиц видно, что квадратурная формула (17) дает второй порядок (по $H$) равномерной сходимости к потенциалу простого слоя и обеспечивает равномерную аппроксимацию потенциала c погрешностью $O({{H}^{2}})$. Погрешность $O({{H}^{2}})$ не зависит от расстояния до $\Gamma $ и справедлива даже на очень малых расстояниях. Улучшенная квадратурная формула из [5], [6] дает первый порядок (по $H$) равномерной сходимости к потенциалу простого слоя и обеспечивает равномерную аппроксимацию потенциала c погрешностью $O(H)$. Тем самым формула (17) превосходит результат, полученный в [5], [6]. Обе эти формулы сохраняют свойство непрерывности потенциала простого слоя при переходе через поверхность $\Gamma $, что также следует из приведенных численных результатов. Из таблиц вытекает, что стандартная квадратурная формула (см. [5], [6]) не дает равномерой сходимости к потенциалу простого слоя и не обеспечивает равномерную аппроксимацию потенциала простого слоя, поскольку расходится вблизи $\Gamma $.
Отметим, что в тесте 2 формула (17) не дает существенного повышения точности по сравнению с улучшенной формулой из [5], [6]. Это объясняется тем, что максимальная погрешность формулы из [5], [6] достигается вблизи полюсов сферы, а в тесте 2 плотность потенциала обращается в нуль на полюсах, за счет чего погрешность формулы из [5], [6] уменьшается. Высокая погрешность формулы из [5], [6] вблизи полюсов сферы вызвана тем, что длина нормали обнуляется на полюсах, но меняется вблизи полюсов наиболее быстро, а формула из [5], [6] не учитывает производные длины нормали. В формуле (17) длина нормали учитывается более точно, чем в формуле из [5], [6], в частности, в (17) учтены производные длины нормали, поэтому в остальных тестах, кроме теста 2, формула (17) показывает более высокую точность.
6. ПРЯМОЕ ЗНАЧЕНИЕ ПОТЕНЦИАЛА ПРОСТОГО СЛОЯ НА ПОВЕРХНОСТИ
Используя полученные результаты, можно построить квадратурную формулу для прямого значения потенциала простого слоя, когда точка $x$ лежит на поверхности $\Gamma $ в одном из узлов. Пусть $x = y({{u}_{{\hat {n}}}},{{v}_{{\hat {m}}}})$, и $y({{u}_{{\hat {n}}}},{{v}_{{\hat {m}}}})$ – один из узлов на поверхности $\Gamma $. Если $(n,m) \ne (\hat {n},\hat {m})$ (т.е. $x = y({{u}_{{\hat {n}}}},{{v}_{{\hat {m}}}}) \ne y({{u}_{n}},{{v}_{m}})$), то интеграл (8) можно считать приближенно равным функции ${{\Theta }_{{nm}}}(x)$, которая найдена в разд. 3 в явном виде. Остается приближенно вычислить интеграл (8), когда $(n,m) = (\hat {n},\hat {m})$ (т.е. $y(u,v)$ принадлежит маленькому кусочку поверхности с центром в точке $x = y({{u}_{{\hat {n}}}},{{v}_{{\hat {m}}}})$). В этом случае, применяя формулу Тейлора с центром в точке $({{u}_{{\hat {n}}}},{{v}_{{\hat {m}}}})$, находим
Из приведенных рассуждений вытекает, что если точка $x$ лежит на поверхности $\Gamma $ в узле $y({{u}_{{\hat {n}}}},{{v}_{{\hat {m}}}})$, то квадратурная формула для прямого значения потенциала простого слоя имеет вид
(20)
$\begin{gathered} {{\left. {{{\mathcal{V}}_{k}}[\mu ](x)} \right|}_{{x = y({{u}_{{\hat {n}}}},{{v}_{{\hat {m}}}}) \in \Gamma }}} \approx \frac{1}{{4\pi }}{\kern 1pt} {{\mu }_{{\hat {n}\hat {m}}}}{\kern 1pt} {\text{|}}\eta (y({{u}_{{\hat {n}}}},{{v}_{{\hat {m}}}})){\kern 1pt} {\text{|}}{{\mathcal{I}}_{{\hat {n},\hat {m}}}} + \\ + \;\frac{1}{{4\pi }}\sum\limits_{\begin{subarray}{l} n = 0,{\kern 1pt} m = 0 \\ (n,m) \ne (\hat {n},\hat {m}) \end{subarray}} ^{n = N - 1,{\kern 1pt} m = M - 1} {{{\mu }_{{nm}}}\exp (ik{\kern 1pt} {\text{|}}y({{u}_{{\hat {n}}}},{{v}_{{\hat {m}}}}) - y({{u}_{n}},{{v}_{m}}){\kern 1pt} {\text{|}}{\kern 1pt} ){{\Theta }_{{nm}}}(y({{u}_{{\hat {n}}}},{{v}_{{\hat {m}}}}))} ,\quad k \geqslant 0. \\ \end{gathered} $Сходимость формулы (20) проверена на тестах 1–4 из разд. 5, полученные результаты приведены в табл. 5, где указана максимальная погрешность вычислений в узловых точках единичной сферы для каждого теста. В таблице по тестам 1 и 3 приводится максимальная относительная погрешность, а по тестам 2 и 4 – максимальная абсолютная погрешность. Для сравнения в табл. 5 приводятся результаты по квадратурной формуле, построенной в [5], [6]. Первое число в табл. 5 – максимальная погрешность квадратурной формулы из [5], [6], а второе число – максимальная погрешность квадратурной формулы (20). Из табл. 5 следует, что квадратурная формула (20) сходится и аппроксимирует прямое значение потенциала простого слоя с погрешностью $O({{H}^{2}})$. Кроме того, квадратурная формула (20) имеет второй порядок сходимости по $H$, тогда как квадратурная формула из [5], [6] дает сходимость первого порядка.
Список литературы
Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: ГИТТЛ, 1951. 659 с.
Гюнтер Н.М. Теория потенциала. М.: ГИТТЛ, 1953. 798 с.
Колтон Д., Кресс Р. Методы интегральных уравнений в теории рассеяния. М.: Мир, 1987. 311 с.
Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов. М.: Мир, 1987.
Крутицкий П.А., Федотова А.Д., Колыбасова В.В. Квадратурная формула для потенциала простого слоя // Дифференц. ур-ния. 2019. Т. 55. № 9. С. 1269–1284.
Крутицкий П.А., Федотова А.Д., Колыбасова В.В. О квадратурной формуле для потенциала простого слоя в трехмерном случае // Препринт ИПМ им. М.В. Келдыша РАН. 2019. № 112. 26 с.
Krutitskii P.A., Kwak D.Y., Hyon Y.K. Numerical treatment of a skew-derivative problem for the Laplace equation in the exterior of an open arc // J. Engineer. Math. 2007. V. 59. P. 25–60.
Крутицкий П.А., Колыбасова В.В. Численный метод решения интегральных уравнений в задаче с наклонной производной для уравнения Лапласа вне разомкнутых кривых // Дифференц. ур-ния. 2016. Т. 52. № 9. С. 1262–1276.
Krutitskii P.A. The 2-dimensional Dirichlet problem in an external domain with cuts // Zeitschrift fur Analysis und ihre Anwendungen. 1998. V. 17. № 2. P. 361–378.
Krutitskii P.A. The skew derivative problem in the exterior of open curves in a plane // Zeitschrift fur Analysis und ihre Anwendungen. 1997. V. 16. № 3. P. 739–747.
Krutitskii P.A. The Dirichlet problem for the dissipative Helmholtz equation in a plane domain bounded by closed and open curves // Hiroshima Math. J. 1998. V. 28. № 1. P. 149–168.
Krutitskii P.A. The Neumann problem for the 2-D Helmholtz equation in a domain, bounded by closed and open curves // Intern. J. Math. and Math. Sci. 1998. V. 21. № 2. P. 209–216.
Krutitskii P.A. On the electric current from electrodes in a magnetized semiconductor film // IMA J. Appl. Math. 1998. V. 60. P. 285–297.
Krutitskii P.A. The Neumann problem on wave propagation in a 2-D external domain with cuts // J. Math. Kyoto Univ. 1998. V. 38. № 3. P. 439–452.
Krutitskii P.A. Wave propagation in a 2-D external domain with cuts // Appl. Anal. 1996. V. 62. № 3–4. P. 297–309.
Krutitskii P.A. The Helmholtz equation in the exterior of slits in a plane with different impedance boundary conditions on opposite sides of the slits // Quarterly of App. Math. 2009. V. 67. № 1. P. 73–92.
Бутузов В.Ф., Крутицкая Н.Ч., Медведев Г.Н., Шишкин А.А. Математический анализ в вопросах и задачах. М.: Физматлит, 2000.
Ильин В.А., Позняк Э.Г. Основы математического анализа. Ч. 2. М.: Физматлит, 1973.
Владимиров В.С. Уравнения математической физики. М.: Физматлит, 1981.
Прудников А.П., Брычков Ю.А., Маричев О.И. Интегралы и ряды. М.: Физматлит, 1981.
Дополнительные материалы отсутствуют.
Инструменты
Журнал вычислительной математики и математической физики