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

Полный текст

Введение

В задачах плазмоники, например при решении нелинейных комплексных дисперсионных уравнений (ДУ), возникает проблема нахождения их комплексных корней. При этом важно построение полной дисперсионной ветви: зависимостей фазовой постоянной (замедления) и постоянной затухания от частоты, например, для поверхностных плазмон-поляритонов (ПП). Мы рассматриваем среды и структуры с диссипацией, поэтому дисперсионные ветви непрерывны и существуют для всех частот, в том числе и в запрещенных для недиссипативных структур зонах. Если в разрешенных зонах при пренебрежимо малой диссипации ПП обычно имеет вид медленной поверхностной волны, а его дисперсия может быть найдена путем решения ДУ в действительной области, то в запрещенной зоне ПП обычно быстрые вытекающие, соответствующие антиповерхностной волне с большим радиационным затуханием. При реальной диссипации ПП классифицируются на втекающие и вытекающие [1‒5], что определяется углом вытекания $\vartheta =\arctan \left({k^{\text{'}}}_{0x}/{k^{\text{'}}}_z\right).$ Здесь мы считаем, что ПП распространяются вдоль оси z, лежащей на поверхности раздела x = 0 структура – вакуум. Сверху в вакууме $k_0^2=k_{0x}^2+k_z^2,$ $k_{0x}={k^{\text{'}}}_{0x}-i{k^{\text{'}\text{'}}}_{0x},$ причем мы ищем решение с зависимостью $\exp \left(i\omega t-i{k}_{z}z\right),$ $k_z={k^{\text{'}}}_z-i{k^{\text{'}\text{'}}}_z.$ При $x<0$ может находиться произвольная структура или среда, поддерживающая указанную волну. Волна прямая, если ${k^{\text{'}}}_z{k^{\text{'}\text{'}}}_z>0$ обратная, если ${k^{\text{'}}}_z{k^{\text{'}\text{'}}}_z<0$ [1; 2]. Случай отсутствия диссипации ${k^{\text{'}\text{'}}}_z=0$ требует для анализа характера волны рассмотрения дисперсии и вычисления групповой скорости. В общем случае определения характера волны следует использовать проинтегрированную по поперечному сечению структуры компоненту вектора Пойнтинга ${\tilde{S}}_z,$ и тогда критерием служит знак величины ${\tilde{S}}_z{k^{\text{'}}}_z$ [5]. Следует отметить, что для вытекающих волн интеграл ${\tilde{S}}_z$ может расходиться, и для определения плотности потока мощности его нужно вычислять в определенных границах, а затем переходить к пределу. Для втекания необходимо условие $\vartheta <0.$ Определяя $S_x,$ видим, что в этом случае энергия втекает из вакуума в структуру, тогда как при $\vartheta >0$ она вытекает из структуры в вакуум. Границей перехода от медленных к быстрым волнам служит условие «отсечки» $n^{\text{'}}={k^{\text{'}}}_z/k_0=1.$ Далее будем решать ДУ вида $k_z=D\left(k_0,k_z\right).$ Простейшее такое ДУ Ценнека имеет вид $k_z=k_0\sqrt{\epsilon \left(\omega \right)/\left(1+\epsilon \left(\omega \right)\right)}$ и описывает ПП над диэлектрическим полупространством. Оно дает явное решение, поэтому уравнением не является. Для него ПП всегда втекающий, причем для быстрой волны Ценнека при ${\epsilon }^{\text{'}}\left(\omega \right)>1$ втекание может быть весьма слабым, как и для весьма медленного ПП, когда ${\epsilon }^{\text{'}}\left(\omega \right)\approx -1$ и диссипация мала: ${\epsilon }^{\text{'}\text{'}}\left(\omega \right)/{\epsilon }^{\text{'}}\left(\omega \right)\ll 1$ [4]. В общем случае правая часть ДУ нелинейно зависит от $k_z$ даже если возможно описание структуры поверхностным импедансом, поэтому нахождение корней ДУ является актуальной задачей.

Одним из распространенных методов поиска корней общего нелинейного уравнения

$f\left(x\right)=0$ (1)

является метод простых итерации (МПИ) [6; 7] или последовательных приближений. Согласно ему осуществляются итерации

 $x_n=\phi \left(x_{n-1}\right), n=\mathrm{0,1,2,...}$ (2)

Невязка уравнения (1) совпадает с самой функцией: $\Delta \left(x\right)=f\left(x\right),$ а для уравнения (2) невязка есть $\Delta \left(x\right)=x-\phi \left(x\right).$ Под невязкой, характеризующей решение, будем понимать квадрат модуля невязки $\delta ={\left|\Delta \left(x\right)\right|}^2={||\Delta \left(x\right)||}^2,$ т. е. квадрат нормы. Функция $\phi \left(x\right)$ может быть получена из (1) бесконечным числом способов: $\phi \left(x\right)=x+g\left(x\right)f\left(x\right),$ где  ‒ $g\left(x\right)$произвольная функция. Обычно для действительного уравнения ее берут непрерывной и знакопостоянной [6]. Если в уравнении (1) функция дифференцируемая p раз, функцию $g\left(x\right)$ также удобно брать таковой. Далее будем предполагать, что все функции дважды, а в некоторых алгоритмах и трижды непрерывно дифференцируемые. МПИ удобен и для комплексного уравнения для определения комплексных корней. Итерационный процесс для (1) и (2) запишем в эквивалентном виде

$x_n=x_{n-1}-{\tau }_n{\Delta }_{n-1},$ (3)

где для уравнения (1) ${\Delta }_n=f\left(x_n\right),$ а для уравнения (2) ${\Delta }_n=x_n-\phi \left(x_n\right).$ Очевидно, классический МПИ (2) реализуется, если параметр итерации равен единице на каждом шаге: ${\tau }_n=1.$ Однако выбор параметра итерации на каждом шаге дает дополнительную возможность улучшить сходимость. Это, например, используется в методе минимальных невязок (ММН) [8; 9]: параметр итерации $\tau {}_n$ выбирается из условия минимума ${\delta }_n$ на следующем шаге. В случае систем линейных алгебраических уравнений (СЛАУ) или для линейных операторных уравнений $\widehat{A}x=b$ или $x=\left(\widehat{I}-\widehat{A}\right)x+b$ в форме (3) $x_n=x_{n-1}-{\tau }_n{\Delta }_{n-1},$ ${\Delta }_{n-1}=\widehat{A}{x}_{n-1}-b$ имеем [8; 9]:

${\tau }_n=\left(\widehat{A}{\Delta }_{n-1},{\Delta }_{n-1}\right)/{||\widehat{A}{\Delta }_{n-1}||}^2.$ (4)

Скалярное произведение понимаем как $\left(x,y\right)=x^{*}\cdot y.$ Для сходимости простых итераций $({\tau }_n\equiv 1)$ достаточно, чтобы для одной из норм выполнялось условие $||\widehat{I}-\widehat{A}||<1.$ Для сходимости ММН достаточно выполнения условия $||\widehat{I}-{\tau }_n\widehat{A}||<1.$ Считаем матрицу эрмитовой, что всегда можно сделать, умножив СЛАУ на матрицу, эрмитово сопряженную исходной и нормированной: $||\widehat{A}||=\mathrm{1,}$ разделив СЛАУ на норму. Рассмотрим квадрат спектральной нормы ${\mu }_i.$ Это наибольшее собственное значение матрицы и по условию ${\mu }_i=1.$ Тогда итерации выполняются с параметром (4), если $0<{\tau }_n<1.$ При выходе из этой области параметр следует корректировать. При ${\tau }_n>1$ следует брать ${\tau }_n=1-\delta ,$ $0<\delta <1.$ При очень малом параметре алгоритм почти зацикливается, поэтому следует брать ${\tau }_n=\delta ,$ не допуская слишком малых значений. Алгоритм метода скорейшего спуска (МСС) аналогичен, но в нем минимизируется энергетическая норма невязки, а параметр итерации приобретает известный вид ${\tau }_n={||{}_{n-1}||}^2/$$\left(\widehat{A}{\Delta }_{n-1},{\Delta }_{n-1}\right).$ ММН и МСС могут не сойтись. Даже в случае одного линейного уравнения итерационный процесс может разойтись, а может сходиться достаточно медленно. Для одного линейного уравнения ММН и МСС дают решение за один шаг в силу скалярной природы задачи: ${\tau }_1=A^{-1}.$ Это также следует из условия ${\left|{\Delta }_1\right|}^2=0.$ В многомерном случае ${||{\Delta }_n||}^2\ne \mathrm{0,}$ и из наложения условия ${||{\Delta }_n||}^2=0$ не следует, вообще говоря, сходящийся алгоритм. Получающееся квадратное уравнение для ${\tau }_n$ может иметь комплексные значения для действительной СЛАУ. Использование матричного параметра итерации ${\widehat{\tau }}_n$ сопряжено с трудностями его определения: такой оптимальный параметр должен быть взят в виде ${\widehat{\tau }}_n={\widehat{A}}^{-1},$ т. е. приводит к той же самой задаче, для решения которой он и выбирается. Матричный параметр целесообразен в специальных случаях, например, когда матрица $\widehat{A}$ имеет большие диагональные элементы. Однако для поиска корня нелинейного уравнения условие ${||{\Delta }_n||}^2=\mathrm{0,}$ весьма продуктивно.

1. Методы коррекции параметра итерации

При решении нелинейного скалярного уравнения соотношение ${\left|{\Delta }_n\right|}^2=0$ может быть наложено, поскольку принципиально может быть достигнуто. Цель работы – модификация МПИ с целью получения сходящихся алгоритмов для решения нелинейных уравнений. В работе [9] рассмотрен ММН в случае расположения начальной точки x близко к корню $\overline{x}.$ Он основан на разложениях Тейлора типа $f\left(x\right)\approx f^{\text{'}}\left(\overline{x}\right)\left(x-\overline{x}\right)+$$f^{\text{'}\text{'}}\left(\overline{x}\right){\left(x-\overline{x}\right)}^2/\mathrm{2,}$ в частном случае – на линеаризации, т. е. на удержании в разложениях такого типа линейного члена. Рассматривая уравнение (2) в форме (3), имеем ${\delta }_n={\left|x_n-x_{n-1}+{\tau }_n{\Delta }_{n-1}\right|}^2.$ Параметр ${\tau }_n$ можно искать из условия, что ${\delta }_n$ имеет минимум, т. е. из условия $\partial {\delta }_n/\partial {\tau }_n^{*}=0$ [8]. Поскольку ${\Delta }_{n-1}^{*}\ne 0$ (иначе итерации уже бы сошлись на шаге $n-1),$ это эквивалентно условию ${\Delta }_n=0.$ Из него имеем ${\tau }_n=\left(x_{n-1}-x_n\right)/{\Delta }_{n-1}.$ Поскольку ${\Delta }_n=\mathrm{0,}$ это означает, что $x_n=\overline{x}=\phi \left(\overline{x}\right).$ С другой стороны, $\phi \left(x_{n-1}\right)=$$\phi \left(\overline{x}+x_{n-1}-\overline{x}\right)\approx$$\overline{x}+{\phi }^{\text{'}}\left(\overline{x}\right)\left(x_{n-1}-\overline{x}\right),$ поэтому имеем соотношение

${\Delta }_n=\overline{x}-x_{n-1}+{\tau }_n\left(x_{n-1}-\overline{x}-{\phi }^{\text{'}}\left(\overline{x}\right)\left(x_{n-1}-\overline{x}\right)\right),$

из которого следует ${\tau }_n={\left[1-{\phi }^{\text{'}}\left(\overline{x}\right)\right]}^{-1}.$ Итерационный процесс приобретает вид

$\begin{array}{l}{x}_n=x_{n-1}\left[1-{\left(1-{\phi }^{\text{'}}\left(x_n\right)\right)}^{-1}\right]+\\ +{\left(1-{\phi }^{\text{'}}\left(x_n\right)\right)}^{-1}\phi \left(x_{n-1}\right).\end{array}$

Поскольку корень  неизвестен, процесс реализуем заменой $x_n\to x_{n-1}$ в форме

$x_n=\tilde{\phi }\left(x_{n-1}\right)=\left(1-{\tau }_n\right)x_{n-1}+{\tau }_n\phi \left(x_{n-1}\right),$ (5)

полагая ${\tau }_n={\left[1-{\phi }^{\text{'}}\left(x_{n-1}\right)\right]}^{-1}$ Итерации (5) имеют форму (2) с переопределением функции $\phi$ на $\tilde{\phi }\left(x\right).$ Рассмотрим, к чему он приводит. Известно, что процесс (2) сходится, если $\left|{\phi }^{\text{'}}\left(x\right)\right|\le q<1$ в области сходимости (принцип сжимающих отображений). Однако если это условие не выполнено, процесс (2) может разойтись. Это требует поиска и подбора вида функций $\phi \left(x\right).$ Удобно строить такие процессы, которые сходятся из любого начального приближения. Конечно, теоретическая сходимость не означает возможности вычисления. Скажем, вполне можно построить процесс, сходящийся за несколько итераций при движении из начальной точки вблизи корня, и сходящийся за 106 или за существенно большее число итераций при движении из удаленной точки. Такие точки нецелесообразно рассматривать и всегда желательно определять области, из которых скорость сходимости высокая. Если рассмотреть форму (5), то для выполнения условия сжимающих отображений достаточно выполнить $\left|1-{\tau }_n\left(1-{\phi }^{\text{'}}\left(x_{n-1}\right)\right)\right|<1/\alpha ,$ где параметр $\alpha >1.$ Выбирая этот параметр на каждом шаге, можно находить параметр итерации и строить сходящийся итерационный процесс. Удобнее, однако, сразу взять $\alpha \to \infty ,$ т. е. положить ${\tau }_n={\left[1-{\phi }^{\text{'}}\left(x_{n-1}\right)\right]}^{-1}.$ В этом случае производная правой части (5) в точке $x_{n-1}$ точно равна нулю, а процесс сходится, даже если ${\phi }^{\text{'}}\left(x_{n-1}\right)>1.$ Такие образом, на каждом шаге процесс (5) подбирает оптимальную функцию $\tilde{\phi }\left(x\right)$ в правой части так, что $\widetilde{{\phi }^{\text{'}}}\left(x_{n-1}\right)=0.$

Рассмотрим ряд частных случаев. Если ${\phi }^{\text{'}}\left(x_{n-1}\right)=\mathrm{0,}$ то ${\tau }_n=\mathrm{1,}$ в окрестности этой точки реализуется хорошо сходящийся процесс (2). Если $\left|{\phi }^{\text{'}}\left(x_{n-1}\right)\right|\gg 1$ (точка вблизи полюса производной), то ${\tau }_n\approx \mathrm{0,}$ поэтому $x_n\approx x_{n-1}.$ Это медленно сходящийся процесс, и для его убыстрения целесообразно сделать одну-две простые итерации (2), а затем перейти к процессу (5) в новой точке, т. е. отойти от полюса производной. Если ${\phi }^{\text{'}}\left(x_{n-1}\right)=x_{n-1},$ то $x_n\left(1-x_{n-1}\right)=-x_{n-1}{x}_{n-1}+\phi \left(x_{n-1}\right).$ Вблизи корня функция отличается от аргумента на малую величину: $\phi \left(x_{n-1}\right)=x_{n-1}+\delta .$ Тогда $x_n=x_{n-1}+\delta /\left(1-x_{n-1}\right).$ Если $\left|x_{n-1}\right|$ существенно меньше единицы или величина $x_{n-1}$ отрицательная, то сходимость хорошая. Если ${\phi }^{\text{'}}\left(x_{n-1}\right)\approx$1 возможна разболтка. Если $\left|{\phi }^{\text{'}}\left(x_{n-1}\right)\right|<1$ то для ее преодоления надо сделать 2‒3 итерации (2), а если $\left|{\phi }^{\text{'}}\left(x_{n-1}\right)\right|>1$ то 2‒3 итерации (5) без пересчета параметра итерации, а затем продолжить процесс (5). Если $\left|1-{\phi }^{\text{'}}\left(x_{n-1}\right)\right|\le \delta ,$ где величина $\delta$ порядка точности представления чисел в ЭВМ (погрешность округления), это означает, что точка вблизи корня, в котором ${\phi }^{\text{'}}\left(\overline{x}\right)=1$ При достижении такой точки целесообразно фиксировать ${\tau }_n$ и выполнять метод простой итерации с ${\tau }_n\approx 1$ Вообще, если ${\phi }^{\text{'}}\left(x_{n-1}\right)=1$ это означает, что $x_{n-1}=\phi \left(x_{n-1}\right),$ т. е. $x_{n-1}$ ‒ это корень. Поэтому рассматриваем случай ${\phi }^{\text{'}}\left(x_{n-1}\right)\approx 1$ Это может быть, только если ${\phi }^{\text{'}}\left(\overline{x}\right)=1$ и процесс подошел к корню. Фактически продолжение итераций (5) означает использование медленно сходящегося процесса (2). В этом случае следует сделать несколько итераций (5) без пересчета параметров, что обычно достаточно для получения результата. При приближенном вычислении производной из-за разболтки сходимость может не достигаться. В этом случае удобно вообще зафиксировать параметр итерации, пока производная слабо отличается от единицы. Если ${\phi }^{\text{'}}\left(\overline{x}\right)=1$ для достижения быстрой сходимости функцию целесообразно заменить на новую $\tilde{\phi }\left(x\right)=x+g\left(x\right)\left(x-\phi \left(x\right)\right).$ На эту функцию наложим условие $\widetilde{{\phi }^{\text{'}}}\left(x_{n-1}\right)=0$ что делает процесс быстро сходящимся. Функцию $g\left(x\right)$ можно взять с точностью до любого множителя, поэтому $g\left(x_{n-1}\right)$ задаем произвольно. Полагаем $\widetilde{{\phi }^{\text{'}}}\left(x\right)=1+g^{\text{'}}\left(x\right)\left(x-\phi \left(x\right)\right)+g\left(x\right)\left(1-{\phi }^{\text{'}}\left(x\right)\right)=a\left(x-x_{n-1}\right).$

Справа вообще можно написать полином, но это усложняет задачу. Имеем $1+g^{\text{'}}\left(x_{n-1}\right)$ $\left(x_{n-1}-\phi \left(x_{n-1}\right)\right)=\mathrm{0,}$ $\tilde{\phi }\left(x\right)=a{\left(x-x_{n-1}\right)}^2/$$2+\tilde{\phi }\left(x_{n-1}\right).$

Параметр $a$ возьмем из условия, что сходимость в точке начального приближения хорошая, например, берем $\widetilde{{\phi }^{\text{'}}}\left(x_0\right)=1/2.$ Таким образом, необходимо интегрировать уравнение

$g^{\text{'}}\left(x\right)=\frac{a\left(x-x_{n-1}\right)-1-g\left(x\right)\left(1-{\phi }^{\text{'}}\left(x\right)\right)}{\left(x-\phi \left(x\right)\right)}$

из начальной точки $x_{n-1}$ с произвольным начальным условием

$g\left(x_{n-1}\right)=\frac{\tilde{\phi }\left(x_{n-1}\right)-x_{n-1}}{x_{n-1}-\phi \left(x_{n-1}\right)}.$

Интегрирование желательно выполнить аналитически, например методом рядов. Если это возможно, к новой функции применимы оба процесса (2) и (5). Начальное условие для $g\left(x\right)$ произвольное и может использоваться для упрощения решения, как и параметр a.

Поскольку при итерации (5) $x_n$ не попадает точно в корень $\overline{x},$ итерации продолжаются, однако погрешности $\left|x_n-x_{n-1}\right|$ не возрастают, а обычно уменьшаются. Численное вычисление производных может приводить к разболтке, т. е. к немонотонному изменению $\left|x_n-x_{n-1}\right|.$ Вблизи корня процесс имеет второй порядок и сходится квадратично, т. е. подобно методу Ньютона [6]. Однако метод Ньютона сходится не всегда. Преимущества процесса (5) в том, что по построению $\widetilde{{\phi }^{\text{'}}}\left(x_{n-1}\right)=\mathrm{0,}$ $n=\mathrm{0,1,2,...,}$  поэтому он обладает бесконечной областью гарантированной сходимости. Однако при его реализации и численном вычислении производных возможно зацикливание. Его определяем условиями, когда невязка $\left|x_n-x_{n-1}\right|$ мала или изменяется немонотонно, а невязка уравнения $\left|{\Delta }_n\right|$ еще не мала для выхода из процесса. Тогда при численном вычислении производных для избегания зацикливания нужно после итерации (5) провести 2‒3 простые итерации (2), если для них $\left|{\phi }^{\text{'}}\right|<\mathrm{1,}$ или по формуле (5) без пересчета параметра итерации, если $\left|{\phi }^{\text{'}}\right|>1.$ С использованием полученного параметра ${\tau }_n$ и функции $f\left(x\right)=x-\phi \left(x\right)$ метод Ньютона в форме $x_n=x_{n-1}-\left(x_{n-1}-\phi \left(x_{n-1}\right)\right)/$$\left(1-{\phi }^{\text{'}}\left(x_{n-1}\right)\right)$ приобретает вид (5), т. е. сходится, если $\left|{\phi }^{\text{'}\text{'}}\left(x_{n-1}\right)\right|<{\left(1-{\phi }^{\text{'}}\left(x_{n-1}\right)\right)}^2/\left|{\Delta }_{n-1}\right|.$ Если вторая производная в корне ограничена, то сходимость имеет место всегда при попадании в малую окрестность корня. Быстрая сходимость в окрестности корня имеет место, если ${\phi }^{\text{'}\text{'}}\left(\overline{x}\right)=0.$

 

Таблица. Сходимость (ход итераций) процессов для вычисления $\mathit{x}\mathbf{=}\sqrt{\mathbf{2}}$

Table. Сходимость (ход итераций) процессов для вычисления $\mathit{x}\mathbf{=}\sqrt{\mathbf{2}}$

n	$x=x/2+1/x$	$x=4{\left(x+2/x\right)}^{-1}$ и алгоритм (5)	Алгоритм (10)	$x=x/2+1/x$	$x=4{\left(x+2/x\right)}^{-1}$	Алгоритм (10)
0	10,0	10,0	10,0	0,1	0,1	0,1
1	5,100000	0,392157	6,733333	10,050000	0,199005	6,733333
2	2,746078	0,728311	4,587899	5,124503	0,390282	4,587899
3	1,737195	1,151281	3,203909	2,757392	0,725323	3,203909
4	1,444238	1,3841813	2,344018	1,741358	1,148529	2,344018
5	1,414526	1,413902	1,847091	1,444943	1,384137	1,847091
6	1,414214	1,414214	1,592322	1,414540	1,413887	1,592322

 

Рассмотрим в качестве примера классическую задачу извлечения квадратного корня $x=\sqrt{a}.$ Как известно, алгоритм $x=a/x$ всегда расходится из любой точки, а алгоритм $x=\left(x+a/x\right)/2$ сходится всегда из любого начального приближения [6]. Причина в том, что в первом случае ${\phi }^{\text{'}}\left(x\right)=-a/x^2,$ $\left|{\phi }^{\text{'}}\left(x\right)\right|\gg 1$ при малых $x^2,$ и в точке корня ${\phi }^{\text{'}}\left(\overline{x}\right)=-\mathrm{1,}$ а во втором случае ${\phi }^{\text{'}}\left(x\right)=\left(1-a/x^2\right)/\mathrm{2,}$ и в точке корня ${\phi }^{\text{'}}\left(\overline{x}\right)=0.$ При этом всегда имеет место попадание в малую окрестность корня. Если в методе Ньютона взять $f\left(x\right)=x^2-a,$ получим сходящийся алгоритм $x=\left(x+a/x\right)/2$ с ${\tau }_n=1/2.$ Если в методе Ньютона взять $f\left(x\right)=\sqrt{x^2-a},$ получим расходящийся алгоритм $x=a/x.$ Вообще для функции $f\left(x\right)={\left(x^2-a\right)}^{\alpha }$ сходимость метода Ньютона будет, если $\left|1-{\alpha }^{-1}\right|<1.$ Следуя вышеизложенному, формируем процесс $x=\tilde{\phi }\left(x\right)=\left(1-\tau \right)x+a\tau /x,$ в котором $\tau ={\left[1+a/x^2\right]}^{-1}.$ Как видно, $\tilde{\phi }\left(x\right)=4{\left(x+a/x\right)}^{-1}.$ Для этой функции $\widetilde{{\phi }^{\text{'}}}\left(\sqrt{a}\right)=0.$ Желательно использовать алгоритм (5), для которого $\widetilde{{\phi }^{\text{'}}}\left(x_n\right)=\mathrm{0,}$ n = 1, 2, … Он эквивалентен алгоритму (2) с приведенной функцией $\tilde{\phi }\left(x\right).$ Однако в общем случае, когда производная неизвестна и вычисляется численно, получить функцию $\tilde{\phi }\left(x\right)$ нельзя. В случае извлечения корня такая формула приобретает вид $x=\left(2x+a/x\right)/3.$ Для нее, как и для процессов $x=x/2+1/x,$ $x=4/\left(x+2/x\right)$ и (5), в таблице приведены сходимости итераций для двух начальных приближений $x_0=10$ и $x_0=\mathrm{0,1}$ при $a=2.$ Еще один процесс с удержанием производной второго порядка следует из равенства $x_n=\phi \left(x_{n-1}-{\tau }_n{\Delta }_{n-1}\right)$ при условии, что величина ${\xi }_n={\tau }_n{\Delta }_{n-1}$ мала по модулю. Получаем

$\begin{array}{l}\phi \left(x_{n-1}-{\xi }_n\right)\approx \phi \left(x_{n-1}\right)-\\ -{\phi }^{\text{'}}\left(x_{n-1}\right){\xi }_n+{\phi }^{\text{'}\text{'}}\left(x_{n-1}\right){\xi }_n^2/\mathrm{2,}\end{array}$

что дает

$\begin{array}{l}{\xi }_n=\xi \left(x_{n-1}\right)\approx -\left(1-{\phi }^{\text{'}}\left(x_{n-1}\right)\right)/{\phi }^{\text{'}\text{'}}\left(x_{n-1}\right)+\\ +\sqrt{{\left[\left(1-{\phi }^{\text{'}}\left(x_{n-1}\right)\right)/{\phi }^{\text{'}\text{'}}\left(x_{n-1}\right)\right]}^2+2{\Delta }_{n-1}/{\phi }^{\text{'}\text{'}}\left(x_{n-1}\right)}.\end{array}$

2. Решение дисперсионных уравнений для поверхностных плазмонов

Поверхностные ПП в настоящее время широко исследуются в нанооптике, фотонике и наноплазмонике [1‒5; 10‒19]. Особенность их моделирования в том, что необходим поиск комплексных корней нелинейных трансцендентных дисперсионных уравнений. Рассмотрим решение комплексного уравнения

$n=\sqrt{1-z^2\left(n\right)}.$ (6)

Это уравнение описывает комплексный коэффициент замедления ПП вдоль поверхности образца с комплексным импедансом $z\left(n\right).$ Поверхностный импеданс характеризует входное сопротивление образца. Металл описывается комплексной зависящей от круговой частоты  диэлектрической проницаемостью (ДП) $\epsilon \left(\omega \right)={\epsilon }_L-{\omega }_p^2/\left({\omega }^2-i\omega {\omega }_c\right),$ где ${\omega }_p,$ ${\omega }_c$ ‒ константы, а ${\epsilon }_L$ на низких частотах также является константой. На высоких частотах мы опишем ее двумя термами дисперсии Лоренца, подбирая параметры так, чтобы удовлетворить экспериментальным данным для металлов в оптическом диапазоне. Далее обозначаем $\epsilon ={\epsilon }^{\text{'}}-i{\epsilon }^{\text{'}\text{'}}.$ Уравнение (6) имеет точное решение $n=\sqrt{\epsilon /\left(\epsilon +1\right)}$ для любого однородного образца с ДП e, занимающего полупространство (J. Zenneck, 1907). Это точное решение получается сшиванием электромагнитных полей или приравниванием входных импедансов, например металла $z\left(n\right)=\sqrt{\epsilon -n^2}/\epsilon$ и характеристического импеданса вакуума $z_0\left(n\right)=\sqrt{1-n^2}$  для E-волны. Однако если образец неоднородный (слоистый или ДП $\epsilon$ зависит от нормальной координаты), входной импеданс $z\left(n\right)$ является сложной функцией, но уравнение (6) все равно справедливо. Решения таких уравнений рассмотрены в ряде работ (см. [1‒4; 14]). Кроме того, в ряде работ показано совместное итерационное решение интегральных уравнений и соответствующих им функционалов, получающиеся из линейной теории рассеяния или линейной электродинамики. Искомый параметр в такие функционалы входит нелинейно, в том числе и через аргументы трансцендентных функций. В нелинейной теории рассеяния итерационные методы, как правило, являются наиболее эффективными для получения решений. Часто в теории рассеяния возникают интегральные уравнения типа Липпмана – Швингера в форме объемных (трехмерных) уравнений Фредгольма второго рода, для которых удобно организовать итерационный процесс типа (2), в котором $\widehat{\phi }\left(x\right)$ ‒ интегральный или интегродифференциальный оператор, а x ‒ функция трехмерных координат и времени. Для таких задач также применимы изложенные итерационные алгоритмы.

Для задач о плазмонах есть области, где модуль производной правой части уравнений типа (6) больше единицы и метод прямой итерации расходится. Решение $n=\sqrt{\epsilon /\left(\epsilon +1\right)}$ для металла (серебро) приведено на рис. 1, кривая 1. Однако если решать уравнение (6) МПИ, то в области ${\epsilon }^{\text{'}}<-1$ он сходится, в области $-1<{\epsilon }^{\text{'}}<1$ расходится, а в области $1<{\epsilon }^{\text{'}}$ снова имеет место сходимость. Действительно, $\phi \left(n\right)=\sqrt{n^2+{\epsilon }^2-\epsilon }/\epsilon ,$ ${\phi }^{\text{'}}\left(n\right)=n/\left(\epsilon \sqrt{n^2+{\epsilon }^2-\epsilon }\right).$ Будем считать потери очень малыми, т. е. ${\epsilon }^{\text{'}\text{'}}\ll \left|{\epsilon }^{\text{'}}\right|.$ В области малых частот $n\approx \mathrm{1,} \left|{\epsilon }^2\right|\gg \mathrm{1,}$ и $\left|{\phi }^{\text{'}}\left(n\right)\right|\approx 0.$ В области плазмонного резонанса ${\epsilon }^{\text{'}}\approx -1$ имеем $n\approx \left(1-i\right)/\sqrt{2{\epsilon }^{\text{'}\text{'}}},$ $\left|n\right|\gg 1$ и $\left|{\phi }^{\text{'}}\left(n\right)\right|\approx 1.$ Если $\left|\epsilon \right|\approx \mathrm{0,}$ то $\left|n\right|\approx \mathrm{0,}$ а $\left|{\phi }^{\text{'}}\left(n\right)\right|\gg 1.$ Если $\epsilon =\mathrm{1,}$ то $n=1$ и ${\phi }^{\text{'}}\left(1\right)=1.$ Нетрудно видеть, что ${\phi }^{\text{'}}\left(n\right)=1/{\epsilon }^2={\left(1-n^2\right)}^2/n^4.$ Поэтому сходимость возникает при $\epsilon >\mathrm{1,}$ когда $1/2<n<1$ и ${\phi }^{\text{'}}\left(n\right)=1/{\epsilon }^2<1.$ Построим процесс, который сходится во всей области комплексных значений ДП. Если взять ${\phi }^{\text{'}}\left(n\right)=1/{\epsilon }^2,$ процесс типа (2) определен функцией

$\tilde{\phi }\left(n\right)=\left[{\epsilon }^2\sqrt{n^2/{\epsilon }^2+1-1/\epsilon }-n\right]/\left({\epsilon }^2-1\right).$ (7)

Однако такой процесс хорошо сходится только в окрестности корня: именно там $\widetilde{{\phi }^{\text{'}}}\left(n\right)=0.$ Процесс с индексом k, сходящийся из любого начального приближения, имеет вид

$n_k=n_{k-1}\left(1-{\tau }_k\right)+{\tau }_k\sqrt{n_{k1}^2/{\epsilon }^2+1-1/\epsilon }.$ (8)

Для него $\partial n_k/\partial n_{k-1}=\mathrm{0,}$ $\partial {\tau }_k/\partial n_{k-1}=\mathrm{0,}$ где

${\tau }_k=\frac{{\epsilon }^2\sqrt{n_k^2/{\epsilon }^2+1-1/\epsilon }}{{\epsilon }^2\sqrt{n_k^2/{\epsilon }^2+1-1/\epsilon }-1}.$

Реально необходимо строить частотную дисперсию ПП. Полагая для низких частот $n_0=\mathrm{1,}$ получаем быстро сходящийся к решению $n\approx 1$ процесс. Это решение на низких частотах определяет быстрый (движущийся почти со скоростью света) ПП. Плавно повышая частоту и используя для начального значения полученное решение на предыдущей частоте, получаем быстро сходящиеся процессы для всех частот. Преимущество таких вычислений в том, что всегда начальное приближение близко к корню, поэтому скорость сходимости высокая, и можно также использовать процесс (7). Результаты расчета приводят к точкам 2 на рис. 1, точно ложащимся на дисперсионную ветвь (кривая 1, рис. 1).

 

Рис. 1. Дисперсия замедления ПП на поверхности серебра $\mathbf{(}{\mathit{\omega }}_{\mathbf{p}}\mathbf{=}\mathbf{1}\mathbf{,}\mathbf{35}\mathbf{\cdot }{\mathbf{10}}^{\mathbf{16}}\mathbf{,}$ ${\mathit{\omega }}_{\mathbf{c}}\mathbf{=}\mathbf{2}\mathbf{,}\mathbf{7}\mathbf{\cdot }{\mathbf{10}}^{\mathbf{13}}\mathbf{,}$ ${\mathit{\epsilon }}_{\mathbf{L}}\mathbf{=}\mathbf{9}\mathbf{,}\mathbf{0}\mathbf{)}$ (кривая 1), на поверхности серебра, покрытой диэлектрической пленкой толщины 250 нм с ДП $\mathbf{4}\mathbf{-}\mathit{i}\mathbf{0}\mathbf{,}\mathbf{004}$ (3), на поверхности серебра, покрытой диэлектрической пленкой и пленкой серебра толщиной 10 нм (4). Точки 2 показывают сходимость итераций (5) для ПП на серебре

Fig. 1. Dispersion of PP deceleration on the silver surface $\mathbf{(}{\mathit{\omega }}_{\mathbf{p}}\mathbf{=}\mathbf{1}\mathbf{,}\mathbf{35}\mathbf{\cdot }{\mathbf{10}}^{\mathbf{16}}\mathbf{,}$ ${\mathit{\omega }}_{\mathbf{c}}\mathbf{=}\mathbf{2}\mathbf{,}\mathbf{7}\mathbf{\cdot }{\mathbf{10}}^{\mathbf{13}}\mathbf{,}$ ${\mathit{\epsilon }}_{\mathbf{L}}\mathbf{=}\mathbf{9}\mathbf{,}\mathbf{0}\mathbf{)}$ (curve 1), on the silver surface covered with a dielectric film 250 nm thick with DP $\mathbf{4}\mathbf{-}\mathit{i}\mathbf{0}\mathbf{,}\mathbf{004}$ (3), on the surface of silver coated with a dielectric film and a silver film 10 nm thick (4). Points 2 show the convergence of iterations (5) for PP on silver

 

Рассмотрим более сложную подобную задачу. Пусть на металле лежит слой диэлектрика с ДП ${\epsilon }_d$ толщины d. Тогда входной импеданс приобретает вид

$z\left(n\right)=z_d\frac{z_0\left(n\right)+i{z}_d\tan \left(k_{0}d\sqrt{{\epsilon }_d-n^2}\right)}{z_d+i{z}_0\left(n\right)\tan \left(k_{0}d\sqrt{{\epsilon }_d-n^2}\right)}.$ (9)

Опять решаем уравнение (6). Результат решения приведен на рис. 1, кривая 3. Еще один пример на рис. 1 (кривая 4) соответствует случаю, когда еще и на слое диэлектрика лежит слой металла толщины $d_m.$ Импеданс строится рекуррентно по формуле (9) подстановкой в нее $z_0\left(n\right)\to z\left(n\right)$ и $z_d\to z_m,$ ${\epsilon }_d\to \epsilon ,$ $d\to d_m.$ Можно получить метод более высокого порядка, оставляя в разложении правой части в точке корня три члена $\phi \left(x\right)=\phi \left(\overline{x}\right)+{\phi }^{\text{'}}\left(\overline{x}\right)\left(x-\overline{x}\right)+{\phi }^{\text{'}\text{'}}\left(\overline{x}\right){\left(x-\overline{x}\right)}^2/2$ и полагая $\partial {\delta }_n/\partial {\tau }_n^{*}=\mathrm{0,}$ что также ведет к условию ${\Delta }_n=0.$ Окончательно имеем

$x=\frac{x\left[2{\phi }^{\text{'}}\left(x\right){\phi }^{\text{'}\text{'}\text{'}}\left(x\right)-{{\phi }^{\text{'}\text{'}}}^2\left(x\right)\right]-\phi \left(x\right){{\phi }^{\text{'}\text{'}}}^2\left(x\right)}{2{\phi }^{\text{'}}\left(x\right){\phi }^{\text{'}\text{'}\text{'}}\left(x\right)}.$ (10)

 

Рис. 2. Дисперсия потерь ПП на поверхности серебра (1), на серебре, покрытом диэлектрической пленкой (3), на серебре, покрытом диэлектрической пленкой и пленкой серебра толщиной 10 нм (4). Точки 2 показывают сходимость итераций (5) для ПП на серебре. Структура соответствует рис. 1

Fig. 2. Dispersion of PP losses on the surface of silver (1), on silver coated with a dielectric film (3), on silver coated with a dielectric film and a silver film with a thickness of 10 nm (4). Points 2 show the convergence of iterations (5) for PP on silver. The structure corresponds to Fig. 1

 

Рис. 3. Дисперсия в слое асимметричного ГММ (вкладка) толщины 50 нм при разных углах a: p/3, p/4 и p/6. Знаки «+» и «‒» в y соответствуют графикам (а) и (б)

Fig. 3. Dispersion in an asymmetric HMM layer (inset) 50 nm thick at different angles a: p/3, p/4, and p/6. The signs «+» and «‒» in y correspond to graphs (a) and (b)

 

Рассмотрим пример более сложной структуры в виде слоя гиперболического метаматериала ГММ второго типа [5] (см. вставка на рис. 3, а) толщины d, образованного тонкими серебряными пленками, разделенными диэлектрическими слоями SiO2. Для упрощения гомогенизации мы считаем слои металла и диэлектрика одинаковой наноразмерной толщины. Тогда гомогенизация справедлива вплоть до оптических частот, и слой, имеющий ось анизотропии вдоль оси z, характеризуется диагональным тензором ДП ${\epsilon }_{xx}={\epsilon }_{yy}=\left({\epsilon }_m+{\epsilon }_d\right)/\mathrm{2,}$ ${\epsilon }_{zz}=2{\epsilon }_m{\epsilon }_d/\left({\epsilon }_m+{\epsilon }_d\right).$  Здесь мы обозначили ДП металла ${\epsilon }_m$ и диэлектрика ${\epsilon }_d.$ Для ДП металла (серебра) мы используем модель Друде – Лоренца с введением двух Лоренцевых термов для описания сложного поведения ${{\epsilon }^{\text{'}}}_m$ вблизи перехода через нуль [20]. Если повернуть ось анизотропии на угол α относительно оси y, получим тензор ДП вида

$\begin{array}{l}\widehat{\epsilon }=\left[\begin{array}{cc}{\cos }^2\left(\alpha \right){\epsilon }_{xx}+{\sin }^2\left(\alpha \right){\epsilon }_{zz}& 0\\ 0& {\epsilon }_{xx}\\ \sin \left(\alpha \right)\cos \left(\alpha \right)\left({\epsilon }_{zz}-{\epsilon }_{xx}\right)& 0\end{array}\right.\\ \left.\begin{array}{c}\sin \left(\alpha \right)\cos \left(\alpha \right)\left({\epsilon }_{zz}-{\epsilon }_{xx}\right)\\ 0\\ {\cos }^2\left(\alpha \right){\epsilon }_{zz}+{\sin }^2\left(\alpha \right){\epsilon }_{xx}\end{array}\right].\end{array}$ (11)

Он соответствует случаю асимметричного ГММ, когда ось анизотропии направлена под углом к границе раздела. Еще более сложным является случай, когда волна распространяется вдоль поверхности в произвольном направлении, т. е. образует некий угол q с осью z. Поворачивая систему координат вдоль оси x так, что новая ось z совпала с направлением распространения, получим в новой системе координат тензор ДП вида

$\tilde{\epsilon }=\left[\begin{array}{ccc}{\tilde{\epsilon }}_{xx}& {\tilde{\epsilon }}_{xy}& {\tilde{\epsilon }}_{xz}\\ {\tilde{\epsilon }}_{xy}& {\tilde{\epsilon }}_{yy}& {\tilde{\epsilon }}_{yz}\\ {\tilde{\epsilon }}_{xz}& {\tilde{\epsilon }}_{yz}& {\tilde{\epsilon }}_{zz}\end{array}\right],$ (12)

с компонентами:

${\tilde{\epsilon }}_{xx}={\cos }^2\left(\alpha \right){\epsilon }_{xx}+{\sin }^2\left(\alpha \right){\epsilon }_{zz},$

$\begin{array}{l}{\tilde{\epsilon }}_{yy}={\cos }^2\left(\theta \right){\epsilon }_{xx}+{\sin }^2\left(\theta \right)\times \\ \times \left({\cos }^2\left(\alpha \right){\epsilon }_{zz}+{\sin }^2\left(\alpha \right){\epsilon }_{xx}\right),\end{array}$

$\begin{array}{l}{\tilde{\epsilon }}_{zz}={\sin }^2\left(\alpha \right){\epsilon }_{xx}+{\cos }^2\left(\theta \right)\times \\ \times \left({\cos }^2\left(\alpha \right){\epsilon }_{zz}+{\sin }^2\left(\alpha \right){\epsilon }_{xx}\right),\end{array}$

${\tilde{\epsilon }}_{xy}=-\sin \left(\theta \right)\sin \left(2\alpha \right)\left({\epsilon }_{zz}-{\epsilon }_{xx}\right)/\mathrm{2,}$

${\tilde{\epsilon }}_{xz}=\cos \left(\theta \right)\sin \left(\alpha \right)\cos \left(\alpha \right)\left({\epsilon }_{zz}-{\epsilon }_{xx}\right),$

$\begin{array}{l}{\tilde{\epsilon }}_{yz}=\sin \left(2\theta \right)\times \\ \times \left[{\epsilon }_{xx}-\left({\cos }^2\left(\alpha \right){\epsilon }_{zz}+{\sin }^2\left(\alpha \right){\epsilon }_{xx}\right)\right]/2.\end{array}$

Мы не приводим довольно сложно получающиеся ДУ. Наиболее простой случай здесь соответствует полупространству с тензором ДП (11), когда волны распадаются на E и H ПП. Также имеется два типа волн для случая конечного слоя (рис. 3). В этом случае ДУ имеет вид

$k_z^{\pm }=\pm k_0{\left[\frac{{\tilde{\epsilon }}_{xx}f\left(k_0,k_z\right)/\Delta -1}{f\left(k_0,k_z\right)/\Delta -1}\right]}^{1/2},$ (13)

где $\Delta ={\widehat{\epsilon }}_{xx}{\widehat{\epsilon }}_{zz}-{\widehat{\epsilon }}_{xz}^2,$ $\psi =\pm d\sqrt{\left(k_0^2{\widehat{\epsilon }}_{xx}-k_z^2\right)\Delta }/{\widehat{\epsilon }}_{xx},$ а входящая в (13) функция представляется в виде

$\begin{array}{l}f\left(k_0,k_z\right)=\\ ={\left[\frac{\left(\rho -{\rho }_0\right)\exp \left(-i\psi \right)-\left(\rho +{\rho }_0\right)\exp \left(i\psi \right)}{\left(\rho -{\rho }_0\right)\exp \left(-i\psi \right)+\left(\rho +{\rho }_0\right)\exp \left(i\psi \right)}\right]}^2.\end{array}$

В случае тензора ДП (12) волны вдоль поверхности слоя гибридные. ДУ может быть получено методом сшивания и на основе матриц передачи (переноса)  – так называемых матриц Берремана. Для рис. 3 также существуют области прямых и обратных волн на дисперсионных кривых. В силу корневых многозначностей в общем случае возможны четыре дисперсионные ветви. Возможно рассмотрение ПП вдоль поверхностей ГММ, построенных на основе графен-диэлектрических структур [5]. ГММ с полупроводниковыми и графеновыми слоями хорошо оптически управляются. Также возможно создание гиротропных структур, управляемых магнитным полем, но получение ДУ в этом случае является весьма сложной задачей.

Заключение

Показано, что алгоритм МПИ с коррекцией сходимости путем подбора правой части, имеющей нулевую производную на каждом шаге, совпадает с ММН первого приближения (при линеаризации) и имеет то преимущество, что всегда сходится (по крайней мере, невязка не возрастает). Такой выбор ${\tau }_n$ приводит к равенству нулю линеаризованной невязки. Вблизи корня сходимость квадратичная. Однако в ряде случаев возможно сильное замедление сходимости. В этих случаях можно предусмотреть специальные методы коррекции шага. Если корней несколько, алгоритм их поиска может состоять из движений из сильно различных начальных приближений.

Рассмотренные адаптивные методы удобны для задач плазмоники и позволяют определять целиком всю комплексную ветвь и даже несколько ветвей (если они есть) дисперсионного уравнения при плавном изменении частоты (или иного параметра), т. е. решать проблему поиска и классификации комплексных корней. На рис. 2 приведены результаты для потерь $n^{\text{'}\text{'}}\left(\omega \right),$ соответствующие рис. 1. Приведенные численные результаты, рис. 1–3 показывают возможность существования обратных ПП, для которых $n^{\text{'}}\left(\omega \right)n^{\text{'}\text{'}}\left(\omega \right)<0.$ Если дисперсионная кривая 4, рис. 1 строится из условия $\mathrm{Re}\left(n\left(\omega \right)\right)>\mathrm{0,}$ то обратным ПП соответствуют «отрицательные» потери (рис. 2, кривая 4).

Об авторах

Михаил Владимирович Давидович

Саратовский национальный исследовательский государственный университет имени Н.Г. Чернышевского

Email: davidovichmv@info.sgu.ru

доктор физико-математических наук, профессор кафедры радиотехники и электродинамики 

Россия, Саратов

Александр Константинович Кобец

ООО «НПФ “Этна плюс”»

Email: kobetzak@info.sgu.ru

директор по качеству 

Россия, Саратов

Кирилл Александрович Саяпин

Саратовский национальный исследовательский государственный университет имени Н.Г. Чернышевского

Автор, ответственный за переписку.
Email: sayapin_kirill@mail.ru

аспирант кафедры радиотехники и электродинамики Института физики

Россия, Саратов

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

Дополнительные файлы