ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ, 2022, том 58, № 8, с.1121-1131
ЧИСЛЕННЫЕ МЕТОДЫ
УДК 517.923+517.925.54
КВАДРАТУРНАЯ ФОРМУЛА
ДЛЯ ПОТЕНЦИАЛА ДВОЙНОГО СЛОЯ
С ДИФФЕРЕНЦИРУЕМОЙ ПЛОТНОСТЬЮ
© 2022 г. П. А. Крутицкий, И. О. Резниченко
Выводится квадратурная формула для потенциала двойного слоя с дифференцируемой
плотностью, заданной на гладкой замкнутой либо незамкнутой поверхности. Проведённые
численные тесты показывают, что эта формула даёт более высокую точность вычислений
вблизи поверхности, где задана плотность потенциала, чем квадратурные формулы, в ко-
торых дифференцируемость плотности не учитывается и плотность предполагается лишь
непрерывной. Преимущество квадратурной формулы в данной работе особенно заметно в
случае, когда плотность потенциала представлена гладкими осциллирующими функциями,
так как она позволяет повысить точность вычисления потенциала без увеличения стоимо-
сти вычислений.
DOI: 10.31857/S037406412208012X, EDN: CGZXBB
Введение. Потенциал двойного слоя для уравнений Лапласа и Гельмгольца использует-
ся при решении краевых задач математической физики методом интегральных уравнений [1,
раздел IV; 2, раздел II; 3, гл. 5]. Для вычисления потенциала двойного слоя с непрерыв-
ной плотностью в прикладных расчётах применяют стандартную квадратурную формулу [4,
гл. 2] либо формулу, основанную на определении телесного угла [5], однако они не позволя-
ют повысить точность вычислений при дифференцируемой плотности, так как не допускают
обобщений на этот случай. В работе [6] предложена новая квадратурная формула для гармо-
нического потенциала двойного слоя с непрерывной плотностью, которую можно обобщить на
случай дифференцируемой плотности, что позволит повысить точность вычислений потенци-
ала двойного слоя, если плотность в потенциале дифференцируема. Указанному обобщению
и посвящена настоящая работа. Улучшение точности вычислений потенциала подтверждается
численными расчётами.
В двумерном случае улучшенная квадратурная формула для потенциала простого слоя
с плотностью, заданной на незамкнутых кривых и имеющей степенные особенности на их
концах, построена в статьях [7, 8] и может применяться для нахождения численных решений
краевых задач для уравнений Лапласа и Гельмгольца вне разрезов и незамкнутых кривых на
плоскости. Такие задачи изучались в работах [9-15].
1. Постановка задачи. Введём декартову систему координат x = (x1, x2, x3) в простран-
стве R3. Пусть Γ - простая гладкая замкнутая поверхность класса C2 либо простая гладкая
ограниченная незамкнутая ориентированная поверхность класса C2, содержащая свои пре-
дельные точки [16, гл. 14, § 1]. Если поверхность Γ замкнутая, то она должна ограничивать
объёмно-односвязную внутреннюю область [17, c. 201]. Предположим, что поверхность Γ па-
раметризована так, что на неё отображается прямоугольник [0, A] × [0, B]:
y = (y1,y2,y3) ∈ Γ, y1 = y1(u,v), y2 = y2(u,v), y3 = y3(u,v); u ∈ [0,A], v ∈ [0,B];
yj(u,v) ∈ C2([0,A] × [0,B]), j = 1,2,3.
(1)
Сферу, поверхность эллипсоида, гладкие поверхности фигур вращения, поверхность тора и
многие другие более сложные поверхности можно параметризовать таким образом.
Введём N точек un с шагом h на отрезке [0, A] и M точек vm с шагом H на отрезке
[0, B] и рассмотрим разбиение прямоугольника [0, A] × [0, B]
A = Nh, B = MH, un = (n + 1/2)h, n = 0,N - 1; vm = (m + 1/2)H, m = 0,M - 1,
на N × M прямоугольников c центрами (un, vm).
8
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
1122
КРУТИЦКИЙ, РЕЗНИЧЕНКО
Известно [16, гл. 14, § 1], что компоненты вектора нормали (не единичного)
η(y) = (η1(y), η2(y), η3(y))
в точке поверхности y = (y1, y2, y3) ∈ Γ выражаются через определители второго порядка
формулами
(y2)u (y3)u
(y3)u (y1)u
(y1)u (y2)u
η1 =
η2 =
η3 =
(y2)v (y3)v,
(y3)v (y1)v,
(y1)v (y2)v,
√
а длина этого вектора равна |η(y)| =
(η1(y))2 + (η2(y))2 + (η3(y))2.
Известно [16, гл. 14, §§ 1, 2], что справедливо равенство
A
B
∫
∫
∫
F (y) dsy = du F (y(u, v))|η(y(u, v))| dv.
Γ
0
0
Потребуем, чтобы
|η(y(u, v))| > 0 для любых точек (u, v) ∈ ((0, A) × (0, B)).
(2)
Из условия (2) следует, что |η(y(u, v))| ∈ C1((0, A) × (0, B)).
Через ny обозначим единичную нормаль в точке y ∈ Γ, т.е. ny = η(y)/|η(y)|. Производная
по нормали ny имеет вид
∂
= |η(y)|-1(η(y), ∇y ).
∂ny
√
Обозначим |x - y(u, v)| =
(x1 - y1(u, v))2 + (x2 - y2(u, v))2 + (x3 - y3(u, v))2 и заметим, что
∑
∂
1
yj - xj
|x - y| =
ηj(y)
∂ny
|η(y)|
|x - y|
j=1
Пусть x ∈ Γ. Рассмотрим потенциал двойного слоя для уравнения Гельмгольца с заданной
на поверхности Γ дифференцируемой плотностью μ(y) ∈ C1(Γ):
∫
1
∂ eik|x-y|
Wk[μ](x) =
μ(y)
dsy =
4π
∂ny |x - y|
Γ
∫
∑
1
1
exp(ik|x - y|)(ik|x - y| - 1)
ηj(y)(yj - xj)
=
μ(y)
dsy =
4π
|η(y)|
|x - y|2
|x - y|
j=1
Γ
∫
A
∫
B
∑
1
ηj(y(u,v))(yj(u,v) - xj)
=
du μ(y(u, v)) exp(ik|x - y(u, v)|)(ik|x - y(u, v)| - 1)
dv =
4π
|x - y(u, v)|3
j=1
0
0
∫
∫
∑
∑
1
=
du
μ(y(u, v))×
4π
n=0 m=0
un-h/2
vm-H/2
∑
ηj(y(u,v))(yj(u,v) - xj)
× exp(ik|x - y(u, v)|)(ik|x - y(u, v)| - 1)
dv,
(3)
|x - y(u, v)|3
j=1
где k ≥ 0. При k = 0 потенциал (3) переходит в потенциал двойного слоя для уравнения
Лапласа.
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
КВАДРАТУРНАЯ ФОРМУЛА ДЛЯ ПОТЕНЦИАЛА ДВОЙНОГО СЛОЯ
1123
При u ∈ [un - h/2, un + h/2] и v ∈ [vm - H/2, vm + H/2] так же, как и в работе [18], можно
показать, что справедливы равенства
|x-y(u, v)| = |x-y(un, vm)|+O(h+H), exp(ik|x-y(u, v)|) = exp(ik|x-y(un, vm)|)+O(h+H)
для любого x ∈ Γ. Константы в оценках функций O(h + H) не зависят от n, m и от
расположения x ∈ Γ. Следовательно,
∑
∑
1
Wk[μ](x) ≈
exp(ik|x - y(un, vm)|)(ik|x - y(un, vm)| - 1) ×
4π
n=0 m=0
∫
∫
∑
ηj(y(u,v))(yj(u,v) - xj)
×
du
μ(y(u, v))
dv.
(4)
|x - y(u, v)|3
j=1
un-h/2
vm
−H/2
Таким образом, чтобы получить квадратурную формулу для потенциала двойного слоя при
x ∈ Γ необходимо вычислить следующий интеграл:
∫
∫
∑
ηj(y(u,v))(yj(u,v) - xj)
du
μ(y(u, v))
dv.
(5)
|x - y(u, v)|3
j=1
un-h/2
vm
−H/2
2. Вычисление интеграла. Пусть точка x не принадлежит части поверхности Γ, на
которой изменяется точка y = y(u, v) при (u - un) ∈ [-h/2, h/2] и (v - vm) ∈ [-H/2, H/2].
Разложим yj(u, v) по формуле Тейлора с центром в точке (un, vm), тогда для j = 1, 2, 3
получим
yj(u,v) = yj(un,vm) + Dj + O(H2 + h2),
где
Dj = (yj)′u(u - un) + (yj)′v(v - vm).
Здесь и далее все производные по u и v берутся в точке (un, vm). Положим
∑
r2 = |x - y(un,vm)|2 =
r2j = 0, rj = yj(un,vm) - xj, j = 1,2,3,
j=1
тогда
yj(u,v) - xj = rj + Dj + O(H2 + h2), j = 1,2,3.
Следовательно,
∑
∑
|x - y(u, v)|2 =
(xj - yj(u, v))2 ≈
(r2j + 2rj Dj + D2j) =
j=1
j=1
= r2 + 2P(u - un) + 2Q(v - vm) + α2(u - un)2 + β2(v - vm)2 + 2δ(u - un)(v - vm) =
= β2(V + δU/β2 + Q/β2)2 - (δU + Q)2/β2 + α2U2 + 2PU + r2,
где
∑
∑
U = u - un, V = v - vm, P = rj(yj)′u, Q = rj(yj)′v,
j=1
j=1
∑
∑
∑
α2 =
((yj )′u)2, β2 =
((yj )′v)2, δ =
(yj )′u(yj)′v.
j=1
j=1
j=1
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
8∗
1124
КРУТИЦКИЙ, РЕЗНИЧЕНКО
Можно показать [16, гл. 14, § 1], что справедливо равенство
α2β2 - δ2 = |η(y(un,vm))|2.
Согласно условию (2) |η(y(un, vm))| > 0 для всех возможных n, m, поэтому
α2β2 - δ2 > 0,
откуда следует, что α2 > 0 и β2 > 0.
Применив формулу Тейлора в точке (un, vm) с остаточным членом в форме Пеано [16,
гл. 10, § 5.3], находим
√
ηj(y(u,v)) = ηj(y(un,vm)) + (ηj)′u(u - un) + (ηj)′v(v - vm) + o(
(u - un)2 + (v - vm)2).
Поскольку μ(y(u, v)) ∈ C1(Γ), аналогично получим
√
μ(y(u, v)) = μ(y(un, vm)) + μ′u(u - un) + μ′v(v - vm) + o(
(u - un)2 + (v - vm)2),
где u ∈ [un - h/2, un + h/2] и v ∈ [vm - H/2, vm + H/2].
∑3
Для вычисления выражения
ηj(y(u,v))(yj(u,v) - xj) с учётом формул
j=1
∑
∑
ηj(y(un,vm))(yj)′u =
ηj(y(un,vm))(yj)′v = 0,
j=1
j=1
отражающих ортогональность вектора нормали и касательных векторов к поверхности (см.
[16, гл. 14, § 1.2]), воспользуемся разложением по формуле Тейлора в точке (un, vm) с оста-
точным членом в форме Пеано
yj(u,v) - xj = rj + (yj)′u(u - un) + (yj)′v(v - vm) +
1
1
+
(yj )′′uu(u - un)2 +
(yj)′′vv(v - vm)2 + (yj)′′uv(u - un)(v - vm) + o((u - un)2 + (v - vm)2),
2
2
тогда
∑
μ(y(u, v))
ηj(y(u,v))(yj(u,v) - xj) ≈ R + ξ4U + ξ5V + ξ1U2 + ξ2V2 + ξ3UV,
j=1
где U = u - un, V = v - vm,
(
)
)
∑
(1
ξ1 =
μ(y(un, vm))
ηj(y(un,vm))(yj)′′uu + (ηj)′u(yj)′
+ μ′u(ηj)′
rj
,
u
u
2
j=1
(
)
)
∑
(1
ξ2 =
μ(y(un, vm))
ηj(y(un,vm))(yj)′′vv + (ηj)′v(yj)′
+ μ′v(ηj)′
rj
,
v
v
2
j=1
∑
ξ3 = (μ(y(un,vm))(ηj(y(un,vm))(yj)′′uv + (ηj)′u(yj)′v + (ηj)′v(yj)′u) + μ′u(ηj)′vrj + μ′v(ηj)′urj),
j=1
∑
ξ4 = (μ(y(un,vm))(ηj)′urj + μ′uηj(y(un,vm))rj),
j=1
∑
∑
ξ5 = (μ(y(un,vm))(ηj)′vrj + μ′vηj(y(un,vm))rj), R = μ(y(un,vm))
ηj(y(un,vm))rj.
j=1
j=1
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
КВАДРАТУРНАЯ ФОРМУЛА ДЛЯ ПОТЕНЦИАЛА ДВОЙНОГО СЛОЯ
1125
Из приведённых соотношений вытекает, что интеграл (5) приближённо равен следующему
интегралу, который обозначим через Knm(x):
∫
∫
∑
μ(y(u, v))
du
ηj(y(u,v))(yj(u,v) - xj)dv ≈
|x - y(u, v)|3
j=1
un-h/2
vm-H/2
∫
∫
R+ξ4U+ξ5V +ξ1U2 +ξ2V2 +ξ3UV
≈ dU
dV =
β3((V + δU/β2 + Q/β2)2 - (δU + Q)2/β4 + (α2U2 + 2PU + r2)/β2)3/2
-h/2
-H/2
= Knm(x).
(6)
Следовательно, чтобы вывести квадратурную формулу для потенциала двойного слоя,
необходимо вычислить интеграл Knm(x) в явном виде. Этот интеграл вычислен в работе [6].
3. Основной результат.
Теорема. Пусть Γ - простая гладкая замкнутая поверхность класса C2, ограничиваю-
щая объёмно-односвязную внутреннюю область, либо простая гладкая ограниченная неза-
мкнутая ориентированная поверхность класса C2, содержащая свои предельные точки.
Пусть Γ допускает параметризацию (1) со свойством (2) и μ(y) ∈ C1(Γ). Тогда для по-
тенциала двойного слоя (3) при x ∈ Γ и k ≥ 0 имеет место квадратурная формула
∑
∑
1
Wk[μ](x) =
exp(ik|x - y(un, vm)|)(ik|x - y(un, vm)| - 1)Knm(x) + εNM (x),
(7)
4π
n=0 m=0
где εNM (x) → 0 при N, M → ∞. Интеграл Knm(x), записанный в (6), вычислен в явном
виде в статье [6], а коэффициенты в Knm(x) берутся из п. 2.
При k = 0 потенциал двойного слоя для уравнения Гельмгольца переходит в потенциал
двойного слоя для уравнения Лапласа, а квадратурная формула (7) принимает вид квадра-
турной формулы для гармонического потенциала двойного слоя.
Если в формуле (7) при вычислении коэффициентов в Knm(x) в п.2 положить μ′u = μ′v = 0,
то она при k = 0 переходит в квадратурную формулу для потенциала двойного слоя из [6], где
плотность μ(y) не предполагается дифференцируемой, а считается всего лишь непрерывной.
Тем самым формула (7) обобщает прежнюю квадратурную формулу для потенциала двойно-
го слоя с непрерывной плотностью на случай дифференцируемой плотности. Ниже эффек-
тивность формулы (7) оценивается на тестовых примерах и сравнивается с эффективностью
других формул.
4. Квадратурная формула, основанная на свойствах телесного угла. Рассмотрим
треугольник с вершинами x1, x2, x3. Известно [19, с. 253], что телесный угол Ω(x; x1, x2, x3),
под которым этот треугольник виден из точки x, с точностью до знака, зависящего от вы-
бранного на треугольнике направления нормали, равен определённому на этом треугольнике
потенциалу двойного слоя с единичной плотностью, умноженному на 4π. С другой стороны,
для телесного угла, под которым треугольник с вершинами x1, x2, x3 виден из точки x,
имеется явная формула [20; 21, результат (359)]:
(√
Ω(x; x1, x2, x3) = 2 arccos
|R1||R2||R3|/2 ×
|R1||R2||R3| + (R1, R2)|R3| + (R2, R3)|R1| + (R3, R1)|R2|
×
√
×
|R1||R2||R3| + (R1, R2)|R3|
)
1
1
×
√
√
,
|R1||R2||R3| + (R2, R3)|R1|
|R1||R2||R3| + (R3, R1)|R2|
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
1126
КРУТИЦКИЙ, РЕЗНИЧЕНКО
где через R1 = x1 - x, R2 = x2 - x, R3 = x3 - x обозначены векторы, а через ( · , · ) -
скалярное произведение. Использовав указанную формулу для телесного угла и параметриза-
цию (1), можно вывести следующую квадратурную формулу для потенциала двойного слоя с
плотностью μ(y) на поверхности Γ:
∑
∑
1
Wk[μ](x) ≈
μ(y(un, vm)) exp(ik|x - y(un, vm)|)(ik|x - y(un, vm)| - 1) ×
4π
n=0 m=0
× [Ω(x; y(un - h/2, vm - H/2), y(un - h/2, vm + H/2), y(un + h/2, vm - H/2)) +
+ Ω(x;y(un + h/2,vm + H/2),y(un - h/2,vm + H/2),y(un + h/2,vm - H/2))] ×
)
(3∑
× sgn
ηj(y(un,vm))(yj(un,vm) - xj)
(8)
j=1
Квадратурная формула (8) используется в инженерных расчётах, однако она не допускает
улучшения за счёт учёта производных в случае дифференцируемой плотности в потенциале
двойного слоя.
5. Численные расчёты для диска. В данных расчётах оценивается эффективность
квадратурных формул вблизи поверхности Γ, когда поверхность Γ является круговым дис-
ком единичного радиуса, расположенным в плоскости x3 = 0 с центром в начале координат.
Тем самым Γ является незамкнутой поверхностью, заданной уравнениями
y1(u,v) = ucos v, y2(u,v) = usin v, y3(u,v) = 0,
(9)
причём (u, v) ∈ [0, 1] × [0, 2π]. Отметим, что в этом случае |η(y(u, v))| = u и направление
нормали η(y) на диске Γ совпадает с направлением оси Ox3. Кроме того, |η(y(0, v))| = 0 для
всех v ∈ [0, 2π], иначе говоря, |η(y)| = 0 в центре диска при такой параметризации. Согласно
[22, § 27.6] плотность потенциала двойного слоя на диске Γ можно найти по формуле
μ(x)|Γ = 2Wk[μ](x)|Γ- .
(10)
Здесь диск Γ рассматривается как двусторонняя поверхность, через Γ- обозначена сторона,
которую мы видим, смотря навстречу вектору нормали ny. Направление единичной нормали
ny совпадает с направлением нормали η, так как вектор ny получается из η в результате
нормировки. В формуле (10) берётся предельное значение потенциала двойного слоя на верх-
ней стороне Γ, т.е. на Γ-. Отметим, что потенциал двойного слоя Wk[μ](x) равен нулю в
плоскости x3 = 0 вне диска Γ.
Эффективность квадратурных формул вблизи поверхности Γ можно оценивать следую-
щим образом. В тестовых примерах известно явное выражение плотности потенциала μ(y) на
диске Γ. Для точек x, расположенных над диском Γ, значение потенциала Wk[μ](x1, x2, x3)
при уменьшении x3 > 0 должно сходиться к значению (1/2)μ(x1, x2) согласно формуле (10)
в силу непрерывности потенциала вплоть до границы. Поэтому приближённые значения по-
тенциала, вычисленные по квадратурным формулам, также должны стремиться к известному
значению (1/2)μ(x1, x2). Однако для совсем малых величин x3 > 0 значения, вычисленные по
квадратурным формулам, перестают приближаться к (1/2)μ(x1, x2) из-за дискретности квад-
ратурных формул. Тесты позволяют оценить расстояния до границы, на которых квадратур-
ные формулы приближаются к предельным значениям потенциала, и расстояние, с которого
они перестают приближаться. Кроме того, тесты показывают, насколько хорошо квадратурные
формулы аппроксимируют предельные значения потенциала вблизи границы.
Во всех тестах приближённое значение потенциала двойного слоя вычислялось по преж-
ней квадратурной формуле из работы [6], по новой квадратурной формуле (7) и по формуле
(8) в некоторых точках на единичных вспомогательных дисках, имеющих центры с коорди-
натами (0, 0, x3) и расположенных параллельно плоскости x3 = 0 над диском Γ. Затем в
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
КВАДРАТУРНАЯ ФОРМУЛА ДЛЯ ПОТЕНЦИАЛА ДВОЙНОГО СЛОЯ
1127
каждой такой точке (x1, x2, x3) на вспомогательном диске был найден модуль разности меж-
ду значением квадратурной формулы в данной точке и явным значением (1/2)μ(x1, x2) при
x3 = 0, x3 - расстояние от вспомогательного диска до диска Γ. На каждом вспомогательном
диске был найден максимум модуля таких разностей по всем указанным точкам для каждой
формулы. Результаты представлены в таблицах.
Координаты точек, которые использовались для оценки максимального модуля разности:
xqlj = yj(uq,vl), j = 1,2,
1
2π
uq =
q, q = 0,2N; vl =
l,
l = 0,1,2,
(11)
2N
2M
где yj(u, v) определяются формулами (9), т.е. эти точки расположены над центрами участ-
ков разбиения диска Γ, серединами границ между такими участками и пересечениями этих
границ. Отметим, что данные точки распределены не по всему вспомогательному диску, а
находятся вблизи радиуса с полярным углом v = H/2.
Вычисления проводились для различных значений M и N. Значения шагов определяются
формулами
h = 1/N, H = 2π/M.
Если N = M/5 = 10, то h ≈ 0.1, H ≈ 0.13; если N = M/5 = 20, то h ≈ 0.05, H ≈ 0.063;
если N = M/5 = 40, то h ≈ 0.025, H ≈ 0.03.
В таблицах для каждого вспомогательного диска с координатой x3 приведены рассчитан-
ные максимальные абсолютные отклонения значений квадратурных формул от явных пре-
дельных значений потенциала двойного слоя на Γ по точкам (11). Первое число в ячейках
таблицы - максимальное отклонение для новой формулы (7) на данном вспомогательном дис-
ке, второе число - максимальное отклонение для прежней квадратурной формулы из статьи [6]
на данном диске, третье число - максимальное отклонение для формулы (8).
В численных тестах квадратурные формулы сравниваются при k = 0, т.е. для потенциала
двойного слоя в случае уравнения Лапласа.
Тест 1. В данном тесте использовалась плотность потенциала
μ(y(u, v)) = 1 - u2.
При этом потенциал двойного слоя непрерывен на кромке диска Γ, но его производные терпят
разрыв, хотя и ограничены. Потенциал в рассматриваемом случае не зависит от полярного
угла v, как и его плотность. В табл. 1 приведены рассчитанные максимальные абсолютные
отклонения квадратурных формул от предельных значений потенциала.
Таблица 1. Максимальные абсолютные отклонения квадратурных формул
в тесте 1
x3
M/5 = N = 10
M/5 = N = 20
M/5 = N = 40
0.03
0.029; 0.034; 0.033
0.029; 0.029; 0.029
0.029; 0.029; 0.029
0.01
0.017; 0.028; 0.22
0.012; 0.018; 0.017
0.012; 0.014; 0.014
0.008
0.019; 0.028; 0.41
0.01; 0.017; 0.016
0.01; 0.012; 0.012
0.006
0.026; 0.031; 0.47
0.0082; 0.016; 0.015
0.0081; 0.011; 0.011
0.004
0.044; 0.051; 0.49
0.0069; 0.015; 0.022
0.0059; 0.0092; 0.0092
0.002
0.1; 0.11; 0.5
0.014; 0.016; 0.49
0.0034 0.0078; 0.0077
Тест 2. В данном тесте плотность потенциала равна μ(y(u, v)) = sin(3πu). При этом по-
тенциал двойного слоя непрерывен на кромке диска Γ, но его производные терпят разрыв, как
и в предыдущем тесте. Плотность потенциала не зависит от полярного угла v, поэтому и сам
потенциал от него не зависит. В табл. 2 приведены рассчитанные максимальные абсолютные
отклонения квадратурных формул от предельных значений потенциала.
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
1128
КРУТИЦКИЙ, РЕЗНИЧЕНКО
Таблица 2. Максимальные абсолютные отклонения квадратурных формул
в тесте 2
x3
M/5 = N = 10
M/5 = N = 20
M/5 = N = 40
0.03
0.11; 0.15; 0.15
0.12; 0.13; 0.13
0.15; 0.15; 0.15
0.01
0.047; 0.12; 0.53
0.05; 0.069; 0.068
0.074; 0.073; 0.073
0.008
0.042; 0.12; 0.53
0.041; 0.067; 0.066
0.062; 0.06; 0.06
0.006
0.036; 0.12; 0.55
0.032; 0.065; 0.063
0.049; 0.047; 0.047
0.004
0.057; 0.13; 0.54
0.023; 0.063; 0.11
0.034; 0.038; 0.037
0.002
0.15; 0.21; 0.52
0.023; 0.062; 0.51
0.018; 0.034; 0.033
Тест 3. В этом тесте использовалась плотность потенциала
μ(y(u, v)) = (u(1 - u))2 cos(10v),
первые производные которого являются непрерывными функциями на кромке диска Γ. В
отличие от предыдущих тестов здесь в формуле (11) вместо l = 0, 1, 2 взяты значения l =
= 0, 2M для учёта зависимости потенциала двойного слоя от полярного угла v. В первых
двух тестах плотность потенциала двойного слоя от v не зависит, поэтому и сам потенциал
не зависит от v.
В табл. 3 приведены данные расчётов максимальных абсолютных отклонений квадратур-
ных формул от предельных значений потенциала.
Таблица 3. Максимальные абсолютные отклонения квадратурных формул
в тесте 3
x3
M/5 = N = 10
M/5 = N = 20
M/5 = N = 40
0.03
0.014; 0.017; 0.017
0.015; 0.015; 0.015
0.015; 0.015; 0.015
0.01
0.0056; 0.01; 0.02
0.0061; 0.007; 0.007
0.0064; 0.0065; 0.0065
0.008
0.0045; 0.01; 0.036
0.0049; 0.006; 0.0059
0.0052; 0.0054; 0.0054
0.006
0.0034; 0.01; 0.034
0.0037; 0.005; 0.0049
0.004; 0.0042; 0.0042
0.004
0.0027; 0.011; 0.034
0.0024; 0.004; 0.0038
0.0027; 0.003; 0.0029
0.002
0.0071; 0.014; 0.032
0.0014; 0.0037; 0.031
0.0014; 0.0017; 0.0017
6. Численные расчёты для сферы. В тестах для сферы оценивается эффективность
квадратурных формул вблизи поверхности Γ, которая является сферой единичного радиуса
и задана параметрически уравнениями
y1(u,v) = cos usin v, y2(u,v) = sin usin v, y3(u,v) = cos v,
(12)
причём (u, v) ∈ [0, 2π]×[0, π]. Отметим, что в данном случае |η(y(u, v))| = sin v и |η(y(u, 0))| =
= |η(y(u, π))| = 0 для всех u ∈ [0, 2π]. Иначе говоря, |η(y)| = 0 на полюсах сферы при такой
параметризации, но условия теоремы выполняются.
В тестах для диска явные значения для потенциала двойного слоя были известны только
на поверхности диска. В отличие от тестов для диска в тестах для сферы известно явное вы-
ражение потенциала двойного слоя во всех точках внутри единичной сферы, поэтому точные
значения потенциала можно сравнить с приближёнными, вычисленными по квадратурным
формулам. Во всех тестах приближённое значение потенциала двойного слоя вычислялось
по новой квадратурной формуле (7), в которой учитывается, что плотность μ(y(u, v)) диф-
ференцируема на Γ, по прежней квадратурной формуле из работы [6], в которой плотность
μ(y(u, v)) считается всего лишь непрерывной на Γ и её дифференцируемость не учитывается,
а также по формуле (8). Вычисления проводились в некоторых точках на вспомогательных
сферах, имеющих центры в начале координат и радиусы, равные 1 - ΔR. Тем самым вспо-
могательные сферы находятся внутри сферы единичного радиуса, на которой задана плот-
ность потенциала на расстоянии ΔR от неё. Затем были рассчитаны значения абсолютных
погрешностей в этих точках, и для каждой вспомогательной сферы определялись максимумы
значений этих погрешностей.
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
КВАДРАТУРНАЯ ФОРМУЛА ДЛЯ ПОТЕНЦИАЛА ДВОЙНОГО СЛОЯ
1129
Координаты точек, которые использовались для оценки максимальной абсолютной погреш-
ности, равны
xqlj = Ryj(uq,vl), j = 1,2,3,
2π
π
uq =
q, q = 0,2N; vl =
l,
l = 1,2M - 1,
2N
2M
где yj (u, v) определяются формулами (12), R - радиус вспомогательной сферы, т.е. эти точки
расположены под центрами участков разбиения единичной сферы, серединами границ между
такими участками и пересечениями этих границ. Отметим, что эти точки распределены по
всей сфере.
Вычисления проводились для различных значений M и N. Значения шагов определяются
формулами h = 2π/N, H = π/M. Если N/2 = M = 16, то h = H ≈ 0.2; если N/2 = M = 32,
то h = H ≈ 0.098; если N/2 = M = 64, то h = H ≈ 0.049.
В таблицах приведены рассчитанные максимальные значения абсолютных погрешностей.
В первом столбце указано отличие радиуса на ΔR вспомогательной сферы от единицы, тем
самым её радиус будет равен 1 - ΔR. В верхней строке указаны значения M, N. Первое
число в ячейках таблицы - максимальная погрешность для новой формулы (7) на данной
вспомогательной сфере, второе число - максимальная погрешность для прежней квадратурной
формулы из [6] на данной сфере, третье число - максимальная погрешность для формулы (8).
В численных тестах квадратурные формулы сравниваются при k = 0, т.е. для потенциала
двойного слоя в случае уравнения Лапласа.
Тест 4. В данном тесте использовалась плотность потенциала μ(y(u, v)) = P10(cos v), тогда
гармонический потенциал двойного слоя внутри единичной сферы Γ имеет вид
10
11|x|
W0[μ](x) =
P10(cos ϑ),
|x| < 1,
21
где
1
P10(cos v) =
(46189 cos10 v - 109395 cos8 v + 90090 cos6 v - 30030 cos4 v + 3465 cos2 v - 63)
256
– полином Лежандра десятой степени от cos v, ϑ - зенитный угол в сферических координа-
тах с центром в начале координат. В табл. 4 приведены максимальные значения абсолютных
погрешностей.
Таблица 4. Максимальная абсолютная погрешность квадратурных формул в
тесте 4
ΔR N/2 = M = 16
N/2 = M = 32
N/2 = M = 64
0.2
0.0077; 0.01; 0.01
0.0022; 0.0026; 0.0024
0.00056; 0.00065; 0.00059
0.1
0.012; 0.025; 0.024
0.0036; 0.0053; 0.0045
0.00094; 0.0014; 0.0011
0.06
0.013; 0.045; 0.041
0.0043; 0.0084; 0.0067
0.0012; 0.0023; 0.0015
0.03
0.019; 0.074; 0.063
0.0048; 0.017; 0.011
0.0018; 0.0046; 0.0022
Тест 5. В этом тесте использовалась плотность потенциала
μ(y(u, v)) = cos(10u) sin10 v.
Гармонический потенциал двойного слоя внутри единичной сферы Γ в данном случае име-
ет вид
11|x|10 cos(10ϕ) sin10 ϑ
W0[μ](x) =
,
|x| < 1,
21
где ϑ и ϕ - зенитный и азимутальный углы соответственно в сферических координатах с
центром в начале координат. В табл. 5 приведены рассчитанные максимальные значения аб-
солютных погрешностей.
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
1130
КРУТИЦКИЙ, РЕЗНИЧЕНКО
Таблица 5. Максимальная абсолютная погрешность квадратурных формул
в тесте 5
ΔR
N/2 = M = 16
N/2 = M = 32
N/2 = M = 64
0.2
0.0046; 0.013; 0.0093
0.0016; 0.0031; 0.002
0.00043; 0.00078; 0.0005
0.1
0.012; 0.058; 0.046
0.0044; 0.011; 0.0068
0.0012; 0.0027; 0.0016
0.06
0.013; 0.11; 0.09
0.0055; 0.021; 0.013
0.0015; 0.0046; 0.0026
0.03
0.061; 0.2; 0.15
0.0047; 0.046; 0.026
0.0011; 0.0085; 0.0042
Заключение. Во всех тестах для диска максимальные отклонения квадратурной фор-
мулы (7) от предельных значений потенциала существенно меньше, чем отклонения у двух
остальных формул. В тестах 2 и 3, когда плотность μ(y(u, v)) задана осциллирующими функ-
циями, на расстояниях от диска x3 ≈ H2 отклонения формулы (7) от предельных значений
потенциала в 2-3 раза меньше, чем у остальных формул. В тесте 1 на таких же расстояниях
от диска отклонения формулы (7) от предельных значений потенциала в 1.5-2 раза меньше,
чем у двух других формул.
В тестах для сферы из табл. 4 и 5 следует, что квадратурная формула для потенциала
двойного слоя из [6], как и квадратурная формула с телесным углом (8), на расстоянии H от
сферы дают погрешность O(H2), а на расстоянии H2 от сферы - погрешность O(H).
Как видно из табл. 4 и 5, в тестах для сферы погрешность новой квадратурной формулы
(7) значительно меньше погрешностей двух других формул для всех приведённых величин M
и N. На расстояниях порядка H от сферы погрешность формулы (7) в 1.3-2 раза меньше, а
на расстояниях порядка H2 от сферы - в 2-3 раза меньше, чем у двух других квадратурных
формул.
Таким образом, сделанные тесты показывают, что новая квадратурная формула (7) обеспе-
чивает более высокую точность вычислений вблизи границы Γ, чем квадратурная формула из
[6] и формула c телесным углом. Эффективность формулы (7) для потенциала двойного слоя с
дифференцируемой плотностью особенно заметна, если плотность в потенциале представлена
гладкими осциллирующими функциями.
СПИСОК ЛИТЕРАТУРЫ
1. Белоцерковский С.М., Лифанов И.К. Численные методы в сингулярных интегральных уравнениях.
М., 1985.
2. Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент. М., 1995.
3. Сетуха А.В. Численные методы в интегральных уравнениях и их приложения. М., 2016.
4. Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов. М., 1987.
5. Гутников В.А., Лифанов И.К., Сетуха А.В. О моделировании аэродинамики зданий и сооружений
методом замкнутых вихревых рамок // Изв. РАН. Механика жидкости и газа. 2006. Т. 2006. № 4.
С. 78-93.
6. Крутицкий П.А., Резниченко И.О. Квадратурная формула для гармонического потенциала двой-
ного слоя // Дифференц. уравнения. 2021. Т. 57. № 7. С. 932-950.
7. 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. of Eng. Math. 2007. V. 59. P. 25-60.
8. Крутицкий П.А., Колыбасова В.В. Численный метод решения интегральных уравнений в задаче с
наклонной производной для уравнения Лапласа вне разомкнутых кривых // Дифференц. уравне-
ния. 2016. Т. 52. № 9. С. 1262-1276.
9. Krutitskii P.A. The Dirichlet problem for dissipative Helmholtz equation in a plane domain bounded by
closed and open curves // Hiroshima Math. J. 1998. V. 28. № 1. P. 149-168.
10. Krutitskii P.A. The Neumann problem for the 2-D Helmholtz equation in a domain bounded by closed
and open curves // Int. J. of Math. and Math. Sci. 1998. V. 21. № 2. P. 209-216.
11. Krutitskii P.A. The skew derivative problem in the exterior of open curves in a plane // Zeitschrift fur
Analysis und Ihre Anwendungen. 1997. Bd. 16. № 3. S. 739-747.
12. Krutitskii P.A. The 2-dimensional Dirichlet problem in an external domain with cuts // Zeitschrift fur
Analysis und Ihre Anwendungen. 1998. Bd. 17. № 2. S. 361-378.
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
КВАДРАТУРНАЯ ФОРМУЛА ДЛЯ ПОТЕНЦИАЛА ДВОЙНОГО СЛОЯ
1131
13. Krutitskii P.A. The Neumann problem in a 2-D exterior domain with cuts and singularities at the tips
// J. of Differ. Equat. 2001. V. 176. № 1. P. 269-289.
14. Krutitskii P.A. The 2-D Neumann problem in a domain with cuts // Rendiconti di Matematica e delle
sue Applicazioni. Ser. VII. 1999. V. 19. № 1. P. 65-88.
15. Krutitskii P.A. The mixed harmonic problem in an exterior cracked domain with Dirichlet condition on
cracks // Comput. and Math. with Appl. 2005. V. 50. P. 769-782.
16. Бутузов В.Ф., Крутицкая Н.Ч., Медведев Г.Н., Шишкин А.А. Математический анализ в вопросах
и задачах. М., 2000.
17. Ильин В.А., Позняк Э.Г. Основы математического анализа. Ч. 2. М., 1973.
18. Крутицкий П.А., Федотова А.Д., Колыбасова В.В. Квадратурная формула для потенциала прос-
того слоя // Дифференц. уравнения. 2019. Т. 55. № 9. С. 1269-1284.
19. Михлин С.Г. Линейные уравнения в частных производных. М., 1977.
20. Van Oosterom A., Strackee J. The solid angle of a plane triangle // IEEE Trans. on Biomedical Eng.
1983. V. 30. P. 125-126.
21. Casey J. A Treatise on Spherical Trigonometry. Dublin, 1889.
22. Владимиров В.С. Уравнения математической физики. М., 1981.
Институт прикладной математики
Поступила в редакцию 16.02.2022 г.
имени М.В. Келдыша РАН, г. Москва,
После доработки 24.04.2022 г.
Московский государственный университет
Принята к публикации 05.07.2022 г.
имени М.В. Ломоносова
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
