ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ, 2021, том 57, № 9, с.1273-1279
ИНТЕГРАЛЬНЫЕ
И ИНТЕГРО-ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ
УДК 517.968.73:534+519.642.7
ОБЪЁМНЫЕ ИНТЕГРАЛЬНЫЕ УРАВНЕНИЯ
С ЗАПАЗДЫВАНИЕМ ПО ВРЕМЕНИ ДЛЯ РЕШЕНИЯ
НЕСТАЦИОНАРНЫХ ЗАДАЧ АКУСТИКИ
© 2021 г. А. Б. Самохин, А. С. Самохина, И. А. Юрченков
Выводятся объёмные интегральные уравнения с запаздыванием по времени, описывающих
нестационарные задачи рассеяния акустического поля на прозрачных трёхмерных струк-
турах. Предлагается эффективный метод численного решения полученных уравнений.
DOI: 10.31857/S0374064121090120
Введение. В работах [1, 2] рассматривались задачи взаимодействия нестационарного элек-
тромагнитного поля с ограниченной материальной средой, окружённой свободным простран-
ством. Используя объёмные сингулярные интегральные уравнения в частотной области [3]
и свойства преобразования Фурье, для решения этих задач получены интегральные уравне-
ния с запаздыванием по времени. В настоящей статье, применяя некоторые результаты работ
[1, 2], рассматриваются нестационарные задачи рассеяния акустического поля на прозрачных
трёхмерных структурах. Выводятся объёмные интегральные уравнения с запаздыванием по
времени, которые описывают указанный класс задач. Предлагается эффективный метод чис-
ленного решения полученных уравнений.
1. Исходное интегральное уравнение. Рассмотрим следующий класс задач нестацио-
нарной акустики. В ограниченной трёхмерной области Q, окружённой свободным простран-
ством, материальная среда характеризуется индексом рефракции n(x), x ∈ Q, а вне Q
индекс рефракции равен единице. Требуется определить акустическое поле в пространст-
ве R3, порождаемое внешним источником f0(x, t), локализованным в конечной области, при-
чём f0(x, t) = 0 при t < 0. В такой постановке соответствующая математическая задача
формулируется следующим образом: найти функцию акустического поля U(x, t) в R3, удо-
влетворяющую волновому уравнению
n(x) ∂2U(x, t)
ΔU(x,t) -
= f0(x,t),
(1)
c2
∂t2
где c - скорость звука в свободном пространстве. К уравнению (1) необходимо добавить ещё
условия, гарантирующие отсутствие волн, приходящих из бесконечности.
Для вывода нестационарных интегральных уравнений будем использовать прямое и обрат-
ное преобразования Фурье:
+∞
∫
1
f (ω) = F [f(t)] =
f (t)eiωt dt, f(t) = F-1[f(ω)] =
f (ω)e-iωt dω.
(2)
2π
−∞
-∞
В формулах (2) исходная функция f(t) и её фурье-образ f(ω) различаются указанием аргу-
ментов t и ω соответственно.
Применяя преобразование Фурье к каждой части уравнения (1), получаем уравнение для
фурье-образов поля
2
n(x)ω
ΔU(x,ω) -
U (x, ω) = f0(x, ω).
(3)
c2
Фурье-образ U(x, ω) должен удовлетворять, кроме уравнения (3), ещё и условию излуче-
ния на бесконечности
(
(
))
√
∂U
ω
lim
r
-i
U
= 0, r = |x| = (x21 + x22 + x23).
(4)
r→∞
∂r
c
1273
1274
САМОХИН и др.
Пусть f(x, ω) - финитная функция в R3 относительно x. Тогда интегральное пред-
ставление
∫
V (x, ω) = - f(y, ω)G(R, ω) dy
(5)
удовлетворяет в R3 уравнению Гельмгольца
ΔV + (ω/c)2V = f
(6)
и условию излучения на бесконечности вида (4). В (5) G(R, ω) - функция Грина уравнения
Гельмгольца, которая в декартовой системе координат имеет вид
G(R, ω) = exp(iωR/c)/(4πR), R = |x - y|, x = (x1, x2, x3), y = (y1, y2, y3).
(7)
Запишем уравнение (3) в следующем виде:
ΔU + (ω/c)2U = f0 - (ω/c)2(n - 1)U.
(8)
Из соотношений (3)-(8) вытекает, что неизвестное поле U(x, ω) имеет следующее интегральное
представление:
∫
(ω)2∫
U (x, ω) = - f0(y, ω)G(R, ω) dy +
(n(y) - 1)U(y, ω)G(R, ω) dy, x ∈ R3.
(9)
c
Q
Первый интеграл в правой части равенства (9) описывает поле U0(x, ω), создаваемое источ-
ником f0(x, ω) в свободном пространстве. Далее, поскольку n(x) = 1 вне Q, из равенства (9)
следует, что неизвестное поле U(x, ω) удовлетворяет в области Q интегральному уравнению
Фредгольма второго рода
)2 ∫
(ω
U (x, ω) -
(n(y) - 1)U(y, ω)G(R, ω) dy = U0(x, ω), x ∈ Q.
c
Q
Зная поле U(x, ω) в области Q, несложно найти поле в любой точке пространства, используя
равенство (9).
Имеют место следующие свойства Фурье-преобразований:
d2f(t)
F-1[-ω2f(ω)] =
,
(10)
dt2
+∞
F-1[f(ω)g(ω)] =
f (τ)g(t - τ) dτ, F-1[eiωΔ] = δ(t - Δ),
(11)
−∞
где δ - дельта-функция Дирака, а Δ - некоторая постоянная. Преобразование Фурье обобщён-
ных функций определяется для класса S′ так называемых обобщённых функций медленного
роста - линейных непрерывных функционалов на пространстве S основных функций, пред-
ставляющем собой линейное подпространство пространства C∞(R), состоящее из функций,
которые при |x| → ∞ убывают вместе со всеми производными быстрее любой степени |x|-1,
с заданной в нём специальной сходимостью (см. [4, гл. II, § 8.1]). Именно, преобразованием
Фурье обобщённой функции f ∈ S′ называется [4, гл. II, § 9.2] обобщённая функция F [f]
такая, что (F [f], ϕ) = (f, F [ϕ]) для всех ϕ ∈ S. Несложно убедиться [4, гл. II, § 9.2], что так
определённое отображение F переводит S′ в S′ и является линейным и непрерывным, а для
обратного отображения F-1 очевидно равенство (F-1[f], ϕ) = (f, F-1[ϕ]) для всех ϕ ∈ S.
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57
№9
2021
ОБЪЁМНЫЕ ИНТЕГРАЛЬНЫЕ УРАВНЕНИЯ С ЗАПАЗДЫВАНИЕМ ПО ВРЕМЕНИ 1275
В силу определений (2) равенство (10) очевидно. Покажем справедливость первого равен-
ства в (11). Имеем
∫
)
∫
∫
)
F
f (τ)g(t - τ) dτ
= eiω(t-τ)
f (τ)g(t - τ) dτ dt =
−∞
-∞
-∞
+∞
∫
)
= f(τ)eiωτ
eiω(t-τ)g(t - τ)d(t - τ) dτ = f(ω)g(ω).
(12)
−∞
-∞
Применяя обратное преобразование Фурье к первому и последнему членам в цепочке равенств
(12), получаем первое равенство в (11). Доказательство второго равенства в (11) следует непо-
средственно из определения преобразования Фурье обобщённых функций и приведено, напри-
мер, в [4, гл. II, § 9.2, пример].
Из представления (7) и второго равенства в (11) вытекает, что
F-1[G(R,ω)] = (4πR)-1δ(t - R/c).
(13)
В силу соотношений (10)-(13) получаем
∫
)
∂2f (y, t - τ)
F-1[G(R,ω)ω2f(y,ω)] = -(4πR)-1
δ(τ - R/c)
dτ
=
∂t2
−∞
∂2f(y,τ)
= -(4πR)-1
,
τ = t - R/c.
(14)
∂τ2
Теперь, применяя обратное преобразование Фурье к уравнению (9) и используя равенство (14),
приходим к следующему объёмному интегральному уравнению с запаздыванием по времени:
∫
1
1
∂2U(y,τ)
U (x, t) = U0(x, t) -
(n(y) - 1)
dy,
4πc2
R
∂τ2
Q
R = |x - y|; τ = t - R/c; U(x,t) = 0 при t < 0.
(15)
Из уравнения (15) очевидно, что акустическое поле в R3 зависит только от значений поля
в области неоднородности Q. Поэтому (15) можно рассматривать как интегральное уравне-
ние в Q.
Из изложенного выше следует
Теорема. Пусть индекс рефракции акустического поля равен n(x) в ограниченной облас-
ти Q ⊂ R3 и единице вне этой области. Тогда акустическое поле U, порождаемое внешним
полем U0, описывается объёмным интегральным уравнением (15), содержащим запаздыва-
ние по времени, зависящее от пространственных переменных.
2. Метод решения. Будем использовать метод коллокации для численного решения урав-
нения (15). Для этого в прямоугольной декартовой системе координат в R3 введём конечное
множество точек - сетку - так, чтобы область Q целиком находилась в прямоугольном парал-
лелепипеде Π со сторонами Nlh, l = 1, 2, 3, где h - шаг сетки по декартовым координатам.
Тогда если x0 = (x01, x02, x03) - вершина этого параллелепипеда с наименьшей суммой координат,
то координаты каждого узла сетки имеют вид x0 + (p1h, p2h, p3h), где pl = 0, Nl, l = 1, 2, 3,
т.е. вид x0 + hp, где p = (p1, p2, p3). Поэтому каждый узел однозначно определяется точкой
p = (p1,p2,p3), pl = 0,Nl, l = 1,2,3, которую назовём сеточными координатами узла и
будем писать p ∈ Q, если узел с сеточными координатами p лежит в области Q. Параллеле-
пипед Π разбивается данной сеткой на ячейки (элементарные кубики) Π(p), p = (p1, p2, p3),
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57
№9
2021
1276
САМОХИН и др.
pl = 0,Nl - 1, где p - вершина кубика с наименьшей суммой сеточных координат. Определим
область
Q как объединение NQ ячеек, центры которых лежат внутри области Q. Центр
кубика Π(p), содержащегося в области
Q, обозначим через x(p). Узловые точки, в которых
определяются значения акустического поля, будем задавать в центрах x(p). Для дискрети-
зации уравнения (15) по времени введём временную сетку t(0), t(1), t(2), . . . с постоянным
шагом δ, где t0 = 0 - начало процесса. Будем обозначать через V (p, t(n)) значение функции
в p-й пространственной узловой точке и во временной узловой точке t(n) = nδ, n = 0,1,2,... ,
а через V (q, τ(n, p, q)) - значение функции в q-й пространственной узловой точке и во вре-
менной точке
τ (n, p, q) = nδ - |x(p) - x(q)|/c.
(16)
Теперь уравнение (15) аппроксимируется следующей системой уравнений, имеющей раз-
мерность ∼NQ:
∑
∂2P(p,t)
1
∂2P(q,τ)
U (p, t) = U0(p, t) + B(0)
+
B(p - q)
,
p ∈ Q, n = 1,2,3,... ,
∂t2
c2
∂τ2
q∈Q
q=p
∫
1
P (q, τ) = (n(q) - 1)U(q, τ), t = t(n), τ = τ(n, p, q), B(p - q) = -
dy.
(17)
4π|x(p) - y|
Π(q)
Значения функции в суммах (17) должны, в связи с запаздыванием по времени, задавать-
ся во временных точках, предшествующих t(n). Это означает, что для всех τ, входящих в
систему (17), должно выполняться условие τ ≤ t(n - 1). Минимальное расстояние между
узловыми точками равно h. Тогда из (16) следует, что между шагами сеток по пространству
и по времени должно выполняться соотношение
h/c ≥ δ.
(18)
Условие (18) совпадает с критерием Куранта-Фридрихса, который является необходимым
условием для устойчивости численного решения гиперболических уравнений, описывающих
распространение волн.
Определим следующую связь между шагами сеток:
h = 2cδ.
(19)
Коэффициент 2 обусловлен тем, чтобы при вычислении производных по разностным форму-
лам не выходить за временные узловые точки, стоящие слева от t(n). Из равенств (16), (19)
вытекает, что для всех значений τ, входящих в (17), выполняется оценка
τ (n, p, q) ≤ (n - 2)δ = t(n - 2), q = p.
Конкретные значения шагов δ и h, которые связаны соотношением (19), зависят от спе-
цифики задачи и требуемой точности решения. Можно предложить следующий способ опре-
деления этих величин. Разложим источник поля f0(x, t), создающий внешнее поле U0(x, t),
в интеграл Фурье по времени. В силу равенства Парсеваля, учитывая, что f0(x, t) = 0 при
t < 0, получаем
∫
∫
∫
∫
I0 =
|f0(x, t)|2 dt dx =
|f0(x, ω)|2 dω dx,
(20)
V
0
V -∞
где V - область локализации источника поля. Определим полосу частот [-ωmax, ωmax], кото-
рая достаточно точно описывает поведение функции f0(x, t), следующей формулой:
∫
∫
I(ωmax)
1-
≪ 1, I(ωmax) =
|f0(x, ω)|2 dω dx.
(21)
I0
V -ωmax
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57
№9
2021
ОБЪЁМНЫЕ ИНТЕГРАЛЬНЫЕ УРАВНЕНИЯ С ЗАПАЗДЫВАНИЕМ ПО ВРЕМЕНИ 1277
Из соотношений (20), (21) несложно найти ωmax. Период, соответствующий частоте ωmax,
равен Tmin = 2π/ωmax. Теперь шаг сетки по времени можно задать формулой
δ = Cπ/ωmax,
(22)
в которой постоянная C зависит от конкретной задачи и, в конечном счёте, определяется
численными экспериментами. Для многих задач C ≈ 1.
Теперь оценим второе слагаемое в правой части равенства (17). Для производной функции
имеем
∂2P (p, t)
ω2max
P (p, t)|.
(23)
≈
∂t2
√
Далее, пусть Ω(ρ) - шар радиуса ρ =
3h/2 с центром в точке x = 0. Очевидно включение
Π(0) ⊂ Ω(ρ). Тогда, используя сферическую систему координат, получаем
∫
∫
1
1
3
|B(0)| =
r dr dS <
r dr dS =
h2.
(24)
4π
4π
8
Π(0)
Ω(ρ)
Теперь из соотношений (19), (22)-(24) вытекает, что
|B(0)|∂2P(p,t)
3π2C2
|n(p) - 1||U(p, t)|.
(25)
<
c2
∂t2
2
Таким образом, если значения постоянной C в оценке (25) достаточно малы, в правой части
равенства (17) можно пренебречь вторым слагаемым по сравнению с U(p, t) из левой части
равенства. Конкретные значения C, при которых это допустимо, определяются численными
экспериментами.
Для максимального расстояния D между пространственными узловыми точками в Q оче-
видно равенство
D = hmax|p - q|, p,q ∈ Q.
(26)
Значит, максимальное время запаздывания между узловыми точками Q равно T = D/c. Тог-
да из (19) следует, что максимальное число узловых точек по времени, в которых необходимо
вычислять функции, входящие в правую часть системы (17), определяется формулой
M = [T/δ] + 1 = [2max|p - q|] + 1, p,q ∈ Q,
(27)
где [·] - целая часть числа.
Теперь рассмотрим аппроксимацию производных функции по времени. Пусть задана два-
жды дифференцируемая функция P (t) и известны её значения в узловых точках t(m) = mδ,
где m - целые числа. Тогда значение производной функции в точке t, находящейся внутри
отрезка [t(m - 1), t(m + 1)], определяется разностной формулой
d2P(t)
P (m - 1) - 2P (m) + P (m + 1)
≈
(28)
dt2
δ2
В качестве центральной точки m целесообразно брать такую узловую точку, которая наи-
более близка к t. Отметим, что формула (28) даёт точное значение производной для любых
квадратичных функций P (t). Можно также использовать более точные формулы для аппрок-
симации, в которые будут входить значения P (m - 2), P (m - 1), P (m), P (m + 1), P (m + 2).
Однако в этом случае шаг по времени должен быть меньше и равняться δ = h/(3c).
Используя формулу (28), значение производных P (t) в точке τ можно приближённо пред-
ставить в виде
{
∑
d2P(τ)
1
[τ/δ],
если τ/δ ≤ [τ/δ] + 1/2,
≈
βjP(m(τ) + j), m(τ/δ) =
(29)
dτ2
δ2
[τ/δ] + 1,
если τ/δ > [τ/δ] + 1/2,
j=-1
где вид βj очевиден из (28).
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57
№9
2021
1278
САМОХИН и др.
Определим безразмерную переменную z формулой y = y(q) + zh. Тогда, учитывая, что
x(p) - y(p) = (p - q)h, из (17) получаем
∫
∫
∫
1
dz1 dz2 dz3
B(p - q) = -h2
= h2B0(p - q).
(30)
√∑3
4π
(pj - qj - zj )2
-1/2 -1/2 -1/2
j=1
Функция B0(p - q) является безразмерной. Если |p - q| ≫ 1, то можно использовать прибли-
жённую формулу B0(p - q) ≈ -1/(4π|p - q|).
Из равенств (16), (19) следует, что τ(n, p, q) = (n - 2|p - q|)δ. Далее, используя представ-
ления (29), (30) для выражений, входящих в сумму (17), получаем
∑
∂2P(q, τ)
B(p - q)
≈ 4B0(p - q)
βjP(q,n - m + j), m = m(2|p - q|).
(31)
∂τ2
j=-1
Вид функции m = m(2|p - q|) вытекает из (29).
Из (17), (31) следует, что последовательные по времени значения поля в узловых точках
области Q могут быть представлены следующей формулой:
∑∑
U (p, n) = U0(p, n) +
B(p - q, l)P (q, n - l),
l=1 q∈Q
P (q, n - l) = (n(q) - 1)U(q, n - l), q ∈ Q, n = 1, 2, . . .
(32)
Функции B(p - q, l) определяются из приведённых выше выражений, а значение M - форму-
лой (27). Очевидно, что B(0, l) = 0. Из (31) следует, что при фиксированных p, q функция
B(p - q, l) имеет только три ненулевых значения по l.
Основные вычислительные затраты в формуле (32) связаны с вычислением сумм. Обо-
значим
∑∑
W (p, n) =
B(p - q, l)P (q, n - l).
(33)
l=1 q∈Q
Для вычисления суммы (33) применим технику на основе быстрого дискретного преобразова-
ния Фурье [5, с. 93; 6]. Введём параллелепипед Π2 со сторонами 2N1h, 2N2h, 2N3h. Доопре-
делим функции B(p1, p2, p3, l) нулями в тех узловых точках параллелепипеда Π2, в которых
они не определены, и продолжим затем эти функции на все целочисленные значения p1, p2, p3
периодически с периодами 2N1,
2N2, 2N3 соответственно. Продолжим функции P (p, n - l)
на все узловые точки p ∈ Π2, полагая их равными нулю в узловых точках, не принадлежащих
области Q. Рассмотрим функцию
∑
∑
W (p, n) =
B(p - q, l)P (q, n - l).
(34)
l=1 q∈Π2
Учитывая изложенное, очевидно, что при p ∈ Q значение W (p1, p2, p3, n) из (34) совпадает
со значением W (p1, p2, p3, n) из (33). Теперь, применяя дискретное преобразование Фурье по
каждой пространственной переменной к обеим частям равенства (34), будем иметь
∑
Wk(k,n) =
BF (k,l)PF (k,n - l), k ∈ Π2.
(35)
l=1
Зная функцию (35), для того чтобы вычислить W (p, n) необходимо выполнить обратное дис-
кретное преобразование Фурье.
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57
№9
2021
ОБЪЁМНЫЕ ИНТЕГРАЛЬНЫЕ УРАВНЕНИЯ С ЗАПАЗДЫВАНИЕМ ПО ВРЕМЕНИ 1279
Оценим требуемое количество арифметических операций и объём памяти для хранения
массивов, которые необходимы для выполнения одного временного шага, определяемого вы-
ражением (32). Массивы B(p, l) и BF (k, l) вычисляются один раз до начала выполнения
алгоритма решения. После вычисления необходимо сохранять в памяти компьютера массив
чисел MB объёма
|MB| ≈ 8MN1, N2, N3, M = [2 max |p - q|] + 1, p, q ∈ Q.
(36)
Из (35) следует, что также необходимо хранить массив комплексных чисел MP объёма
|MP | ≈ 8MN1, N2, N3.
(37)
При вычислении акустического поля для очередного временного шага требуется, согласно (34)
и (35), выполнить дискретное преобразование Фурье только для функции P (p, n - 1) в парал-
лелепипеде Π2, поскольку преобразования Фурье функций в предыдущие моменты времени
известны. Таким образом, число Tt арифметических операций, требуемое для выполнения
одного временного шага, оценивается формулами
Tt ≈ 2TFF + TS, TS = 10MN1N2N3,
TFF ≈ 10N1N2N3(LOG (N1) + LOG(N2) + LOG(N3)).
(38)
В выражении (38) через LOG (N) обозначен целочисленный логарифм, т.е. сумма всех прос-
тых делителей числа N.
Обозначим N = N1N2N3. Тогда из (36)-(38) следует, что объём требуемой памяти и число
арифметических операций, необходимое для выполнения одного временного шага, практиче-
ски пропорциональны MN ≈ N4/3.
Важной характеристикой алгоритма является количество временных шагов, которое нуж-
но выполнить для получения полного решения поставленной задачи. Если внешнее поле огра-
ничено во времени, то критерием остановки алгоритма может быть следующий: значения аку-
стического поля в области Q близки к нулю.
Заключение. Рассмотренные в работе методы могут быть использованы для решения
многих важных классов задач нестационарной акустики. Отметим, что описанный выше ме-
тод решения с соответствующими изменениями может быть применён и к нестационарным
интегральным уравнениям, описывающим линейные материальные среды с временной и про-
странственной дисперсией.
Работа выполнена при финансовой поддержке Российского научного фонда (проект 20-11-
20087).
СПИСОК ЛИТЕРАТУРЫ
1. Самохин А.Б. Интегральные уравнения для нестационарных задач электродинамики в материаль-
ных средах // Дифференц. уравнения. 2002. Т. 38. № 9. С. 1288-1290.
2. Самохин А.Б., Самохина А.С., Кобаяси К. Численные методы решения нестационарного объёмного
сингулярного интегрального уравнения электродинамики // Дифференц. уравнения. 2019. Т. 55.
№ 9. С. 1293-1300.
3. Самохин А.Б. Объемные сингулярные интегральные уравнения для задач рассеяния на трехмерных
диэлектрических структурах // Дифференц. уравнения. 2014. Т. 50. № 9. С. 1215-1230.
4. Владимиров В.С. Уравнения математической физики. М., 1988.
5. Воеводин В.В., Тыртышников Е.Е. Вычислительные процессы с теплицевыми матрицами. М., 1987.
6. Сигов А.С., Андрианова Е.Г., Жуков Д.О., Зыков С.В., Тарасов И.Е. Квантовая информатика:
обзор основных достижений // Рос. технол. журн. 2019. Т. 7. № 1. С. 5-37.
Российский технологический университет (МИРЭА),
Поступила в редакцию 25.03.2021 г.
г. Москва,
После доработки 25.03.2021 г.
Институт проблем управления
Принята к публикации 08.06.2021 г.
им. В.А. Трапезникова РАН, г. Москва
9
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57
№9
2021
