Введение
Численные методы решения волновых уравнений иг-
рают важную роль не только в конкретных технических
приложениях, но и в фундаментальной науке в целом.
К таким методам относится и FDTD (Finite-Difference Time-
Domain) [1], некоторые особенности которого и являются
предметом данной работы.
Основное достоинство метода FDTD — простота реа-
лизации расчетного алгоритма. Именно это обуславлива-
ет широкое применение FDTD в самых разнообразных при-
ложениях: биологии и медицине [2–5], экологии, геологии
и минералогии [6, 7], оптике, фотонике, электронике, связи
и телекоммуникациях [8–13]. Кроме множества статей, так
или иначе связанных с методом FDTD, на эту тему имеется
и обширная учебная литература [14–17].
Несмотря на давнюю историю развития, внимание ис-
следователя по-прежнему привлекают фундаментальные
основы метода FDTD. К числу этих основ относится и во-
прос оценки корректности решений, полученных методом
FDTD в разнообразных постановках задачи для сигналов
с различной формой спектра [14, 17, 18].
Хорошо известно (см., например, учебную литерату-
ру [14–17]), что основным параметром, регламентирующим
Известия Коми научного центра Уральского отделения Российской академии наук № 5 (71), 2024
Серия «Физико-математические науки»
www.izvestia.komisc.ru
73
точность вычислений методом FDTD, является число Ку-
ранта, которое для 2D-случая (1 пространственное + 1 вре-
менное измерения) имеет вид
Sc =
cΔt
Δx
(1)
и связывает вместе основной физический параметр c —
скорость света в вакууме с численными параметрами за-
дачи Δx и Δt, определяющими шаг дискретизации про-
странства-времени.
Замечание 1. Качество численного решения определяется
не только выбором значения числа Куранта Sc, но и требо-
ванием к спектральному составу сигнала, а также уровнем
начального тока его источника (см. Утверждение 2 в нашей
предыдущей работе [18]). Именно влиянию последних двух
факторов и была посвящена статья [18], в которой число
Куранта было выбрано оптимальным образом для модели-
рования электромагнитных процессов в вакууме, а имен-
но — было положено, что Sc = 1. Данное исследование яв-
ляется логическим продолжением работы [18].
1. Мотивация и цель работы
Выбор Sc = 1 в общем случае не обеспечивает ка-
чество получаемого численного решения в однородных
недиспергирующих средах, оптически отличных от ваку-
ума. Это можно легко увидеть, изучая рис. 1 и 2.
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
P1
P2
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
Ez[m]
m
q = 130; "r = 1
q = 230; "r = 4
J q
q
w = 10; d = 30
Рисунок 1. Моделирование распространения импульсов гауссова вида
в вакууме P1 и диэлектрике P2 с относительной диэлектрической про-
ницаемостью εr = 4.
Figure 1. Simulating of the propagation of Gaussian form pulses in vacuum
P1 and dielectric P2 with relative permittivity εr = 4.
Рисунки 1 и 2 построены в соответствии с алгоритмом
Йи, подробно обсуждавшимся нами в предыдущей рабо-
те [18] (см. там формулы (13), (19)–(21)), с той разницей, что
теперь уравнения обновления (19)–(20) в явном виде учи-
тывают материальные параметры среды (а именно — ее
диэлектрическую εr и магнитную μr проницаемости), в ко-
торой распространяются сигналы
Hq+1
2
y

m +
1
2

= Hq−12
y

m +
1
2

+
+
Sc
ημr
(Eq
z [m + 1] − Eq
z [m]) ,
(2)
Eq+1
z [m] = Eq
z [m] − J q+12
[m]+
+
Scη
εr

Hq+12
y

m +
1
2

− Hq+12
y

m − 1
2

,
(3)
а простейшие поглощающие граничные условия (21) из [18]
заменены нами на более общие уравнения, построенные на
основе дифференциальных уравнений адвекции электро-
магнитного поля второго порядка точности. Все детали ре-
ализации данной схемы подробно изложены в статье [18]
и пособии [14].
Замечание 2. Моделирование работы направленного ис-
точника тока в формализме TF/SF в данной статье так-
же претерпело некоторые изменения по сравнению с [18]
и сводится к вычислению полей по схеме
Hq+12
y

s − 1
2

= Hq−1
2
y

s − 1
2

− Sc
ημr
Einc
z [0, q] , (4)
Eq+1
z [s] = Eq
z [s] +
√Sc
εrμr
Einc
z

−1
2
, q +
1
2

. (5)
Здесь (как и в работе [18]) s = 50 — это фиксированный
во всех численных экспериментах номер узла сетки, за-
дающий расположение в пространстве точечной антен-
ны, формирующей падающее поле Einc
z . Вычисление элек-
трического Ez[s] и магнитного Hy

s − 1
2

полей, соглас-
но (4), (5) по заданному типу падающей волны Einc
z , поз-
воляет имитировать работу источника тока J q, формиру-
ющего волну, излучаемую в правую область сетки m ⩾ s.
В дальнейшем везде под работой источника тока J q будем
иметь в виду именно эту схему.
􀀀0:5
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
P3
P4
􀀀0:5
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
Ez[m]
m
q = 150; "r = 1
q = 250; "r = 4
J q
q
lm = 50; d = 1
Рисунок 2. Моделирование распространения импульсов в форме вейвле-
та Рикера в вакууме P3 и диэлектрике P4 с относительной диэлектри-
ческой проницаемостью εr = 4.
Figure 2. Simulating of pulse propagation in the form of a Ricker wavelet in
vacuum P3 and dielectric P4 with relative permittivity εr = 4.
Параметры источников тока J q, формирующих сигна-
лы, представленные на рис. 1 и 2, также подробно описа-
ны нами в работе [18] (см. там формулы (32)–(34) и соот-
ветствующий текст) и выбраны такими, чтобы импульсы P1
и P3, распространяющиеся в вакууме, можно было считать
корректным численным решением задачи в смысле Опре-
деления 2, данного в [18]. Исходя из этого, импульсы P1
и P3 можно считать «опорными», сравнивая с которыми
74
Известия Коми научного центра Уральского отделения Российской академии наук № 5 (71), 2024
Серия «Физико-математические науки»
www.izvestia.komisc.ru
импульсы P2 и P4 соответственно, легко оценить влияние
среды на корректность численного решения задачи.
Рисунки 1 и 2 построены при одном и том же значении
числа Куранта Sc = 1, для которого в источнике [14] утвер-
ждается, что оно минимизирует численные ошибки FDTD-
расчета, а в пособии [17] заявлено, что такой выбор позво-
ляет получить точное решение задачи.
Вместе с тем, совершенно очевидно, что FDTD-решения,
изображенные импульсами P2 на рис. 1 и P4 на рис. 2, не
являются корректными в смысле Определения 2, данно-
го нами в работе [18], что находится в явном противоре-
чии с отмеченным в предыдущем абзаце. Эта некоррект-
ность есть результат явления, повсеместно возникающего
при численных расчетах методом FDTD и известного в ли-
тературе [14, 17] под названием «численная дисперсия».
Действительно, сигналы P2 на рис. 1 и P4 на рис. 2 не
представляют собой ни исходный гауссов импульс, ни вей-
влет Рикера соответственно. При распространении данных
волновых пакетов в диэлектрике с εr = 4 за достаточно
длительное время (Δq = 230 и 250 соответственно) их
форма существенно исказилась (отметим, что этот эффект
при использованных параметрах достаточно слаб и отчет-
ливо начинает проявляться только к концу моделирова-
ния). Это противоречит изначальной математической мо-
дели симулируемого нами явления, так как среда полага-
лась нами недиспергирующей. Это явное противоречие по-
рождает закономерный вопрос: может ли метод FDTD вооб-
ще использоваться для получения корректных результатов
моделирования распространения сигналов в недисперги-
рующих материалах? И если это возможно, то как на кор-
ректность получаемых решений влияет выбор числа Ку-
ранта Sc?
Несмотря на отмеченное выше противоречие, метод
FDTD был использован нами ранее в ряде работ, в кото-
рых либо рассматривалось распространение волн в слу-
чайно-неоднородных средах [12, 13], либо исследовались
особенности решения однородной и неоднородной задач
Коши методом FDTD [18]. И в первом, и во втором случаях
все проблемы, связанные с численной дисперсией, полно-
стью игнорировались, что отчасти обусловлено тем, что во
всех этих работах существенную часть трассы, вдоль ко-
торой рассматривалось распространение электромагнит-
ных сигналов, представляло собой воздушное простран-
ство (моделируемое нами неотличимым от вакуума). Вме-
сте с тем, подобное игнорирование данной стороны дела
не может быть полностью оправданным, что требует, если
не полной корректировки численной дисперсии, то хотя бы
уточнения эффектов, связанных с ней.
Далее отметим, что в книге [14] без каких-либо дока-
зательств постулируется, что значение Sc = 1 — макси-
мально возможное, что, вообще говоря, абсолютно не соот-
ветствует смыслу определения (1), которое не накладыва-
ет никаких ограничений на произвольно выбираемые про-
странственно-временные шаги сетки Δx и Δt и их соот-
ношение между собой. В источнике [17] более осторожно
утверждается, что выбор Sc > 1 приводит к экспоненци-
альному нарастанию шума, вызванному ошибками округле-
ния, всегда имеющими место при численном моделирова-
нии, что в конечном счете полностью разрушает получен-
ное решение.
Замечание 3. Следует также отметить, что в работах [14]
и [17] числами Куранта называются разные величины, что
создает дополнительные трудности для понимания всех
отмеченных вопросов. Это связано с тем, что в [14] в опре-
делении (1) величина c — это скорость света в вакууме (как
принято и нами в данной работе), в то время как в [17] счи-
тается, что c — это скорость распространения волны в дан-
ной конкретной среде.
Кроме уже отмеченных проблем, также возникает во-
прос о возможности применения метода FDTD для моде-
лирования электродинамики в «не самых обыкновенных
средах» (пусть и в пренебрежении явлением дисперсии).
Здесь имеются в виду следующие примеры: распростране-
ние радиоволн в ионосфере Земли [19], электромагнитные
волны терагерцового или рентгеновского диапазона в про-
водниках, полупроводниках или диэлектриках [20,21], сре-
ды с отрицательными значениями относительных прони-
цаемостей (в частности, «левые» среды [22–25]). Во всех
этих случаях решающую роль играет дисперсия электро-
магнитных волн, аномальный характер которой с матема-
тической стороны дела состоит в том, что проницаемости εr
и μr могут принимать значения, меньшие единицы, и даже
более того, — отрицательные.
При этом обычные алгоритмы, встречающиеся в лите-
ратуре по методу FDTD [1, 12–18], не позволяют моделиро-
вать распространение электромагнитных волн в таких сре-
дах, просто задавая значения 0 < εr, μr < 1 и εr, μr < 0.
Решению отмеченных вопросов и посвящена настоя-
щая статья, цель которой состоит, таким образом, в рас-
смотрии численной дисперсии метода FDTD и оценке вли-
яния выбора значений числа Куранта на качество модели-
рования распространения сигналов в однородных недис-
пергирующих средах.
2. Численная дисперсия в методе FDTD
Для получения выражения, описывающего численную
дисперсию в сетке Йи, вернемся к исходным дискретным
аналогам уравнений Ампера (3) и Фарадея (2), которые
в более компактной и удобной форме можно записать с по-
мощью операторов сдвига в пространственно-временной
сетке b Sχ
x и b Sτ
t (здесь χ и τ — параметры сдвига, имеющие
форму χ, τ = p/2, ∀p ∈ Z). Действие этих операторов по
определению имеет вид
b Sχ
x : b Sχ
x ψq[m] = ψq [m + χ] , (6)
b Sτ
t : b Sτ
t ψq[m] = ψq+τ [m], (7)
где ψ обозначает произвольную компоненту электромаг-
нитного поля (в рассматриваемом нами 2D-случае — Ez
или Hy).
Несложно проверить, что с помощью операторов (6) и (7)
уравнение Ампера (3) в отсутствии сторонних источников
(J = 0) может быть записано в форме
Известия Коми научного центра Уральского отделения Российской академии наук № 5 (71), 2024
Серия «Физико-математические науки»
www.izvestia.komisc.ru
75
b S 12
t ε
 
b S 12
t
− b S−12
t
Δt
!
Eq
z [m] =
= b S 12
t
 
b S 12
x − b S−12
x
Δx
!
Hq
y [m].
(8)
Здесь учтено определение числа Куранта (1), связь
скорости света в вакууме с его диэлектрической ε0
и магнитной μ0 проницаемостями c = 1/
√
ε0μ0 и вы-
ражение для характеристического импеданса вакуума
η =
p
μ0/ε0 ≈ 120π (для справки см., например, учеб-
ники [20, 21]). Кроме того, при записи (8) используется обо-
значение для абсолютной диэлектрической проницаемо-
сти среды ε = εrε0.
Определим для удобства также конечно-разностные
операторы ˜∂i согласно
˜∂i =
b S 12
i
− b S−1
2
i
Δi
, (9)
где i — это x или t. С помощью (9) закон Ампера в стацио-
нарной среде без сторонних токов (8) окончательно может
быть записан в форме Йи [14]
ε b S 1
2
t
˜∂tEq
z [m] = b S 12
t
˜∂xHq
y [m]. (10)
Таким же образом, с помощью (6), (7) и (9), дискретный
аналог закона Фарадея (2) приводится к форме Йи
μ b S 1
2
x ˜∂tHq
y [m] = b S 12
x ˜∂xEq
z [m]. (11)
Теперь рассмотрим плоскую монохроматическую вол-
ну частоты ω, распространяющуюся в сетке Йи вправо
(m ⩾ s)
Eq
z [m] = E0 ei(ωqΔt−˜βmΔx), (12)
Hq
y [m] = H0 ei(ωqΔt−˜βmΔx), (13)
где ˜ β — волновое число, т. е. постоянная распространения
плоской монохроматической волны в FDTD-сетке (отлича-
ющаяся от соответствующей константы β в непрерывном
пространстве), а E0 и H0 — комплексные амплитуды на-
пряженности электрического и магнитного полей.
Лемма 1. Действие операторов (9) на произвольную компо-
ненту ψ плоской монохроматической волны (12), (13) сво-
дится к
˜∂tψ = i
2
Δt
sin

ωΔt
2

ψ, (14)
˜∂xψ = −i
2
Δx
sin
 
˜ βΔx
2
!
ψ. (15)
Доказательство. Проверяется непосредственной подста-
новкой (12), (13) в (9) с учетом (6), (7). Выпишем здесь в явном
виде только действие на ψ операторов сдвига (6), (7) при
параметрах χ, τ = ±1
2
b S±12
t ψ = e±iωΔt/2ψ, b S±12
x ψ = e∓i ˜βΔx/2ψ, (16)
с помощью которых равенства (14) и (15) получаются эле-
ментарно.
Утверждение 1. Дисперсионное соотношение для сетки Йи
может быть представлено в форме
sin

ωΔt
2

=
√Δt
εμΔx
sin
 
˜ βΔx
2
!
. (17)
Доказательство. Используя результаты (14) и (15) Леммы 1,
а также (16) в законе Ампера (10) для плоской монохрома-
тической волны, запишем
iε
2
Δt
sin

ωΔt
2

eiωΔt/2Eq
z [m] =
= −i
2
Δx
sin
 
˜ βΔx
2
!
eiωΔt/2Hq
y [m].
(18)
Подставляя в последнее равенство явные выражения
полей (12), (13) для плоской монохроматической волны и со-
кращая общие множители, получаем
ε
1
Δt
sin

ωΔt
2

E0 = − 1
Δx
sin
 
˜ βΔx
2
!
H0. (19)
Отсюда приходим к выражению для численного импе-
данса в сетке Йи
E0
H0
= − Δt
εΔx
·
sin
˜ βΔx
2
sin
ωΔt
2
. (20)
Выполняя аналогичные действия применительно к за-
кону Фарадея (11), последовательно получаем
iμ
2
Δt
sin

ωΔt
2

e−i ˜βΔx/2Hq
y [m] =
= −i
2
Δx
sin
 
˜ βΔx
2
!
e−i ˜βΔx/2Eq
z [m],
(21)
μ
1
Δt
sin

ωΔt
2

H0 = − 1
Δx
sin
 
˜ βΔx
2
!
E0 (22)
и соответствующий импеданс в форме
E0
H0
= −μΔx
Δt
·
sin
ωΔt
2
sin
˜ βΔx
2
. (23)
Приравнивая правые части (20) и (23) и выполняя в по-
лучившемся равенстве перекрестное произведение со-
множителей, получаем
sin2

ωΔt
2

=
Δ2t
εμΔ2
x
sin2
 
˜ βΔx
2
!
. (24)
Извлекая квадратный корень в последнем равенстве,
окончательно приходим к (17), что и завершает доказатель-
ство.
76
Известия Коми научного центра Уральского отделения Российской академии наук № 5 (71), 2024
Серия «Физико-математические науки»
www.izvestia.komisc.ru
Замечание 4. Дисперсионное соотношение (17) метода
FDTD существенно отличается от своего непрерывного
аналога, имеющего вид [20, 21]
β = ω
√
εμ. (25)
Вместе с тем отметим, что это отличие становится исчезаю-
ще малым при достаточно малой дискретизации простран-
ства-времени. Действительно, сохраняя в рядах Тейлора
разложения синуса при Δt и Δx → 0 только слагаемые
первого порядка малости, из (17) легко получить
˜ β = ω
√
εμ. (26)
Подчеркнем, однако, что (26) справедливо только в пре-
деле бесконечно малой дискретизацииΔt,Δx → 0. В об-
щем же случае конечной дискретизации пространства-
времени фазовые скорости волны в FDTD-сетке Йи ˜cp
и в непрерывном пространстве cp будут отличны.
Утверждение 2. Отклонение фазовой скорости волны
в сетке Йи от соответствующего значения в непрерывном
случае может быть описано равенством
˜cp
cp
=
π
√
εrμr
Nλ arcsin
√
εrμr
Sc
sin

πSc
Nλ
, (27)
где параметр Nλ — есть число узлов пространственной
сетки на длину волны в свободном пространстве
λ = NλΔx. (28)
Доказательство. Используем связь фазовой скорости
с волновым числом в непрерывном пространстве [20, 21]
и FDTD-сетке
cp =
ω
β
, ˜cp =
ω
˜ β
, (29)
что позволяет привести левую часть равенства (27) к виду
˜cp
cp
=
β
˜ β
=
βΔx
2
˜βΔx
2
. (30)
Далее применим (25) к числителю последнего равенства
β = ω
√
εμ = 2π
c
λ
√
ε0εrμ0μr =
2π
λ
√
εrμr, (31)
где длина волны в вакууме λ дискретизуется согласно (28),
что дает
βΔx
2
=
π
√
εrμr
Nλ
. (32)
Применяя теперь к знаменателю правой части (30) дис-
персионное соотношение (17), в котором используются пре-
образования аналогичные (31), а также определение числа
Куранта (1) совместно с (28), получаем
˜ βΔx
2
= arcsin
√
εrμr
Sc
sin

πSc
Nλ

. (33)
Подстановка двух последних равенств в (30) приводит
к результату (27), что завершает доказательство.
Cоотношение (27) в рамках ограничений, рассматрива-
емых в настоящей статье (согласно которым исследуются
недиспергирующие однородные среды, для которых и εr
и μr — есть некоторые действительные константы), име-
ет очевидный смысл на следующей области определения
параметров:
Nλ ∈ N\1, εr, μr, Sc ∈ (0,+∞) ⊂ R, (34)
где ни множество натуральных чисел, больших едини-
цы N\1, ни множество положительных вещественных чи-
сел считаются не содержащими непосредственно эле-
мент +∞. Выход за рамки области определения (34) тре-
бует отдельного подробного исследования, и мы вернемся
к нему в последнем разделе данной статьи.
Прежде чем перейти к обсуждению особенностей вы-
ражения (27) отметим, что диэлектрическая и магнитная
проницаемости среды входят в него только в комбинации
nr =
√
εrμr, (35)
которая имеет физический смысл относительного показа-
теля преломления среды.
Пример 1. Рассмотрим распространение плоской монохро-
матической волны с Nλ = 10 в стекле с показателем пре-
ломления nr = 1.5 в случае числа Куранта Sc = 1. От-
ношение (27) при этом равно ˜cp/cp ≈ 0.9777, что в про-
центном отношении составляет численную ошибку пример-
но 2.23%. Таким образом, в данной ситуации на каждую
единицу пути, равную длине волны, FDTD-расчет накапли-
вает существенную фазовую ошибку порядка 8.03◦. Заме-
тим также, что улучшение дискретизации длины волны в
два раза (Nλ = 20) сокращает соответствующие ошибки
до значений 0.48% и 1.89◦ (т. е. примерно четырехкратно),
как это и должно происходить для вычислительного мето-
да второго порядка точности.
0:70
0:75
0:80
0:85
0:90
0:95
1:00
1:05
1:10
2 4 6 8 10 12 14 16 18 20
1
2
3
4
~cp/cp
N
Рисунок 3. Зависимость фазовой скорости, определяемой согласно мето-
ду FDTD по отношению к ее точному значению (27), от параметра дискре-
тизации длины волныNλ для некоторых сред с относительными показа-
телями преломления nr. Число Куранта Sc = 1.
Figure 3. Dependence of the phase velocity ratio, determined according to
the FDTD method with respect to its exact value (27), on the wavelength discretization
parameterNλ for some media with relative refractive indicesnr.
The Courant number Sc = 1.
Известия Коми научного центра Уральского отделения Российской академии наук № 5 (71), 2024
Серия «Физико-математические науки»
www.izvestia.komisc.ru
77
Особенности численной дисперсии метода FDTD, опре-
деляемые выражением (27), и дополняющие приведенный
выше пример, представлены на рис. 3, где показано семей-
ство кривых, для которых nr > Sc (а именно, nr =
√
2 для
кривой 1, nr =
√
3 — для линии 2 и nr = 2 — для 3 соот-
ветственно). Видно, что увеличение показателя преломле-
ния среды при плохой дискретизации длины волны приво-
дит к существенному запаздыванию волны, моделируемой
методом FDTD, что выражается в значительном отклонении
отношения ˜cp/cp от единицы.
Кроме того, на рис. 3 изображена кривая 4, отвеча-
ющая случаю nr < Sc (в данном конкретном случае вы-
брано значение nr =
p
1/2). Этот пример демонстрирует,
что в средах оптически менее плотных, чем вакуум, FDTD-
расчет приводит к распространению моделируемых волн
с опережением по сравнению с истинной скоростью.
Следствие 1. Из рис. 3 видно, что точность FDTD-расчета
падает с уменьшением Nλ, а кроме того — по мере откло-
нения проницаемости среды от единицы (что имеет место
и в диэлектриках, и в магнетиках, и в магнитных диэлек-
триках).
Следствие 2. Легко вычислить предел отношения (27), ко-
торый независимо от величин nr и Sc равен
lim
Nλ
→+∞
˜cp
cp
= 1. (36)
Это означает, что точность FDTD-расчета улучшается с ро-
стом Nλ.
Замечание 5. Ухудшение точности FDTD-расчета в среде
с nr > 1 связано с тем, что длина волны в такой среде
короче, чем в вакууме, и малой дискретизации длины вол-
ны Nλ оказывается недостаточно для корректного расче-
та. Однако простой перенос этого утверждения на случай
среды с εr < 1 не так очевиден и требует уточнения.
Следствие 3. Кроме того, рис. 3 указывает на то, что
высокочастотные компоненты волнового пакета в средах
с εr > 1 имеют тенденцию к запаздыванию, в то время
как в средах с εr < 1 эта особенность меняется на про-
тивоположную — высокочастотные компоненты волново-
го пакета в сетке Йи распространяются с большей фазо-
вой скоростью, чем это имеет место в действительности.
Другими словами, в средах оптически более плотных, по
сравнению с вакуумом, FDTD-расчет приводит к возникно-
вению численной дисперсии, имеющей характер аномаль-
ной дисперсии [26]. И наоборот, — при моделировании оп-
тически менее плотных сред, влияние FDTD-дискретиза-
ции соответствует поведению нормальной дисперсии (при
этом «красные» компоненты волнового пакета «обгоняют
синие»).
Замечание 6. Структура знаменателя (27) имеет явное
сходство с формой дисперсионных уравнений, получаемых
в задачах о распространении волн в средах с периодиче-
скими неоднородностями [19,27–29]. А именно, — здесь фи-
гурирует функция
φ(Nλ, Sc, nr) =
nr
Sc
sin
πSc
Nλ
, (37)
характер которой определяет полосы пропускания и не-
пропускания сетки Йи. Действительно, легко понять, что
область значений параметров (Nλ, Sc, nr), на которой вы-
полняется условие |φ(Nλ, Sc, nr)| > 1, отвечает непро-
пусканию волн, так как выражение arcsin φ при этом не
имеет смысла в вещественных значениях.
Подобная картина явления имеет место в модели Кро-
нига-Пенни, а также в уравнениях Хилла и Матье, описы-
вающих параметрические колебания [19, 27–29], а также
распространение волн в системах с пространственной пе-
риодичностью в расположении неоднородностей. В каче-
стве периодической структуры в нашем случае (в методе
FDTD) выступает непосредственно сама расчетная сетка
Йи. При этом пропускание и непропускание волн (кванто-
во-механический аналог — разрешенные и запрещенные
зоны) определяется как раз величинами Nλ, Sc и nr.
0
20
40
60
80
100
50 100 150 200 250 300 350
Sc
N
􀀀1
0
1
A1
A2
B
Рисунок 4. (Nλ, Sc)-диаграмма полос пропускания сетки Йи
при nr = 100.
Figure 4. (Nλ, Sc)-diagram of Yee grid bandwidths when nr = 100.
Рисунок 4 иллюстрирует (Nλ, Sc)-диаграмму по-
лос пропускания сетки Йи, построенную на основе
функции (37) для сравнительно большого значения по-
казателя преломления среды nr = 100. Области по-
лос непропускания {Ai}N
i=1, соответствующие условию
|φ(Nλ, Sc, nr)| > 1, изображены на этой диаграмме од-
нотонной заливкой серого цвета. Общее число таких об-
ластей N (на рис. 4 отмечены первые две из них) зависит
от величины показателя преломления среды nr и растет
с увеличением последнего. Полосы пропускания представ-
лены на данной диаграмме градиентной заливкой, меня-
ющейся от черного (φ = +1) к белому (φ = −1) цвету.
Именно в рамках последних областей имеет смысл дис-
персионное соотношение в форме (27).
Также на рис. 4 серым цветом выделена область B,
определяемая неравенством Sc > nc, в рамках которой
метод FDTD также не может обеспечивать корректного чис-
ленного решения задачи о распространении волны в одно-
родной недиспергирующей среде. Аргументы, поддтвер-
ждающие справедливость данного утверждения, приведе-
ны в следующем разделе статьи при обсуждении рис. 6.
3. Связь численной дисперсии с числом Ку-
ранта
Теперь, после проведенного предварительного иссле-
дования, обсудим алгоритм коррекции численной диспер-
78
Известия Коми научного центра Уральского отделения Российской академии наук № 5 (71), 2024
Серия «Физико-математические науки»
www.izvestia.komisc.ru
сии с помощью специального выбора числа Куранта. Его
возможность опирается на соотношение (27), из которого
непосредственно вытекает следующий результат.
Утверждение 3. При любых заданных εr и μr, принадле-
жащих области определения (34), всегда можно устранить
вычислительные ошибки, связанные с численной диспер-
сией, задавая число Куранта равным
Sc =
√
εrμr. (38)
Доказательство. Проверяется непосредственной подста-
новкой (38) в (27).
Замечание 7. Значение числа Куранта (38) можно назвать
«магическим», так как в этом случае фазовая скорость вол-
ны в FDTD-сетке ˜cp точно совпадает с реальным значени-
ем фазовой скорости cp, независимо от выбранной дискре-
тизации длины волны Nλ. В то же время в источнике [14]
«магическим» называется значение Sc = 1, что справед-
ливо только для вакуума, но никак не в общем случае.
Здесь же отметим, что наш выбор (38) эквивалентен зада-
нию Sc = 1, используемому в книге [17].
Несколько примеров, иллюстрирующих идею приме-
нения «магического» числа Куранта (38) для коррекции
программной дисперсии в недиспергирующих однородных
средах, ограниченных выбором параметров (34), приведе-
ны на рис. 5.
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
P5 P1 P6
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
q = 130; nr = 1:0
q = 80; nr = 2:0
Ez[m] q = 155; nr  0:7
m
w = 10; d = 30
J q
q
Рисунок 5. Моделирование распространения импульсов гауссова вида
в средах, оптически более P5 (nr = 2) и менее P6 (nr =
p
1/2) плот-
ных, чем вакуум, с применением коррекции программной дисперсии (38).
Figure 5. Simulating of Gaussian pulse propagation in media optically
more P5 (nr = 2) and less dense P6 (nr =
p
1/2) than vacuum, using
correction of numeric dispersion (38).
Рисунок 5 подтверждает идею, сформулированную
в виде Утверждения 3 — численная дисперсия здесь дей-
ствительно не имеет места ни в случае среды оптически
более плотной по сравнению с вакуумом (импульс P5), ни
в случае оптически менее плотной среды (импульсP6). Ри-
сунок 5 следует сопоставить с рис. 1, на котором численная
дисперсия отчетливо проявляется на форме импульса P2.
Замечание 8. Отметим также, что распространение элек-
тромагнитного поля в пространстве со временем (опреде-
ляемым, как и всегда, величиной дискретного индекса q),
представленное на рис. 5 для импульса P5, происходит
с видимым опережением аналогичного процесса для им-
пульса P2 на рис. 1 в два раза. Следует подчеркнуть, что
это опережение — только кажущееся, так как его причи-
ной является именно двукратное отличие в числах Куран-
та, выбранных в случае рис. 1 и 5. В действительности же
(при учете поправки на Sc), все временны́е характристи-
ки сигналов на рис. 1 и 5 идентичны. Другими словами, ес-
ли считать пространственный шаг решетки Δx фиксиро-
ванным параметром, не меняющимся при переходе от мо-
делирования в случае Sc = 1, представленного на рис. 1,
к моделированию с числом Куранта (38), то «цена» каждо-
го временного шага Δt (а вместе с ним, и полная продол-
жительность моделирования) изменится в
√
εrμr раз. Тот
же эффект (с точностью до замены слов «опережение» на
«запаздывание») имеет место и в случае оптически менее
плотных сред (см. импульс P6) на рис. 5.
Замечание 9. Также отдельно укажем на то, что модели-
рование распространения импульса P6 в среде оптически
менее плотной, чем вакуум, в соответствии с расчетным ал-
горитмом (2) – (5) при выборе числа Куранта Sc = 1 прин-
ципиально невозможно. Несложно убедиться в том, что вы-
числения в этом случае очень быстро приводят к расходи-
мостям, не давая никакой полезной информации, и прин-
ципиально не описывают данный частный случай.
Последнее замечание прекрасно иллюстрирует рис. 6,
на котором показан лавинообразный процесс накопления
численных ошибок в процессе расчета по алгориму (2) –
(5) в случае превышения числа Куранта Sc, по сравнению
с nr, всего на 0.1%. Данная иллюстрация построена при
тех же параметрах источника сигнала J q, что и на рис. 1
и 5, для случая вакуума, хотя принципиально такая же кар-
тина имеет место и в случае любых других однородных
недиспергирующих сред из области определения (34).
􀀀1
0
1
0 25 50 75 100 125 150 175 200
􀀀1
0
1
􀀀1
0
1
q = 50
Eq
z [m]
m
q = 100
Eq
z [m]
q = 130
Eq
z [m]
Рисунок 6. Временная динамика разрушения численного решения, по-
лучаемого в процессе вычислений согласно алгоритму (2) – (5) при
Sc − nr = 10−3.
Figure 6. Time dynamics of the numerical solution destruction obtained
during computation according to the algorithm (2) – (5) when
Sc − nr = 10−3.
Как видно из рис. 6, к моменту времени q = 130 уро-
вень шума, произвольно возникающего в результате вы-
числений по схеме (2) – (5) при Sc − nr = 10−3, двухкрат-
но превышает величину полезного сигнала по абсолютному
Известия Коми научного центра Уральского отделения Российской академии наук № 5 (71), 2024
Серия «Физико-математические науки»
www.izvestia.komisc.ru
79
значению. При этом пространственная протяженность ча-
сти сетки Йи, затронутой данным шумом, составляет 3/4 от
всей длины расчетной области. Последнее обстоятельство
существенно искажает форму заднего фронта моделиру-
емого полезного сигнала. Далее с последующим развити-
ем процесса при q > 130 полезное решение разрушается
полностью.
Дополнительные численные эксперименты, выполнен-
ные нами, также показывают, что эффекты, связанные
с накоплением ошибок, вследствие самовозбуждения сет-
ки, имеют место всегда при выполнении неравенства
Sc > nr. Так, при Sc − nr = 10−4 самовозбуждение сет-
ки становится заметным уже после прохождения полезного
сигнала, однако его лавинообразный характер наблюдает-
ся и в этом случае. Все эти наблюдения можно обобщить
в форме следующего утверждения.
Утверждение 4. Расчетный алгоритм (2) – (5) быстро рас-
ходится и не может быть использован для получения кор-
ректного численного решения задачи о распространении
сигналов в недиспергирующих однородных средах при
выборе Sc > nr.
Доказательство. Исчерпывающая аргументация справед-
ливости данного утверждения на «физическом уровне
строгости» приведена выше при обсуждении рис. 6.
Применением данного утверждения как раз и обясня-
ется наличие области B, показанной на рис. 4, хотя с фор-
мальной точки зрения функция (37) при этом не превышает
единицу по модулю.
Замечание 10. В то же время отметим, что выбор значений
Sc ⩽ nn, согласованный с условием |φ(Nλ, Sc, nr)| ⩽ 1,
не приводит к столь драматическим последствиям, кото-
рые имеют место в Утверждении 4.
Важным следствием проведенного таким образом по-
дробного исследования, основные результаты которого
сдержатся в Утверждениях 3 и 4, является еще один факт.
Следствие 4. Задание числа Куранта в форме (38) — един-
ственно возможный оптимальный выбор с точки зрения
применения численного FDTD-алгоритма (2) – (5) для мо-
делирования распространения сигналов в недиспергиру-
ющих однородных средах.
Доказательство. Оптимальность такого выбора обуслов-
лена тем, что он устраняет численную дисперсию при
моделировании сигналов любой формы и спектрального
состава, не противоречащих Утверждению 2 работы [18],
распространяющихся в широком классе недиспергирую-
щих однородных сред, описываемых областью определе-
ния (34). Единственность вытекает из Утверждения 4.
Замечание 11. Проблема, оставшаяся нерешенной до конца
на данный момент (возникающая при внимательном изуче-
нии применения алгоритма (2) – (5) к однородным недис-
пергирующим средам, оптически отличным от вакуума),
состоит в том, что в данном случае помимо основного сиг-
нала с избранной направленностью всегда наблюдается
сигнал сравнительно малой амплитуды обратной направ-
ленности. Назовем для краткости этот последний сигнал —
обратным Pb, в противоположность исходному — прямо-
му Pf . Иллюстрация этой проблемы приведена на рис. 7,
который построен при тех же параметрах источника сиг-
нала, что и рис. 1 и 5.
Рисунок 7 демонстрирует формирование обратного им-
пульса Pb при моделировании распространения сигнала
гауссовой формы в средах, оптически менее плотных, чем
вакуум при оптимальном выборе числа Куранта. Видно, что
характер этого явления существенно зависит от величи-
ны относительного показателя преломления nr среды. При
стремлении последнего к нулю этим явлением уже нель-
зя пренебречь, так как оно приводит к искажению харак-
теристик и полезного прямого сигнала Pf . В то же время,
в случае nr > 10−1 соответствующая численная ошибка
сравнительно мала, и может быть оценена не превышаю-
щей уровня в 1%.
􀀀1
0
1
0 25 50 75 100 125
Pb Pf
􀀀1
0
1
􀀀1
0
1
1 2
E20
z [m]
m
1 2
E40
z [m]
1 2
E60
z [m]
Рисунок 7. Динамика формирования прямого Pf и обратного Pb им-
пульсов при моделировании распространения сигнала гауссовой формы
в средах, оптически менее плотных, чем вакуум при оптимальном вы-
боре числа Куранта. Штрих-пунктирная кривая 1 соответствует случаю
Sc = nr = 10−1, сплошная линия 2 — Sc = nr = 10−2.
Figure 7. Dynamics of forward Pf and backward Pb impulses formed in
simulating the propagation of a Gaussian-shaped signal in media optically
less dense than vacuum under optimal choice of the Courant number. The
dashed-dotted curve 1 corresponds to the case of Sc = nr = 10−1, the
solid line 2 — Sc = nr = 10−2.
4. Границы применимости основных резуль-
татов
В завершении данной работы обсудим вопрос о приме-
нимости полученных здесь результатов за рамками обла-
сти определения (34).
Во-первых, укажем на то, что интервал значений
0 < εr, μr < 1, входящий в область определения (34), ис-
пользуемый при выполнений условий Утверждения 3, уже
сам по себе расширяет диапазон применимости развитого
в данной работе расчетного алгоритма на область экзоти-
ческих сред, обычно не рассматриваемых в литературе.
Во-вторых, отметим, что наш численный алгоритм (2) –
(5) остается справедливым и при расширении области
опредления (34) на интервал εr, μr ∈ (−∞,+∞) при
условии одновременной отрицательности проницаемостей
среды εrμr > 0.
Среды с одновременно отрицательными проницаемо-
стями, как известно [23–25], называются левыми. В таких
средах могут распространяться обратные волны [30], опре-
деляемые тем фактом, что для них скалярное произведе-
ние волнового вектора k (в нашей работе всюду k = βex)
80
Известия Коми научного центра Уральского отделения Российской академии наук № 5 (71), 2024
Серия «Физико-математические науки»
www.izvestia.komisc.ru
и вектора Умова-Пойнтинга S — отрицательно
(k · S) < 0, (39)
где поток энергии, переносимый волной, и определяемый
вектором S, равен [20–23, 26]
S =
c
4π
[E ×H] . (40)
Справедливость такого обобщения демонстрирует
рис. 8, на котором представлены результаты сравнения
волновых пакетов, формируемых тем же источником сиг-
нала, что и на рис. 1, 5 – 7, зафиксированные в одинаковый
момент времени q = 100, при распространении в правой
среде с εr = μr = +1 (нижняя половина рисунка) и левой
среде с εr = μr = −1 (верхняя половина).
􀀀1
0
1
0 25 50 75 100 125 150 175 200
􀀀1
0
1
Ez
Hy
E80
z [m]; H80
y [m]
m
Ez
Hy
E80
z [m]; H80
y [m]
Рисунок 8. Мгновенные снимки прямых (нижняя половина) и обратных
(верхняя часть) волн, распространяющихся в сетке Йи согласно расчет-
ному алгоритму (2) – (5).
Figure 8. Instantaneous snapshots of forward (bottom half) and backward
(top half) waves propagating in the Yee grid according to the computational
algorithm (2) – (5).
Легко убедиться в том, что для волнового пакета, по-
казанного в верхней части рис. 8 вектор Умова-Пойнтин-
га (40) S ↑↓ ex, что как раз и определяет обратную волну
согласно условию (39).
И наконец, поясним, что справедливость алгоритма (2) –
(5) при расширении области опредления (34) на интер-
вал εr, μr ∈ (−∞,+∞) при условии неодновременной
отрицательности проницаемостей среды εrμr < 0 в рам-
ках данной работы нами не исследовалась. Это связано
с тем, что при εrμr < 0 постоянная распространения (25)
оказывается мнимой величиной, что соответствует силь-
ному затуханию волн в таких средах, которые вследствии
этого оказываются сильно диспергирующими (для них рас-
пространение плоских волн оказывается невозможным, а
вместе с этим теряет простой смысл Утверждение 2, кото-
рое требует иной формулировки в данном случае). Все это
выходит за рамки темы данной работы.
Заключение
Таким образом, в данной работе исследовано явле-
ние численной дисперсии при FDTD-моделировании рас-
пространения электромагнитных сигналов в недисперги-
рующих однородных средах. Сформулированы несколько
утверждений, определяющих характер этой дисперсии, а
также описано влияние на нее числа Куранта. Определено
оптимальное значение числа Куранта, устраняющее чис-
ленную дисперсию волновых пакетов. Исследованы гра-
ницы применимости разработанного метода моделирова-
ния, и впервые указано на возможность его применения
для сред, оптически менее плотных, чем вакуум, а также
для левых сред.
Вместе с тем, тема исследования все еще остается до-
статочно интересной, так как в данной статье не исследо-
вались многие интересные вопросы, часть из которых при-
ведена ниже.
Некоторые открытые вопросы
1. Чем может быть объяснено с физической точки зре-
ния ухудшение точности FDTD-расчета в средах
с nr < 1? Объяснение аналогичное ситуации в сре-
дах, для которых nr > 1, и основанное на уменьше-
нии длины волны в них (и, соответственно, плохой ее
дискретизации) здесь не годится.
2. В чем причина формирования обратного импуль-
са Pb и как его устранить? Возможно, он является
следствием каких-то аналитических ошибок, допу-
щенных при записи (2) – (5) либо неустранимых оши-
бок машинного представления чисел с плавающей
запятой при компьютерных вычислениях.
3. Возможно ли корректировать численную диспер-
сию при моделировании распространения сигналов
в неоднородных средах? Можно ли этого добиться,
изменяя число Куранта динамически в процессе мо-
делирования?
Этот список интересных вопросов ни в коей мере не
претендует на полноту и вполне может быть расширен.
Авторы заявляют об отсутствии конфликта интересов.

1. Yee, K. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media / K. Yee // IEEE Trans. on Ant. and Prop. – 1966. – Vol. 14, № 3. – P. 302–307.

2. Miyazaki, Y. FDTD analysis of spatial filtering of scattered waves for optical CT of medical diagnosis / Y. Miyazaki, K. Kouno // IEEJ Trans. FM. – 2009. – Vol. 129, № 10. – P. 693–698.

3. Tan, T. Single realization stochastic FDTD for weak scattering waves in biological random media / T. Tan, A. Taflove, V. Backman // IEEE Trans. AP. – 2013. – Vol. 61, № 2. – P. 818–828.

4. Stark, J. Light scattering microscopy measurements of single nuclei compared with GPU-accelerated FDTD simulations / J. Stark [et al.] // Phys. Med. Biol. – 2016. – Vol. 61, № 7. – P. 2749–2761.

5. Nzao, A. B. S. Analysis and FDTD Modeling of the Influences of Microwave Electromagnetic Waves on Human Biological Systems / A. B. S. Nzao // Open Journal of Applied Sciences. – 2022. – Vol. 12. – P. 912–929.

6. Glubokovskikh, S. Seismic monitoring of CO2 geosequestration: CO2CRC Otway case study using full 4D FDTD approach / S. Glubokovskikh [et al.] // International Journal of Greenhouse Gas Control. – 2016. – Vol. 49. – P. 201–216.

7. Yu, J. Modeling of Whole-Space Transient Electromagnetic Responses Based on FDTD and its Application in the Mining Industry / J. Yu, R. Malekian, J. Chang, B. Su // IEEE Trans. Indust. Inform. – 2017. – Vol. 13, № 6. – P. 2974–2982.

8. Fantoni, A. A model for the refractive index of amorphous silicon for FDTD simulation of photonics waveguides / A. Fantoni, P. Loureniço, M. Vieira // International Conference on Numerical Simulation of Optoelectronic Devices (NUSOD), Copenhagen, Denmark. – 2017. – P. 167–168.

9. Mishra, C. S. FDTD approach to photonic based angular waveguide for wide range of sensing application / C. S. Mishra [et al.] // Optik. – 2019. – Vol. 176. – P. 56–59.

10. Mohanty, S. P. FDTD method to photonic waveguides for application of optical demultiplexer at 3-communication windows / S. P. Mohanty, S. K. Sahoo, A. Panda, G. Palai // Optik. – 2019. – Vol. 185. – P. 146–150.

11. Bakirtzis, S. FDTD-Based Diffuse Scattering and Transmission Models for Ray Tracing of Millimeter-Wave Communication Systems / S. Bakirtzis, T. Hashimoto, C. D. Sarris // IEEE Trans. AP. – 2021. – Vol. 69, № 6. – P. 3389–3398.

12. Makarov, P. Simulation of Electromagnetic Wave Propagation in Magnetic Randomly Inhomogeneous Magnetic Media / P. Makarov [et al.] // IEEE Magnetics Letters. – 2022. – Vol. 13. – P. 1–5.

13. Макаров, П. А. Моделирование распространения электромагнитных волн в магнитно-неоднородных средах / П. А. Макаров, В. А. Устюгов, В. И. Щеглов // Известия Коми научного центра Уральского отделения Российской академии наук. Серия «Физико-математические науки». – 2022. – № 5 (57). – C. 100–105.

14. Schneider, J. B. Understanding the Finite-Difference Time-Domain Method / J. B. Schneider. – www.eecs.wsu.edu/~schneidj/ufdtd, 2010. – 403 p.

15. Inan, U. S. Numerical electromagnetics. The FDTD method / U. S. Inan, R. A. Marshall. – Cambridge: Cambridge University Press, 2011. – 406 p.

16. Taflove, A. Advances in FDTD computational electrodynamics photonics and nanotechnology / A. Taflove, A. Oskooi, S. G. Johnson. – Boston : Artech House, 2013. – 639 p.

17. Langtangen, H. P. Finite Difference Computing with PDEs: A Modern Software Approach / H. P. Langtangen, S. Linge. – Springer Cham, 2017. – XXIII. – 507 p.

18. Макаров, П. А. Особенности численного моделирования уравнений Максвелла методом FDTD в однородной и неоднородной формулировках / П. А. Макаров, В. А. Устюгов, В. И. Щеглов // Известия Коми научного центра Уральского отделения Российской академии наук. Серия «Физико-математические науки». — 2023. — № 4 (62). – C. 96–107.

19. Виноградова, М. Б. Теория волн / М. Б. Виноградова, О. В. Руденко, А. П. Сухоруков. – Москва : Наука, 1979. – 384 с.

20. Бредов, М. М. Классическая электродинамика / М. М. Бредов, В. В. Румянцев, И. Н. Топтыгин. – Москва : Наука, 1985. – 400 с.

21. Кугушев, А. М. Основы радиоэлектроники. Электродинамика и распространение радиоволн / А. М. Кугушев, Н. С. Голубева, В. Н. Митрохин. – Москва : Изд-во МГТУ им. Н. Э. Баумана, 2001. – 368 с.

22. Шустер, А. Введение в теоретическую оптику / А. Шустер. – Ленинград, Москва : ОНТИ, гл. ред. общетех. лит., 1935. – 376 с.

23. Веселаго, В. Г. Электродинамика веществ с одновременно отрицательными значениями ε и μ / В. Г. Веселаго // УФН. – 1967. – T. 92, № 3. – C. 517–526.

24. Pendry, J. Negative refraction / J. Pendry // Contemporary Physics – 2004. – V. 45, № 3. – C. 191–202.

25. Агранович, В. М. Пространственная дисперсия и отрицательное преломление света / В. М. Агранович, Ю. Н. Гартштейн // УФН – 2006. – T. 176, № 10. – C. 1052–1068.

26. Ландсберг, Г. С. Оптика / Г. С. Ландсберг. – Москва : ФИЗМАТЛИТ, 2010. – 848 с.

27. Коткин, Г. Л. Лекции по аналитической механике / Г. Л. Коткин, В. Г. Сербо, А. И. Черных. – Москва, Ижевск : НИЦ РХД, 2017. – 236 с.

28. Карлов, Н. В. Колебания, волны, структуры / Н. В. Карлов, Н. А. Кириченко. – Москва : Физматлит, 2008. – 498 с.

29. Флюгге, З. Задачи по квантовой механике, Т. I / З. Флюгге. – Москва : Мир, 1974. – 342 с.

30. Шевченко, В. В. Прямые и обратные волны: три определения, их взаимосвязь и условия применимости / В. В. Шевченко // УФН. – 2007. – T. 177, № 3. – C. 301–306.