1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 62, № 1, 2022

Журнал вычислительной математики и математической физики, 2022, T. 62, № 1, стр. 90-104

Прецизионные методы решения одномерных уравнений Максвелла в слоистых средах

А. А. Белов 1, 2, *, Ж. О. Домбровская 1, **

1 МГУ
119991 Москва, Ленинские горы, 1, стр. 2, Россия

2 РУДН
117198 Москва, ул. Миклухо-Маклая, 6, Россия

* E-mail: aa.belov@physics.msu.ru
** E-mail: dombrovskaya@physics.msu.ru

Поступила в редакцию 15.04.2021
После доработки 15.04.2021
Принята к публикации 17.09.2021

Полный текст (PDF)

Аннотация

Для системы стационарных и нестационарных одномерных уравнений Максвелла ранее была построена бикомпактная разностная схема. Ее шаблон включает только один шаг пространственной сетки, причем границы слоев выбираются в качестве узлов сетки. Схема явно учитывает условия сопряжения на границе раздела сред. Это позволяет проводить расчеты обобщенных решений с разрывом решения и его производной. Для решения нестационарных задач применяется оригинальный метод спектрального разложения, который позволяет учитывать произвольный закон дисперсии среды. В данной работе построена новая форма записи данной схемы. Это позволило снизить трудоемкость расчетов в 4 раза. Такой выигрыш является существенным. Впервые дано строгое обоснование предложенной схемы. Библ. 93.

Ключевые слова: уравнения Максвелла, бикомпактные схемы, слоистые среды, условия сопряжения, дисперсия вещества.

1. ВВЕДЕНИЕ

Задачи для системы уравнений Максвелла встречаются во многих физических и технических приложениях. Например, это задачи магнитной гидродинамики в токамаках и астрофизических объектах, распространение излучения по волноводам, расчеты элементов радиотехнических цепей и микросхем. В последнее время особенное внимание уделяется задачам оптики диэлектрических наноструктур. На их основе разрабатывают элементную базу для оптической вычислительной техники. Многие реальные задачи рассматриваются в нестационарной постановке.

Для таких задач предложено большое количество методов расчета. Им посвящены тысячи статей. Однако точность и надежность этих подходов зачастую оказываются недостаточными. Не претендуя на библиографическую полноту, в данной работе мы привели анализ трудностей, с которыми сталкиваются наиболее известные методы.

Эти трудности удалось преодолеть в недавно предложенных подходах [1], основанных на построении бикомпактных схем и на методе спектрального разложения. Эти методы позволяют проводить расчеты обобщенных решений даже с сильным разрывом на границе раздела сред. Они реализованы пока только для системы одномерных уравнений Максвелла. Расчеты представительных тестовых задач, приведенные в [1], показали, что бикомпактные схемы кардинально превосходят другие методы по точности и надежности.

Данная работа посвящена дальнейшему развитию бикомпактных схем для системы одномерных уравнений Максвелла. Построена новая форма записи этих схем, которая позволила снизить трудоемкость расчетов в 4 раза. Такой выигрыш по трудоемкости является существенным. Впервые дано строгое обоснование предложенных схем. Доказаны теоремы о сходимости, причем для произвольных неоднородных сред и неравномерных сеток.

2. ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ МАКСВЕЛЛА

2.1. Методы

Приведем обзор наиболее популярных алгоритмов решения уравнений Максвелла. Многие из них (матричные методы [2]–[4], модальные методы [5], псевдоспектральные методы [6], методы мультипольной аппроксимации [7]–[10], методы конечных разностей и конечных элементов в частотной области [11]–[15]) дают хорошие результаты для достаточно узких классов задач.

Кудрявцев и Трашкеев предложили нетривиальный подход, который заключается в сведении электродинамической задачи к системе уравнений газодинамического типа с помощью формализма потенциалов [16]. Это позволяет использовать хорошо разработанные почти монотонные разностные схемы высокого порядка точности. Среди них, например, методы, предложенные в работах Тишкина и Фаворского [17], схемы классов Essential non-oscillatory (ENO) [18], Weighted essential non-oscillatory (WENO) [19] и др.

Наиболее универсальны методы конечных разностей (FDTD) и конечных элементов (FETD) во временной области [20]–[23]. Практически все они сводятся к схеме типа “крест” [24], в которой электрическое и магнитное поля относятся к целым и полуцелым пространственным узлам и временным слоям. Такие схемы называются схемами “крест”. Они впервые появились в газодинамике и были исследованы Курантом. Аналогичный подход применяется в так называемом методе конечных объемов во временной области (FVTD), в котором схема строится с помощью разностной аппроксимации интегральных уравнений Максвелла [22].

Явные схемы метода FDTD являются лишь условно устойчивыми, поэтому для них необходимо соблюдать выполнение условия Куранта. Для неявных схем этого условия нет, однако они оказываются гораздо более трудоемкими. В литературе описан ряд таких схем [22]. В одномерном случае работоспособны чисто неявная схема и схема Кранка-Николсон. В многомерном случае на каждом временном слое нужно решать линейную систему с сильно разреженной матрицей огромной размерности. Для экономичного решения таких задач проводят их факторизацию, т.е. расщепляют многомерную задачу на произведение одномерных. Для прямоугольной геометрии наиболее эффективной оказывается эволюционная факторизация, предложенная Калиткиным [25].

Методы типа FDTD, FETD, FVTD чрезвычайно популярны. К настоящему времени опубликовано более 2000 статей, посвященных развитию этих методов и их применению в различных задачах. Однако несмотря на такое обилие работ, принципиальные недостатки этих методов до сих пор не устранены. Эти недостатки перечислены далее.

2.2. Методы конечных разностей, конечных объемов и конечных элементов во временной области

Эти методы имеют следующие принципиальные недостатки.

Резкое ухудшение точности на границах раздела, понижающее общую точность расчета. Схемы методов FDTD, FETD, FVTD не являются компактными, поскольку ячейка Йе занимает 1.5 шага по пространству (ячейки магнитного поля сдвинуты на полшага относительно узлов электрического поля). Поэтому как бы ни были выбраны узлы пространственной сетки, граница раздела сред (и, значит, разрыв поля) попадает внутрь ячейки. Поэтому условие аппроксимации схемы Йе фактически не выполняется. Заметим, что схема обладает некоторым “запасом прочности” и в случае слабых разрывов обеспечивает сходимость (но медленную).

Поэтому на практике погрешность расчета слоистых диэлектрических и магнитных сред оказывается очень большой [26], [27]. При приемлемой густоте сеток добиться удовлетворительно точности зачастую не получается. Такие расчеты занимают много часов на мощных компьютерах, однако нередко результаты согласуются с экспериментом лишь с визуальной графической точностью.

В [28]–[31] предлагалось усреднение показателя преломления в ячейке, внутрь которой попадает граница раздела. В [32], [33] разрыв показателя преломления сглаживали, т.е. заменяли границу раздела средой с непрерывно изменяющимся показателем преломления. В работе [34] проводили искусственную сшивку решений до и после границы раздела с помощью фиктивных ячеек. Однако перечисленные подходы являются эвристическими. Они вносят физическую погрешность в решение. Оценить ее количественно и исключить практически невозможно. Известно, что размывание или усреднение границ раздела может сильно ослаблять расчетное отражение. Нетрудно заметить, что усреднение показателя преломления вносит погрешность $O(1)$ вблизи границы раздела. Однако такая оценка недостаточна для экстраполяционного уточнения.

Для преодоления этой трудности рассматривались схемы FDTD с явным выделением особенности [35]. Эти схемы применялись к задачам с филаментными токами. Такой подход работоспособен, если имеется только одна граница раздела. Однако наличие уже двух границ раздела приводит к многократным переотражениям. В этом случае схемы с явным выделением особенностей становятся крайне громоздкими. В то же время однородные разностные схемы намного удобнее [36].

Таким образом, методы типа FDTD не позволяют явно учесть условия сопряжения на границе раздела сред, вытекающие из уравнений Максвелла. Это принципиально ограничивает их точность.

В методах FETD реализация правильных условий на границе раздела представляет проблему. Конечные элементы Неделека [37] обеспечивают тангенциальную непрерывность сеточной функции. Элементы Равьяра-Томаса (см., например, [38] и цитированную литературу) дают нормальную непрерывность сеточной функции. Однако эти условия являются частными случаями физических граничных условий (а именно, граничное условие для вектора напряженности магнитного поля ${\mathbf{H}}$ включает разрыв, пропорциональный величине поверхностных токов). В многочисленных реализациях разрывного метода Галеркина [39]–[45] условия на границе раздела также нарушаются. Вместо них вводятся фиктивные потоки между ячейками. Эти потоки рассматриваются как степени свободы. В результате в расчете возникают нефизичные решения [46].

Численная дисперсия. Схемы методов FDTD и FETD обладают численной дисперсией, причем последняя анизотропна [23]. Иначе говоря, скорости распространения плоских волн по оси координатной сетки или по диагонали неодинаковы. Искривленный волновой фронт (например, цилиндрический или сферический) при распространении “самопроизвольно” деформируется. Это принципиальный дефект, присущий всем методам данного класса.

Отсюда следует, что схема Йе не может быть строго консервативной. В каждой ячейке дисбаланс отличен от нуля, но его величина стремится к нулю при уменьшении шага сетки. Такие схемы называют почти консервативными. Как известно, наличие дисбаланса приводит к ухудшению количественной точности.

Немонотонность. Схема Йе является немонотонной [27]. Это значит, что для неоднородных сред в численном решении возникают нефизичные пилообразные осцилляции, которые могут быть сильны. В результате точность катастрофически ухудшается, особенно вблизи границ раздела сред.

Фиктивные граничные условия. Существующие варианты метода FDTD, FETD, FVTD не позволяют задать непосредственно волну, падающую на рассеиватель из бесконечности. Падение волны задают с помощью фиктивного источника поля (метод разделения полного и рассеянного полей TF/SF) [23]. Этот прием представляется искусcтвенным. Уход волны на бесконечность описывают с помощью условий Мура [47] либо с помощью поглощающих стенок (т.н. идеально согласованного слоя PML) [48]. Условия Мура более физичны, однако они имеют лишь первый порядок точности. При использовании PML неизбежно возникают переотражения, которые могут ощутимо ухудшить точность расчета.

Погрешность учета дисперсии материалов. Плазма является сильно диспергирующей средой. Многие твердые тела (например, Si, Ta и их оксиды, Ge и др.) также обладают сильной дисперсией. Чтобы это учесть, в уравнение для вектора электрической индукции добавляют интегральный член, содержащий свертку напряженности электричесокого поля и диэлектрической восприимчивости. Система уравнений Максвелла становится интегродифференциальной. Такие задачи намного более трудны, чем чисто дифференциальные.

Для некоторых законов дисперсии (полиномиальный, отношение двух полиномов, модели Дебая и Лоренца) эту свертку удается вычислить точно [49]–[59]. Однако реальные законы дисперсии приходится приближать комбинацией таких модельных законов. Другой подход заключается в прямом сеточном вычислении этой свертки, для которой фактически записывается некоторая явная схема [60]–[72]. Оба способа учета дисперсии в методах FDTD, FETD, FVTD вносят погрешность в результаты расчетов, причем эта погрешность нередко оказывается весьма существенной [73], [74].

2.3. Контроль точности

Единственным надежным способом контроля фактической точности для разностных методов являются расчет на сгущающихся сетках и применение оценок точности по методу Ричардсона.

В работах [75], [76] использовались оценки точности и экстраполяция по Ричардсону, но сгущение проводилось только для сетки по времени, причем не на каждом ее шаге. Результат такой процедуры не дает представления о действительной точности. Чтобы применить метод Ричардсона, необходимо проводить глобальное (а не локальное!) сгущение сетки, причем одновременно по всем переменным.

В [77] исследована аппроксимация и выписан остаточный член схемы Йе. Однако авторы не приводят ни априорных, ни апостериорных оценок точности.

В работах [28], [31] анонсируется применение экстраполяции по методу Ричардсона. Однако результат экстраполяции, полученный в этих работах, не соответствует теоретически ожидаемому [78]. Кроме того, из-за усреднения показателя преломления локальная ошибка вблизи границы раздела есть $O(1)$. В этом случае оценки точности по Ричардсону неправомерны, а экстраполяция может ухудшать точность.

3. СЛОИСТЫЕ СРЕДЫ

Чтобы проводить расчеты обобщенных решений в слоистых средах, необходимо выполнение двух специфических требований: схемы должны быть бикомпактными и консервативными.

3.1. Бикомпактность

Границы раздела сред не должны попадать внутрь шаблона разностной схемы. В противном случае внутри шаблона нарушается гладкость, и невязка оказывается огромной. Для этого, во-первых, нужно выбирать такие сетки, в которых во всех границах раздела расположены узлы, и ни одна граница не попадает внутрь ячейки. Такие сетки называют специальными [79]. Во-вторых, необходимо использовать двухточечные схемы, шаблон которых включает только один шаг по пространству. Такие схемы называют бикомпактными.

Если выбраны специальные сетки, то внутри каждого шага решение является гладким. Поэтому бикомпактные схемы на специальных сетках обеспечивают теоретический порядок точности в любой норме даже при наличии сильных разрывов, и их погрешность быстро убывает при уменьшении шага сетки. Таким образом, использование бикомпактных схем для задач в слоистых средах является принципиально необходимым. Впервые понятие бикомпактности было сформулировано в [80]–[82] и применено к уравнению теплопроводности в слоистой среде.

3.2. Консервативность

Уравнения математической физики выводятся из интегральных соотношений (уравнений баланса), выражающих законы сохранения для малого объема. Разностная схема должна отражать основные свойства сплошной среды. Важнейшим свойством является выполнение разностных аналогов основных законов сохранения. Такие схемы называются консервативными [36], [83]–[85]. При этом уравнение баланса для всей расчетной области должно быть алгебраическим следствием балансов в отдельных ячейках. Для построения консервативных схем в [36], [83], [84] предложен интегроинтерполяционный метод.

В сложных задачах может быть несколько законов сохранения. Примером такой задачи является газодинамика: одновременно должны выполняться законы сохранения массы, импульса и энергии. Схемы, выражающие все нужные законы сохранения, называются полностью консервативными. Для задач электродинамики исходными физическими законами сохранения являются интегральные уравнения Максвелла.

Дифференциальные формы основных уравнений применимы только если решение является достаточно гладким. Если же решение имеет разрыв или излом, то в точке разрыва дифференциальное уравнение неприменимо. В этом случае записывают интегральные уравнения баланса для малого объема, захватывающего области по обе стороны от разрыва. Далее этот объем стягивают в точку, совпадающую с разрывом. Это дает условия на разрыве. Например, в газодинамике таковыми являются условия Гюгонио, определяющие скачок температуры и давления в ударной волне.

В электродинамике применение интегральных уравнений Максвелла к указанному объему дает известные условия сопряжения полей на границе раздела сред. Эти условия являются прямым следствием исходных уравнений баланса. Тем самым, чтобы схема для расчета электродинамики слоистых сред была полностью консервативной, она обязательно должна учитывать не только уравнения Максвелла внутри слоев, но и условия сопряжения на границах раздела. В противном случае возникает сеточный дисбаланс, который влияет на решение во всей области.

Например, в схеме Йе метода FDTD это приводит к появлению нефизичных осцилляций и резкому падению точности (порядок точности ухудшается до первого). Для этой схемы доказано, что для однородной среды из разностных роторных уравнений следует разностный аналог дивергентных [22]. Однако эта схема не учитывает физически корректных условий сопряжения на границах раздела. Поэтому она неконсервативна согласно указанному выше определению по Тихонову-Самарскому. То же относится и к другим конечно-разностным и конечно-элементным методам, не учитывающим условия сопряжения.

4. СТАЦИОНАРНАЯ ЗАДАЧА

4.1. Постановка задачи

1. Рассмотрим распространение плоской электромагнитной волны частоты $\omega $ перпендикулярно набору из $Q$ плоско-параллельных пластин общей толщиной $a$. Обозначим координаты границ слоев через ${{\xi }_{1}} \leqslant {{\xi }_{2}} \leqslant \ldots \leqslant {{\xi }_{{Q - 1}}}$. Пусть все объемные и поверхностные токи зависят от времени как $ \sim {\kern 1pt} {{e}^{{ - i\omega t}}}$.

Ориентируем координатную ось $z$ перпендикулярно пластинам. Тогда в одномерной задаче электрическое поле ${\mathbf{E}}$ направлено по оси $x$, магнитное поле ${\mathbf{H}}$ – по оси $y$, волновой вектор ${\mathbf{k}}$ – по оси $z$. Здесь $k = \omega {\text{/}}c$, где $c$ – скорость света. Диэлектрическая проницаемость ${{\varepsilon }_{q}}$, магнитная восприимчивость ${{\mu }_{q}}$, проводимость ${{\sigma }_{q}}$ объемная плотность токов ${{{\mathbf{J}}}_{q}}$ $q$-й пластины, а также поверхностная плотность токов ${{{\mathbf{j}}}_{q}}$ на $q$-й границе не зависят от $x$ и $y$. Векторы ${{{\mathbf{J}}}_{q}}$ и ${{{\mathbf{j}}}_{q}}$ направлены по оси $x$. Амплитуды полей зависят только от $z$. Одномерная постановка актуальна во многих приложениях (например, в интегральной фотонике).

2. Чтобы получить консервативную схему, будем исходить из интегральных уравнений Максвелла

(1)
$\int\limits_\Gamma {{{{\mathbf{H}}}_{q}}d{\mathbf{l}}} = \frac{{4\pi }}{c}\int\limits_S {({{{\mathbf{J}}}_{q}} + {{\sigma }_{q}}{{{\mathbf{E}}}_{q}})d{\mathbf{s}}} - \frac{{i\omega }}{c}\int\limits_S {{{{\mathbf{D}}}_{q}}d{\mathbf{s}}} ,\quad {{{\mathbf{D}}}_{q}} = {{\varepsilon }_{q}}{{{\mathbf{E}}}_{q}},\quad 1 \leqslant q \leqslant Q;$
(2)
$\int\limits_\Gamma {{{{\mathbf{E}}}_{q}}d{\mathbf{l}}} = \frac{{i\omega }}{c}\int\limits_S {{{{\mathbf{B}}}_{q}}d{\mathbf{s}}} ,\quad {{{\mathbf{B}}}_{q}} = {{\mu }_{q}}{{{\mathbf{H}}}_{q}},\quad 1 \leqslant q \leqslant Q.$
Здесь $S$ – произвольная поверхность, ограниченная контуром $\Gamma $. Используется система единиц СГС.

3. На внешних границах поставим условия излучения

(3)
$\frac{{\partial {{{\mathbf{E}}}_{1}}}}{{\partial z}} + ik{{{\mathbf{E}}}_{1}} = 2ik{{{\mathbf{E}}}^{0}},\quad z = 0;\quad \frac{{\partial {{{\mathbf{E}}}_{Q}}}}{{\partial z}} - ik{{{\mathbf{E}}}_{Q}} = 2ik{{{\mathbf{E}}}^{a}}{{e}^{{ - ika}}},\quad z = a.$
Здесь ${{{\mathbf{E}}}^{0}}$, ${{{\mathbf{E}}}^{a}}$ – заданные амплитуды волн, падающих справа и слева.

Поясним условия (3). Волны, распространяющиеся в положительном и отрицательном направлении оси $z$, имеют вид ${{{\mathbf{E}}}_{ + }} \sim {{e}^{{i\omega t - ikz}}}$ и ${{{\mathbf{E}}}_{ - }} \sim {{e}^{{i\omega t + ikz}}}$ соответственно. В общем случае по обе стороны от рассеивателя поле представимо в виде суперпозиции ${{{\mathbf{E}}}_{ + }}$ и ${{{\mathbf{E}}}_{ - }}$. При $z < {{\xi }_{1}}$ падающей является волна ${{{\mathbf{E}}}_{ + }}$, а амплитуда волны ${{{\mathbf{E}}}_{ - }}$ неизвестна. Нетрудно проверить, что для ${{{\mathbf{E}}}_{ - }}$ выполнено ${{\partial }_{z}}{{{\mathbf{E}}}_{ - }} + i\omega {{c}^{{ - 1}}}{{{\mathbf{E}}}_{ - }} = 0$ независимо от ее амплитуды. Для ${{{\mathbf{E}}}_{ + }}$ справедливо ${{\partial }_{z}}{{{\mathbf{E}}}_{ - }} + i\omega {{c}^{{ - 1}}}{{{\mathbf{E}}}_{ - }} = 2ik{{{\mathbf{E}}}^{0}}$, это позволяет явно задать амплитуду ${{{\mathbf{E}}}_{ + }}$, не затрагивая волну ${{{\mathbf{E}}}_{ - }}$. При $z > {{\xi }_{{Q - 1}}}$ волны “меняются ролями”: задана падающая волна ${{{\mathbf{E}}}_{ - }}$, а неизвестной – волна ${{{\mathbf{E}}}_{ + }}$. Поэтому при $z = a$ граничное условие записывается аналогично таковому при $z = 0$, но со сменой знака при $ik$ в левой части.

Условия (3) являются реализацией парциальных условий излучения Свешникова [86] для непрерывного спектра (а не дискретного спектра волноводных мод). Они позволяют естественно учитывать падение и отражение лазерного луча. Во многих коммерческих пакетах программ (Lumeical, Optiwave и т.д.) падающую волну задают моделью фиктивного источника [22].

4. На границах слоев поставим условия сопряжения

(4)
$\begin{gathered} {{{\mathbf{e}}}_{z}} \times ({{{\mathbf{E}}}_{q}} - {{{\mathbf{E}}}_{{q - 1}}}) = 0,\quad {{{\mathbf{e}}}_{z}} \cdot ({{{\mathbf{D}}}_{q}} - {{{\mathbf{D}}}_{{q - 1}}}) = 0, \\ {{{\mathbf{e}}}_{z}} \times ({{{\mathbf{H}}}_{q}} - {{{\mathbf{H}}}_{{q - 1}}}) = \frac{{4\pi }}{c}{{{\mathbf{j}}}_{q}},\quad {{{\mathbf{e}}}_{z}} \cdot ({{{\mathbf{B}}}_{q}} - {{{\mathbf{B}}}_{{q - 1}}}) = 0. \\ \end{gathered} $
Равенства (1)–(4) есть постановка стационарной задачи.

4.2. Разностная схема

1. Введем сетку $0 = {{z}_{0}} < {{z}_{1}} < \ldots < {{z}_{N}} = a$, ${{z}_{{n + 1}}} - {{z}_{n}} = \Delta {{z}_{n}}$ так, чтобы границы слоев были узлами. Введем сеточные значения полей. В общем случае поля могут быть разрывными, поэтому в каждом внутреннем узле нужно определить левое и правое предельные значения. Для каждого шага $\Delta {{z}_{{n - 1/2}}}$, $1 \leqslant n \leqslant N$, введем значения полей ${{E}_{{2n - 2}}}$, ${{H}_{{2n - 2}}}$, относящиеся к левой границе ${{z}_{{n - 1}}}$, и значения полей ${{E}_{{2n - 1}}}$, ${{H}_{{2n - 1}}}$, относящиеся к правой границе ${{z}_{n}}$. Таким образом, в каждом внутреннем узле правое и левое предельные значения каждого поля могут быть неодинаковы.

2. В качестве поверхности $S$ для (1) возьмем прямоугольник ${{z}_{{n - 1}}} \leqslant z \leqslant {{z}_{n}}$, $0 \leqslant y \leqslant 1$ в плоскости $x = 0$. Контурный интеграл вычисляется точно

(5)
$\int {{{{\mathbf{H}}}_{q}}d{\mathbf{l}}} = \int\limits_0^1 {{{{\left. {{{{({{H}_{y}})}}_{q}}} \right|}}_{{{{z}_{{n - 1}}}}}}dy} - \int\limits_0^1 {{{{\left. {{{{({{H}_{y}})}}_{q}}} \right|}}_{{{{z}_{n}}}}}dy} = {{H}_{{2n - 2}}} - {{H}_{{2n - 1}}}.$
Поверхностный интеграл в правой части (1) вычислим по гибридной квадратурной формуле, комбинируя правила трапеций и средних
(6)
$\int\limits_S {A{{{\mathbf{E}}}_{q}}d{\mathbf{s}}} = \int\limits_{{{z}_{{n - 1}}}}^{{{z}_{n}}} {dz} \int\limits_0^1 {dyA{{{({{E}_{x}})}}_{q}}} \approx \frac{1}{2}\Delta {{z}_{{n - 1/2}}}{{A}_{{n - 1/2}}}({{E}_{{2n - 1}}} + {{E}_{{2n - 2}}}),\quad A = - \frac{{i\omega }}{c} + \frac{{4\pi }}{c}\sigma .$
Интеграл от плотности тока $J$ вычислим по формуле средних. При этом значения диэлектрической проницаемости ${{\varepsilon }_{{n - 1/2}}} \equiv {{\varepsilon }_{q}}\left( {0.5({{z}_{{n - 1}}} + {{z}_{n}})} \right)$, магнитной восприимчивости ${{\mu }_{{n - 1/2}}} \equiv {{\mu }_{q}}\left( {0.5({{z}_{{n - 1}}} + {{z}_{n}})} \right)$, проводимости ${{\sigma }_{{n - 1/2}}} \equiv {{\sigma }_{q}}\left( {0.5({{z}_{{n - 1}}} + {{z}_{n}})} \right)$ и объемного тока ${{J}_{{n - 1/2}}} \equiv {{J}_{q}}\left( {0.5({{z}_{{n - 1}}} + {{z}_{n}})} \right)$ припишем к середине ячейки. Номер пластины $q$ ради компактности записи будем опускать. В результате получим

(7)
${{H}_{{2n - 1}}} - {{H}_{{2n - 2}}} = \frac{{i\omega }}{{2c}}{{A}_{{n - 1/2}}}\Delta {{z}_{{n - 1/2}}}({{E}_{{2n - 1}}} + {{E}_{{2n - 2}}}) + \frac{{4\pi }}{c}{{J}_{{n - 1/2}}}\Delta {{z}_{{n - 1/2}}},\quad 1 \leqslant n \leqslant N.$

В качестве поверхности $S$ (2) выберем прямоугольник ${{z}_{{n - 1}}} \leqslant z \leqslant {{z}_{n}}$, $0 \leqslant x \leqslant 1$ в плоскости $y = 0$. Аналогично предыдущему получим

(8)
${{E}_{{2n - 1}}} - {{E}_{{2n - 2}}} = \frac{{i\omega }}{{2c}}{{\mu }_{{n - 1/2}}}\Delta {{z}_{{n - 1/2}}}({{H}_{{2n - 1}}} + {{H}_{{2n - 2}}}),\quad 1 \leqslant n \leqslant N.$

3. Запишем условия сопряжения в каждом внутреннем узле

(9)
${{E}_{{2n}}} = {{E}_{{2n - 1}}},\quad {{H}_{{2n}}} - {{H}_{{2n - 1}}} = - \frac{{4\pi }}{c}{{j}_{n}},\quad 1 \leqslant n \leqslant N - 1.$
Здесь поверхностный ток ${{j}_{n}} = j({{z}_{n}})$ относится ко внутреннему узлу. Таким образом, каждый узел сетки может быть границей раздела сред.

4. Наконец, построим аппроксимацию условий излучения (3). Выразим производную ${{\partial }_{z}}{{E}_{q}}$ из дифференциального уравнения Максвелла ${\text{rot}}{{{\mathbf{E}}}_{q}} = i\omega {{c}^{{ - 1}}}{{\mu }_{q}}{{{\mathbf{H}}}_{q}}$ и подставим значение ${{H}_{0}}$ на левой границе и ${{H}_{{2N - 1}}}$ на правой. После сокращения общих множителей получим

(10)
$\begin{array}{*{20}{c}} {\sqrt {{{\varepsilon }_{0}}} {{E}_{0}} + \sqrt {{{\mu }_{0}}} {{H}_{0}} = \sqrt {{{\varepsilon }_{0}}} {{E}^{0}},\quad z = 0} \\ {\sqrt {{{\varepsilon }_{N}}} {{E}_{{2N - 1}}} - \sqrt {{{\mu }_{N}}} {{H}_{{2N - 1}}} = \sqrt {{{\varepsilon }_{N}}} {{E}^{a}},\quad z = a.} \end{array}$
В первом из условий (10) значения ${{\varepsilon }_{0}} = \varepsilon (0)$ и ${{\mu }_{0}} = \mu (0)$ целесообразно взять при $z = 0$. Тогда граничное условие будет передано не приближенно, а точно. Аналогично, во втором условии мы полагаем ${{\varepsilon }_{N}} = \varepsilon (a)$, ${{\mu }_{N}} = \mu (a)$.

4.3. Алгебраическая система

1. Равенства (7)–(10) образуют систему $4N$ алгебраических уравнений относительно $4N$ неизвестных ${{E}_{{2n - 2}}}$, ${{E}_{{2n - 1}}}$, ${{H}_{{2n - 2}}}$, ${{H}_{{2n - 1}}}$. Она аналогична системе, приведенной в [1]. Эта система имеет 7-диагональную матрицу. Трудоемкость метода Гаусса для нее составляет $ \approx {\kern 1pt} 200N$ арифметических действий.

Порядок системы (7)–(10) можно значительно понизить. Для этого исключим из (9) величины ${{E}_{{2n - 1}}} = {{E}_{{2n}}}$ и ${{H}_{{2n - 1}}} = {{H}_{{2n}}} + 4\pi {{c}^{{ - 1}}}{{j}_{n}}$ при $1 \leqslant n \leqslant N - 1$. Тогда система (7)–(10) преобразуется к виду

(11)
$\sqrt {{{\varepsilon }_{0}}} {{E}_{0}} + \sqrt {{{\mu }_{0}}} {{H}_{0}} = \sqrt {{{\varepsilon }_{0}}} {{E}^{0}},$
(12)
${{H}_{{2n}}} - {{H}_{{2n - 2}}} - \frac{{i\omega }}{{2c}}{{A}_{{n - 1/2}}}\Delta {{z}_{{n - 1/2}}}({{E}_{{2n}}} + {{E}_{{2n - 2}}}) = \frac{{4\pi }}{c}({{J}_{{n - 1/2}}}\Delta {{z}_{{n - 1/2}}} - {{j}_{n}}),$
(13)
${{E}_{{2n}}} - {{E}_{{2n - 2}}} - \frac{{i\omega }}{{2c}}{{\mu }_{{n - 1/2}}}\Delta {{z}_{{n - 1/2}}}({{H}_{{2n - 1}}} + {{H}_{{2n - 2}}}) = \frac{{4\pi i\omega }}{{2{{c}^{2}}}}{{\mu }_{{n - 1/2}}}\Delta {{z}_{{n - 1/2}}}{{j}_{n}},$
(14)
${{H}_{{2N - 1}}} - {{H}_{{2N - 2}}} - \frac{{i\omega }}{{2c}}{{A}_{{N - 1/2}}}\Delta {{z}_{{N - 1/2}}}({{E}_{{2N - 1}}} + {{E}_{{2N - 2}}}) = \frac{{4\pi }}{c}{{J}_{{N - 1/2}}}\Delta {{z}_{{N - 1/2}}},$
(15)
${{E}_{{2N - 1}}} - {{E}_{{2N - 2}}} - \frac{{i\omega }}{{2c}}{{\mu }_{{N - 1/2}}}\Delta {{z}_{{N - 1/2}}}({{H}_{{2N - 1}}} + {{H}_{{2N - 2}}}) = 0,$
(16)
$\sqrt {{{\varepsilon }_{N}}} {{E}_{{2N - 1}}} - \sqrt {{{\mu }_{N}}} {{H}_{{2N - 1}}} = \sqrt {{{\varepsilon }_{N}}} {{E}^{a}}.$
В уравнениях (12), (13) индекс $n$ пробегает значения от 1 до $N - 1$. Таким образом, система (11)–(16) содержит $2N + 2$ уравнения. Она и является разностной схемой для задачи (1)–(4).

2. Составим вектор неизвестных $\{ {{E}_{0}},{{H}_{0}},{{E}_{2}},{{H}_{2}}, \ldots ,{{E}_{{2N - 1}}},{{H}_{{2N - 1}}}\} $. Тогда система (11)–(16) имеет 5-диагональную матрицу. Трудоемкость метода Гаусса составляет $ \approx {\kern 1pt} 50N$ арифметических действий, что в 4 раза меньше, чем для системы из [1]. Такой выигрыш по трудоемкости является существенным.

Определитель системы (11)–(16) содержит различные комбинации степеней шагов $\Delta {{z}_{{n - 1/2}}}$. Нам удалось вычислить главный член этого определителя для произвольной неравномерной сетки и неоднородной среды

(17)
${\text{Det}} = \sqrt {{{\varepsilon }_{0}}{{\mu }_{N}}} + \sqrt {{{\varepsilon }_{N}}{{\mu }_{0}}} + ia + O(\Delta {{z}^{2}}).$
Если $\Delta z$ достаточно мало, то определитель отличен от нуля. Поэтому система (11)–(16) имеет и при том единственное решение.

4.4. Аппроксимация

Будем считать, что внутри $q$-й пластины функции ${{\varepsilon }_{q}}(z)$, ${{\mu }_{q}}(z)$, ${{J}_{q}}(z)$ являются дважды непрерывно дифференцируемыми. Тогда точное решение может иметь особенности (разрывы) только на границах раздела, а внутри пластин также имеет вторые непрерывные производные.

Напомним, что разностная схема строится на специальной сетке, узлами которой являются все границы раздела. Шаблон схемы использует только один шаг по пространству. Такие схемы называются бикомпактными. Для системы уравнений Максвелла такие схемы ранее не предлагались.

Таким образом, внутри шаблона решение является гладким, и для него все квадратурные формулы имеют второй порядок аппроксимации. Ошибка $O(\Delta {{z}^{2}})$ вносится лишь при вычислении поверхностных интегралов типа (6). Контурные интегралы вычислены точно. Сеточные условия сопряжения и граничные условия также являются точными.

Введем класс функций, которые а) дважды непрерывно дифференцируемы всюду, кроме границ раздела ${{\xi }_{q}}$, и б) могут иметь разрывы в точках ${{\xi }_{q}}$, причем величины этих разрывов подчиняются условиям сопряжения (9). Из сказанного выше следует, что построенная схема имеет в данном классе функций аппроксимацию $O(\Delta {{z}^{2}})$.

4.5. Устойчивость

Для простоты сначала проведем доказательство для случая непроводящей среды $\sigma = 0$. Далее обобщим результат на случай $\sigma \ne 0$.

1. Обоснуем устойчивость по граничным условиям. Пусть сначала ${{J}_{{n + 1/2}}} = 0$, ${{j}_{n}} = 0$, ${{E}^{a}} = 0$, ${{E}^{0}} \ne 0$. Рассмотрим первый шаг сетки ${{h}_{{1/2}}}$. Очевидно, амплитуда поля $E$ в узле ${{z}_{0}}$ равна ${{E}_{0}} = {{E}^{0}}$. Найдем, во сколько раз поле увеличивается за шаг $\Delta {{z}_{{1/2}}}$, предполагая, что на правой границе этого шага отражение отсутствует. В этом случае $N = 1$, и разностная схема состоит из уравнений (11), (14)–(16). Ее явное решение имеет вид

(18)
$\begin{array}{*{20}{c}} {{{E}_{0}} = {{E}^{0}},\quad {{E}_{2}} = \frac{{1 - 0.5\omega \Delta {{z}_{{1/2}}}{{c}^{{ - 1}}}({\text{n}}_{{1/2}}^{{''}} - i{\text{n}}_{{1/2}}^{'})}}{{1 + 0.5\omega {{\Delta }_{{1/2}}}{{c}^{{ - 1}}}({\text{n}}_{{1/2}}^{{''}} - i{\text{n}}_{{1/2}}^{'})}}{{E}_{0}} \equiv {{A}_{{1/2}}}{{E}_{0}},} \\ {{{H}_{{0,2}}} = \sqrt {\frac{{{{\varepsilon }_{{1/2}}}}}{{{{\mu }_{{1/2}}}}}} {{E}_{{0,2}}}.} \end{array}$
Здесь ${\text{n}} = {\text{n}}{\kern 1pt} '\; + i{\text{n}}{\kern 1pt} '{\kern 1pt} ' = \sqrt {\varepsilon \mu } $ – показатель преломления.

Если отражение в точке ${{z}_{1}}$ есть, то амплитуда прошедшей волны заведомо не превосходит ${\text{|}}{{E}_{2}}{\text{|}}$. Аналогично можно оценить нарастание амплитуды поля за шаг $\Delta {{z}_{{3/2}}}$:

(19)
$\left| {{{E}_{4}}} \right| \leqslant \left| {{{A}_{{3/2}}}{{E}_{2}}} \right| \leqslant \left| {{{A}_{{3/2}}}{{A}_{{1/2}}}{{E}^{0}}} \right|,\quad \left| {{{H}_{4}}} \right| = \left| {\sqrt {\frac{{{{\varepsilon }_{{3/2}}}}}{{{{\mu }_{{3/2}}}}}} {{E}_{4}}} \right| \leqslant \left| {\sqrt {\frac{{{{\varepsilon }_{{3/2}}}}}{{{{\mu }_{{3/2}}}}}} {{A}_{{3/2}}}{{A}_{{1/2}}}{{E}^{0}}} \right|.$
Продолжая эти рассуждения, найдем оценку величины поля ${{E}_{{2n}}}$ во всех узлах ${{z}_{n}}$
(20)
$\left| {{{E}_{{2n}}}} \right| \leqslant \prod\limits_n \,\left| {{{A}_{{n + 1/2}}}} \right|\left| {{{E}^{0}}} \right|,\quad \left| {{{H}_{{2n}}}} \right| \leqslant \mathop {max}\limits_n \left| {\sqrt {\frac{{{{\varepsilon }_{{n + 1/2}}}}}{{{{\mu }_{{n + 1/2}}}}}} } \right|\prod\limits_n \,\left| {{{A}_{{n + 1/2}}}} \right|\left| {{{E}^{0}}} \right|.$
Реальные среды являются поглощающими, и для них ${\text{n}}{\kern 1pt} '{\kern 1pt} ' > 0$. Легко видеть, что в этом случае величина $\left| {{{A}_{{n + 1/2}}}} \right|$ не превосходит единицы. Поэтому
(21)
$\left| {{{E}_{{2n}}}} \right| \leqslant \left| {{{E}^{0}}} \right|w$
при всех n. Отсюда следует устойчивость по граничному условию при $z = 0$. Аналогично доказывается устойчивость по граничному условию при $z = a$.

2. Докажем устойчивость по правой части. Рассмотрим объемные токи. Пусть сначала ${{E}^{0}} = {{E}^{a}} = 0$, $j = 0$, ${{J}_{{n + 1/2}}} \ne 0$ при некотором $n = {{n}_{0}}$, а все остальные значения ${{J}_{{n + 1/2}}} = 0$. Этот ток возбуждает поля ${{E}_{{2{{n}_{0}}}}}$, ${{H}_{{2{{n}_{0}}}}}$ в узле ${{n}_{0}}$ и ${{E}_{{2{{n}_{0}} + 2}}}$, ${{H}_{{2{{n}_{0}} + 2}}}$ в узле ${{n}_{0}} + 1$, которые далее распространяются в положительном и отрицательном направлении оси $z$. Поэтому устойчивость решения относительно малых возмущений ${{J}_{{{{n}_{0}} + 1/2}}}$ сводится к а) устойчивости значений полей ${{E}_{{2{{n}_{0}}}}}$, ${{H}_{{2{{n}_{0}}}}}$, ${{E}_{{2{{n}_{0}} + 2}}}$, ${{H}_{{2{{n}_{0}} + 2}}}$ и б) устойчивости волн, создаваемых этими полями и распространяющихся в положительном направлении оси $z$.

Выразим ${{E}_{{2{{n}_{0}}}}}$, ${{H}_{{2{{n}_{0}}}}}$, ${{E}_{{2{{n}_{0}} + 2}}}$, ${{H}_{{2{{n}_{0}} + 2}}}$ через ${{J}_{{{{n}_{0}} + 1/2}}}$, считая, что за пределы данной ячейки волна уходит без отражения. Как и в предыдущем случае, разностная схема состоит из уравнений (11), (14)–(16). Приведем ее явное решение

(22)
$\begin{array}{*{20}{c}} {{{H}_{{2{{n}_{0}}}}} = \frac{{2\pi }}{c}\frac{{{{J}_{{{{n}_{0}} + 1/2}}}\Delta {{z}_{{{{n}_{0}} + 1/2}}}}}{{1 - 0.5i\omega \Delta {{z}_{{{{n}_{0}} + 1/2}}}{{c}^{{ - 1}}}{{{\text{n}}}_{{{{n}_{0}} + 1/2}}}}},\quad {{H}_{{2{{n}_{0}} + 2}}} = - {{H}_{{2{{n}_{0}}}}},} \\ {{{E}_{{2{{n}_{0}}}}} = {{E}_{{2{{n}_{0}} + 2}}} = {{H}_{{2{{n}_{0}} + 2}}}\sqrt {\frac{{{{\mu }_{{{{n}_{0}} + 1/2}}}}}{{{{\varepsilon }_{{{{n}_{0}} + 1/2}}}}}} .} \end{array}$
Эти значения можно рассматривать как источники волн, распространяющихся в положительном и отрицательном направлении оси $z$. Для этих волн справедлива оценка $\left| {{{E}_{{2n}}}} \right| \leqslant \left| {{{E}_{{2{{n}_{0}}}}}} \right|$. Тем самым, для всех $n$ имеют место неравенства

(23)
$\left| {{{H}_{{2n}}}} \right| \leqslant \frac{{2\pi }}{c}\frac{{\left| {{{J}_{{{{n}_{0}} + 1/2}}}} \right|\Delta {{z}_{{{{n}_{0}} + 1/2}}}}}{{\left| {1 - 0.5i\omega \Delta {{z}_{{{{n}_{0}} + 1/2}}}{{c}^{{ - 1}}}{{{\text{n}}}_{{{{n}_{0}} + 1/2}}}} \right|}} \equiv B,\quad \left| {{{E}_{{2n}}}} \right| \leqslant B\mathop {max}\limits_n \left| {\sqrt {\frac{{{{\mu }_{{n + 1/2}}}}}{{{{\varepsilon }_{{n + 1/2}}}}}} } \right|.$

Если объемные токи присутствуют в нескольких ячейках, то в выражениях (23) нужно провести суммирование по соответствующим ${{n}_{0}}$. Это дает

(24)
$\left| {{{H}_{{2n}}}} \right| \leqslant \frac{{2\pi a}}{c}\mathop {max}\limits_n \left| {{{J}_{{n + 1/2}}}} \right| \equiv C,\quad \left| {{{E}_{{2n}}}} \right| \leqslant C\mathop {max}\limits_n \left| {\sqrt {\frac{{{{\mu }_{{n + 1/2}}}}}{{{{\varepsilon }_{{n + 1/2}}}}}} } \right|.$
Здесь учтено то, что для реальных физических сред ${\text{n}}{\kern 1pt} {\text{'}} \leqslant 1$, ${\text{n}}{\kern 1pt} {\text{''}} \leqslant 0$; поэтому знаменатель (23) можно оценить снизу единицей. Из (24) следует устойчивость решения относительно малых возмущений объемных токов.

3. Полностью аналогично обосновывается устойчивость решения относительно малых возмущений поверхностных токов. Пусть сначала ${{E}^{0}} = {{E}^{a}} = 0$, $J = 0$, ${{j}_{{{{n}_{0}}}}} \ne 0$, ${{j}_{n}} = 0$ при $n \ne {{n}_{0}}$. Этот ток возбуждает поля ${{E}_{{2{{n}_{0}}}}}$, ${{H}_{{2{{n}_{0}}}}}$ в узле ${{n}_{0}}$, которые далее распространяются в положительном и отрицательном направлении оси $z$. Для установления связи между указанными значениями полей и током ${{j}_{{{{n}_{0}}}}}$ рассмотрим среду, состоящую из двух ячеек: $\Delta {{z}_{{{{n}_{0}} - 1/2}}}$ и $\Delta {{z}_{{{{n}_{0}} + 1/2}}}$. Подчеркнем, что левое и правое предельные значения поля $H$ в узле ${{n}_{0}}$ будут различны: ${{H}_{{2{{n}_{0}} - 1}}} \ne {{H}_{{2{{n}_{0}}}}}$. При этом левое и правое предельные значения поля $E = {{E}_{{2{{n}_{0}}}}}$ одинаковы. Будем считать, что за пределы ячеек волна уходит без отражения. Разностная схема содержит граничные условия (11), (16), одну пару уравнений (12), (13) и пару уравнений (14), (15). После несложных, но громоздких вычислений получим

(25)
${{H}_{{2{{n}_{0}}}}} = \frac{{2\pi }}{c}{{\mu }_{{{{n}_{0}} + 1/2}}}{{j}_{{{{n}_{0}}}}},\quad {{H}_{{2{{n}_{0}} - 1}}} = \frac{{2\pi }}{c}({{\mu }_{{{{n}_{0}} + 1/2}}} + 2){{j}_{{{{n}_{0}}}}},\quad {{E}_{{2{{n}_{0}}}}} = \sqrt {\frac{{{{\mu }_{{{{n}_{0}} + 1/2}}}}}{{{{\varepsilon }_{{{{n}_{0}} + 1/2}}}}}} {{H}_{{2{{n}_{0}}}}}.$
Созданные этими источниками волны распространяются в положительном и отрицательном направлении оси $z$. Поэтому для правого и левого предельных значений полей во всех узлах справедливы оценки
(26)
$\left| {{{H}_{{2n}}}} \right|,\;\left| {{{H}_{{2n - 1}}}} \right| \leqslant \frac{{2\pi }}{c}(\mathop {max}\limits_n {\text{|}}{{\mu }_{{n + 1/2}}}{\kern 1pt} {\text{|}} + 2){\kern 1pt} {\text{|}}{\kern 1pt} {{j}_{{{{n}_{0}}}}}{\kern 1pt} {\text{|}} \equiv D,\quad {\text{|}}{{E}_{{2n}}}{\kern 1pt} {\text{|}} \leqslant \mathop {max}\limits_n \left| {\sqrt {\frac{{{{\mu }_{{{{n}_{0}} + 1/2}}}}}{{{{\varepsilon }_{{{{n}_{0}} + 1/2}}}}}} } \right|D.$
Если поверхностный ток имеется на нескольких границах раздела, то оценки (26) нужно просуммировать по соответствующим ${{n}_{0}}$. Заметим, что эта сумма будет включать фиксированное число слагаемых, не зависящее от числа шагов сетки. Из оценки (26) следует устойчивость решения относительно малых возмущений поверхностных токов $j$.

4. Чтобы учесть проводимость, во всех предыдущих формулах достаточно сделать замену $ - i\omega \varepsilon {\text{/}}c \to - i\omega \varepsilon {\text{/}}c + 4\pi \sigma {\text{/}}c$. При этом доказательство устойчивости проводится полностью аналогично, но выкладки оказываются более громоздкими.

5. Таким образом, проведено полное обоснование устойчивости разностной схемы по граничным условиям и правой части. Оно справедливо для произвольных неравномерных сеток и неоднородных сред. Отсюда и из однозначной разрешимости системы уравнений (11)–(16) следует корректность предложенной разностной схемы.

4.6. Сходимость

1. Согласно известным теоремам теории разностных схем [36], из корректности разностной схемы и установленной выше аппроксимации вытекает сходимость. Напомним, что аппроксимация имеет место на классе функций, имеющих разрывы на границах раздела сред, причем величина разрыва должна подчиняться условию сопряжения. При этом необходимо использовать специальные сетки. Поэтому справедлива

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

2. Для предложенной схемы применим метод сгущения сеток и оценки точности по правилу Ричардсона [87]–[89]. В литературе предлагались разные многосеточные методы (т.н. sub-gridding techniques), см., например, [90]–[92]. Эти методы используют локальное сгущение для построения адаптивных сеток. Однако эти методы не дают гарантированных оценок точности. Напротив, глобальное сгущение сеток по методу Ричардсона дает асимптотически точное значение погрешности. Этот метод существенно повышает надежность расчета. Он имеет строгое обоснование [93].

3. Схема фактически построена интегроинтерполяционным методом. Поэтому для сред с непрерывно изменяющимися свойствами консервативность схемы следует из аддитивности интегралов. При наличии границ раздела консервативность обеспечивается явным учетом условий сопряжения. Таким образом, схема является полностью консервативной в смысле определения Тихонова-Самарского.

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

5. НЕСТАЦИОНАРНАЯ ЗАДАЧА

5.1. Постановка задачи

В постановке (1)–(4) следует заменить множитель $( - i\omega )$ на производную ${{\partial }_{t}}$. Объемные ${{J}_{q}}(z,t)$ и поверхностные ${{j}_{q}}(z,t)$ токи могут произвольно зависеть от времени. Слева и справа на рассеиватель падают волновые пакеты ${{{\mathbf{E}}}^{0}}(t - z{\text{/}}c)exp( - i{{\omega }^{0}}t + i{{k}^{0}}z)$, ${{{\mathbf{E}}}^{a}}(t + z{\text{/}}c)exp( - i{{\omega }^{0}}t - i{{k}^{0}}z)$ с несущей частотой ${{\omega }^{0}}$ и заданными огибающими ${{E}^{{0,a}}}$. Будем считать огибающие финитными функциями, а поля в начальный момент времени – нулевыми. Таким образом, постановка нестационарной задачи имеет следующий вид:

(27)
$\int\limits_\Gamma {{{{\mathbf{H}}}_{q}}d{\mathbf{l}}} = \frac{{4\pi }}{c}\int\limits_S {{{{\mathbf{J}}}_{q}}d{\mathbf{s}}} + \frac{1}{c}\frac{\partial }{{\partial t}}\int\limits_S {{{{\mathbf{D}}}_{q}}d{\mathbf{s}}} ,\quad 1 \leqslant q \leqslant Q;$
(28)
$\int\limits_\Gamma {{{{\mathbf{E}}}_{q}}d{\mathbf{l}}} = - \frac{1}{c}\frac{\partial }{{\partial t}}\int\limits_S {{{{\mathbf{B}}}_{q}}d{\mathbf{s}}} ,\quad 1 \leqslant q \leqslant Q;$
(29)
$\begin{array}{*{20}{c}} {\frac{\partial }{{\partial z}}{{{\mathbf{E}}}_{1}} - \frac{1}{c}\frac{\partial }{{\partial t}}{{{\mathbf{E}}}_{1}} = 2ik{{{\mathbf{E}}}^{0}}(t - z{\text{/}}c)exp( - i{{\omega }^{0}}t + i{{k}^{0}}z),\quad z = 0;} \\ {\frac{\partial }{{\partial z}}{{{\mathbf{E}}}_{Q}} + \frac{1}{c}\frac{\partial }{{\partial t}}{{{\mathbf{E}}}_{Q}} = 2ik{{{\mathbf{E}}}^{a}}(t + z{\text{/}}c)exp( - i{{\omega }^{0}}t - i{{k}^{0}}z),\quad z = a;} \end{array}$
(30)
$\begin{gathered} {{{\mathbf{e}}}_{z}} \times ({{{\mathbf{E}}}_{q}} - {{{\mathbf{E}}}_{{q - 1}}}) = 0,\quad {{{\mathbf{e}}}_{z}}({{{\mathbf{D}}}_{q}} - {{{\mathbf{D}}}_{{q - 1}}}) = 0, \\ {{{\mathbf{e}}}_{z}} \times ({{{\mathbf{H}}}_{q}} - {{{\mathbf{H}}}_{{q - 1}}}) = \frac{{4\pi }}{c}{{{\mathbf{j}}}_{q}},\quad {{{\mathbf{e}}}_{z}}({{{\mathbf{B}}}_{q}} - {{{\mathbf{B}}}_{{q - 1}}}) = 0; \\ \end{gathered} $
${{{\mathbf{E}}}_{q}} = 0,\quad {{{\mathbf{H}}}_{q}} = 0,\quad t = 0,\quad 1 \leqslant q \leqslant Q,\quad t = 0.$
Условия излучения (29) интерпретируются аналогично стационарному случаю. При этом учитывается, что решения, соответствующие распространяющимся волнам, являются автомодельными ${{{\mathbf{E}}}_{ \pm }} = {{{\mathbf{E}}}_{ \pm }}(z \mp ct)$. Они обобщают известные условия Мура [47]
(32)
$\frac{\partial }{{\partial t}}{\mathbf{E}} - c\frac{\partial }{{\partial z}}{\mathbf{E}} = 0,\quad z = 0;\quad \frac{\partial }{{\partial t}}{\mathbf{E}} + c\frac{\partial }{{\partial z}}{\mathbf{E}} = 0,\quad z = a.$
Условия Мура также исходят из автомодельного вида решения, но они описывают уход волны на бесконечность. Они применимы, если источник электромагнитных волн является локализованным. Однако для задач, в которых падающая волна приходит из бесконечности, такие условия нефизичны. В этом случае в литературе применяют модель фиктивного источника. Условия (3) и (29) записываются более естественно. Они в литературе не описаны.

При распространении волнового пакета в линейной диспергирующей среде для каждой гармоники реализуются свои $\varepsilon $ и $\mu $ и своя скорость распространения. Поэтому пакет следует разбить на монохроматические компоненты, решить задачу (1)–(4) для каждой из них в отдельности и просуммировать полученные решения. Разбиение пакета на компоненты сводится к прямому преобразованию Фурье, суммирование решений – к обратному.

5.2. Метод спектрального разложения

1. При распространении волнового пакета в линейной диспергирующей среде для разных спектральных компонент решения реализуются разные значения $\varepsilon $ и $\mu $ и, соответственно, разные скорости распространения. Разложим пакет на монохроматические компоненты, решим стационарную задачу для каждой компоненты и просуммируем полученные решения. Спектральное разложение исходного пакета сводится к прямому преобразованию Фурье, суммирование полученных решений – к обратному преобразованию Фурье.

2. Выполним численное преобразование Фурье падающих волновых пакетов, объемных и поверхностных токов по формуле трапеций, используя одинаковые сетки по времени и наборы частот $\{ {{\omega }_{m}}\} $, ${{\omega }_{{m + 1}}} - {{\omega }_{m}} = \Delta {{\omega }_{m}}$. Для каждой частоты получим амплитуды падающих волн $E_{m}^{{0,a}}$ и источников поля ${{({{J}_{q}})}_{m}}$, ${{({{j}_{q}})}_{m}}$. По описанному выше методу решим набор стационарных задач, соответствующих различным частотам ${{\omega }_{m}}$. Затем просуммируем полученные спектральные амплитуды решения по всем частотам ${{\omega }_{m}}$. Это дает решение $E(z,t)$, $H(z,t)$ исходной нестационарной задачи.

Назовем предложенный алгоритм методом спектрального разложения.

5.3. Сходимость

1. Если функции ${{E}^{{0,a}}}$, ${{j}_{q}}(z,t)$, ${{J}_{q}}(z,t)$ имеют вторые непрерывные производные, то квадратуры, выражающие прямое и обратное преобразования Фурье, сходятся со вторым порядком точности: $O(\Delta {{\omega }^{2}})$ и $O(\Delta {{t}^{2}})$ соответственно. Как показано выше, аппроксимация разностных схем для стационарных задач есть $O(\Delta {{z}^{2}})$. Поэтому аппроксимация схемы есть $O(\Delta {{z}^{2}} + \Delta {{\omega }^{2}} + \Delta {{t}^{2}})$.

2. Покажем, что метод спектрального разложения устойчив относительно малых возмущений входных данных. Предварительно заметим, что в реальных задачах падающие импульсы и токи являются локализованными функциями времени: отрезок $T$, на котором они существенно отличны от нуля, является конечным. Соответственно, интервал частот $\Omega $, на котором отличны от нуля их фурье-образы, также конечен. Чаще всего используется временная модуляция гауссовой огибающей. Хотя она формально отлична от нуля на всей прямой $ - \infty < t < + \infty $, однако такая огибающая очень быстро убывает по мере удаления от максимума. Для получения физически разумной точности 0.1% достаточно взять $T \sim 3\sigma $, где $\sigma $ – ширина огибающей на половине высоты. Машинная точность ${{10}^{{ - 16}}}$ достигается при $T \sim 6\sigma $.

Пусть в одну из функций ${{E}^{{0,a}}}$, ${{j}_{q}}$, ${{J}_{q}}$ внесена ошибка $\delta $, зависящая от $t$ и, возможно, $z$. Погрешность соответствующего фурье-образа равна

(33)
$\tilde {\delta }(\omega ,z) = \frac{1}{{\sqrt {2\pi } }}\int\limits_{ - \infty }^{ + \infty } {\delta (t,z){{e}^{{ - i\omega t}}}dt} \approx \sum\limits_{k = 0}^K \,\delta ({{t}_{k}},z){{e}^{{ - i\omega {{t}_{k}}}}}{{g}_{k}}\Delta {{t}_{k}}.$
Здесь ${{g}_{0}} = {{g}_{K}} = 0.5$, ${{g}_{k}} = 1$, $k \ne 1,K$ – веса квадратурной формулы трапеций. Легко видеть, что погрешность (33) не превосходит
(34)
$\left| {\tilde {\delta }} \right| \leqslant \frac{1}{{\sqrt {2\pi } }}max\left| \delta \right|T.$
Такое возмущение входных данных приводит к возмущению $\delta \tilde {E}$, $\delta \tilde {H}$ решений стационарных задач, для которого справедливы оценки (21), (24), (26), в которые вместо исходных входных данных следует подставить величину (34).

Выполним обратное преобразование Фурье погрешностей $\delta \tilde {E}$, $\delta \tilde {H}$:

(35)
$\delta E(z,t) = \frac{1}{{\sqrt {2\pi } }}\int\limits_{ - \infty }^{ + \infty } {\delta \tilde {E}(\omega ,z){{e}^{{i\omega t}}}dt} \approx \sum\limits_{m = 0}^M \,\delta \tilde {E}({{\omega }_{m}},z){{g}_{m}}\Delta {{\omega }_{m}}.$
Справедлива следующая оценка:
(36)
$\left| {\delta E} \right| \leqslant \frac{1}{{2\pi }}{{C}_{{max}}}max\left| \delta \right|T\Omega .$
Здесь ${{C}_{{max}}}$ – наибольшая из мажорирующих констант в оценках (21), (24), (26). Аналогично выписывается оценка для погрешности $\delta H$. Это завершает обоснование устойчивости предложенной схемы. Оно проведено для произвольных неравномерных сеток по всем переменным и неоднородных сред (в том числе, диспергирующих).

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

Теорема 2. Схема метода спектрального разложения сходится со вторым порядком точности.

4. Для предложенного метода применимы оценки точности по правилу Ричардсона и экстраполяционное уточнение. Заметим, что сетки по переменным $z$, $t$, $\omega $ необходимо сгущать одновременно и в одинаковое число раз.

Несмотря на простоту, метод спектрального разложения является принципиально новым. Он позволяет учитывать произвольный закон дисперсии, включая таблично заданный. Этот метод рекомендуется использовать как общий подход построения разностных схем для гиперболических задач в линейной диспергирующей среде.

Авторы искренне благодарны Н.Н. Калиткину за ценные замечания и обсуждения и А.Н. Боголюбову за внимание к работе.

Список литературы

  1. Белов А.А., Домбровская Ж.О. Бикомпактная разностная схема для уравнений Максвелла в слоистых средах // Докл. АН. 2020. Т. 492. С. 15–19.

  2. Hansen W.N. Electric fields produced by the propagation of plane coherent electromagnetic radiation in a stratified medium // J. Opt. Soc. Am. 1968. V. 58. № 9. P. 380–390.

  3. Berreman D.W. Optics in stratified and anisotropic media: 4 × 4-matrix formulation // J. Opt. Soc. Am. 1972. V. 62. № 9. P. 502–510.

  4. Cotter N.P.K., Preist T.W., Sambles J.R. Scattering-matrix approach to multilayer diffraction // J. Opt. Soc. Am. A. May 1995. V. 12. № 9. P. 1097–1103.

  5. Tishchenko A.V., Hamdoun M., Parriaux O. Two-dimensional coupled mode equation for grating waveguide excitation by a focused beam // Opt. Quant. Electron. May 2003. V. 35. № 1. P. 475–491.

  6. Fornberg B. A Practical Guide to Pseudospectral Methods. Cambridge University Press, Cambridge. 1996.

  7. Wriedt T. Generalized Multipole Techniques for Electromagnetic and Light Scatterin. Elsevier Science, Amsterdam. 1999.

  8. Yurkin M.A., Maltsev V.P., Hoekstra A.G. Convergence of the discrete dipole approximation. ii. an extrapolation technique to increase the accuracy // J. Opt. Soc. Am. A. May 2006. V. 23. № 1. P. 2592–2601.

  9. Yurkin M.A., Huntermann M. Rigorous and fast discrete dipole approximation for particles near a plane interface // J. Phys. Chem. C. May 2015. V. 119. № 1. P. 29088–29094.

  10. Shcherbakov A.A., Tishchenko A.V. Generalized source method in curvilinear coordinates for 2d grating diffraction simulation // JQSRT. May 2017. V. 187. № 1. P. 76–96.

  11. Albani M., Bernardi P. A numerical method based on the discretization of maxwell equations in integral form // IEEE Trans. Microwave Theory Techn. May 1974. V. 22. № 4. P. 446–450.

  12. Christ A., Hartnagel H.L. Three-dimensional finite-difference method for the analysis of microwave-device embedding // IEEE Trans. Microwave Theory Techn. May 1987. V. 35. № 8. P. 688–696.

  13. Beilenhoff K., Heinrich W., Hartnagel H.L. Improved finite-difference formulation in frequency domain for three-dimensional scattering problems // IEEE Trans. Microwave Theory Techniques. May 1992. V. 40. № 3. P. 540–546.

  14. Wang J.-S., Mittra R. A finite element cavity resonance method for waveguide and microstrip line discontinuity problems // IEEE Trans. Microwave Theory Techniques. May 1994. V. 42. № 3. P. 433–440.

  15. Ivinskaya A.M., Lavrinenko A.V., Shyroki D.M. Modeling of nanophotonic resonators with the finite-difference frequency-domain method // IEEE Trans. Antennas Propag. 2011. V. 59. № 3. P. 4155–4161.

  16. Кудрявцев А.Н., Трашкеев С.И. Формализм двух потенциалов для численного решения уравнений Максвелла // Ж. вычисл. матем. и матем. физ. 2013. Т. 153. № 11. С. 1823–1834.

  17. Вязников К.В., Тишкин В.Ф., Фаворский А.П. Построение монотонных разностных схем повышенного порядка аппроксимации для систем уравнений гиперболического типа // Матем. моделирование. 1989. Т. 1. № 5. С. 95–120.

  18. Harten A., Osher S. Uniformly high-order accurate nonoscillatory schemes. i // SIAM J. Numer. Anal. 1987. V. 24. № 2. P. 279–309.

  19. Liu X.-D., Osher S., Chan T. Weighted essentially non-oscillatory schemes // J. Comput. Phys. 1994. V. 115. № 1. P. 200–212.

  20. Shlager K.L., Schneider J.B. A selective survey of the finite-difference timedomain literature // IEEE Antennas Propagat. Mag. 1995. V. 37. № 4. P. 39–57.

  21. Sullivan D.M. Electromagnetic simulation using the FDTD method. IEEE Press. 2000.

  22. Inan U.S., Marshall R.A. Numerical electromagnetics. The FDTD method. Cambridge University Press, Cambridge. 2011.

  23. Taflove A., Johnson S.G., Oskooi A. Advances in FDTD Computational Electromagnetics: Photonics and Nanotechnology. Artech House, London. 2013.

  24. Yee K. Numerical solution of inital boundary value problems involving maxwell’s equations in isotropic media // IEEE Trans. Antennas. Propag. 1966. V. 14. № 3. P. 302–307.

  25. Калиткин Н.Н. Улучшенная факторизация параболических схем // Докл. АН. Информатика. 2005. Т. 402. № 4. С. 467–471.

  26. Cangellaris A.C., Wright D.B. Analysis of the numerical error caused by the stair-stepped approximation of a conducting boundary in fdtd simulations of electromagnetic phenomena // IEEE Trans. Antennas. Propag. 1991. V. 39. № 10. P. 1518–1525.

  27. Домбровская Ж.О., Боголюбов А.Н. Немонотонность схемы fdtd при моделировании границ раздела между диэлектриками // Уч. зап. физ. факультета МГУ. 2017. № 4. С. 1740302.

  28. Werner G.R., Cary J.R. A stable fdtd algorithm for non-diagonal, anisotropic dielectrics // J. Comp. Phys. 2007. V. 226. P. 1085–1101.

  29. Oskooi A.F., Kottke C., Johnson S.G. Accurate finite-difference time-domain simulation of anisotropic media by subpixel smoothing // Opt. Lett. 2009. V. 34. № 18. P. 2778–2780.

  30. Bauer C.A., Werner G.R., Cary J.R. A second-order 3d electromagnetic algorithm for curved interfaces between anisotropic dielectrics on a yee mesh // J. Comput. Phys. 2011. V. 230. № 5. P. 2060–2075.

  31. Werner G.R., Bauer C.A., Cary J.R. A more accurate, stable, fdtd algorithm for electromagnetics in anisotropic dielectrics // J. Comput. Phys. 2013. V. 255. P. 436–455.

  32. Hirono T., Shibata Y., Lui W.W., Seki S., Yoshikuni Y. The second-order condition for the dielectric interface orthogonal to the yee-lattice axis in the fdtd scheme // IEEE Microwave Guided Wave Lett. 2000. V. 10. № 9. P. 359–361.

  33. Hwang K.P., Cangellaris A.C. Effective permittivities for second-order accurate fdtd equations at dielectric interfaces // IEEE Microw. Wireless Compon. Lett. 2001. V. 11. № 4. P. 158–160.

  34. Armenta R.B., Sarris C.D. A second-order domain-decomposition method for modeling material interfaces in finite-difference discretizations // Proc. of the IEEE/MTT-S International Symposium. 2012. P. 502–505.

  35. Taflove A., Hagness S.C. Computational Electrodynamics: The Finite-Difference Time-Domain Method. Artech House, London. 2005.

  36. Самарский А.А. Теория разностных схем. М.: Наука, 1989.

  37. Nédélec J.C. Mixed finite elements in R3 // Numer. Math. 1980. V. 35. P. 315–341.

  38. Kirby R.C., Logg A., Terrel A.R. Automated Solution of Differential Equations by the Finite Element Method. Lecture Notes in Computational Science and Engineering. Common and Unusual Finite Elements. Vol. 84. 2012.

  39. Hesthaven J.S., Gottlieb D. Stable spectral methods for conservation laws on triangles with unstructured grids // Comput. Methods Appl. Mech. Engng. 1999. V. 175. P. 361–381.

  40. Hesthaven J.S. Spectral penalty methods // Appl. Numer. Math. 2000. V. 33. P. 23–41.

  41. Hesthaven J.S., Teng C.H. Stable spectral methods on tetrahedral elements // SIAM J. Sci. Comput. 2000. V. 21. P. 2352–2380.

  42. Dridi K., Hesthaven J.S., Ditkowski A. Staircase-free finite-difference timedomain formulation for general materials in complex geometries // IEEE Trans. Antennas Propag. 2001. V. 49. № 5. P. 749–756.

  43. Hesthaven J.S., Warburton T. High-order nodal methods on unstructured grids. i. time-domain solution of maxwell’s equations // J. Comput. Phys. 2002. V. 181. P. 1–34.

  44. Piperno S., Fezoui L. A discontinuous galerkin fvtd method for 3d maxwell equations // CERMICS research report. 2003. P. 2003–4733.

  45. Hesthaven J.S., Warburton T. High order nodal discontinuous galerkin methods for the maxwell eigenvalue problem // Phil. Trans. R. Soc. Lond. A. 2004. V. 362. P. 493–524.

  46. Hesthaven J.S. High-order accurate methods in time-domain computational electromagnetics: A review // Advances in imaging and electron physics. 2003. V. 127. P. 59–123.

  47. Mur G. Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations // IEEE Trans. Electromagn. Comp., EMC. 1981. V. 23. № 4. P. 377–382.

  48. Berenger J.P. Three-dimensional perfectly matched layer for the absorption of electromagnetic waves // J. Comput. Phys. 1996. V. 127. № 2. P. 363–367.

  49. Sullivan D.M. Frequency-dependent fdtd methods using z transforms // IEEE Trans. Antennas Propagat. 1992. V. 40. № 10. P. 1223–1230.

  50. Sullivan D.M. Z-transform theory and the fdtd method // IEEE Trans. Antennas Propagat. 1996. V. 44. № 1. P. 28–34.

  51. Abdijalilov K., Grebel H. Z-transform theory and fdtd stability // IEEE Trans. Antennas Propagat. 2004. V. 52. № 11. P. 2950–2954.

  52. Kelley D.F., Destan T.J., Luebbers R.J. Debye function expansions of complex permittivity using a hybrid particle swarm-least squares optimization approach // IEEE Trans. Antennas Propagat. 2007. V. 55. № 7. P. 1999–2005.

  53. Lin Z., Fang Y., Hu J., Zhang C. On the fdtd formulations for modeling wideband lorentzian media // IEEE Trans. Antennas Propagat. 2011. V. 59. № 4. P. 1338–1346.

  54. Dong X.T., Venkatarayalu N.V., Guo B., Yin W.Y., Gan Y.B. General formulation of unconditionally stable adi-fdtd method in linear dispersive media // IEEE Trans. Microwave Theory Techn. 2004. V. 52. № 1. P. 170–174.

  55. Young J.L., Nelson R.O. A summary and systematic analysis of fdtd algorithms for linearly dispersive media // Antennas Propag. Magazine IEEE. 2001. V. 43. № 1. P. 61–126.

  56. Cai H., Hu X., Xiong B., Zhdanov M.S. Finite-element time-domain modeling of electromagnetic data in general dispersive medium using adaptive pade series // Antennas Propag. Magazine IEEE. 2017. V. 109. № 12. P. 194–205.

  57. Lin Z., Zhang C., Ou P., Jia Y., Feng L. A generally optimized fdtd model for simulating arbitrary dispersion based on the maclaurin series expansion // J. Lightwave Technol. 2010. V. 28. № 19. P. 2843–2850.

  58. Weedon W.H., Rappaport C.M. A general method for fdtd modeling of wave propagation in arbitrary frequency-dispersive media // IEEE Trans. Antennas Propagat. 1997. V. 45. № 3. P. 401–410.

  59. Pereda J.A., Vegas A., Prieto A. Fdtd modeling of wave propagation in dispersive media by using the mobius transformation technique // IEEE Trans. Microwave Theory Techn. 2002. V. 50. № 7. P. 1689–1695.

  60. Bolomey J.C., Durix C., Lesselier D. Time domain integral equation approch for inhomogeneous and dispersive slab problems // IEEE Trans. Antennas Propagat. 1978. V. AP-26. № 9. P. 658–667.

  61. Luebbers R., Hunsberger F.P., Kunz K.S., Standler R.B., Schneider M. A frequency-dependent finite-difference time-domain formulation for dispersive materials // IEEE Trans. Electromagn. Compatibility. 1990. V. 32. № 3. P. 222–227.

  62. Luebbers R., Hunsberger F.P., Kunz K.S. A frequency-dependent finitedifference time-domain formulation for transient propagation in plasmas // IEEE Trans. Antennas Propagat. 1991. V. 39. № 1. P. 29–43.

  63. Luebbers R., Hunsberger F.P. Fdtd for nth-order dispersive media // IEEE Trans. Antennas Propagat. 1992. V. 40. № 11. P. 1297–1301.

  64. Giannakis I., Giannopoulos A. A novel piecewise linear recursive convolution approach for dispersive media using the finite-difference time-domain method // IEEE Trans. Antennas Propagat. 2014. V. 62. № 5. P. 2669–2678.

  65. Young J.L. Propagation in linear dispersive media: Finite difference time-domain methodologies // IEEE Trans. Antennas Propagat. 1995. V. 43. № 4. P. 422–426.

  66. Kashiwa T., Yoshida N., Fukai I. A treatment by the finite-difference timedomain method of the dispersive characteristics associated with orientation polarization // Inst. Electron. Inform. Communicat. Eng. Trans. 1990. V. E73. № 8. P. 1326–1328.

  67. Kashiwa T., Fukai I. A treatment by the fd-td method for the dispersive characteristics associated with electronic polarization // Microwave Opt. Tech. Lett. 1990. V. 3. № 6. P. 203–205.

  68. Gandhi O.P. A frequency-dependent finite-difference time-domain formulation for general dispersive media // IEEE Trans. Microwave Theory Tech. 1993. V. 41. № 4. P. 658–665.

  69. Joseph R.M., Hagness S.C., Taflove A. Direct time integration of maxwell’s equations in linear dispersive media with absorption for scattering and propagation of femtosecond electromagnetic pulses // Opt. Lett. 1991. V. 16. № 18. P. 1412–1414.

  70. Maradei F. A frequency-dependent wetd formulation for dispersive materials // IEEE Trans. Magnet. 2001. V. 37. № 5. P. 3303–3306.

  71. Kobidze G., Gao J., Shanker B., Michielssen E. A fast time domain integral equation based scheme for analyzing scattering from dispersive objects // IEEE Trans. Magnet. 2005. V. 53. № 3. P. 1215–1226.

  72. Zhuansun X., Ma X. Integral-based exponential time differencing algorithms for general dispersive media and the cfs-pml // IEEE Trans. Antennas Propagat. 2012. V. 60. № 7. P. 3257–3264.

  73. Petropoulos P.G. Stability and phase error analysis of fdtd in dispersive dielectrics // IEEE Trans. Antennas Propagat. 1994. V. 42. № 1. P. 62–69.

  74. Young J.L., Kittichartphayak A., Kwok Y.M., Sullivan D. On the dispersion errors related to (fd)2/td type schemes // IEEE Trans. Microwave Theory Techn. 1995. V. 43. № 8. P. 1902–1910.

  75. Fomberg B. Some numerical techniques for maxwell’s equations in different types of geometries // Topics in Comput. Wave Propagat. Lecture Notes in Comput. Sci. and Engng, vol 31. 2003. P. 265–299.

  76. Fornberg B., Zuev J., Lee J. Stability and accuracy of time-extrapolated adi-fdtd methods for solving wave equations // J. Comput. Appl. Math. 2007. V. 200. P. 178–192.

  77. Garcia S.G., Lee T.W., Hagness S.C. On the accuracy of the adi-fdtd method // IEEE Antennas Wirel. Propag. Lett. 2002. V. 1. № 1. P. 31–34.

  78. Домбровская Ж.О., Боголюбов А.Н. Повышение точности одномерной схемы Йе методом сгущения сеток // Известия Российской академии наук. Серия физическая. 2017. Т. 81. № 1. С. 117–120.

  79. Толстых А.И. Компактные разностные схемы и их применение в задачах аэрогидродинамики. М.: Наука, 1990.

  80. Калиткин Н.Н., Корякин П.В. Бикомпактные схемы и слоистые среды // Докл. АН. 2008. Т. 419. № 6. С. 744–748.

  81. Калиткин Н.Н., Корякин П.В. Одномерные и двумерные бикомпактные схемы в слоистых средах // Матем. моделирование. 2009. Т. 21. № 8. С. 44–62.

  82. Калиткин Н.Н., Корякин П.В. Численные методы. Т. 2. Методы математической физики. М.: Академия, 2013.

  83. Тихонов А.Н., Самарский А.А. О сходимости разностных схем в классе разрывных коэффициентов // Докл. АН. 1959. Т. 8. № 3. С. 529–532.

  84. Самарский А.А. Введение в теорию разностных схем. М.: Наука, 1971.

  85. Самарский А.А., Попов Ю.П. Разностные методы решения задач газовой динамики. М.: Наука, 1992.

  86. Свешников А.Г. Принципы излучения // Докл. АН. 1950. Т. 3. № 5. С. 517–520.

  87. Richardson L.F., Gaunt J.A. The deferred approach to the limit // Phil. Trans. A. 1927. V. 226. P. 299–349.

  88. Марчук Г.И. Методы расщепления. М.: Наука, 1988.

  89. Калиткин Н.Н., Альшин А.Б., Альшина Е.А., Рогов Б.В. Вычисления на квазиравномерных сетках. М.: Физматлит, 2005.

  90. Meglicki Z., Gray S.K., Norris B. Multigrid fdtd with chombo // Comp. Phys. Commun. 2007. V. 176. P. 109–120.

  91. Van Londersele A., De Zutter D., Ginste D.V. An in-depth stability analysis of nonuniform fdtd combined with novel local implicitization techniques // J. Comp. Phys. 2017. V. 342. P. 177–193.

  92. Balsara D.S., Simpson J.J. Making a synthesis of fdtd and dgtd schemes for computational electromagnetics // IEEE J. Multiscale Multiphys. Comp. Techniques. 2020. V. 5. P. 99.

  93. Рябенький В.С., Филлипов А.Ф. Об устойчивости разностных уравнений. М.: Гостехиздат, 1956.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал