Решение методом Монте–Карло нелинейного уравнения Шредингера Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

АПВПМ-2019

РЕШЕНИЕ МЕТОДОМ МОНТЕ-КАРЛО НЕЛИНЕЙНОГО УРАВНЕНИЯ ШРЕДИНГЕРА

В, Л, Лукинов

Институт вычислительной математики и математической геофизики СО РАН, 630090, Новосибирск Новосибирский государственный университет, 630090, Новосибирск Сибирский государственный университет телекоммуникаций и информатики, 630102, Новосибирск

УДК 519.633

Б01: 10.24411/9999-016А-2019-10046

Проведено сравнение численных решений нелинейного уравнения Шредингера путем применения нелинейных аналогов прямого (РИРТ) и обратного преобразования (ЕШРТ) Фурье и прямой схемы Эйлера. При РИРТ происходит разложение спектра поступающего сигнала путём решения уравнения Захарова-Шабата. Распространение нелинейной части спектра описывается известным уравнением. Для нахождения получаемого сигнала при ВХКТ применяются разработанные автором методы Монте-Карло для численного решения интегрального уравнения Гельфанда-Левитана-Марченко. Для оценки влияния аддитивного гауссовского белого шума использовалось численное решение на основе схемы Эйлера. В обоих численных схемах моделировался случайный гауссовский белый шум источника. Проведено моделирование потери мощности сигнала при прохождении случайной гетерогенной среды и изучено влияние шумов, возникающих при максимальных нагрузках оптического волокна. Проведено моделирование взаимодействия солитонов, распространяющихся навстречу друг другу с разных концов оптоволоконной среды. Для моделирования использовалось сочетание численной схемы на основе РИРТ и ВИРТ до момента взаимодействия и неявной безусловно-устойчивой схемы Кранка-Николсона для решения дискретной дифференциальной задачи на равномерных сетках по пространству и времени. Ключевые слова: статистическое моделирование, нелинейное уравнений Шредингера.

Введение

Данная работа посвящена изучению влияния случайных шумов при распространении и взаимодействии оптических сигналов в оптоволоконных каналах. На пропускную способность траффика нелинейных каналов оказывают влияние шумы с различного рода источниками: обратное рэлеевское рассеяние, спонтанное ра-мановское рассеяние и в усилителях вследствие спонтанной эмиссии фотонов. Оптический шум совместно с дисперсией и нелинейностью это три ключевых физических эффекта, влияющих на распространение оптических сигналов в оптоволоконных каналах, которые точно описываются нелинейным уравнением Шредингера (]МЬ8Е), учитывающим непрерывное взаимодействие между дисперсией и нелинейностью [1]. Известно, что КЪБЕ (без возмущения) принадлежит классу интегрируемых нелинейных систем [2]. В частности, это обосновывает применение нелинейных аналогов преобразования Фурье ^И?) для перехода от пространственно-временных к спектрально-временным переменным. В настоящее время можно выделить два подхода: нелинейное частотное разделение каналов ^ГИМ) и нелинейный обратный синтез (N18), использующие для передачи сигналов дискретную (солитонную) и непрерывные части нелинейного спектра [3,4].

Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (код проекта 17-01-00698 А) и Института вычислительной математики и математической геофизики СО РАН (гос. задание 0315-2016-0002).

ISBN 978-5-901548-42-4

1 Основные нелинейные преобразования в оптической связи

Рис. 1: Схема прямого ОТТ (а) и обратного ЮТТ (Ь) для входного потенциала ц(Ь), стремящегося к 0 при

В передатчике передаваемая последовательность двоичных данных сначала кодируется цифровым способом в комплексную форму волны в(1), используя произвольный формат модуляции и метод кодирования. Линейный спектр Фурье этой кодированной комплексной формы сигнала определяется следующим образом

5 И

в^) ехр

(1)

После кодирования линейный спектр Фурье кодированной комплексной формы волны ^(ш) отображается на непрерывную часть нелинейного спектра комплексного сигнала ц(Ь) с использованием обратного нелинейного преобразования Фурье BNFT (то-есть путем решения уравнения СЬМЕ). Обозначим через г(£) непрерывную часть нелинейной спектра д^). Тогда операция отображения в блоке BNFT может быть выражена как:

•К£)15=-ш/2 = -ЭД

Комплексный сигнал, затем подается в Щ модулятор для прямого преобразования в оптическую среду и передается в оптоволокно.

В>№Т блок отображает закодированную информацию на непрерывную часть нелинейного спектра сложного сигнала. Комплексный передаваемый сигнал (выход ВХГТ блока) в этом случае не содержит солитонов, иначе говоря, нелинейный спектр не содержит дискретной части. На приеме реальные и мнимые части передаваемого сигнала обнаруживаются с помощью когерентного приемника. Нелинейный спектр принятого сигнала получается прямым нелинейным преобразованием Фурье ХГТ (т. е. путем решения 2БР).

При распространении по каналу оптического волокна без потерь (например, с распределенным Раманов-скнм усилением), эволюция г(£) тривиальна [1]:

г(Ь,0= г(0 ■

где Ь — расстояние передачи.

Восстановление линейного спектра Фурье исходного кодированного комплексного сигнала после декодера возможно путем применения одноступенчатого линейного дисперсионного удаления:

§(и)= -г(Ь,0

2 Ц2Ь

После комплексный сигнал я^) можно восстановить используя операцию ПТТ, наконец, сигнал можно подавать в стандартный декодер для обнаружения данных.

2 Нелинейное уравнение Шредингера (N1(8)

Распространение комплексной медленно изменяющейся огибающей оптического поля ц(г,Ь) вдоль оптоволокна без потерь можно описать нелинейным уравнением Шредингера

Щг - у Яы + 1я\я\2

0

(2)

где г — расстояние та которое распространяется сигнал, £ — время распространения сигнала. Для определенности будем рассматривать КЪБЕ фокусирующего типа, при котором постоянный коэффициент хроматической дисперсии < 0 в уравнении (2). Дисперсионные члены более высокого порядка здесь не рассматриваются. Мгновенный коэффициент нелинейности Керра 7 равен

7 = П2^0/сАет, (3)

где п2 — показатель преломления, Ает — эффективная облаеть моды, с — скорость света в вакууме, — угловая несущая частота. Здесь рассматривается случай, когда постоянный коэффициент хроматической дисперсии < 0 в уравнении (2), то-есть так называемого фокусирующиюго типа КЪБЕ. Дисперсионные члены более высокого порядка здесь не рассматриваются. Мгновенный коэффициент нелинейности Керра 7 равен

7 = п2Шо/сАет, (4)

где п2 — показатель преломления, Ает — эффективная облаеть моды, с — вакуумная скорость света иш0 — угловая несущая частота.

В дальнейшем КЪБЕ будет рассматриваться в нормализованном виде:

гЧг - Чгг + Я^2 = 0, (5)

который может быть получен с помощью следующей замены:

Тр--> ^ --> г л/Ч^ —> Ч (6)

где Т3 — свободный параметр (например, характерный временной масштаб входного сигнала) и связанная пространственная шкала равна Zs = Т2 Заметим, что все величины д, £ и г в нормированном уравнении теперь безразмерны.

2.1 Прямое нелинейное преобразование Фурье

Подобно прямому И?, целью FNFT является разложение сигнала в ЮТ спектральные данные. Это преобразование производится путем решения задачи Захарова-Шабата (2БР) [2]. Последнее соответствует задаче рассеяния для неэрмитовой (в случае аномальной дисперсии) системы уравнений Дирака для двух вспомогательных функций '^2 (£) с профилем входного импульса КЪБЕ д(0, £) = д^), являющимся эффективным потенциалом в уравнениях:

= <1(Ф2 - г&1, = -аОМ + ■ (7)

Здесь £ является (как правило, комплексным) собственным значением, £ = £ + ц(Ь) является комплексным сопряжением потенциала ц^), который считается распадающимся при Ь —> (при выполнении определенных условий, наложенных на скорость затухания, см. [2]).

Для определения непрерывной части нелинейного спектра (вещественных £ = £), фиксируются два линейно-независимых решения уравнения Поста (7):

Ф(*,€) = [Ф1,Ф2]Т, Ф^,0 = [4>2, -ФТ]Т, (8)

с "начальным" состоянием слева:

Ф(Ш = [е"*4, 0]Т (9)

Таким же образом, справа фиксируются два других решения Поста

Ф(и) = [^2]Т, Ф(г,0 = [ф2, -а]Т, (ю)

= М«г ]Т (И)

Эти два набора решений линейно зависимы и могут быть выражены через коэффициенты рассеяния Поста а(£) и Ь(^):

Ф(Ь,® = а(£)Ф (г,0+ ь(€)*$€) (и

Ф (¿,0 = -а(£)*(Ъ£) + Ь(£)Ф (1,0 1 ;

Коэффициент отражения, производящий непрерывную часть нелинейного спектра, определяется как:

т = (13)

Солитоны соответствуют комплексным собственным значениям (¡п, где а(^п) = 0. Прямое NFT отображает начальное поле д(0, ^ па множество данных рассеяния ^ = [(г(£), £ — вещественное); ((п, уп =

(Ъ((п)а£(Сп))-1)^) где индекс п пробегает все дискретные собственные значения ЕБР (дискретная недиспер-

( )

( )

играет роль частоты. Нелинейная спектральная (N8) функция, определяемая как:

N (ш) = — (01г=-о,/2 (14)

(0, )

в линейном пределе.

2.2 Обратное нелинейное преобразование Фурье

( )

пня уравнения Гельфанда-Левитана-Марченко (СЬМЕ) для неизвестной функции К^ ,х) [5—7,9]. СЬМЕ без солитонов можно записать в виде:

г г

х)

К(t,x)+F(t + x)+ J ds J dr К(t, r)F(s + y)F(r + s)=0. (15)

— TO —TO

Здесь F(x) является линейным FT для r(^):

+ TO

F(x) = ±- ( r(d)е—*х<% (16)

-TO

Решая уравнение (26) для К (t,x), обратное можно получить в виде

q(t) = 2 lim К(t,x) (17)

3 Численные методы

3.1 Вычисление непрерывного спектра

Непрерывный спектр (т.е. а(£), Ь(£) и г(£)) может быть вычислен путем непосредственного интегрирования ZSP (7), а затем оценивая пределы для соответствующих компонентов функции Поста как:

а(0= ШтоМ*, О

b(0 = Jirf^2(t, Ое-^ (18)

t ^ + то

Для решения уравнения ZSP можно использовать легко распараллеливающийся метод PC А (7). При этом потенциал q(t) отсекается гае диапазона (— Т0; Т0). Внутри рассматриваемой области q(t) берется равным константе qn = q(tn) та каждом элементарном подинтервале (tn — At/2; tn + At/2), где tn = —Т0 + nAt, At = Т0/М — это шаг по времени и 2М + 1 — общее число разбиений рассматриваемой области. Идея метода РСА основана на том, что уравнение (7) может быть решено точно внутри каждого элементарного

Ф^n + At/2, i) = Т(qn, n — At/2,0, (19)

где матрица перехода Т(qn, £) задается в виде

" ' — ^

Т(qn, О = exp

A — n

V — In К J

(20)

Задачу рассеяния можно решить путем "продолжения" решения итеративно, начиная от —То до правой границы Т0, используя набор матриц переноса Т(qn, £) заданных уравнением (20). Итоговый результат может быть выражен в виде

Ф(То — At/2,0 = Ц(0$(—То + At/2,0

2м С21)

П(0=ПТ (qn, О [¿i)

п=1

Начальное условие, определенное на правом конце временного интервала, может быть представлено в ввиде

Ф(— Т0 — At/2, Í) = (1,0)Te-í«(-T°-Aí/2), (22)

на левом конце интервала имеет место следующее равенство

а(П e—«(To-Aí/2)

Ф(Т0 — A/2 0= ^(To-At/2), (23)

Коэффициенты Иоста могут быть выражены в виде

К0 = П 2i(0e-*At, ( j

где Пи(0 и П21 (О — соответствующие компоненты матрицы П(0- ® общем случае, когда потенциал q(t) обрывается вне интервала (Тт^,Ттах) с произвольными границами, выражение (24) может быть представлено в виде

а(0 = Пii (Ое^—-T—) Ъ(0 = П 2i(0e

11 ^ = -2i^(Tmax +Tmín-At) (25)

Из (21) видно, что можно вычислить матрицы переноса Т(хг) независимо друг от друга. В результате алгоритм РСА может быть легко реализован параллельно, что позволяет существенно сократить время вычислений.

В качестве примера рассмотрим прямоугольный импульс

„^Aíe ¡Ъ.Щ

Непрерывный спектр этого прямоугольного импульса определяется следующим образом:

m = (i —ve+ И|2 cot (уё2^(Т2 — Т1)) j

Определим нормированную среднеквадратичную ошибку (RMSE) для РСА следующим образом:

RMSE =

1 N

I Г-хаа (£к) — Гпишеггс(£ к )|2

к=1

3.2 Вычисление ВМЕТ

Для удобства в численных расчетах (I [,.\1К можно переписать в виде двух связанных интегральных уравнений:

X

Ат(х, г)+ / Р(г + х)^2(х, у)Лу = 0

(26)

А2(х,г) - / р(г + х)А1(х,у)йу = р(х + г),х>г,

где Р^) это обратное линейное преобразование Фурье для г(£)- Для численного анализа далее используется замена

и(х, в) = Ат(х,х — в), у(х, т)=А2(х,т — х), (27)

так, что (Н.МК (26) можно переписать в виде

2х

и(х, в) + § Р(1 — з)у(х, т)йт = 0

8Т (28)

у(х, т) — § Р(Ь — з)и(х, я)с1з = Р(т)

0

Функции и(х, в) и у(х, в) определены внутри интервала 0 < т < 2х < 2То- ^^^Т отображение г(£) на временную область задается формулой

д(х) = 2у(х, 2х — 0) (29)

Следуя процедуре дискретизации, представленной в [10], разделим интервал 0 <т< 2Т0, где функция Р(т) известна, в отрезки длины к = 2Т0/К. Дискретные перемеиные тп, вк и хт определяются по формулам

вк = к(к — 1)/2, к =1, 2,...,т

тп = к(п — 1)/2, п = 1, 2,... ,т (30)

хт = тк/2, т = 1, 2, . . . , N

Определим сеточные функции

и^ =и(хп, Тт), =у(хп, тт), Рп = Р(пк) (31)

Используя прямоугольную квадратурную схему для приближения интегралов в уравнении (28), получаем следующую дискретную форму (Н.МК:

4т) + к Е Рп-Аг] = 0

п=

--к

(т) Ь. Ц! (т) 17>

Vп - к £ Рп-кик = Рп

(32)

к = 1 т

ч(т) = у™ (33)

Уравнения (32) могут быть записаны в матричной форме

ет( ) = ь(т) (34)

где Ь(т) формируется го нулевого вектора размерности т и вектора размерности т с компонентами Рп; Ст в (34) является квадратной матрицей размерностей 2т х 2т, которая имеет следующий вид:

ш ( Е(т) №(т) \

ст = (ш) Е(т) (35)

4 Результаты

Целью работы является изучение нелинейных эффектов при распространении и взаимодействии лазерных импульсов в длинном оптическом волокне, возникающих из-за модуляции неустойчивости. Распространение длинных импульсов в оптическом волокне описывается КЪБ [1]:

дА 1 дА Ш д2 . . „12 „

где А^, г) — амплитуда электромагнитного поля на несущей частоте ш0 = 2-пп0/Х0, амплитуда нормируется степенью Р = |А|2; с — скорость света в оптическом волокне с показателем преломления п0; = — 10кт-1 • пт-2 — дисперсия па длине волны 1550пт; 7 = 3\¥-1 • кт-1 — коэффициент нелинейности;

Impulse range transmitted over S km Impulse range transmitted over S km

-1-1-1-1-1-1-r-^ ^-1-1-1-1-1-1

-300 -200 -100 0 100 200 300 -300 -200 -100 0 100 200 300

Time [ns) Time (ns)

а = 0.1МВ/кт — коэффициент оптических потерь. Граничное условие на входе в оптическое волокно дается как сумма гауссовского импульса с шириной в половине максимума То = 100пв и амплитудой шума Ап (1):

¿2 , \ /4\и2

A(t, 0) = ^PÖexpf^-2Ы2^ ^J^^+Mt)-

Здесь Ро = £о/То это средняя мощность энергии импульса ео- Случайный шум имеет ненулевую спектральную плотность мощности Рп(ш) = 1Ап(ш)|2 (Рп(ш) = еп/А) в спектральной области для А = 1 пт относительно центральной несущей частоты £ — средняя энергия шума за период т = 2гш Отношение энергий импульса и шума ео/еп принималось равным 0.7 [3]

Для численного решения уравнения (36) использовалась описанная выше модифицированная N18 схема со временем Т = 1тв, при этом количество точек разбиения бралось равным N = 220. Пространственный шаг полагался равным А г = 0.1тт. Начальное распределение было задано гауссовским импульсом с длительность в половине максимальной интенсивности 100ш и переменной пиковой мощности Рр, а именно А(0,1) = \[Р~Р ехр(—12/2Т2), и к нему добавлен белый гауссовский шум со спектральной полосой 1пт. Количество смоделированных импульсов составляло 106.

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

fl] G. Agrawal, Nonlinear fiber optics. Academic Press. New York. 1996

[2] В. E. Захаров и А.Б. Шабат. Точная теория двумерной самофокусировки и одномерной самомодуляции волн в нелинейных средах // Сов. Физ. 34. 6269. 1972.

[3] S.T. Lee. J.E. Prilepsky. S.K. Tnritsyn. Nonlinear inverse synthesis for high spectral efficiency transmission in optical fibers // OPTICS EXPRESS. Vol. 22, No 22, 2014.

[4] Z. Dong, S. Hari, G. Tao, Z. Ivangping, M. I. Yousefi, L. Chao, et al., Nonlinear Frequency Division Multiplexed Transmissions Based on NFT // PTL, IEEE, vol. 27, pp. 1621 1623, 2015.

[5] M. J. Ablowitz, D. J. Каир, A. C. Newell, and H. Segur, The inverse scattering transform-Fourier analysis for nonlinear problems // Stud. Appl. Math. 53, 249315, 1974.

[6] V. E. Zakharov, S. V. Manakov, S. P. Novikov, and L. P. Pitaevskii, Theory of Solitons. The Inverse Scattering Method. Colsnltants Bureau, New York, 1984.

[7] М. J. Ablowitz and Н. Segur, Solitons and the Inverse Scattering Transform. SIAM, Philadelphia, 1981.

[8] Crank, J.; Nicolson, P., "A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type" // Proc. Camb. Phil. Soc. 43 (1): 50-67, 1947. doi:10.1017/S0305004100023197

[9] A. C. Newell, Solitons in mathematics and physics. SIAM, Philadelphia, 1985.

[10] О. V. Belai, L. L. Frumin, E. V. Podivilov, and D. A. Shapiro, "Efficient numerical method of the fiber Bragg grating synthesis", //J. Opt. Soc. Am. В 24(7), 1451-1457, 2007.

Лукинов Виталий Леонидович — к.ф.-м.н., ст. науч. сотр. Института вычислительной математики и математической геофизики СО РАН;

e-mail: vitaliy.lukinov@sscc.ru. Дата поступления — 30 апреля 2019 г.