ВведениеЧисленные методы решения волновых уравнений иг-рают важную роль не только в конкретных техническихприложениях, но и в фундаментальной науке в целом.К таким методам относится и 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.ru73точность вычислений методом 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:00:51:00 25 50 75 100 125 150 175 200P1P20:00:51:00 25 50 75 100 125 150 175 200Ez[m]mq = 130; "r = 1q = 230; "r = 4J qqw = 10; d = 30Рисунок 1. Моделирование распространения импульсов гауссова видав вакууме P1 и диэлектрике P2 с относительной диэлектрической про-ницаемостью εr = 4.Figure 1. Simulating of the propagation of Gaussian form pulses in vacuumP1 and dielectric P2 with relative permittivity εr = 4.Рисунки 1 и 2 построены в соответствии с алгоритмомЙи, подробно обсуждавшимся нами в предыдущей рабо-те [18] (см. там формулы (13), (19)–(21)), с той разницей, чтотеперь уравнения обновления (19)–(20) в явном виде учи-тывают материальные параметры среды (а именно — еедиэлектрическую εr и магнитную μr проницаемости), в ко-торой распространяются сигналыHq+12ym +12= Hq−12ym +12++Scημr(Eqz [m + 1] − Eqz [m]) ,(2)Eq+1z [m] = Eqz [m] − J q+12[m]++ScηεrHq+12ym +12− Hq+12ym − 12,(3)а простейшие поглощающие граничные условия (21) из [18]заменены нами на более общие уравнения, построенные наоснове дифференциальных уравнений адвекции электро-магнитного поля второго порядка точности. Все детали ре-ализации данной схемы подробно изложены в статье [18]и пособии [14].Замечание 2. Моделирование работы направленного ис-точника тока в формализме TF/SF в данной статье так-же претерпело некоторые изменения по сравнению с [18]и сводится к вычислению полей по схемеHq+12ys − 12= Hq−12ys − 12− ScημrEincz [0, q] , (4)Eq+1z [s] = Eqz [s] +√ScεrμrEincz−12, q +12. (5)Здесь (как и в работе [18]) s = 50 — это фиксированныйво всех численных экспериментах номер узла сетки, за-дающий расположение в пространстве точечной антен-ны, формирующей падающее поле Eincz . Вычисление элек-трического Ez[s] и магнитного Hys − 12полей, соглас-но (4), (5) по заданному типу падающей волны Eincz , поз-воляет имитировать работу источника тока J q, формиру-ющего волну, излучаемую в правую область сетки m ⩾ s.В дальнейшем везде под работой источника тока J q будемиметь в виду именно эту схему.􀀀0:50:00:51:00 25 50 75 100 125 150 175 200P3P4􀀀0:50:00:51:00 25 50 75 100 125 150 175 200Ez[m]mq = 150; "r = 1q = 250; "r = 4J qqlm = 50; d = 1Рисунок 2. Моделирование распространения импульсов в форме вейвле-та Рикера в вакууме P3 и диэлектрике P4 с относительной диэлектри-ческой проницаемостью εr = 4.Figure 2. Simulating of pulse propagation in the form of a Ricker wavelet invacuum 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.ru75b S 12t ε b S 12t− b S−12tΔt!Eqz [m] == b S 12t b S 12x − b S−12xΔx!Hqy [m].(8)Здесь учтено определение числа Куранта (1), связьскорости света в вакууме с его диэлектрической ε0и магнитной μ0 проницаемостями c = 1/√ε0μ0 и вы-ражение для характеристического импеданса вакуумаη =pμ0/ε0 ≈ 120π (для справки см., например, учеб-ники [20, 21]). Кроме того, при записи (8) используется обо-значение для абсолютной диэлектрической проницаемо-сти среды ε = εrε0.Определим для удобства также конечно-разностныеоператоры ˜∂i согласно˜∂i =b S 12i− b S−12iΔi, (9)где i — это x или t. С помощью (9) закон Ампера в стацио-нарной среде без сторонних токов (8) окончательно можетбыть записан в форме Йи [14]ε b S 12t˜∂tEqz [m] = b S 12t˜∂xHqy [m]. (10)Таким же образом, с помощью (6), (7) и (9), дискретныйаналог закона Фарадея (2) приводится к форме Йиμ b S 12x ˜∂tHqy [m] = b S 12x ˜∂xEqz [m]. (11)Теперь рассмотрим плоскую монохроматическую вол-ну частоты ω, распространяющуюся в сетке Йи вправо(m ⩾ s)Eqz [m] = E0 ei(ωqΔt−˜βmΔx), (12)Hqy [m] = H0 ei(ωqΔt−˜βmΔx), (13)где ˜ β — волновое число, т. е. постоянная распространенияплоской монохроматической волны в FDTD-сетке (отлича-ющаяся от соответствующей константы β в непрерывномпространстве), а E0 и H0 — комплексные амплитуды на-пряженности электрического и магнитного полей.Лемма 1. Действие операторов (9) на произвольную компо-ненту ψ плоской монохроматической волны (12), (13) сво-дится к˜∂tψ = i2ΔtsinωΔt2ψ, (14)˜∂xψ = −i2Δxsin ˜ βΔx2!ψ. (15)Доказательство. Проверяется непосредственной подста-новкой (12), (13) в (9) с учетом (6), (7). Выпишем здесь в явномвиде только действие на ψ операторов сдвига (6), (7) припараметрах χ, τ = ±12b S±12t ψ = e±iωΔt/2ψ, b S±12x ψ = e∓i ˜βΔx/2ψ, (16)с помощью которых равенства (14) и (15) получаются эле-ментарно.Утверждение 1. Дисперсионное соотношение для сетки Йиможет быть представлено в формеsinωΔt2=√ΔtεμΔxsin ˜ βΔx2!. (17)Доказательство. Используя результаты (14) и (15) Леммы 1,а также (16) в законе Ампера (10) для плоской монохрома-тической волны, запишемiε2ΔtsinωΔt2eiωΔt/2Eqz [m] == −i2Δxsin ˜ βΔx2!eiωΔt/2Hqy [m].(18)Подставляя в последнее равенство явные выраженияполей (12), (13) для плоской монохроматической волны и со-кращая общие множители, получаемε1ΔtsinωΔt2E0 = − 1Δxsin ˜ βΔx2!H0. (19)Отсюда приходим к выражению для численного импе-данса в сетке ЙиE0H0= − ΔtεΔx·sin˜ βΔx2sinωΔt2. (20)Выполняя аналогичные действия применительно к за-кону Фарадея (11), последовательно получаемiμ2ΔtsinωΔt2e−i ˜βΔx/2Hqy [m] == −i2Δxsin ˜ βΔx2!e−i ˜βΔx/2Eqz [m],(21)μ1ΔtsinωΔt2H0 = − 1Δxsin ˜ βΔx2!E0 (22)и соответствующий импеданс в формеE0H0= −μΔxΔt·sinωΔt2sin˜ βΔx2. (23)Приравнивая правые части (20) и (23) и выполняя в по-лучившемся равенстве перекрестное произведение со-множителей, получаемsin2ωΔt2=Δ2tεμΔ2xsin2 ˜ βΔx2!. (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. Отклонение фазовой скорости волныв сетке Йи от соответствующего значения в непрерывномслучае может быть описано равенством˜cpcp=π√εrμrNλ arcsin√εrμrScsinπScNλ, (27)где параметр Nλ — есть число узлов пространственнойсетки на длину волны в свободном пространствеλ = NλΔx. (28)Доказательство. Используем связь фазовой скоростис волновым числом в непрерывном пространстве [20, 21]и FDTD-сеткеcp =ωβ, ˜cp =ω˜ β, (29)что позволяет привести левую часть равенства (27) к виду˜cpcp=β˜ β=βΔx2˜βΔx2. (30)Далее применим (25) к числителю последнего равенстваβ = ω√εμ = 2πcλ√ε0εrμ0μr =2πλ√εrμr, (31)где длина волны в вакууме λ дискретизуется согласно (28),что даетβΔx2=π√εrμrNλ. (32)Применяя теперь к знаменателю правой части (30) дис-персионное соотношение (17), в котором используются пре-образования аналогичные (31), а также определение числаКуранта (1) совместно с (28), получаем˜ βΔx2= arcsin√εrμrScsinπScNλ. (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:700:750:800:850:900:951:001:051:102 4 6 8 10 12 14 16 18 201234~cp/cpNРисунок 3. Зависимость фазовой скорости, определяемой согласно мето-ду FDTD по отношению к ее точному значению (27), от параметра дискре-тизации длины волныNλ для некоторых сред с относительными показа-телями преломления nr. Число Куранта Sc = 1.Figure 3. Dependence of the phase velocity ratio, determined according tothe FDTD method with respect to its exact value (27), on the wavelength discretizationparameterNλ for some media with relative refractive indicesnr.The Courant number Sc = 1.Известия Коми научного центра Уральского отделения Российской академии наук № 5 (71), 2024Серия «Физико-математические науки»www.izvestia.komisc.ru77Особенности численной дисперсии метода FDTD, опре-деляемые выражением (27), и дополняющие приведенныйвыше пример, представлены на рис. 3, где показано семей-ство кривых, для которых nr > Sc (а именно, nr =√2 длякривой 1, nr =√3 — для линии 2 и nr = 2 — для 3 соот-ветственно). Видно, что увеличение показателя преломле-ния среды при плохой дискретизации длины волны приво-дит к существенному запаздыванию волны, моделируемойметодом FDTD, что выражается в значительном отклоненииотношения ˜cp/cp от единицы.Кроме того, на рис. 3 изображена кривая 4, отвеча-ющая случаю nr < Sc (в данном конкретном случае вы-брано значение nr =p1/2). Этот пример демонстрирует,что в средах оптически менее плотных, чем вакуум, FDTD-расчет приводит к распространению моделируемых волнс опережением по сравнению с истинной скоростью.Следствие 1. Из рис. 3 видно, что точность FDTD-расчетападает с уменьшением Nλ, а кроме того — по мере откло-нения проницаемости среды от единицы (что имеет местои в диэлектриках, и в магнетиках, и в магнитных диэлек-триках).Следствие 2. Легко вычислить предел отношения (27), ко-торый независимо от величин nr и Sc равенlimNλ→+∞˜cpcp= 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) =nrScsinπScNλ, (37)характер которой определяет полосы пропускания и не-пропускания сетки Йи. Действительно, легко понять, чтообласть значений параметров (Nλ, Sc, nr), на которой вы-полняется условие |φ(Nλ, Sc, nr)| > 1, отвечает непро-пусканию волн, так как выражение arcsin φ при этом неимеет смысла в вещественных значениях.Подобная картина явления имеет место в модели Кро-нига-Пенни, а также в уравнениях Хилла и Матье, описы-вающих параметрические колебания [19, 27–29], а такжераспространение волн в системах с пространственной пе-риодичностью в расположении неоднородностей. В каче-стве периодической структуры в нашем случае (в методеFDTD) выступает непосредственно сама расчетная сеткаЙи. При этом пропускание и непропускание волн (кванто-во-механический аналог — разрешенные и запрещенныезоны) определяется как раз величинами Nλ, Sc и nr.02040608010050 100 150 200 250 300 350ScN􀀀101A1A2BРисунок 4. (Nλ, Sc)-диаграмма полос пропускания сетки Йипри nr = 100.Figure 4. (Nλ, Sc)-diagram of Yee grid bandwidths when nr = 100.Рисунок 4 иллюстрирует (Nλ, Sc)-диаграмму по-лос пропускания сетки Йи, построенную на основефункции (37) для сравнительно большого значения по-казателя преломления среды nr = 100. Области по-лос непропускания {Ai}Ni=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:00:51:00 25 50 75 100 125 150 175 200P5 P1 P60:00:51:00 25 50 75 100 125 150 175 200q = 130; nr = 1:0q = 80; nr = 2:0Ez[m] q = 155; nr 0:7mw = 10; d = 30J qqРисунок 5. Моделирование распространения импульсов гауссова видав средах, оптически более P5 (nr = 2) и менее P6 (nr =p1/2) плот-ных, чем вакуум, с применением коррекции программной дисперсии (38).Figure 5. Simulating of Gaussian pulse propagation in media opticallymore P5 (nr = 2) and less dense P6 (nr =p1/2) than vacuum, usingcorrection 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).􀀀1010 25 50 75 100 125 150 175 200􀀀101􀀀101q = 50Eqz [m]mq = 100Eqz [m]q = 130Eqz [m]Рисунок 6. Временная динамика разрушения численного решения, по-лучаемого в процессе вычислений согласно алгоритму (2) – (5) приSc − nr = 10−3.Figure 6. Time dynamics of the numerical solution destruction obtainedduring computation according to the algorithm (2) – (5) whenSc − nr = 10−3.Как видно из рис. 6, к моменту времени q = 130 уро-вень шума, произвольно возникающего в результате вы-числений по схеме (2) – (5) при Sc − nr = 10−3, двухкрат-но превышает величину полезного сигнала по абсолютномуИзвестия Коми научного центра Уральского отделения Российской академии наук № 5 (71), 2024Серия «Физико-математические науки»www.izvestia.komisc.ru79значению. При этом пространственная протяженность ча-сти сетки Йи, затронутой данным шумом, составляет 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%.􀀀1010 25 50 75 100 125Pb Pf􀀀101􀀀1011 2E20z [m]m1 2E40z [m]1 2E60z [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 insimulating the propagation of a Gaussian-shaped signal in media opticallyless dense than vacuum under optimal choice of the Courant number. Thedashed-dotted curve 1 corresponds to the case of Sc = nr = 10−1, thesolid 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 =c4π[E ×H] . (40)Справедливость такого обобщения демонстрируетрис. 8, на котором представлены результаты сравненияволновых пакетов, формируемых тем же источником сиг-нала, что и на рис. 1, 5 – 7, зафиксированные в одинаковыймомент времени q = 100, при распространении в правойсреде с εr = μr = +1 (нижняя половина рисунка) и левойсреде с εr = μr = −1 (верхняя половина).􀀀1010 25 50 75 100 125 150 175 200􀀀101EzHyE80z [m]; H80y [m]mEzHyE80z [m]; H80y [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 computationalalgorithm (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. Возможно ли корректировать численную диспер-сию при моделировании распространения сигналовв неоднородных средах? Можно ли этого добиться,изменяя число Куранта динамически в процессе мо-делирования?Этот список интересных вопросов ни в коей мере непретендует на полноту и вполне может быть расширен.Авторы заявляют об отсутствии конфликта интересов.