ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ, 2021, том 57, № 7, с.932-950

ЧИСЛЕННЫЕ МЕТОДЫ

УДК 519.644.5+517.956.224

КВАДРАТУРНАЯ ФОРМУЛА ДЛЯ ГАРМОНИЧЕСКОГО

ПОТЕНЦИАЛА ДВОЙНОГО СЛОЯ

Выводится квадратурная формула для гармонического потенциала двойного слоя с непре-

рывной плотностью, заданной на гладкой поверхности (замкнутой или разомкнутой).

Полученная квадратурная формула даёт более высокую точность вычислений, чем стан-

дартная квадратурная формула, что подтверждается численными тестами. Преимущество

новой квадратурной формулы особенно заметно вблизи поверхности, где значения, под-

считанные по стандартной квадратурной формуле, стремятся к бесконечности, тогда как

новая формула обеспечивает приемлемую точность вычислений.

DOI: 10.31857/S0374064121070074

Введение. Гармонический потенциал двойного слоя используется при решении краевых

задач для уравнения Лапласа методом интегральных уравнений. Такие задачи возникают в

различных областях математической физики, например, в теории обтекания препятствий по-

током идеальной жидкости, в электростатике, в стационарной теплопроводности, в теории

фильтрации и т.д. Численное решение краевых задач с помощью потенциала двойного слоя

обсуждалось в монографиях [1-3] и состоит из двух этапов. На первом этапе, численно решая

граничное интегральное уравнение, находят плотность потенциала. На втором этапе, подстав-

ляя численное значение плотности в квадратурную формулу, находят решение краевой задачи

в любой точке области. Однако квадратурные формулы для потенциалов, используемые в ин-

женерных расчётах [4, гл. 2], не дают равномерной аппроксимации потенциала в области и

не сохраняют свойство непрерывности потенциалов вплоть до границы области. Более того,

вблизи некоторых точек на границе области потенциалы, подсчитанные по этим квадратурным

формулам, стремятся к бесконечности, хотя сами потенциалы ограничены вблизи границы.

При использовании стандартных квадратурных формул для повышения точности приходит-

ся либо уменьшать шаг, либо проводить дополнительные построения вблизи границы, что

увеличивает стоимость вычислений. Актуальной является задача по получению улучшенных

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

В двумерном случае улучшенная квадратурная формула для потенциала простого слоя

с плотностью, заданной на разомкнутых кривых и имеющей степенные особенности на их

концах, построена в [5, 6]. Эта формула может применяться при нахождении численных реше-

ний краевых задач для уравнения Лапласа вне разрезов и разомкнутых кривых на плоскости

с использованием метода потенциалов и граничных интегральных уравнений. Такие задачи

изучались указанным методом в работах [7-9]. В [10] предложена улучшенная квадратурная

формула для потенциала простого слоя в трёхмерном случае, обеспечивающая равномерную

аппроксимацию и равномерную сходимость в области. Кроме того, эта формула сохраняет

свойство непрерывности потенциала простого слоя при переходе через границу области.

В настоящей работе построена улучшенная квадратурная формула для гармонического

потенциала двойного слоя. Формула обеспечивает более высокую точность вычислений, чем

стандартная квадратурная формула. Преимущество улучшенной формулы по сравнению со

стандартной особенно заметно на небольших расстояниях от границы.

1. Постановка задачи. Введём в пространстве трёхмерную декартову систему координат

x = (x₁,x₂,x₃) ∈ R³. Пусть Γ - либо простая гладкая замкнутая поверхность класса C²,

либо простая гладкая ограниченная разомкнутая ориентированная поверхность класса C²,

содержащая свои предельные точки [11, гл. 14, § 1]. Если поверхность Γ замкнутая, то она

932

КВАДРАТУРНАЯ ФОРМУЛА

933

должна ограничивать объёмно-односвязную внутреннюю область [12, c. 201]. Предположим,

что поверхность Γ параметризована так, что на неё отображается прямоугольник:

y = (y₁,y₂,y₃) ∈ Γ, y₁ = y₁(u,v), y₂ = y₂(u,v), y₃ = y₃(u,v); u ∈ [0,A], v ∈ [0,B];

y_j(u,v) ∈ C²([0,A] × [0,B]), j = 1,2,3.

(1)

Сферу, поверхность эллипсоида, гладкие поверхности фигур вращения, поверхность тора и

многие другие более сложные поверхности можно параметризовать таким образом.

Введём N точек u_n с шагом h на отрезке [0, A] и M точек v_m с шагом H на отрезке

[0, B] равенствами

u_n = (n + 1/2)h, n = 0,N - 1 и v_m = (m + 1/2)H, m = 0,M - 1; A = Nh, B = MH.

Прямоугольник [0, A] × [0, B], который отображается на поверхность Γ, разобьём на NM

маленьких прямоугольничков размера h × H, тогда точки (u_n, v_m) являются серединами

этих прямоугольничков.

Известно [11, гл. 14, § 1], что компоненты вектора нормали (не обязательно единичного)

η(y) = (η₁(y), η₂(y), η₃(y)) в точке поверхности y = (y₁, y₂, y₃) ∈ Γ выражаются через опреде-

лители второго порядка по формулам



y₂)_u (y₃)_u

(y₃)_u (y₁)_u

(y₁)_u (y₂)_u

η₁ =(



η₂ =



η₃ =



(2)

(y₂)_v (y₃)_v,

(y3)v (y1)v,

(y1)v (y2)v.

√

Положим |η(y)| =

(η₁(y))² + (η₂(y))² + (η₃(y))². Известно [11, гл. 14, §§ 1, 2], что

∫

F (y) ds_y = du dv F (y(u, v))|η(y(u, v))|.

Потребуем, чтобы выполнялось неравенство

|η(y(u, v))| > 0 для всех (u, v) ∈ (0, A) × (0, B).

(3)

Из условия (3) следует, что |η(y(u, v))| ∈ C¹((0, A) × (0, B)).

Через n_y обозначим единичную нормаль в точке y ∈ Γ, т.е. n_y = η(y)/|η(y)|. Производная

вдоль нормали n_y имеет вид

∂

= |η(y)|-1(η(y), ∇_y ).

∂n_y

√

Обозначим |x - y(u, v)| =

(x₁ - y₁(u, v))² + (x₂ - y₂(u, v))² + (x₃ - y₃(u, v))² и заметим, что

∑

∂

y_j - x_j

(η(y), y - x)

|x - y| =

η_j(y)

∂n_y

|η(y)|

|x - y|

|η(y)| |x - y|

j=1

∑

∂

(η(y), y - x)

-1

η_j(y)(y_j - x_j).

∂n_y |x - y|

|η(y)||x - y|³

j=1

Гармонический потенциал двойного слоя используется при решении краевых задач для

уравнения Лапласа методом интегральных уравнений. Пусть x ∈ Γ. Рассмотрим гармониче-

ский потенциал двойного слоя с заданной на поверхности Γ плотностью μ(y) ∈ C⁰(Γ):

∫

∂

W₀[μ](x) =

μ(y)

ds_y =

du dv μ(y(u, v))|η(y(u, v))|

4π

∂n_y |x - y|

4π

∂n_y |x - y(u,v)|

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

934

КРУТИЦКИЙ, РЕЗНИЧЕНКО

∫

∑

μ(y(u, v))

η_j(y(u,v))(y_j(u,v) - x_j).

(4)

4π

|x - y(u, v)|³

n=0 m=0

j=1

un-h/2

vm-H/2

Пусть μ_nm = μ(y(u_n, v_m)), тогда

μ(y(u, v)) = μ_nm + o(1)

(5)

для u ∈ [u_n - h/2, u_n + h/2] и v ∈ [v_m - H/2, v_m + H/2]. Следовательно,

∫

∑

W₀[μ](x) ≈ -

μ_nm

η_j(y(u,v))(y_j(u,v)-x_j). (6)

4π

|x - y(u, v)|³

n=0 m=0

j=1

un-h/2

vm-H/2

Таким образом, чтобы получить квадратурную формулу для гармонического потенциала двой-

ного слоя при x ∈ Γ, необходимо вычислить интеграл

∫

∑

η_j(y(u,v))(y_j(u,v) - x_j),

(7)

|x - y(u, v)|³

j=1

un-h/2

vm-H/2

который будем называть каноническим интегралом.

2. Вычисление канонического интеграла. Пусть точка x не принадлежит тому ку-

сочку поверхности Γ, на котором изменяется точка y = y(u, v), когда (u - u_n) ∈ [-h/2, h/2]

и (v - v_m) ∈ [-H/2, H/2]. Разложим функцию y_j(u, v) по формуле Тейлора с центром разло-

жения в точке (u_n, v_m), тогда для j = 1, 2, 3 получим

y_j(u,v) = y_j(u_n,v_m) + D_j + O(H² + h²),

где D_j = (y_j)^′u(u - u_n) + (y_j)^′v(v - v_m). Здесь и далее все производные по u и v берутся в

точке (u_n, v_m). Положим

∑

r² = |x - y(u_n,v_m)|² =

r2j = 0, r_j = y_j(u_n,v_m) - x_j, j = 1,2,3,

j=1

тогда

y_j(u,v) - x_j = r_j + D_j + O(H² + h²), j = 1,2,3.

Следовательно,

∑

|x - y(u, v)|² =

(x_j - y_j(u, v))² ≈

(r2j + 2r_j D_j + D2j) =

j=1

= r² + 2P(u - u_n) + 2Q(v - v_m) + α²(u - u_n)² + β²(v - v_m)² + 2δ(u - u_n)(v - v_m) =

= β²(V + δU/β² + Q/β²)² - (δU + Q)²/β² + α²U² + 2PU + r²,

где U = u - u_n, V = v - v_m,

∑

P = r_j(y_j)^′u, Q = r_j(y_j)^′v, α² =

((y_j )^′u)², β² =

((y_j )^′v)², δ =

(y_j )^′u(y_j )^′v.

j=1

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

КВАДРАТУРНАЯ ФОРМУЛА

935

Производные по u и v берутся в точке u = u_n, v = v_m. Несложно показывается [11,

гл. 14, § 1], что справедливо равенство

α²β² - δ² = |η(y(u_n,v_m))|².

Так как, согласно условию (3), |η(y(u_n, v_m))| > 0 для всех возможных n и m, то

α²β² - δ² > 0.

(8)

Отсюда следует, что α² > 0 и β² > 0.

Применяя формулу Тейлора с центром разложения в точке (u_n, v_m) с остаточным членом

в форме Пеано [11, гл. 10, § 5.3], находим

√

η_j(y(u,v)) = η_j(y(u_n,v_m)) + (η_j)^′u(u - u_n) + (η_j)^′v(v - v_m) + o(

(u - u_n)² + (v - v_m)²).

Производные по u и v берутся в точке (u_n, v_m).

Для вычисления выражения

∑

η_j(y(u,v))(y_j(u,v) - x_j)

j=1

с учётом формул

∑

η_j(y(u_n,v_m))(y_j)^′u =

η_j(y(u_n,v_m))(y_j)^′v = 0,

j=1

означающих ортогональность вектора нормали и касательных векторов к поверхности (см.

[11, гл. 14, § 1.2]), воспользуемся разложением по формуле Тейлора с центром разложения в

точке (u_n, v_m) с остаточным членом в форме Пеано

y_j(u,v) - x_j = r_j + (y_j)^′u(u - u_n) + (y_j)^′v(v - v_m)+

(y_j )^′′uu(u - u_n)² +

(y_j)^′′vv(v - v_m)² + (y_j)^′′uv(u - u_n)(v - v_m) + o((u - u_n)² + (v - v_m)²),

тогда

∑

η_j(y(u,v))(y_j(u,v) - x_j) ≈ R + ξ₄U + ξ₅V + ξ₁U² + ξ₂V² + ξ₃UV,

j=1

где U = u - u_n, V = v - v_m,

)

∑

⁽1

ξ₁ =

η_j(y(u_n,v_m))(y_j)^′′uu + (η_j)^′u(y_j)^′

ξ₂ =

η_j(y(u_n,v_m))(y_j)^′′vv + (η_j)^′v(y_j)^′

j=1

(

)

∑

ξ₃ =

η_j(y(u_n,v_m))(y_j)^′′uv + (η_j)^′u(y_j)^′v + (η_j)^′v(y_j)^′

j=1

∑

ξ₄ = (η_j)^′ur_j, ξ₅ =

(η_j )^′vr_j, R =

η_j(y(u_n,v_m))r_j.

j=1

Все производные по u, v берутся в точке (u_n, v_m).

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

936

КРУТИЦКИЙ, РЕЗНИЧЕНКО

Из приведённых соотношений вытекает, что канонический интеграл (7) приближённо равен

следующему интегралу, который обозначим через K_nm(x):

∫

∑

η_j(y(u,v))(y_j(u,v) - x_j) ≈

|x - y(u, v)|³

j=1

un-h/2

vm-H/2

∫

R+ξ₄U+ξ₅V +ξ₁U² +ξ₂V² +ξ₃UV

≈ dU

β³((V + δU/β² + Q/β²)² - (δU + Q)²/β⁴ + (α²U² + 2PU + r²)/β²)3/2

−h/2

-H/2

= K_nm(x).

Следовательно, чтобы вывести квадратурную формулу для потенциала двойного слоя,

необходимо вычислить интеграл K_nm(x) в явном виде.

2.1. Вычисление интегралов по dV. Введём обозначения

z = V + δU/β² + Q/β² = V + c, c = δU/β² + Q/β²,

(9)

a = -c² + (α²U² + 2PU + r²)/β² = -(δU + Q)²/β⁴ + (α²U² + 2PU + r²)/β² =

((α²β² - δ²)U² + 2(P β² - δQ)U + r²β² - Q²).

β⁴

Как показано выше, из неравенства (8) вытекает, что α² > 0 и β² > 0. Кроме того, из

неравенства (8) следует, что a ≡ 0, поскольку a является квадратным трёхчленом по U, в

котором коэффициент при U² положителен: (α²β² - δ²)/β⁴ > 0.

Покажем, что a ≥ 0. Положим

D_j = (y_j)^′uU - (y_j)^′vc, где c определено в (9), и проведём

следующие преобразования с учётом введённых обозначений:

∑

(r_j +Dj )² =

(r2j + 2r_j D_j +D2j) = r² + 2

r_j D_j +

j=1

= r² + 2PU - 2Qc + α²U² + β²c² - 2δUc =

(

δU + Q

⁽δU + Q)2)

(δU + Q)²

=β2 c2 - 2

+ α²U² + 2PU + r² =

β²

(

)₂

δU + Q

(δU + Q)²

=β2 -c +

+ α²U² + 2PU + r² =

β²

(δU + Q)

+ α²U² + 2PU + r² = aβ² ≥ 0.

(10)

β²

Так как β² > 0, то, поделив полученное неравенство на β², заключаем, что a ≥ 0. Следова-

тельно, квадратный трёхчлен a неотрицательный.

Обозначим

b=ξ₂c² -ξ₅c-ξ₃Uc+R+ξ₄U+ξ₁U² =

(

) (

)

δ²

δQ

ξ₅Q

Q²

=U²

+ U 2ξ₂

+ξ₄ -ξ₃

-ξ₅

+R-

+ξ₂

ξ₂ β4+ξ1 -ξ3 β2

β⁴

β²

β⁴

z_± = ±H/2 + c = ±H/2 + (δU + Q)/β².

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

КВАДРАТУРНАЯ ФОРМУЛА

937

Применяя введённые обозначения, вычислим в K_nm(x) интеграл по dV, переходя к перемен-

ной z:

∫

R+ξ₄U+ξ₅V +ξ₁U² +ξ₂V² +ξ₃UV

((V + δU/β² + Q/β²)² - (δU + Q)²/β⁴ + (α²U² + 2P U + r²)/β²)3/2

−H/2

∫

ξ₂(V + c - c)² + (ξ₃U + ξ₅)(V + c - c) + R + ξ₄U + ξ₁U²

((V + c)² + a)3/2

−H/2

∫

)



ξ₂z² + (ξ₃U + ξ₅ - 2ξ₂c)z + b

⁽⁽b

z+



= dz

-ξ₂

z - ξ₃U - ξ₅ + 2ξ₂c

√



(z² + a)3/2

z² + a

z_-



√

z+

+ ξ2 ln|z + z² + a|



z_-

где использованы интегралы 1.2.43.17-1.2.43.19 из справочника [13]. Заметим, что

∫

√

ξ₂

dU ln |z +

z² + a||z+z_- = ξ2βθnm(x),

−h/2

где функция θ_nm(x) найдена в явном виде в работе [10]. Тогда интеграл K_nm(x) можно

представить в виде

K_nm(x) =

(ξ₂βθ_nm(x) + J(H) - J(-H)).

β³

Так как функция θ_nm(x) найдена в [10], то задача сводится к вычислению интеграла

∫

(b/a - ξ₂)z_± - ξ₃U - ξ₅ + 2ξ₂c

J (±H) =

√

z2± + a

−h/2

∫

(b/a - ξ₂)(±H/2 + c) - ξ₃U - ξ₅ + 2ξ₂c

= dU

√

(±H/2 + c)² + a

-h/2

Достаточно вычислить интеграл J(H). Интеграл J(-H) вычисляется по тем же формулам,

что и интеграл J(H), с заменой в них H на -H.

Вычислим интеграл J(H). Распишем величины, входящие в подынтегральную функцию,

в виде многочленов по U:

a = C₂U²+C₁U+C₀, C₂ = (α²-δ²/β²)/β², C₁ = (2P -2δQ/β²)/β², C₀ = (r²-Q²/β²)/β²;

z²⁺ + a = B₂U² + B₁U + B₀, B₂ = α²/β², B₁ = (Hδ + 2P)/β², B₀ = H²/4 + (HQ + r²)/β²;

b = A₂U² + A₁U + A₀, A₂ = ξ₁ - ξ₃δ/β² + ξ₂δ²/β⁴, A₁ = ξ₄ - ξ₅δ/β² - ξ₃Q/β² + 2ξ₂δQ/β⁴,

A₀ = R - ξ₅Q/β² + ξ₂Q²/β⁴;

bz₊ = (A₂U² + A₁U + A₀)(δU/β² + H/2 + Q/β²) = E₃U³ + E₂U² + E₁U + E₀, E₃ = A₂δ/β²,

E₂ = A₂(H/2 + Q/β²) + A₁δ/β², E₁ = A₁(H/2 + Q/β²) + A₀δ/β², E₀ = A₀(H/2 + Q/β²);

ξ₂z₊ + ξ₃U + ξ₅ - 2ξ₂c = F₁U + F₀, F₁ = ξ₃ - ξ₂δ/β², F₀ = ξ₅ - ξ₂Q/β² + ξ₂H/2.

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

938

КРУТИЦКИЙ, РЕЗНИЧЕНКО

Применяя введённые обозначения, запишем интеграл J(H) в виде

J (H) = J₁ - J₂,

∫

E₃U³ + E₂U² + E₁U + E₀

J₁ =

√

(C₂U² + C₁U + C₀)

B₂U² + B₁U + B₀

−h/2

∫

F₁U + F₀

J₂ =

√

(11)

B₂U² + B₁U + B₀

-h/2

2.2. Вычисление интегралов (11). Используя деление многочленов и учитывая, что

C₂ > 0 в силу (8), приведём интеграл J₁ к виду

J₁ = J₁₁ + J₁₂,

где

∫

L₁U + L₀

J₁₁ =

√

B₂U² + B₁U + B₀

-h/2

∫

S₁U + S₀

J₁₂ =

√

(12)

(C₂U² + C₁U + C₀)

B₂U² + B₁U + B₀

−h/2

здесь

E₃

E₂C₂ - E₃C₁

L₁ =

L₀ =

S₁ = E₁ - C₀L₁ - C₁L₀, S₀ = E₀ - C₀L₀.

C₂

C²

Если B²¹ - 4B₂B₀ = 0 и -B₁/(2B₂) ∈ [-h/2, h/2], то в силу [10, п. 2] точка x лежит на

маленьком кусочке проходящей через y(u_n, v_m) касательной плоскости, по которому ведёт-

ся интегрирование в каноническом интеграле после линеаризации y(u, v) вблизи y(u_n, v_m).

В данном случае в знаменателе канонического интеграла с такой линеаризацией возникает

особенность в точке x. Однако если в числителе провести такую же линеаризацию, а нор-

маль приближённо заменить нормалью в точке y(u_n, v_m), то числитель будет тождественно

равен нулю, так как x лежит в касательной плоскости, проходящей через y(u_n, v_m). С другой

стороны, исходный канонический интеграл (без линеаризации) не имеет особенности, так как

x ∈ Γ, и может быть оценён как O(hH). Поэтому если выполняется условие B²¹ - 4B₂B₀ =

= 0 и -B₁/(2B₂) ∈ [-h/2,h/2], то будем считать, что K_nm(x) ≈ 0. Далее предполагаем, что

указанное условие не выполняется.

Так как B₂ > 0, интегралы J₁₁ и J₂ находятся с помощью табличных интегралов 2.261

и 2.264 из справочника [14]:

L₁ √

J₁₁ =

B₂U² + B₁U + B₀ +

B₂

(

)







√

U=h/2

L₁B₁





+ L₀ -

B₂U + B₁ + 2

B₂

B₂U² + B₁U + B₀

²



2B₂

√B₂ ln

U=-h/2

F₁ √

J₂ =

B₂U² + B₁U + B₀ +

B₂

(

)







√

U=h/2

F₁B₁





+ F₀ -

B₂U + B₁ + 2

B₂

B₂U² + B₁U + B₀

²



2B₂

√B₂ ln

U=-h/2

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

КВАДРАТУРНАЯ ФОРМУЛА

939

Остаётся вычислить интеграл J₁₂. Способ вычисления интеграла зависит от знака дискри-

минанта квадратного трёхчлена C₂U² + C₁U + C₀, стоящего в знаменателе подынтегральной

функции.

Первый случай: C²¹ - 4C₂C₀ > 0. В п. 2.1 показано, что квадратный трёхчлен, обо-

значенный через a, неотрицательный; следовательно, его дискриминант неположительный:

C²¹ - 4C₂C₀ ≤ 0, поэтому первый случай не реализуется.

Второй случай: C²¹ - 4C₂C₀ = 0. В этом случае C₂U² + C₁U + C₀ = C₂(U - U₁)²,

где U₁ = -C₁/(2C₂) - корень этого трёхчлена. Если U₁ ∈ [-h/2, h/2], то в силу (10) точка

x лежит на маленьком кусочке проходящей через y(u_n, v_m) касательной плоскости, по кото-

рому ведётся интегрирование в каноническом интеграле после линеаризации y(u, v) вблизи

y(u_n, v_m). В данном случае в знаменателе канонического интеграла с такой линеаризацией

возникает особенность в точке x. Однако если в числителе провести такую же линеаризацию,

а нормаль приближённо заменить нормалью в точке y(u_n, v_m), то числитель будет тождест-

венно равен нулю, так как x лежит в касательной плоскости, проходящей через y(u_n, v_m).

С другой стороны, исходный канонический интеграл (без линеаризации) не имеет особенно-

сти, так как x ∈ Γ, и может быть оценён как O(hH). Поэтому если U₁ ∈ [-h/2, h/2], то

будем считать, что K_nm(x) ≈ 0.

Пусть U₁ ∈ [-h/2, h/2]. Воспользовавшись равенством

S₁U + S₀

S₁(U - U₁) + S₀ + U₁S₁

S₁

S₀ + U₁S₁

(U - U₁)²

U-U₁

(U - U₁)²

получим

∫

S₁ dU

J₁₂ =

√

C₂(U - U₁)

B₂U² + B₁U + B₀

−h/2

∫

S₀ + U₁S₁

√

C₂

(U - U₁)²

B₂U² + B₁U + B₀

−h/2

Сделав в интеграле J₁₂ замену t = 1/(U - U₁), будем иметь [12, ч. 1, гл. 7, § 10, п. 5]:

t₁ =

t₂ = -

- пределы интегрирования,

h - 2U₁

h + 2U₁

t₁

∫

S₁

sgn (t) dt

J₁₂ = -

√

C₂

(U²¹B₂ + U₁B₁ + B₀)t² + (2B₂U₁ + B₁)t + B₂

t₂

∫

t₁

S₀ + U₁S₁

sgn (t)t dt

√

C₂

(U²¹B₂ + U₁B₁ + B₀)t² + (2B₂U₁ + B₁)t + B₂

t₂

∫t1

∫

t₁

S₁

S₀ + U₁S₁

t dt

sgn (C₁)

√

sgn (C₁)

√

C₂

ω₂t² + ω₁t + B₂

C₂

ω₂t² + ω₁t + B₂

t₂

где ω₁ = 2B₂U₁ + B₁, ω₂ = U²¹B₂ + U₁B₁ + B₀. Как показано выше, a ≥ 0, поэтому и ω₂ ≥ 0.

Если ω₂ > 0, то с помощью табличных интегралов 2.261 и 2.264 из справочника [14] находим



t₁

-S₁ sgn (C₁)



J₁₂ =

ln |2ω₂t + ω₁ + 2√ω2√ω2t2 + ω1t + B2|

C₂

√ω₂

t₂

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

940

КРУТИЦКИЙ, РЕЗНИЧЕНКО

(

S₀ + U₁S₁

1^√

sgn (C₁)

ω₂t² + ω₁t + B₂ -

C₂

ω₂

)

t₁

ω₁



√ω₂√ω2t2 + ω1t + B2|



2ω₂

√ω₂ ln|2ω2t+ω1 +2

t₂

Если ω₂ = 0, а ω₁ = 0, то, пользуясь интегралами 1.2.18.5, 1.2.18.6 из справочника [13],

получаем



t₁

S₁ 2

√



S₀ + U₁S₁

2(ω₁t - 2B₂)√

t1



J₁₂ = -

sgn (C₁)

ω₁t + B₂

sgn (C₁)

ω₁t + B₂



C₂ ω₁

C₂

3ω²¹

t₂

Если ω₂ = 0 и ω₁ = 0, то



t₁

S₁

sgn (C₁)



S₀ + U₁S₁ sgn (C₁)

t1

J₁₂ = -

t^



C₂

√B₂

C₂

2√B2 t2

t₂

Тем самым во втором случае интеграл J₁₂ вычислен явно.

Третий случай: C²¹-4C₂C₀ < 0. В этом случае многочлен C₂U²+C₁U +C₀ неприводим

над R. Рассмотрим различные варианты вычисления интеграла J₁₂ (см. (12)), воспользовав-

шись методом, предложенным в [12, ч. 1, гл. 7, § 10, п. 5, (7.75)] или в [14, раздел 2.25].

Вариант 1. Если B₁ = B₂C₁/C₂, то в интеграле J₁₂ достаточно сделать замену перемен-

ных U = t - C₁/(2C₂). Тогда

C₁

t₁ =

t₂ = -

t₂ < t₁,

2C₂

∫t1

[S₁t + (-S₁C₁/(2C₂) + S₀)] dt

J₁₂ =

√

C₂[t² + C₀/C₂ - C²¹/(4C²²)]

B₂t² + B₀ - B₂C²¹/(4C²²)

t₂

∫

(

)∫t1

)

⁽S₁

dt²

S₁C₁

√

+ -

+S₀

√

C₂

[t² + σ₁]

B₂t² + σ₂

2C₂

[t² + σ₁]

B₂t² + σ₂

t=t₂

t₂

где

C₀

C²¹

B₂C²¹

σ₁ =

> 0, σ₂ = B₀ -

≥ 0.

C₂

4C²²

4C²

Вводя обозначения

∫

S₁

dt²

J₁₂₁ =

√

2C₂

[t² + σ₁]

B₂t² + σ₂

t=t₂

(

)∫t1

S₁C₁

J₁₂₂ =

+S₀

√

C₂

√B₂

2C₂

[t² + σ₁]

t² + σ₂/B₂

t₂

получаем

J₁₂ = J₁₂₁ + J₁₂₂.

(13)

Воспользовавшись табличным интегралом 2.246 из справочника [14], находим интеграл

J121 в явном виде:

1₁) если σ₁B₂ - σ₂ > 0, то

)

t₁

S₁

( √B₂t2 + σ₂



J₁₂₁ =

;



C₂

√σ₁B₂ - σ₂ arctg

√σ₁B₂ - σ₂

t₂

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

КВАДРАТУРНАЯ ФОРМУЛА

941

2₁) если σ₁B₂ - σ₂ < 0, то

√



S₁



B₂t² + σ₂ -

√σ₂ - σ₁B₂t1

J₁₂₁ =

√



;



2C₂

√σ₂ - σ₁B₂ ln

B₂t² + σ₂ +

√σ₂ - σ₁B₂

t₂

3₁) если σ₁B₂ - σ₂ = 0, то



t₁

S₁



J₁₂₁ = -

√



C₂

B₂t² + σ₂

t₂

Используя табличные интегралы 1.2.43.17, 1.2.45.10, 1.2.45.13 и 1.2.11.10 из справочни-

ка [13], находим интеграл J₁₂₂ в явном виде:

1₂) если σ₂ > 0 и σ₂ - B₂σ₁ < 0, то

(

)

√



S₁C₁

σ₁ - σ₂/B₂ +

σ₁(t² + σ₂/B₂)|t1

J₁₂₂ =

+S₀

√



;



C₂

√B₂

2C₂

σ₁(σ₁ - σ₂/B₂)

t² + σ₁

t₂

2₂) если σ₂ - B₂σ₁ > 0, то

(

)

√



S₁C₁

σ₂/B₂ - σ₁

t1



J₁₂₂ =

+S₀

√

arctg

√

;



C₂

√B₂

2C₂

σ₁(σ₂/B₂ - σ₁)

σ₁(t² + σ₂/B₂)

t₂

3₂) если σ₁ - σ₂/B₂ = 0, то

(

)



S₁C₁

t1

J₁₂₂ =

+S₀

√



;



C₂

√B₂

2C₂

σ₁

t² + σ₁

t₂

4₂) если σ₂ = 0 и |B₁/(2B₂)| > h/2, то

(

)



S₁C₁

t²

t1

J₁₂₂ =

+S₀

sgn (B₁)



C₂

√B₂

2C₂

2σ₁

|t² + σ₁|

t₂

Итак, в этом варианте интеграл J₁₂ вычисляется явно по формуле (13).

Вариант 2. Пусть B₁ = B₂C₁/C₂. В [10, п. 2] показано, что в (12) квадратный трёхчлен

под корнем квадратным неотрицателен (этот результат вытекает также и из приведённых

в п. 2.1 выкладок), поэтому его дискриминант неположителен, т.е. B²¹/(4B²²) - B₀/B₂ ≤ 0.

Положим χ²¹ = B₀/B₂ - B²¹/(4B²²) ≥ 0, где без нарушения общности χ₁ ≥ 0. Поэтому могут

представиться только две возможности: χ₁ > 0 и χ₁ = 0. Рассмотрим их по отдельности.

Вариант 2а: χ₁ > 0. Запишем квадратный трёхчлен в виде

[(

)₂

]

[(

)₂

]

B₁

B₂U² + B₁U + B₀ = B₂ U +

+χ²

=B₂χ²

U+

/χ²¹ + 1 = B2χ21(sh2 t + 1),

2B₂

где сделана гиперболическая замена переменных

(

)

[(

)

]

B₁

sh t = U +

/χ₁, U = χ₁ sh t -

t = arcsh U +

/χ₁ .

2B₂

Теперь рассмотрим в (12) второй квадратный трёхчлен в знаменателе (учитывая условие,

которым выделяется третий случай) и линейную функцию в числителе - соответственно

(

)

χ₁B₁

B²¹

B₁C₁

C₂U² + C₁U + C₀ = C₂ χ²¹ sh² t -

sh t +

+ C₁χ1 sh t + C₀ -

B₂

4B²²

2B₂

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

942

КРУТИЦКИЙ, РЕЗНИЧЕНКО

(

)

χ₁B₁C₂

B²¹C₂

B₁C₁

= C₂χ21 sh2 t + C₁χ₁ -

sht +

+ C₀ = ν2 sh2 t + ν1 sht + ν₀ > 0

B₂

4B²²

2B₂

B₁S₁

S₁U + S₀ = S₁χ₁ sht -

+ S₀ = ε1 sht + ε₀,

2B₂

где ν₂ = C₂χ²¹ > 0, ν₁ = C₁χ₁ - χ₁B₁C₂/B₂, ν₀ = B²¹C₂/(4B²²) - B₁C₁/(2B₂) + C₀, ε₁ = S₁χ₁,

ε₀ = S₀ - B₁S₁/(2B₂). Так как ν₂ sh² t + ν₁ sh t + ν₀ > 0, то дискриминант этого квадратного

относительно sh t трёхчлена отрицательный, т.е.

ν²¹ - 4ν₂ν₀ < 0.

(14)

Кроме того, ν₁ = 0 в силу условий, принятых в варианте 2 и в варианте 2а.

Учитывая, что dU = χ₁ ch t dt и sh² t + 1 = ch² t, получаем

∫

[(

)

]

ε₁ sht + ε₀

B₁

J₁₂ =

ch t dt

t_± = arcsh

/χ₁ .

√B₂

(ν₂ sh² t + ν₁ sh t + ν₀) ch t

2B₂

Разобьём знаменатель на произведение двух линейных функций от sh t с комплексными ко-

эффициентами

∫

(

)

ε₁ sht + ε₀

λ_-

λ₊

J₁₂ =

√

√B₂

ν₂(sh t - G_-)(sh t - G₊)

ν₂

B₂

sh t - G_-

sh t - G₊

где G_-, G₊ - комплексные корни уравнения ν₂ sh² t + ν₁ sh t + ν₀ = 0 относительно sh t:

√

-ν₁ -

ν²¹ - 4ν₂ν₀

-ν₁ +

ν²¹ - 4ν₂ν₀

G- =

= g₁ - ig₂, G₊ =

= g₁ + ig₂,

2ν₂

√

ν₁

|ν²¹ - 4ν₂ν₀|

g₁ = -

= 0, g₂ =

> 0, g₁, g₂ ∈ R; G₊ = G_-.

(15)

2ν₂

Отметим, что g₂ = 0 в силу (14). Коэффициенты λ_± зависят от G_±:

ε₁G_- + ε₀

ε₁G₊ + ε₀

λ_- =

λ₊ =

(16)

G- - G+

-2g₂i

G+ - G-

2g₂i

Так как справедливы равенства

λ_±

e^t

2λ_±e^t

sht - G_±

(e^t - e^-t)/2 - G_± e^t

e2t - 2G_±e^t - 1

то, перейдя к замене ζ = e^t, dζ = e^t dt, dt = dζ/ζ,

(

[(

)

])

(

)

])

B₁

^[(h

B₁

ζ₁ = exp arcsh

/χ₁

ζ₂ = exp arcsh

/χ₁

(17)

2B₂

будем иметь

∫

(

)

∫

ζ₂

(

2λ_-e^t

2λ₊e^t

dζ

2λ_-ζ

J₁₂ =

√

ν₂

√B₂

e2t - 2G_-e^t - 1

e2t - 2G₊e^t - 1

ν₂

B₂

ζ ζ² - 2G_-ζ - 1

ζ₁

)

∫

ζ₂

(

)

2λ₊ζ

λ_-

λ₊

√

dζ

ζ² - 2G₊ζ - 1

ν₂

B₂

ζ² - 2G_-ζ - 1

ζ² - 2G₊ζ - 1

ζ₁

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

КВАДРАТУРНАЯ ФОРМУЛА

943

В данном случае нельзя применить стандартные формулы для нахождения первообразной

от вещественных подынтегральных функций, поскольку в нашем случае коэффициенты ком-_√√

плексные. Введём обозначения: Z1- = G_- - G2- + 1, Z₁₊ = G_- + G2- + 1,

√

Z2- = G₊ - G²⁺ + 1, Z₂₊ = G₊ + G²⁺ + 1,

(18)

γ_- =

√

γ₊ =

√

(19)

G2- + 1

G²⁺ + 1

Знаменатель в дробях (19) не обращается в нуль, поскольку g₁ = 0, как отмечено выше.

В результате

∫

ζ₂

(

)

λ_-

λ₊

J₁₂ =

dζ

ν₂

√B₂

(ζ - Z1-)(ζ - Z₁₊)

(ζ - Z2-)(ζ - Z₂₊)

ζ₁

[

∫

ζ₂

(

)

∫

ζ₂

(

)]

-2

√ λ_-γ_-

dζ

+λ₊γ₊

dζ

ν₂

B₂

ζ-Z1-

ζ-Z₁₊

ζ-Z2-

ζ-Z₂₊

ζ₁

Принимая во внимание выражения (15), (16), (18), (19), получаем λ_- = λ₊, Z1- = Z2-,

Z₁₊ = Z₂₊, γ₊ = γ_-. Тогда

∫

ζ₂

[

(

)

(

)]

-2

J₁₂ =

dζ λ₊γ₊

+λ₊γ₊

ν₂

√B₂

ζ-Z2-

ζ-Z₂₊

ζ-Z2-

ζ-Z₂₊

ζ₁

Так как ζ₁, ζ₂ ∈ R, то под интегралом находится сумма двух комплексно сопряжённых функ-

ций, поэтому интеграл J₁₂ вещественный:

[

∫

ζ₂

(

)]

-4

J₁₂ =

√ Re λ₊γ₊

dζ

ν₂

B₂

ζ-Z2-

ζ-Z₂₊

ζ₁

В итоге приходим к равенству

-4

J₁₂ =

√ Re [λ₊γ₊(ln(Z2- - ζ₂) - ln(Z2- - ζ₁) - ln(Z₂₊ - ζ₂) + ln(Z₂₊ - ζ₁))].

(20)

ν₂

B₂

Логарифмы с комплексными аргументами Z2± - ζ_l, где l = 1, 2, преобразуются по формуле

ln(Z2± - ζ_l) = ln |Z2± - ζ_l| + i arg(Z2± - ζ_l), l = 1, 2.

(21)

Рассмотрим величины, входящие в обозначения (18), (19). Так как g₂ = Im G₊ > 0, то

можно считать, что arg G₊ ∈ (0, π), тогда arg(G²⁺ + 1) ∈ (0, 2π) и

√

arg G²⁺ + 1 =

arg(G²⁺ + 1) ∈ (0, π),

√

поэтому Im G²⁺ + 1 > 0. Более того, если arg G₊ ∈ (0, π/2), то arg G²⁺ + 1 ∈ (0, π/2), а

√

если arg G₊ ∈ (π/2, π), то arg G²⁺ + 1 ∈ (π/2, π). Следовательно,

(

)

√

sgn Re G²⁺ + 1

= sgn Re G₊ = sgn g₁.

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

944

КРУТИЦКИЙ, РЕЗНИЧЕНКО

Положим

G²⁺ + 1 = g²¹ - g²² + 1 + 2g₁g₂i = |G²⁺ + 1|exp (iΨ), Ψ = arg(G²⁺ + 1),

√

g²¹ - g²² + 1

cos Ψ =

|G²⁺ + 1| = (g²¹ - g²² + 1)² + 4g²¹g²².

(22)

|G²⁺ + 1|

Используя тригонометрические формулы [15, с. 109], получаем

(

)

√

G²⁺ + 1 =

|G²⁺ + 1| exp (iΨ/2) =

|G²⁺ + 1| cos

+ isin

√

1 + cosΨ

|G²⁺ + 1| + g²¹ - g²² + 1

cos

= sgn (g₁)

2|G²⁺ + 1|

√

1 - cosΨ

|G²⁺ + 1| - (g²¹ - g²² + 1)

sin

(23)

2|G²⁺ + 1|

√

Поскольку arg G²⁺ + 1 = Ψ/2 ∈ (0, π), то (-Ψ/2) ∈ (-π, 0). Следовательно,

exp(-iΨ/2)

γ₊ =

√

G²⁺ + 1

|G²⁺ + 1|

√

|G²⁺ + 1| + g²¹ - g²²

cos(Ψ/2)

Re γ₊ =

√

= sgn (g₁)

√

2|G²⁺ + 1|

|G²⁺ + 1|

√

|G²⁺ + 1| - (g²¹ - g²²

+ 1)

sin(Ψ/2)

Im γ₊ = -

√

(24)

2|G²⁺ + 1|

|G²⁺ + 1|

где величины g₁ и g₂ определены в (15), а |G²⁺ + 1| - в (22).

Лемма. Пусть G₊ = g₁ + ig₂, где g₁, g₂ - вещественные числа и g₂ > 0, тогда

√

g₂ = ImG₊ > |Im G²⁺ + 1|.

Доказательство. Используя обозначения и соотношения из (22) и (23), находим

√

Im G²⁺ + 1 =

|G²⁺ + 1| sin

√

|G²⁺ + 1| - (g²¹ - g²² + 1)

|G²⁺ + 1|

(25)

2|G²⁺ + 1|

Так как g₂ > 0, то несложно проверить, что

√

|G²⁺ + 1| = (g²¹ - g²² + 1)² + 4g²¹g²² < g²¹ + g²² + 1,

поэтому справедлива оценка

(|G²⁺ + 1| - (g²¹ - g²² + 1)) < g²².

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

КВАДРАТУРНАЯ ФОРМУЛА

945

Из (25) вытекает равенство

(

√

)₂

Im G²⁺ + 1

(|G²⁺ + 1| - (g²¹ - g²² + 1)).

√

Следовательно, (Im G²⁺ + 1)² < g²², откуда |Im G²⁺ + 1| < g₂. Лемма доказана.

Следствие. Если выполнены условия леммы, то

(

)

√

Im (Z2± - ζ_l) = Im G₊ ± G²⁺ + 1

> 0, l = 1, 2.

√

Так как Z2± = G₊ ± G²⁺ + 1 в силу (18), а ζ₁ и ζ₂ - вещественные числа, то

Z2± - ζ_l = Re(Z2± - ζ_l) + iIm (Z2± - ζ_l),

√

Re (Z2± - ζ_l) = g₁ ±

|G²⁺ + 1| cos

- ζ_l, Im (Z2± - ζ_l) = g₂ ±

|G²⁺ + 1| sin

√(

√

)₂

(

√

)₂

|Z2± - ζ_l| =

g₁ ±

|G²⁺ + 1| cos

-ζ_l

+ g₂ ±

|G²⁺ + 1| sin

l = 1,2.

(26)

Согласно следствию к лемме справедливо неравенство Im (Z2± - ζ_l) > 0, поэтому можно

считать, что arg(Z2± - ζ_l) ∈ (0, π), где l = 1, 2. Следовательно,

√

g₁ ±

|G²⁺

+ 1| cos(Ψ/2) - ζ_l

Re(Z2± - ζ_l)

arg(Z2± - ζ_l) = arccos

= arccos

l = 1,2.

(27)

|Z2± - ζ_l|

В соотношениях (26), (27) используются обозначения из (15), (17), (22), (23).

Из равенств (16) вытекает, что

ε₁G₊ + ε₀

ε₁(g₁ + ig₂) + ε₀

ε₁

ε₁g₁ + ε₀

λ₊ =

G+ - G-

2ig₂

2g₂

поэтому

ε₁

ε₁g₁ + ε₀

Reλ₊ =

Imλ₊ = -

2g₂

Положим

ε₁

ε₁g₁ + ε₀

ε₁

ε₁g₁ + ε₀

λ₊γ₊ = Λ₁ + iΛ₂, Λ₁ =

Reγ₊ +

Im γ₊, Λ₂ =

Imγ₊ -

Re γ₊,

2g₂

где Re γ₊ и Im γ₊ однозначно определяются в (24) с использованием (15), (22). Применяя

формулу (26), находим сумму логарифмов

s₁ = ln |Z2- - ζ₂| - ln |Z2- - ζ₁| - ln |Z₂₊ - ζ₂| + ln |Z₂₊ - ζ₁| =

[(

√

)₂

(

√

)₂]

∑∑

(-1)q+l ln g₁

+ (-1)^q

|G²⁺ + 1| cos

-ζ_l

+ g₂

+ (-1)^q

|G²⁺ + 1| sin

q=1 l=1

Используя формулу (27), вычисляем сумму аргументов

s₂ = arg(Z2- - ζ₂) - arg(Z2- - ζ₁) - arg(Z₂₊ - ζ₂) + arg(Z₂₊ - ζ₁) =

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

946

КРУТИЦКИЙ, РЕЗНИЧЕНКО

√

∑∑

g₁ + (-1)^q

|G²⁺ + 1| cos(Ψ/2) - ζ_l

(-1)q+l arccos

√(

)₂

(

)₂

√

q=1 l=1

g₁+ (-1)^q

|G²⁺ +1| cos

−ζ_l

+ g₂+ (-1)^q

|G²⁺ +1| sin

В формулах для s₁, s₂ использованы обозначения из (15), (17), (22), (23).

Согласно (20) и (21) получаем

-4

J₁₂ =

√ Re [(Λ₁ + iΛ₂)(s₁ + is₂)].

ν₂

√B₂ Re[λ+γ⁺(s1 +is2)]=

ν₂

B₂

В итоге приходим к равенству

-4

J₁₂ =

√ (Λ₁s₁ - Λ₂s₂).

ν₂

B₂

Таким образом, в варианте 2а интеграл J₁₂ вычислен явно.

Вариант 2б: χ₁ = 0. В этом случае

[(

)₂

]

(

)₂

B₁

B₂U² + B₁U + B₀ = B₂ U +

+χ²

=B₂ U+

= B₂(U + Ω)²,

2B₂

где Ω = B₁/(2B₂). Как отмечено выше, считаем, что Ω ∈ [-h/2, h/2]. Тогда

∫

S₁U + S₀

J₁₂ =

√

(C₂U² + C₁U + C₀)

B₂U² + B₁U + B₀

−h/2

∫

(S₁U + S₀) sgn (U + Ω) dU

sgn (Ω)

(S₁U + S₀) dU

√

√B₂

(C₂U² + C₁U + C₀)(U + Ω)

B₂

(C₂U² + C₁U + C₀)(U + Ω)

-h/2

[

∫

]

sgn (Ω)

U dU

√ p₁

+p₂

+p₃

(28)

B₂

C₂U² + C₁U + C₀

U+Ω

−h/2

-h/2

где

C₂(S₁Ω - S₀)

S₀

C₀(S₁Ω - S₀)

S₁Ω - S₀

p₁ = -

p₂ =

p₃ =

C₁Ω - C₀ - C₂Ω²

Ω(C₁Ω - C₀ - C₂Ω²)

C₁Ω - C₀ - C₂Ω

Интегралы в (28) табличные. Воспользуемся формулами 1.2.8.19 и 1.2.8.13 из [13]. Учитывая,

что C²¹ - 4C₂C₀ < 0, находим

(

sgn (Ω)

p₁

J₁₂ =

p₃ ln |U + Ω| +

ln |C₂U² + C₁U + C₀| -

√B₂

2C₂

(

)

(

))

p₁C₁

2C₂U + C₁

2p₂

2C₂U + C₁

U=h/2

√

arctg

√

arctg

√



C₂

4C₂C₀ - C²¹

U=-h/2

Таким образом, в варианте 2б интеграл J₁₂ также вычислен явно.

3. Основной результат. Сформулируем основной результат этой работы.

Теорема. Пусть Γ - простая гладкая замкнутая поверхность класса C², ограничиваю-

щая объёмно-односвязную внутреннюю область, или простая гладкая ограниченная разомкну-

тая ориентированная поверхность класса C², содержащая свои предельные точки. Пусть Γ

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

КВАДРАТУРНАЯ ФОРМУЛА

947

допускает параметризацию (1) со свойством (3) и μ(y) ∈ C⁰(Γ). Тогда для гармонического

потенциала двойного слоя (4) при x ∈ Γ имеет место квадратурная формула

∑

W₀[μ](x) ≈ -

μ_nmK_nm(x),

(29)

4π

n=0 m=0

где интегралы K_nm(x) вычислены в явном виде в п. 2.

4. Стандартная квадратурная формула. Квадратурная формула (29) является аль-

тернативой стандартной квадратурной формуле для потенциала двойного слоя вне поверхно-

сти Γ, используемой в инженерных расчётах [4, гл. 2]:

∑

W₀[μ](x) ≈ -

μ_nm

η_j(y(u_n,v_m))(y_j(u_n,v_m) - x_j).

(30)

4π

|x - y(un,vm)|3

n=0 m=0

j=1

Стандартная квадратурная формула получается из формулы (6) заменой канонического ин-

теграла из (7) на следующее его приближённое значение:

∑

η_j(y(u_n,v_m))(y_j(u_n,v_m) - x_j).

|x - y(u_n, v_m)|³

j=1

Пусть для упрощения Γ - замкнутая поверхность. Значение, вычисленное по стандартной

квадратурной формуле (30), стремится к бесконечности, если точка x стремится изнутри или

извне Γ к одной из точек y(u_n, v_m) ∈ Γ, хотя сам потенциал двойного слоя в наших предпо-

ложениях ограничен на Γ и непрерывно продолжим на Γ как изнутри, так и извне. Другими

словами, если точка x стремится к точке на Γ изнутри или извне, то потенциал двойного

слоя стремится к конечному пределу. Стандартная квадратурная формула не сохраняет важ-

нейшие свойства потенциала двойного слоя, а именно, ограниченность на Γ и непрерывную

продолжимость на Γ извне и изнутри. Стандартная квадратурная формула расходится вблизи

границы.

5. Численные тесты. Тестирование улучшенной (29) и стандартной (30) квадратурных

формул проведено в случае, когда поверхность Γ является сферой единичного радиуса, кото-

рая задана параметрически уравнениями:

y₁(u,v) = cos usin v, y₂(u,v) = sin usin v, y₃(u,v) = cos v,

(31)

причём (u, v) ∈ [0, 2π]×[0, π]. Отметим, что в данном случае |η(y(u, v))| = sin v и |η(y(u, 0))| =

= |η(y(u, π))| = 0 для всех u ∈ [0, 2π]. Иначе говоря, |η(y)| = 0 на полюсах сферы при такой

параметризации, но условия теоремы выполняются.

Согласно [16, § 27.6] плотность потенциала двойного слоя на поверхности Γ можно найти

по формуле

μ(x)|_Γ = W₀[μ](x)|_Γ- - W₀[μ](x)|_Γ+ .

Здесь поверхность Γ рассматривается как двусторонняя, через Γ^- обозначена сторона, кото-

рую мы видим, глядя навстречу вектору нормали n_y, а через Γ⁺ - противоположная сторона.

В формуле берутся предельные значения потенциала двойного слоя на разных сторонах Γ. От-

метим, что направление единичной нормали n_y совпадает с направлением нормали η, так

как вектор n_y получается из η в результате нормировки. Пусть теперь Γ - единичная сфера,

заданная параметризацией (31), тогда формулы (2) для нормали η определяют внутреннюю

нормаль на сфере, а значит, Γ^- - внутренняя сторона единичной сферы, а Γ⁺ - её внешняя

сторона.

В рассматриваемых тестовых примерах для потенциала двойного слоя с заданной на еди-

ничной сфере плотностью известно явное выражение во всём пространстве, поэтому точные

значения потенциала можно сравнить с приближёнными, вычисленными по квадратурным

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

6^∗

948

КРУТИЦКИЙ, РЕЗНИЧЕНКО

формулам. Во всех тестах приближённое значение потенциала двойного слоя вычислялось по

стандартной квадратурной формуле (30) и по улучшенной квадратурной формуле (29) в неко-

торых точках на вспомогательных сферах, имеющих центры в начале координат и радиусы,

равные 1 ± ΔR. Тем самым вспомогательные сферы находятся либо внутри, либо снаружи

сферы единичного радиуса, на которой задана плотность потенциала, на расстоянии ΔR от

неё. Затем были рассчитаны значения абсолютных погрешностей в этих точках, и для каждой

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

Координаты точек, которые использовались для оценки максимальной абсолютной погреш-

ности, выбирались следующими:

x^qlj = Ry_j(u_q,v_l), j = 1,2,3,

2π

u_q =

q, q = 0,2N; v_l =

l = 1,2M - 1,

(32)

где величина y_j(u, v) определяется формулами (31), R - радиус вспомогательной сферы.

В частности, эти точки расположены над и под центрами участков разбиения единичной сфе-

ры, серединами границ между такими участками и пересечениями этих границ. Отметим, что

эти точки распределены по всей сфере.

Вычисления проводились для различных значений M и N. Величины шагов выбирались

по формулам h = 2π/N, H = π/M. Если N = M = 25, то h ≈ 0.25, H ≈ 0.13; если

N = M = 50, то h ≈ 0.126, H ≈ 0.063; если N = M = 100, то h ≈ 0.063, H ≈ 0.031.

В таблицах приведены рассчитанные максимальные значения абсолютных погрешностей.

В левом столбце указано отличие радиуса вспомогательной сферы от единицы: для внутренних

сфер радиус равен 1 - ΔR, для внешних 1 + ΔR. В верхней строке указаны значения M, N.

Первое число в ячейках таблицы - максимальная погрешность для стандартной квадратурной

формулы на данной вспомогательной сфере, а число после точки с запятой - максимальная

погрешность на данной сфере для улучшенной формулы.

Тест 1. В данном тесте использовалась плотность потенциала μ(y(u, v)) = 1; тогда гар-

монический потенциал двойного слоя имеет вид [16, § 27.2]

{

при |x| < 1,

W₀[μ](x) =

при |x| > 1.

В табл. 1 приведены рассчитанные максимальные значения абсолютных погрешностей.

Таблица 1. Максимальная абсолютная погрешность квад-

ратурных формул в тесте 1

ΔR M = N = 25 M = N = 50 M = N = 100

Внутренние сферы

0.1

0.096; 0.021

0.0079; 0.0045

0.0020; 0.0011

0.06

0.44; 0.041

0.054; 0.0083

0.0055; 0.0018

0.03

2.42; 0.094

0.45; 0.022

0.056; 0.0041

Внешние сферы

0.1

0.11; 0.021

0.0081; 0.0039

0.0020; 0.00092

0.06

0.47; 0.048

0.061; 0.0080

0.0055; 0.0016

0.03

2.46; 0.15

0.47; 0.023

0.060; 0.0041

Тест 2. В данном тесте использовалась плотность потенциала μ(y(u, v)) = cos v; тогда

гармонический потенциал двойного слоя имеет вид

⎧

⎪2|x| cos ϑ

⎨

при |x| < 1,

W₀[μ](x) =

⎪

cos ϑ

⎩-

при |x| > 1,

3|x|²

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

КВАДРАТУРНАЯ ФОРМУЛА

949

где ϑ - зенитный угол в сферических координатах с центром в начале координат. В табл. 2

приведены рассчитанные максимальные значения абсолютных погрешностей.

Таблица 2. Максимальная абсолютная погрешность квад-

ратурных формул в тесте 2

ΔR M = N = 25 M = N = 50 M = N = 100

Внутренние сферы

0.1

0.032; 0.011

0.0078; 0.0030

0.0020; 0.00076

0.06

0.18; 0.017

0.020; 0.0044

0.0055; 0.0012

0.03

1.14; 0.043

0.19; 0.0086

0.020; 0.0021

Внешние сферы

0.1

0.040; 0.0076

0.0077; 0.0016

0.0020; 0.00042

0.06

0.20; 0.017

0.020; 0.0030

0.0055; 0.00080

0.03

1.16; 0.050

0.20; 0.0086

0.020; 0.0017

Тест 3. В данном тесте использовалась плотность потенциала μ(y(u, v)) = cos u sin v; тогда

гармонический потенциал двойного слоя имеет вид

⎧

⎪2|x| cos ϕ sin ϑ

⎨

при |x| < 1,

W₀[μ](x) =

⎪

cos ϕ sin ϑ

⎩-

при |x| > 1,

3|x|²

где ϑ и ϕ - зенитный и азимутальный углы в сферических координатах с центром в начале

координат. В табл. 3 приведены рассчитанные максимальные значения абсолютных погреш-

ностей.

Таблица 3. Максимальная абсолютная погрешность квадра-

турных формул в тесте 3

ΔR M = N = 25 M = N = 50

M = N = 100

Внутренние сферы

0.1

0.096; 0.023

0.0055; 0.0049

0.000028; 0.0012

0.06

0.44; 0.043

0.054; 0.0088

0.0021; 0.0019

0.03

2.42; 0.096

0.45; 0.022

0.056; 0.0043

Внешние сферы

0.1

0.12; 0.019

0.0082; 0.0036

0.000069; 0.00085

0.06

0.48; 0.045

0.061; 0.0077

0.0029; 0.0016

0.03

2.46; 0.15

0.47; 0.023

0.060; 0.0040

Тест 4. В данном тесте использовалась плотность потенциала μ(y(u, v)) = (3 cos² v - 1)/2;

тогда гармонический потенциал двойного слоя имеет вид

⎧

3|x|²(3 cos² ϑ - 1)

⎪

при |x| < 1,

⎨

W₀[μ](x) =

⎪

(3 cos² ϑ - 1)

⎩-

при |x| > 1,

5|x|³

где ϑ - зенитный угол в сферических координатах с центром в начале координат. В табл. 4

приведены рассчитанные максимальные значения абсолютных погрешностей.

Заключение. Из таблиц видно, что улучшенная квадратурная формула обеспечивает бо-

лее высокую точность вычислений вблизи границы Γ, чем стандартная. Кроме того, значения,

вычисленные по стандартной формуле, стремятся к бесконечности при приближении к гра-

нице. Тесты показывают, что улучшенная формула даёт хорошую точность вычислений для

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

950

КРУТИЦКИЙ, РЕЗНИЧЕНКО

всех точек, расположенных на расстоянии H и более от границы Γ. В этом случае улучшен-

ная квадратурная формула имеет второй порядок сходимости и обеспечивает погрешность

вычислений порядка O(hH). На расстояниях порядка hH до границы улучшенная формула

даёт погрешность O(H), а стандартная формула расходится.

Таблица 4. Максимальная абсолютная погрешность квад-

ратурных формул в тесте 4

ΔR M = N = 25 M = N = 50 M = N = 100

Внутренние сферы

0.1

0.047; 0.012

0.0079; 0.0026

0.0020; 0.00063

0.06

0.22; 0.022

0.027; 0.0045

0.0055; 0.0011

0.03

1.21; 0.048

0.22; 0.011

0.028; 0.0022

Внешние сферы

0.1

0.057; 0.010

0.0077; 0.0019

0.0020; 0.00044

0.06

0.24; 0.023

0.031; 0.0039

0.0055; 0.00080

0.03

1.23; 0.074

0.23; 0.012

0.03; 0.0020

СПИСОК ЛИТЕРАТУРЫ

1. Белоцерковский С.М., Лифанов И.К. Численные методы в сингулярных интегральных уравнениях.

М., 1985.

2. Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент. М., 1995.

3. Сетуха А.В. Численные методы в интегральных уравнениях и их приложения. М., 2016.

4. Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов. М., 1987.

5. 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.

6. Крутицкий П.А., Колыбасова В.В. Численный метод решения интегральных уравнений в задаче с

наклонной производной для уравнения Лапласа вне разомкнутых кривых // Дифференц. уравне-

ния. 2016. Т. 52. № 9. С. 1262-1276.

7. Крутицкий П.А. Смешанная задача для уравнения Лапласа вне разрезов на плоскости // Диффе-

ренц. уравнения. 1997. Т. 33. № 9. С. 1181-1190.

8. Krutitskii P.A. The 2-dimensional Dirichlet problem in an external domain with cuts // Zeitschr. für

Analysis und ihre Anwend. 1998. Bd. 17. № 2. S. 361-378.

9. Krutitskii P.A. The Dirichlet problem for the two-dimensional Laplace equation in a multiply connected

domain with cuts // Proc. of the Edinburgh Math. Soc. 2000. V. 43. № 2. P. 325-341.

10. Крутицкий П.А., Федотова А.Д., Колыбасова В.В. Квадратурная формула для потенциала прос-

того слоя // Дифференц. уравнения. 2019. Т. 55. № 9. С. 1269-1284.

11. Бутузов В.Ф., Крутицкая Н.Ч., Медведев Г.Н., Шишкин А.А. Математический анализ в вопросах

и задачах. М., 2000.

12. Ильин В.А., Позняк Э.Г. Основы математического анализа. Ч. 2. М., 1973.

13. Прудников А.П., Брычков Ю.А., Маричев О.И. Интегралы и ряды. М., 1981.

14. Градштейн И.С., Рыжик И.М. Таблицы интегралов, сумм, рядов и произведений. М., 1963.

15. Цыпкин А.Г., Цыпкин Г.Г. Математические формулы. М., 1985.

16. Владимиров В.С. Уравнения математической физики. М., 1981.

Институт прикладной математики

Поступила в редакцию 03.03.2021 г.

им. М.В. Келдыша РАН, г. Москва,

После доработки 03.03.2021 г.

Московский государственный университет

Принята к публикации 27.04.2021 г.

им. М.В. Ломоносова

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021