Введение
Ни для кого не является секретом актуальность и вос-
требованность в науке и технических приложениях совре-
менных численных методов, позволяющих просто, нагляд-
но и информативно моделировать процессы распростра-
нения электромагнитных сигналов в различных условиях.
Одним из таких методов является техника FDTD (Finite-
Difference Time-Domain) — метод численного решения вол-
новых уравнений, разработанный К. Йи в 1966 г. [1], осно-
ванный на аппроксимации производных отношением ко-
нечных разностей
df
dx
≈ Δf
Δx
.
Данная аппроксимация выполняется путем некоторо-
го «огрубления задачи», состоящего в дискретизации
пространства-времени, при котором дифференциальные
уравнения для непрерывных функций заменяются конеч-
но-разностным уравнением для функций дискретных пе-
ременных.
Эта схема в настоящее время нашла широкое примене-
ние в различных областях науки, техники и самых разно-
образных приложениях. Так, известны примеры использо-
вания метода для решения биологических и медицинских
задач [2–5], а также проблем в области экологии, геологии,
минералогии и геологоразведки [6, 7]. Достаточно очевид-
ным, однако нисколько не потерявшим своей актуально-
сти, является применение метода FDTD в оптике, фотони-
ке, электронике, связи и телекоммуникациях [8–11]. Помимо
непосредственно научной литературы, существует множе-
ство учебных изданий [12–14], посвященных методу FDTD.
96
Известия Коми научного центра Уральского отделения Российской академии наук № 4 (62), 2023
Серия «Физико-математические науки»
www.izvestia.komisc.ru
Такое широкое использование метода, очевидно, обос-
новывается большим числом его достоинств, немалое зна-
чение среди которых имеет простота реализации расчет-
ного алгоритма. Конечно, ценой этой простоты являет-
ся некоторая ресурсоемкость выполняемых вычислений,
определяемая заметными затратами машинного времени
и памяти. На первых порах это существенно ограничива-
ло применение метода FDTD, однако в настоящее время,
в связи со значительным прогрессом вычислительной тех-
ники, эти ограничения потеряли свое значение. Еще од-
ним важным преимуществом рассматриваемого метода яв-
ляется возможность исследования сигналов со сложной
структурой спектра в ходе однократного численного экс-
перимента. В связи с этим в электродинамических зада-
чах FDTD обычно используют для моделирования поведе-
ния импульсных сигналов, однако часто рассматривается
и квазигармонический режим нагрузки передающей ли-
нии [12].
Приведенный выше краткий обзор источников может
создать ошибочное впечатление о том, что метод FDTD яв-
ляется очень хорошо разработанным инструментом с дос-
конально изученной научной базой. Последнее не совсем
корректно, так как, несмотря на обширную библиографию
по методу FDTD, его тщательное обоснование и всесторон-
нее исследование по сей день не исчерпало себя. Во мно-
гом это обуславливается тем, что в имеющейся литературе
практически не рассматриваются вопросы строгой оценки
корректности решений, полученных методом FDTD в раз-
личных постановках задачи.
Во всех известных авторам данной статьи работах ме-
тод FDTD используется для решения уравнений Максвел-
ла с источниками (т.е. в форме неоднородных дифферен-
циальных уравнений). Вместе с тем имеется возможность
применять этот метод и в иной постановке задачи, пере-
формулируя ее в виде задачи Коши для однородной систе-
мы уравнений. При этом естественным образом возникают
вопросы эквивалентности данных формулировок и оцен-
ки корректности решений, полученных в их рамках. Кроме
того, численное решение уравнений Максвелла в однород-
ной формулировке требует строгого отбора корректных на-
чальных условий для соответствующей задачи Коши и от-
сев недопустимых конфигураций электромагнитного поля.
Более того, даже при традиционной схеме решения
в неоднородной формулировке часто возникают вопросы
строгой оценки корректности применения метода FDTD для
сигналов со сложной формой спектра. Этот факт имеет ме-
сто быть даже для гармонических сигналов, однако в ли-
тературе (включая наши предыдущие работы [18,19] по этой
теме) он освещается недостаточно подробно.
Анализу отмеченных вопросов и посвящена настоя-
щая работа. Основные ее результаты представляют собой
две гипотезы, сформулированные в форме утверждений
в основном тексте статьи. Строгое математическое дока-
зательство этих утверждений не приводится, однако об-
суждаются их мотивировка и обоснование справедливости
на «физическом уровне строгости». Кроме того, приводит-
ся аргументация, подкрепленная численными примерами,
подтверждающими сделанные выводы.
1. Основные формулировки
Микроскопические уравнения Максвелла в вакууме,
записанные в форме дифференциальных уравнений для
трехмерных вектор-функций в системе единиц СИ, имеют
следующий вид [15–17]:
∇ × E = −μ0
∂H
∂t
, ∇ ×H = Jext + ε0
∂E
∂t
, (1)
∇ · E =
ρext
ε0
, ∇ ·H = 0. (2)
Здесь E и H — это напряженности электрического и маг-
нитного полей, ε0 и μ0 — диэлектрическая и магнитная
проницаемости вакуума, а плотность заряда ρext и плот-
ность тока Jext играют роль сторонних источников элек-
тромагнитного поля.
Будем решать задачу в отсутствии сторонних зарядов
ρext = 0. (3)
Отбрасывая уравнения (2), дающие в принятых нами
условиях (3) статическое решение для полей, и вводя обо-
значения (здесь использована краткая запись ∂t = ∂/∂t)
ψ =

E
H

, b L =

∇× μ0∂t
∇× −ε0∂t

, φ =

0
Jext

, (4)
законы Фарадея и Максвелла-Ампера (1) можно записать в
компактной матрично-операторной форме
b Lψ = φ. (5)
Следуя общепринятой классификации [17], решениями
уравнений (5) могут быть как нестационарные, так и ква-
зистационарные или стационарные поля. С математиче-
ской точки зрения выражение (5) представляет собой со-
кращенную запись неоднородной системы линейных диф-
ференциальных уравнений в частных производных перво-
го порядка. Очевидно, что в отсутствии сторонних токов
Jext = 0 (6)
выражение (5) обращается в однородную систему
b Lψ = 0, (7)
решениями которой теперь могут быть исключительно
нестационарные поля.
В дальнейшем везде, как для однородных (7), так и для
неоднородных (5) формулировок, будем решать задачи
Коши с начальными условиями ψ0 = ψ(r, 0). Обозначим
символом ¯ ψ(r, t) точное решение соответствующей зада-
чи Коши для неоднородной системы (5), которое формаль-
но можно записать в виде
¯ ψ(r, t) = b L−1
ψ0
φ. (8)
Определение 1. Будем называть произвольную функ-
цию f(r, t) функцией ограниченной длительности, если
существует такой момент времени τ , после достижения ко-
торого она тождественно обращается в нуль
f(r, t) ≡ 0, ∀ r, t > τ.
Известия Коми научного центра Уральского отделения Российской академии наук № 4 (62), 2023
Серия «Физико-математические науки»
www.izvestia.komisc.ru
97
Утверждение 1. Всякая задача Коши для неоднородной
системы (5) с источниками φ(r, t) ограниченной длитель-
ности τ эквивалентна задаче Коши для однородной систе-
мы (7) с начальным условием
ψ0 = ¯ ψ(r, t − τ ).
Эквивалентность этих задач состоит в том, что их ре-
шения идентичны.
Обсуждение данного утверждения отложим вплоть до
раздела 4, так как для осмысления заложенного в него со-
держания предварительно требуется обсудить множество
деталей, изложенных в последующих разделах.
2. Общая схема численного FDTD-решения
В основных чертах материал данного раздела следу-
ет нашим предыдущим работам [18, 19] и классическим ру-
ководствам [12–14], однако для полноты картины и ясно-
сти понимания воспроизведем здесь все необходимые вы-
кладки.
С целью упрощения изложения и максимальной на-
глядности полученных результатов будем рассматривать
одномерный случай. Это нисколько не умаляет общности
обсуждаемых результатов, так как их масштабирование на
случаи большей размерности не вносит ничего принципи-
ально нового.
Итак, пусть электромагнитное поле распространяется
вдоль декартовой осиOx, а Jext = J(x, t)ez. В этих усло-
виях E = Ez(x, t)ez и H = Hy(x, t)ey, а основная си-
стема уравнений (5) принимает вид
∂Ez
∂x
− μ0
∂Hy
∂t
= 0,
∂Hy
∂x
− ε0
∂Ez
∂t
= J. (9)
Дискретизируем пространство-время и выполним пе-
реход от непрерывных функций к их точечным аналогам
согласно преобразованию
f(x, t) → fq[m] = bD f(x, t) = f(mΔx, qΔt), (10)
где Δx, Δt — фиксированные шаги пространственно-
временной сетки, а m, q ∈ N0 — соответствующие ин-
дексы, определяющие ее узел. Обратное преобразование
fq[m] → f(x, t) будем обозначать символом bD −1.
Далее, следуя классической работе [1] и используя пре-
образование (10), введем две сетки для электрическо-
го Eq
z [m] и магнитного Hq±1
2
y

m ± 1
2

полей, смещенные
по отношению друг к другу в шахматном порядке. Это поз-
воляет записать конечно-разностный аналог уравнений (9)
в форме
Hq+1
2
y

m +
1
2

= Hq−12
y

m +
1
2

+
+
Δt
μ0Δx
(Eq
z [m + 1] − Eq
z [m]) ,
(11)
Eq+1
z [m] = Eq
z [m] − Δt
ε0
Jq+12
[m]+
+
Δt
ε0Δx

Hq+12
y

m +
1
2

− Hq+12
y

m − 1
2

.
(12)
Равенства (11) и (12) отличает от аналогов, полученных
нами ранее [18, 19], то, что они записаны для вакуума. В
англоязычной литературе [12–14] за уравнениями (11) и (12)
закрепилось название Update Equations, отражающее тот
факт, что при заданных начальных условиях

E0
z
H
−1
2
y

= bD ψ0, (13)
данные уравнения позволяют вычислять временную дина-
мику электромагнитного поля рекуррентно.
На практике систему уравнений (11), (12) используют
в эквивалентной, но несколько более общей, форме. Для
этого вместо одновременного задания и Δx, и Δt регла-
ментируют величину их отношения числом Куранта [12]
Sc =
cΔt
Δx
, (14)
где скорость света в вакууме c, как известно [15–17], свя-
зана с параметрами ε0 и μ0 следующим образом
c =
√ 1
ε0μ0
. (15)
Замечание 1. Выбор Sc = 1 определяет временной шаг
сетки Δt точно соответствующим времени распростране-
ния света в вакууме между двумя соседними узлами сет-
ки, расположенными на расстоянии Δx друг относитель-
но друга. Несмотря на кажущуюся оптимальность, такой
выбор заведомо не является корректным даже в вакуу-
ме при моделировании достаточно высокочастотных сиг-
налов, период которых T ⩽ Δt.
С помощью (14), (15) несложно привести коэффициенты
уравнений (11), (12) к виду
Δt
μ0Δx
=
Sc
η
,
Δt
ε0Δx
= Scη, (16)
где использовано обозначение характеристического импе-
данса вакуума
η =
r
μ0
ε0
≈ 377 Ом. (17)
Замечание 2. Введение импеданса вакуума связано с вы-
бором системы физических единиц и теории измерений
в широком смысле. Как известно, в системе СГС подоб-
ной сущности в принципе не возникает. Поэтому исполь-
зование η — есть следствие того, что мы выполняем ис-
следование, используя систему СИ. Величина η влияет ис-
ключительно на взаимное отношение между электриче-
ской и магнитной компонентами электромагнитного поля,
однако на динамику (и физику в принципе) это, конечно,
не влияет. Переход к единицам СГС выполнить достаточно
просто. При этом нет необходимости вводить импеданс ва-
куума, а число Куранта (14) трансформируется в параметр,
имеющий размерность обратной скорости и определяемый
как отношение шагов сетки
Sc =
Δt
Δx
.
98
Известия Коми научного центра Уральского отделения Российской академии наук № 4 (62), 2023
Серия «Физико-математические науки»
www.izvestia.komisc.ru
В подавляющем большинстве случаев амплитуда и знак
функции источника J не имеют значения. Это связано
с тем, что при вычислении сечения рассеяния, коэффи-
циента отражения и других подобных величин всегда ис-
пользуется нормировка по падающему полю. Поэтому нет
необходимости явно задавать коэффициентΔt/ε0 в (12) —
достаточно считать его содержащимся в самой функции
источника. Таким образом, можно сделать переход
Δt
ε0
Jq+1
2 → J q+12
. (18)
Окончательно, с учетом всех принятых нами соглаше-
ний, дискретные аналоги уравнений (5) принимают вид
Hq+12
y

m +
1
2

= Hq−12
y

m +
1
2

+
+
Sc
η
(Eq
z [m + 1] − Eq
z [m]) ,
(19)
Eq+1
z [m] = Eq
z [m] − J q+12
[m]+
+ Scη

Hq+1
2
y

m +
1
2

− Hq+12
y

m − 1
2

.
(20)
Далее будем рассматривать электромагнитное поле
в некоторой активной области — ограниченной обла-
сти пространстваX ⊂ R → m ∈ [0,N − 1]. Здесь и да-
лееN — это число узлов электрической и магнитной сеток,
отвечающих данной активной области.
Для простоты моделирования, чтобы не учитывать мно-
гократные переотражения электромагнитного поля от гра-
ниц активной области, будем использовать поглощающие
граничные условия (ABC — Absorbing Boundary Conditions),
которые в рассматриваемом нами одномерном случае мо-
гут быть записаны в виде:
Eq+1
z [0] = Eq
z [1],
Hq+1
2
y

N − 1
2

= Hq−12
y

N − 3
2

.
(21)
Таким образом, состояния полей на границах рассмат-
риваемой нами области не фиксированы жестко, а динами-
чески изменяются в соответствии с тем, что «к ним подхо-
дит изнутри». За счет этого реализуется «убегание» поля
из активной области за ее пределы, что моделирует пове-
дение электромагнитного поля в открытом пространстве.
Уравнения (13), (19)–(21) позволяют реализовать итера-
ционный алгоритм Йи [1, 12–14], определяющий всю дина-
мику электромагнитного поля ψqY
[m]. Зададимся теперь
вопросом определения корректности такого решения. Для
этого требуется оценить справедливость равенства
ψqY
[m] = bD ¯ ψ(x, t), (22)
где правая часть понимается как результат дискретиза-
ции (10) решения задачи, точного в смысле (8).
Очевидно, что строгое выполнение равенства (22)
принципиально недостижимо в силу неизбежных оши-
бок численного счета, поэтому требуется более реали-
стичная процедура оценки корректности решения ψqY
[m].
Здесь следует упомянуть следующее. В наших более ран-
них работах [18, 19] отмечалось, что «успех реализации»
сформулированного метода во многом определяется «пра-
вильно подобранными» параметрами пространственно-
временной сетки, для которых приводились следующие
оценки
Δx ≪ λ, Δt ≪ 1/ν, (23)
где λ — длина волны, а ν — частота исследуемого сигнала.
На текущий момент нами установлено, что выполнение
неравенств (23) не является ни достаточным, ни необходи-
мым условием «успешности» численного решения уравне-
ний (5). Примеры, подкрепляющие этот тезис, показаны на
рис. 1 и 2, подробное описание и анализ которых приводит-
ся нами в следующем разделе. Данные примеры отчетливо
демонстрируют, что требуется как более точная формули-
ровка самого понятия «успешности» решения ψqY
[m], по-
лученного в результате численного расчета по схеме (13),
(19)-(21), так и более корректный критерий его достижения.
3. Задача Коши для неоднородного уравнения
Как следует из вышесказанного, традиционно уравне-
ния Максвелла (5) численно решают следуя схеме Йи (13),
(19)-(21). С точки зрения математики данная методика
представляет собой итерационный алгоритм решения дис-
кретного аналога задачи Коши для системы неоднородных
уравнений. При этом на практике интересуются именно ре-
акцией системы на действие источников J → ψqY
[m], по-
этому во всех известных нам публикациях полагают на-
чальное состояние системы невозбужденным, т.е. выбира-
ют
ψ0 ≡ 0. (24)
Вообще говоря, задание именно такого тривиально-
го начального условия не является обязательным. Вме-
сте с тем нами установлено, что далеко не любой выбор
ψ0 ̸= 0 приводит к корректному решению задачи. В связи
с этим отложим обсуждение данного вопроса до раздела 4,
а в рамках этого пункта будем исследовать именно тради-
ционную схему решения с начальным условием (24).
Определение 2. Будем называть решениеψqY
[m], получен-
ное в ходе итерационного алгоритма Йи (13), (19)-(21), кор-
ректным с точностью ϵ > 0 в заданной области X × T =
{(x, t)|x ∈ X, t ∈ T}, где X ⊂ R, а T ⊆ R+ ∪ 0, в слу-
чае выполнения условия

bD
−1ψqY
[m] − ¯ ψ(x, t)
¯ ψ(x, t)

⩽ ϵ. (25)
Менее строго, будем называть решение ψqY
[m] корректным
тогда, когда

bD
−1ψqY
[m] − ¯ ψ(x, t)

≪ ¯ ψ(x, t). (26)
Теперь, прежде чем формулировать критерий кор-
ректности решения ψY[m], приведем несколько конкрет-
ных примеров, демонстрирующих особенности поведе-
ния ψY[m] при различных источниках сигнала J . Во всех
Известия Коми научного центра Уральского отделения Российской академии наук № 4 (62), 2023
Серия «Физико-математические науки»
www.izvestia.komisc.ru
99
наших примерах будем полагать, что источник поля — то-
чечный, он расположен в узле сетки с номером m = 50
и генерирует сигналы ограниченной длительности τ . Кроме
того, будем считать, что формируемые источником сигна-
лы обладают определенной направленностью и для опре-
деленности распространяются вправо (m ⩾ 50). Для про-
граммной реализации такой схемы возбуждения необхо-
димо использовать технику, известную в литературе под
названием TF/SF — Total-Field/Scattered-Field (см., напри-
мер, [12–14]). Мы не будем останавливаться на ее описа-
нии, поскольку для основного содержания статьи это не
так важно, а все необходимые детали можно найти в ука-
занных источниках.
3.1. Квазигармонические сигналы. Рассмотрим сначала
сигналы с наиболее простым спектром, а именно: пусть ис-
точник «включается» в момент времени t = 0 и ток в нем
представляет собой синусоидальные функции времени
Js(x, t) = δ(x − 50Δx)Asin

2πct
λ0

(27)
и меняющиеся со временем по закону косинуса
Jc(x, t) = δ(x − 50Δx)Acos

2πct
λ0

. (28)
Здесь δ(x) — дельта-функция Дирака, указывающая по-
ложение источника, A — амплитуда создаваемого им
сигнала. Во всех наших численных примерах будем по-
лагать A = 1. Далее, при достижении времени τ , ис-
точник «выключается», и ток в нем обращается в нуль
J(x, t)|t>τ ≡ 0.
Процессы включения/выключения источника матема-
тически можно смоделировать, умножая функции (27)
и (28) на комбинацию функций Хевисайда вида θ(t)[1 −
θ(t − τ )]. Здесь первый множитель регулирует момент
включения источника, второй — определяет время его вы-
ключения. В связи с тем, что источник, по нашему предпо-
ложению, работает ограниченное время, сигналы, созда-
ваемые им, не являются строго монохроматическими, од-
нако всегда можно выбрать τ ≫ Δt. Это позволит прене-
бречь влиянием такого рода немонохроматичности и счи-
тать всю систему работающей в квазигармоническом ре-
жиме. Для достижения этого зададим τ = 100Δt.
Теперь, после перечисления всех предварительных
деталей реализации схемы (13), (19)-(21), обсудим ее
непосредственные результаты, полученные для источни-
ков (27) и (28).
Так, на рис. 1 показан результат численного расчета для
синусоидального источника (27), генерирующего сигнал
с основной длиной волны в спектре λ0 = 3Δx. Очевидно,
что для этого случая сильное условие (23) заведомо не вы-
полняется. Вместе с тем из рисунка следует, что результаты
моделирования в этой ситуации можно считать корректны-
ми с точностью не хуже, чем ϵ = 0.14. Важно отметить, что
результаты дискретизации J q[m] = bD J(x, t) тока (27) в
этих условиях корректны с той же степенью точности, что
и сами результаты моделирования, представляющие собой
состояние электромагнитного поля.
􀀀1:0
􀀀0:5
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
􀀀1:0
􀀀0:5
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
Ez[m]
m
q = 110
J q
q
Рисунок 1. Результаты моделирования для синусоидального точечного
источника (27) при λ0 = 3Δx (все остальные параметры расчета по-
дробно описаны в основном тексте).
Figure 1. Simulation results for a sinusoidal point source (27) at λ0 = 3Δx
(all other calculation parameters are described in detail in the main text).
Подтверждениями описанных выше фактов являются
следующие результаты. На верхней части рис. 1 хорошо
виден сигнал, имеющий синусоидальную форму и рас-
пространяющийся вправо так, что его передний фронт
к моменту времени q = 110 успел преодолеть дистанцию
Δm = 110. При этом область пространства, занимаемая
волновым пакетом, соответствует ровно ℓ = 33λ0, как это
и должно быть при избранных нами параметрах моделиро-
вания.
Кривые, изображенные на рис. 1, строго говоря, не яв-
ляются синусоидальными сигналами. Причины этого легко
понять, если учесть, что при выбранной нами длине вол-
ны λ0 = 3Δx и числе Куранта Sc = 1 временные отсче-
ты дискретизируют задающий сигнал всего три раза за его
период T . Вкупе с неизбежными ошибками округления при
численном счете это как раз и приводит к отмеченному по-
ведению J q. Стоит, однако, подчеркнуть, что на характе-
ре самого поля ψqY
[m] указанные ошибки дискретизации
тока (27) сказываются не столь существенно, как в слу-
чае источника (28). Решающим здесь является именно про-
цесс начального развития сигнала, который хоть и являет-
ся достаточно резким (см. начальный участок нижней кри-
вой J q, изображенный на рис. 1, в сравнении с передним
фронтом поля Ez вблизи узлаm = 160 на верхней кривой
того же рисунка), но эта резкость не столь принципиальна,
так как является не характеристикой самого сигнала (см.
обсуждения, связанные с рис. 2, и утверждение 2), а лишь
следствием не самой удачной его дискретизации.
Выполненные нами дополнительные расчеты показы-
вают, что с увеличением λ0 корректность (в смысле опре-
деления 2) схемы решения (13), (19)-(21) для синусоидаль-
ного источника (27) растет (а ϵ, соответственно, уменьша-
ется), так как при этом все более справедливыми становят-
ся условия (23). Кроме того, твердо установлено, что при
λ0 ⩽ Δx ни о какой корректности решения говорить не
приходится.
После столь подробного анализа результатов модели-
рования для синусоидального источника (27) перейдем
к обсуждению особенностей работы источника (28). На
100
Известия Коми научного центра Уральского отделения Российской академии наук № 4 (62), 2023
Серия «Физико-математические науки»
www.izvestia.komisc.ru
рис. 2 показаны итоги расчетов (13), (19)-(21) в случае, когда
λ0 = 100Δx. При этом условия (23) заведомо выполняют-
ся, и результаты расчета должны быть существенно лучше
тех, что представлены на рис. 1.
􀀀2:0
􀀀1:0
0:0
1:0
2:0
0 25 50 75 100 125 150 175 200
􀀀1:0
􀀀0:5
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
Ez[m]
m
q = 110
J q
q
Рисунок 2. Результаты моделирования для точечного источника (28) при
λ0 = 100Δx.
Figure 2. Simulation results for a point source (28) at λ0 = 100Δx.
В действительности все обстоит ровно наоборот —
мгновенные значения решения ψqY
[m], представленные
на рис. 2, не могут быть признаны корректными в смыс-
ле определения 2, так как представляют собой некоторый
«пилообразный» сигнал (как мы увидим впоследствии, это
достаточно типичная картина, возникающая в ходе FDTD-
моделирования). Более того, выделяется и то, что сигнал
«проникает» в область m < 50, в которой его не должно
быть согласно принятой нами схеме TF/SF.
Вместе с тем заметим, в данном конкретном случае ре-
зультат усреднения мгновенных значений ψqY
[m], кото-
рый можно получить как полусумму его верхней и нижней
огибающих, вполне соответствует поведению ожидаемого
точного решения ¯ ψq[m] и может быть признан корректным.
Так ли это в случае произвольного источника — вопрос на
данный момент открытый и требует более подробного изу-
чения.
Попытаемся понять причину столь разительного отли-
чия результатов моделирования в рассмотренных нами
случаях. Для этого отметим, что в некоторых источниках
(см., например, [12]) упоминается влияние на результаты
моделирования постоянной составляющей в спектре ис-
следуемого сигнала. Результатом Фурье-преобразования
˜ f(ω) = b Ff(t) =
√1
2π
∞ Z
−∞
f(t) e−iωt dt (29)
временно́й части сигналов (27) и (28) являются следующие
обобщенные функции
˜ Js(ω) = iA
r
π
2
[δ(ω + ω0) − δ(ω − ω0)] , (30)
˜ Jc(ω) = A
r
π
2
[δ(ω + ω0) + δ(ω − ω0)] . (31)
Обращает на себя внимание тот факт, что в правых ча-
стях выражений (30) и (31) фигурируют дельта-функции
Дирака. Это является следствием того, что мы вычисляли
преобразование Фурье (29) для неограниченных во вре-
мени (и потому строго монохроматических) функций (27)
и (28). Учет ограниченной длительности данных сигна-
лов (другими словами, их квазигармоничности) приводит
к трансформации дельта-функций в функции более слож-
ного вида. В частности, первое слагаемое выражения (30)
преобразуется к виду
iτ
2
√
2π
Aexp

−i
(ω + ω0)τ
2

sinc
(ω + ω0)τ
2
.
Здесь использовано достаточно стандартное обозначе-
ние функции sinc x = sin(x)/x. Для остальных слагае-
мых в (30) и (31) имеют место аналогичные замены, которые
не выписаны полностью в виду их громоздкости.
Для нас существенно то, что как для Фурье-трансфор-
мант (30) и (31), так и для более точных их аналогов, мож-
но с уверенностью утверждать отсутствие постоянной со-
ставляющей сигнала, т. е. b F[Js,c(t)](0) = 0.
3.2. Импульсы и негармонические сигналы. Теперь об-
судим особенности моделирования сигналов с более слож-
ным спектральным составом. Начнем со случая отдельных
импульсов, в котором ограничимся только двумя примера-
ми: импульсами гауссова вида
Jg(t) = Aexp
"
−(t − t0)2
w2
#
(32)
и импульсами в форме вейвлета Рикера
Jr(t) = A

1 − 2 [πνm(t − t0)]2

×
× exp

−[πνm(t − t0)]2

.
(33)
Здесь, как и всегда, амплитуда A = 1, а t0, w и νm — есть
параметры сигналов (32) и (33), определяющие их форму.
Величина t0 в обоих случаях задает основную временную
задержку, т. е. тот момент времени, к которому достигает-
ся J(t0) = A. Значение w регламентирует основную про-
должительность гауссова сигнала, так что его подавляю-
щая мощность источника тока сосредоточена во времен-
ном интервале [t0 − w, t0 + w]. Нам будет удобно вы-
ражать этот параметр в единицах шагов временной сетки
w = wΔt. Величина νm задает ту частоту в спектре вей-
влета Рикера, на которую приходится его максимум. Для
целей сопоставления всех результатов этого раздела нам
удобнее задавать не νm, а отвечающую ей длину волны
λm = c/νm, которую можно выразить в единицах шагов
пространственной сетки λm = lmΔx.
Для гауссова импульса временную задержку t0 можно
непосредственно задавать в форме t0 = dgΔt, в то время
как для вейвлета Рикера ее удобнее представлять в еди-
ницах обратной частоты, множителем dr, определенным
согласно:
t0 = dr
1
νm
= dr
λm
c
= dr
lm
Sc
Δt. (34)
Известия Коми научного центра Уральского отделения Российской академии наук № 4 (62), 2023
Серия «Физико-математические науки»
www.izvestia.komisc.ru
101
Последнее равенство здесь записано с помощью опреде-
ления числа Куранта (14). В дальнейшем нижний индекс
у параметра d в явном виде указывать не будем, так как
его смысл виден из контекста.
На рис. 3 показаны результаты моделирования соглас-
но схеме (13), (19)-(21) для источника (32), формирующего
направленные вправо гауссовы импульсы с разной шири-
ной w. На верхней части рисунка показаны мгновенные
распределения поля Ez[m] в пространстве для обоих им-
пульсов, вычисленные для удобства в такие моменты вре-
мени q, чтобы дистанции Δm, пройденные данными сиг-
налами, отличались в два раза. Нижняя часть того же ри-
сунка иллюстрирует временную динамику токаJ q, форми-
рующего эти импульсы.
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
P2 P1
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
Ez[m]
m
w = 10; q = 130
w = 20; q = 80
J q
q
w = 10
w = 20
Рисунок 3. Результаты моделирования для точечного источника (32) при
d = 30.
Figure 3. Simulation results for a point source (32) at d = 30.
Беглого взгляда, обращенного на данный рисунок,
вполне достаточно для утверждения, что результаты чис-
ленного расчета электромагнитного поля для сигнала
меньшей ширины w вполне корректны. В то же время ре-
зультат моделирования поля, созданного более широким
импульсом тока, проявляет те же особенности, что были от-
мечены нами при анализе работы источника (28). А именно,
мгновенные значения поля ψqY
[m], рассчитанные в этом
случае, демонстрируют «пилообразное поведение», так что
истинную форму сигнала можно восстановить только как
полусумму его огибающих.
Вместе с тем, как было отмечено нами ранее, форма
гауссова сигнала определяется не только его шириной w,
но и задержкой d. Влияние этого параметра можно оценить
изучая рис. 4, на котором представлены результаты рас-
чета для гауссова источника (32), формирующего сигналы
одинаковой шириной w = 10, но с разной задержкой d.
На рис. 4 показано, что уменьшение временной задерж-
ки d гауссова импульса при том же значении его шириныw
ухудшает корректность результатов расчета. Это ухудше-
ние восприниматся не столь значительным, поскольку из-
менение значения d от 30 до 20 приw = 10 пусть и увели-
чивает начальный ток J 0
g примерно в 100 раз, но сама его
величина остается сравнительно небольшой и составляет
≈ 2 % от максимального значения. Заметим, что даль-
нейшее сокращение d, например до значения 10, ухуд-
шит корректность результатов расчета сигнала P3 силь-
нее, чем для более широкого импульса P2, представлен-
ного на рис. 3.
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
P3 P1
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
Ez[m]
m
d = 30; q = 130
d = 20; q = 70
J q
q
d = 30
d = 20
Рисунок 4. Результаты моделирования для точечного источника (32) при
w = 10.
Figure 4. Simulation results for a point source (32) at w = 10.
Таким образом, рассмотренные примеры приводят к
выводу о том, что на корректность результатов численного
расчета по алгоритму Йи (13), (19)-(21) самым существенным
образом влияет именно величина начального тока J 0 ис-
точника (на рис. 3 этот уровень отмечен небольшой верти-
кальной стрелкой при нулевой абсциссе).
Заметим, что спектр гауссова сигнала
˜ Jg(ω) =
w √
2
Aexp

−
ωw
2
2
− iωt0

, (35)
вообще говоря, имеет достаточно существенный вклад по-
стоянной составляющей. Поэтому для подтверждения вы-
водов, сформулированных в предыдущем абзаце, необхо-
димо исключить влияние данного фактора. При этом удоб-
но исследовать импульсы как раз в форме вейвлета Рике-
ра, поскольку их спектральный состав
˜ Jr(ω) =
ω2
2
√
2π3ν3m
A×
× exp
"
−

ω
2πνm
2
− iωt0
# (36)
не имеет постоянной составляющей.
На рис. 5 представлены результаты численного расче-
та динамики электромагнитного поля согласно алгоритму
Йи (13), (19)-(21) для источника (33), формирующего направ-
ленные вправо импульсы в форме вейвлета Рикера с раз-
ной длиной волны, отвечающей максимуму в спектре lw.
Как и всегда до этого, верхняя часть рисунка представля-
ет собой распределение электрического поля в простран-
стве, тогда как нижняя — дает представление о развитии
со временем тока, являющегося источником сигнала.
102
Известия Коми научного центра Уральского отделения Российской академии наук № 4 (62), 2023
Серия «Физико-математические науки»
www.izvestia.komisc.ru
􀀀0:5
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
P5 P4
􀀀0:5
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
Ez[m]
m
lm = 50; q = 150
lm = 3; q = 53
J q
q
lm = 50
lm = 3
Рисунок 5. Результаты моделирования для точечного источника (33) при
d = 1.
Figure 5. Simulation results for a point source (33) at d = 1.
Параметры lw, использованные для двух источников,
действие которых показано на рис. 5, подобраны так, что-
бы формируемые ими импульсы существенно отличались.
Так, при малом значении lw = 3 (при этом λm ≳ Δx,
но λm ̸≫ Δx) пространственная протяженность соответ-
ствующего волнового пакета настолько узка, что оказыва-
ется сопоставимой с пространственным шагом сетки, одна-
ко его спектральный состав достаточно широк. В противо-
положность этому, для большого значения lw = 50 все об-
стоит ровно наоборот. Существенно, что для обоих импуль-
сов временная задержка выбрана одинаковой и единич-
ной d = 1, так что изначально (в момент времени q = 0)
оба импульса практически неотличимы от нуля, а макси-
мума достигают к моментам времени q = 3 и q = 50 соот-
ветственно.
􀀀1:0
0:0
1:0
2:0
0 25 50 75 100 125 150 175 200
􀀀0:5
0:0
0:5
1:0
0 25 50 75 100 125 150 175 200
Ez[m]
m
d = 1; q = 150
d = 0; q = 150
J q
q
d = 1
d = 0
Рисунок 6. Результаты моделирования для точечного источника (33) при
lm = 50.
Figure 6. Simulation results for a point source (33) at lm = 50.
Как и для рис. 3 и 4, мгновенные значения электри-
ческого поля на верхней части рис. 5 приведены для мо-
ментов времени, когда импульсы P4 и P5 успели пройти
дистанции, отличающиеся друг от друга в два раза. При
этом очень хорошо видно, что результаты моделирования
являются корректными в обоих случаях, и это обстоятель-
ство никак не связано со спектральным составом обоих
импульсов. Совершенно иначе обстоит дело тогда, когда
существенным оказывается начальное значение тока ис-
точника J 0
r
̸= 0. Добиться этого можно уменьшая значе-
ние задержки d. Сопоставление результатов моделирова-
ния с критическим случам d = 0 показано на рис. 6, все
остальные параметры на котором заданы такими же, как
и для импульса P4 на рис. 5.
Рис. 6 прекрасно иллюстрирует те выводы, которые по-
степенно формулировались при обсуждении результатов
моделирования для всех предыдущих источников сигнала.
Можно утверждать, что рис. 6 в некотором смысле комби-
нирует результаты, представленные ранее на рис. 1 и 2.
3.3. Результаты раздела. Дополнительно нами установ-
лено, что все особенности численного решения задачи Ко-
ши для неоднородного уравнения (5) с начальным услови-
ем (24) согласно схеме (13), (19)-(21), рассмотренные в дан-
ном пункте, остаются справедливыми и для других часто
встречающихся сигналов, таких как треугольные импуль-
сы, меандр и т. д. Таким образом, можно сформулировать
основной результат данного раздела в следующей форме.
Утверждение 2. Необходимым и достаточным условием
корректности численного решения ψqY
[m] задачи Коши
для неоднородного уравнения Максвелла в вакууме (5)
с некоторым начальным условием ψ0, полученного в ходе
итерационного алгоритма Йи (13), (19)-(21), являются следу-
ющие ограничения:
1. Шаг пространственной сеткиΔx⩽λmin (где λmin —
это минимальная длина волны в спектре сигнала, ге-
нерируемого источником);
2. Начальный ток источника не имеет существенных
скачков по сравнению с заданной начальной конфи-
гурацией электрического поля J 0[m] ≪ E0
z [m].
Замечание 3. Требование Δt ⩽ 1/νmax, вообще говоря,
излишне, так как оно обеспечивается автоматически при
соблюдении первого условия Утверждения 2 и разумном
выборе числа Куранта (Sc ⩽ 1).
Замечание 4. Характер спектра сигналов, генерируемых
источником тока, не влияет на корректность моделиро-
вания в рассматриваемом нами случае распространения
в вакууме. Вместе с тем спектральный состав (и в особен-
ности его постоянная составляющая b F [φ(t)] (0)) сигна-
лов — один из ключевых факторов, влияющих на коррект-
ность моделирования в случае их генерации и распростра-
нения в той или иной среде. Этот вопрос очень обширен
и требует отдельного рассмотрения.
4. Задача Коши для однородного уравнения
До этого момента с помощью алгоритма Йи мы исследо-
вали поведение электромагнитного поля ψqY
[m], создава-
емого источником тока J q[m] при нулевой начальной кон-
фигурации поля. С математической точки зрения решалась
задача Коши для неоднородного уравнения Максвелла (5)
с тривиальным начальным условием (24). Теперь посмот-
рим на задачу с иной точки зрения. А именно, будем счи-
тать то решение ψqY
[m], которое нами получено для неод-
нородной задачи Коши к моменту завершения действия
ограниченного по времени источника, начальным состоя-
нием поля ψ0[m] = ψqY
[m]

q⩾τ/Δt.
Известия Коми научного центра Уральского отделения Российской академии наук № 4 (62), 2023
Серия «Физико-математические науки»
www.izvestia.komisc.ru
103
Поскольку оператор b L, определенный согласно (4), аб-
солютно одинаков как для однородного, так для неод-
нородного уравнения, то совершенно ясно, что решение
сформулированной таким образом задачи Коши для одно-
родного уравнения полностью эквивалентно решению со-
ответствующей задачи Коши для неоднородного уравне-
ния. Именно это и имелось в виду при формулировании
утверждения 1, причем единственным необходимым усло-
вием здесь является только ограниченность во времени
источника поля.
4.1. Конфигурации полей, соответствующие бегущим
волнам. Наиболее просто продемонстрировать эквива-
лентность задач Коши, о которых говорилось выше, можно
в случае конфигураций электромагнитного поля с той или
иной избранной направленностью. В частности, в преды-
дущем разделе с помощью техники TF/SF мы везде моде-
лировали волны, распространяющиеся в вакууме вправо от
источника, расположенного в точкеm = 50. При этом, по-
скольку векторы E,Hи k составляют в электромагнитной
плоской волне правую тройку [15–17], всегда справедливо
Hq
y [m] = −Eq+1
2
z

m − 1
2

η
. (37)
В случае распространения волны влево предыдущее
равенство должно быть заменено соотношением:
Hq
y [m] =
Eq+12
z

m + 1
2

η
. (38)
Пользуясь этим, зададим начальную конфигурацию
электромагнитного поля вида:
ψ0[m] =
0
@ E∗
z [m]
E∗
z

m + 1
2

/η
1
A, (39)
где E∗
z [m] — распределение напряженности электриче-
ского поля, показанное на рис. 4 для импульса P1. Да-
лее, следуя описанной ранее схеме Йи (13), (19)-(21), будем
решать задачу без источников J = 0. Результат решения
в данном случае приведен на рис. 7.
0
5
10
15
20
25
30
35
40
45
0 25 50 75 100 125 150 175 200
q
m
Рисунок 7. Решение задачи Коши для однородного уравнения Максвел-
ла с начальным условием в виде гауссова импульса (39), направленного
влево.
Figure 7. Solution of the Cauchy problem for the homogeneous Maxwell equation
with initial condition in the form of Gaussian pulse (39) directed to the
left.
Как и следовало ожидать, картина, изображенная на
рис. 7, представляет собой гауссов импульс, распростра-
няющийся из начального состояния влево без искаже-
ния формы. Результаты моделирования вполне можно при-
знать корректными с точки зрения определения 2. На-
ми установлено, что тоже самое относится к импульсам
и ограниченным сигналам любой другой формы — коль
скоро начальная конфигурация поля представляет со-
бой (37) или (38) (для любого наперед заданного распре-
деления E∗
z [m]), решение задачи Коши для однородного
уравнения всегда будет корректным.
4.2. Конфигурации поля, соответствующие излуче-
нию ненаправленных источников. В случае, когда гово-
рить о конкретной направленности электромагнитного по-
ля нельзя, начальное состояние электромагнитного поля
должно представлять собой суперпозицию обоих конфи-
гураций (37) и (38) одновременно.
􀀀1
0
1
0 25 50 75 100 125 150 175 200
􀀀1
0
1
􀀀1
0
1
E0
z ; H0
y
m
E0
z H0
y
Hq
y [m]
q = 15
Eq
z [m]
q = 30
Рисунок 8. Некорректно заданная начальная конфигурация электромаг-
нитного поля вида (40) для ненаправленного гауссовского импульса и его
развитие во времени.
Figure 8. Incorrectly specified initial configuration of the electromagnetic
field of the form (40) for a non-directional Gaussian pulse and its evolution
over time.
Здесь очень важно подчеркнуть, что распределение
напряженности электрического поля E∗
z [m], показанное
на рис. 4 для импульса P1 и уже использованное нами ра-
нее в предыдущем подпункте, не может рассматриваться
теперь в качестве «отправной точки» для задания ψ0[m].
Это связано с тем, что для данной конфигурации электри-
ческого поля импульсы, разбегающиеся в противополож-
ные от источника стороны, еще «не разделились оконча-
тельно» (другими словами, прошло лишь время τ/2, но ни-
как не τ ). Это вызывает существенные численные ошибки
и неправильные результаты моделирования в целом, что
очень легко понять, анализируя рис. 8.
Нижняя часть данного рисунка демонстрирует началь-
ную конфигурацию электромагнитного поля, для которой
принято заданным распределение напряженности элек-
трического поля E∗
z [m], соответствующее импульсу P1
рис. 4. Распределение напряженности магнитного поля за-
дается с помощью комбинации (38) для левой полови-
ны пакета и (37) для правой соответственно. Очевидно,
что здесь имеется существенная «нефизическая» особен-
104
Известия Коми научного центра Уральского отделения Российской академии наук № 4 (62), 2023
Серия «Физико-математические науки»
www.izvestia.komisc.ru
ность — конечный разрыв напряженности магнитного поля
в центре волнового пакета.
Разрывное поведение магнитного поля в начальной
конфигурации волнового пакета ответственно за дальней-
шее некорректное поведение решения. А именно, верхняя
часть рис. 8 иллюстрирует, что в последующие моменты
времени как магнитное, так и электрическое поле демон-
стрируют «пилообразное поведение» в центральной зоне
пространства, расположенной между двумя «половинка-
ми начального гауссова импульса, расходящимися влево
и вправо».
Остается нераскрытым лишь следующий вопрос. Дей-
ствительно ли указанная начальная конфигурация поля
некорректна? Если попытаться придать ей некий физиче-
ский смысл, то такое поведение решения можно интерпре-
тировать следующим образом. В начальный момент време-
ни источник тока оказался «резко выключен» именно в той
ситуации, когда формируемый им сигнал достиг максиму-
ма. Очевидно, что при этом неизбежно должно проявить се-
бя явление самоиндукции, действие которого как раз и на-
блюдается в виде «бесконечного пилообразного сигнала»
в постоянно расширяющейся центральной области. Впро-
чем, если вспомнить развиваемые нами ранее соображе-
ния об усреднении некорректно ведущих себя сигналов по
их огибающим, то данным явлением, вообще говоря, мож-
но пренебречь, так как средние значения напряженностей
полей в центральной области всегда остаются нулевыми,
и две половинки начального гауссова импульса действи-
тельно расходятся влево и вправо. Насколько можно поль-
зоваться такой интерпретацией решения — вопрос, по всей
видимости, еще открытый, но не лишенный смысла.
4.3. Иные конфигурации поля. Метод FDTD в том виде,
в котором он сформулирован в данной работе, вообще го-
воря, некорректен для задач статики, поскольку уравнения
Максвелла (2) нами были отброшены в самом начале и ни-
как не принимались в расчет при алгоритмизации схемы
решения (13), (19)-(21). В связи с этим очевидно, что конфи-
гурации электромагнитного поля вида
ψe0
[m] =

Ez[m]
0

(40)
или
ψm
0 [m] =

0
Hy[m]

(41)
в общем случае не могут рассматриваться нами в качестве
корректно заданных начальных условий для однородной
(а равно и неоднородной) задачи Коши. Вместе с тем до-
статочно интересным представляется исследование пове-
дения решений, определяемых алгоритмом Йи (13), (19)-(21),
и в этих случаях.
В качестве примеров достачно посмотреть на рис. 9
и 10, которые иллюстрируют поведение численного реше-
ния однородной задачи Коши для начальных конфигура-
ций вида (40). На рис. 9 исходное распределение Ez[m]
задано равным конфигурации электрического поля им-
пульса P1 (см. рис. 3), а на рис. 10 — соответствующим им-
пульсу P4 (см. рис. 5).
􀀀1
0
1
0 25 50 75 100 125 150 175 200
􀀀1
0
1
􀀀1
0
1
E0
z ; H0
y
m
E0
z H0
y
Hq
y [m]
q = 15
Eq
z [m]
q = 30
Рисунок 9. Некорректно заданная начальная конфигурация электромаг-
нитного поля вида (40) для ненаправленного гауссова импульса и ее раз-
витие со временем.
Figure 9. Incorrectly specified initial configuration of the electromagnetic
field of the form (40) for a non-directional Gaussian pulse and its evolution
over time.
Очевидно, что результаты моделирования ни в том, ни
в другом случае, строго говоря, не могут быть признаны
корректными, хотя для гауссова импульса решение вы-
глядит «несколько более плавным». Понятно, что руковод-
ствоваться такими соображениями при оценке корректно-
сти решения нельзя, тем более, что здесь важнейшим яв-
ляется вопрос правильной интерпретации полученных ре-
зультатов. Данный вопрос на текущий момент нами еще не
закрыт, а его более подробное изучение может быть про-
должено в последующих работах.
􀀀1
0
1
0 25 50 75 100 125 150 175 200
􀀀1
0
1
􀀀1
0
1
E0
z ; H0
y
m
E0
z H0
y
Hq
y [m]
q = 15
Eq
z [m]
q = 30
Рисунок 10. Некорректно заданная начальная конфигурация электромаг-
нитного поля вида (40) для ненаправленного импульса в форме вейвлета
Рикера и ее развитие со временем.
Figure 10. Incorrectly specified initial configuration of the electromagnetic
field of the form (40) for a non-directional impulse in the form of a Ricker
wavelet and its evolution over time.
Заключение
Таким образом, в данной работе исследованы особен-
ности численного решения уравнений Максвелла методом
FDTD в различных формулировках. Показано, что для слу-
чая ограниченных во времени источников поля задачу Ко-
ши для неоднородной системы уравнений можно сформу-
лировать в форме эквивалентной задачи Коши для систе-
мы однородных уравнений, в качестве начального условия
Известия Коми научного центра Уральского отделения Российской академии наук № 4 (62), 2023
Серия «Физико-математические науки»
www.izvestia.komisc.ru
105
которой выступает решение исходной задачи, полученное
к моменту окончания действия источника. Указано, что для
бегущих волн с избранной направленностью такая пере-
формулировка возможна всегда, в то время как для расхо-
дящихся волн при определении начальной конфигурации
поля имеются некоторые особенности, которые обязатель-
но следует учитывать.
Кроме того, нами определен критерий для оценки сте-
пени корректности численного решения, полученного ме-
тодом FDTD согласно алгоритму Йи. С помощью данного
критерия проанализированы особенности численного ре-
шения однородных и неоднородных задач Коши для раз-
личных форм начальных конфигураций электромагнитных
полей и задающих импульсов. Сформулированы необходи-
мые и достаточные условия корректности получаемых ре-
шений.
Авторы благодарны Р.Н. Скандакову (ФМИ ФИЦ Коми
НЦ УрО РАН), В.С. Власову (Сыктывкарский государствен-
ный университет имени Питирима Сорокина) и В.В. Темно-
ву (LSI, Ecole Polytechnique, CEA/DRF/IRAMIS, CNRS, Institut
Polytechnique de Paris) за конструктивную критику, полез-
ные замечания и внимание, стимулировавшие данное ис-
следование.

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, 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 / 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. Schneider, J.B. Understanding the finite-difference time-domain method / J.B. Schneider. - www.eecs.wsu.edu/~schneidj/ufdtd, 2010. - 403 p.

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

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

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

16. Ландау, Л.Д. Теоретическая физика: Т. II. Теория поля / Л.Д. Ландау, Е.М. Лифшиц. - Москва: ФИЗМАТЛИТ, 2006. - 536 с.

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

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

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