CC BY
16
УДК 111.111
А. В. Гутор, Е. С. Мануйлович
Московский физико-технический институт (государственный университет)
Выбор оптимального метода моделирования распространения импульсов в волоконных световодах
Правильный выбор метода моделирования распространения импульсов в волоконных световодах позволяет существенно сократить время моделирования с заданной точностью. В данной работе выполняется сравление и анализ наиболее перспективных модификаций фурье-метода расщепления по физическим процессам. В первой части работы приводится обобщённое уравнение распространения импульсов в волоконном световоде, выполняются замены переменных для удобного расчёта распространения импульсов в солитонном режиме. Численное решение сравнивается с точным аналитическим решением с заданной точностью для солитонов НО порядков. Выполняется сравнительный анализ четырёх методов. Показано, что выигрыш в скорости одного метода по сравнению с другим достигает до 327%.
Ключевые слова: оптический световод, распространение импульсов в волоконном световоде, фурье-метод расщепления по физическим процессам, нелинейное уравление Шредингера.
А. V. Gutor, Е. S. Manuylovich Moscow Institute of Physics and Technology (State University)
Selection of the optimal method for modelling of pulse propagation in optical fibers
The proper choice of a method for modeling pulse propagation in optical fibers can significantly reduce simulation time with given accuracy. Various modifications of the split step Fourier method (SSFM) are analyzed and compared in this paper. The first part contains the generalized Schrodinger Equation describing pulse propagation in optical fibers. A change of variables for convenient calculation of pulse propagation in the soliton mode is performed. The numerical solution is compared with the exact analytical solution with given accuracy for solitons of 1-10 orders. A comparative analysis of four methods is carried out. Significant performance gain up to 327% of one method in comparison with others is demonstrated.
Key words: optical fiber, optical pulse propagation, split-step Fourier method, nonlinear Schrodinger Equation.
1. Введение
Волоконно-оптические линии передачи являются на сегодняшний момент основными каналами передачи информации. В волоконном световоде посредством модуляции лазерного излучения инфракрасного диапазона (длина волны излучения вблизи 1550 нм) работают наиболее широкие пропускные каналы информации. Промышленные волоконно-оптические лазеры (длина волны излучения вблизи 1064 нм) с КПД, приближающимся к 50% «от розетки», и не требующие значительного обслуживания, заняли лидирующее
© Гутор А. В., Мануйлович Е. С., 2018
(с) Федеральное государственное автономное образовательное учреждение высшего образования «Московский физико-технический институт (государственный университет)», 2018
место в обработке материалов. Волоконно-оптические датчики по некоторым своим характеристикам сильно превосходят аналоги, спроектированные на базе предшествующих технологий.
Проектирование компонент и оптических систем, макетирование волоконно-оптических схем затрудняется без теоретического апробирования и сопровождения. В первом приближении линейной среды, когда показатель преломления считается постоянной величиной, возможно использовать аналитические оценки в расчётах. По мере приближения к границам возможностей технологий, эти оценки не дают достаточного представления о физических процессах. Появляется спрос на программные продукты, позволяющие учитывать большую часть физических процессов в случае распространения импульсов в волоконно-оптических световодах. В большинстве случаев такие световоды имеют ступенчатый профиль показателя преломления [1]. Наиболее значимые эффекты, которые учитываются в расчётах: затухание, дисперсия, фазовая само- и кросс-модуляция, комбинационное рассеяние, нелинейные эффекты высших порядков.
Начиная с конца 90-х годов для системного решения описанных выше задач ведутся исследования и разрабатываются эффективные методы решения обобщённого управления распространения [1]:
- 2А+г £ Чт ££+1+~ £ > |A(z, *) / , * -i01 w I • (1)
n=l V 0
Многие проекты переросли в коммерческие продукты, когда, помимо самого решения уравнения, добавлены инструменты удобного задания начальных условий, визуализации физических процессов, подобно измерениям на измерительных приборах, таких как осциллограф, оптический анализатор спектра, анализатор ошибок, измеритель шумов и др. К таким продуктам относится программное обеспечение OptSim, VPItransmissionMaker, Optiwave Photonic. В большинстве случаев оно направлено на решение инженерных задач. Однако в применении к научному поиску такие программы не дают полного понимания механизма и критериев расчётов. Поэтому задача по эффективному решению уравнения (1) и анализа результатов актуальна.
В работах [3] 1984 г. и [4] 1999 г. показано преимущество фурье-метода расщепления по физическим процессам по сравнению с традиционными сеточными численными методами. Выбор модификаций и оптимизация этого метода позволяют существенно выиграть в скорости получения решения с заданной точностью [6-11]. Кроме того, в литературе встречаются работы, в которых предлагается использовать другие методы. Так, конкуретными по скорости являются вычисления методами на базе алгоритма Кренка-Николсона, сравнительный анализ и суть которого раскрываются в статьях [9], [11]. В другой работе [12] вместо гармонических функций предлагается использование вейвлетов.
2. Постановка задачи и описание численного метода решения
Целью настоящего исследования является определение оптимального шага по переменкой точности. Для описания распространения электромагнитных импульсов в волоконном световоде необходимо численно решить обобщенное уравнение распространения для медленно меняющейся амплитуды A(z, t) (1).
Произведем стандартную замену переменных:
TT А A z Т n2 LD ч 1 . & m
Перепишем уравнение (1) в новых переменных, оставив дисперсионные члены до 3 порядка включительно. В новых обозначениях получаем
Щ = - 2и - 12*9п№ ^ - + гМ2 + г) ] Ж^и^т-т')^ .
° (3)
Для расчетов будем пользоваться симметризованпой схемой второго порядка точности по шагу. Для этого представим уравнение (3) в операторном виде:
ди
Ж
Ь + М(и) и. (4)
Здесь
х % д2 д3
Ь = -2 - 281ап(Р2)- 5дЙ3 (5)
и
(те
и (а, т)^К(1')р ^ ,т -т')|^т'| . (6)
Функция отклика среды в нелинейном операторе может быть представлена в виде (1):
Я^) = (1 - /ц)5(1) + /КНК, (7)
где первое слагаемое (с дельта-функцией) отвечает за мгновенный отклик, а второе - за рамановский (запаздывающий) отклик, функция рамановского отклика определяется из эксперимента по спектру ВКР-усиления, здесь будем пользоваться моделью [2], в которой
Ък = (т-2 + т-2)п ехр(-—) ехр(-), 1 2 2
где п = 12,1 фс, т2 = 32 фс. Примем также /я = 0,18 [1]. После подстановки (7) в (6) получаем
М(и) = 1(1 - /к)М2\и|2 + г/КМ\ош(НК')-
N2 д М2 д
-(1 - Ы—)—[|и|2и] - ^Ы—)—\и^ш(кК'|и|2)].
те
Здесь сопу(йд'|и|2) = / Ък(т')р(т - т'^йт'.
°
Решать уравнение (4) будем фурье-методом расщепления по физическим параметрам, то есть будем считать, что на каждом шаге дисперсия и нелинейность действуют независимо:
ц+н
Ъ ~ Г - Ъ -
и(а + к,Т)=ехр(2Ь)ехр(! ,и)с!()ехр(2Ь)и(а,Т). (8)
?
В этой схеме оператор нелинейности действует в середине шага, что позволяет сделать
Ъ
Интеграл в (8) в некоторых случаях приблизим трапецией:
?+н
I М((, и УК = 2[М(0 + М(а + ЪЪ)]. (9)
?
Так как в начале каждого шага мы не знаем, чему равно N(0, + Ъ) будем вычислять его по упрощенной схеме, в которой интеграл заменим просто на Ъ N(£,)'■
ЪЪ и' (а + Ъ,Т)= ехр( 2Ь) ■ еМЪЩа)) ■ ехр(-Ь) ■ и (а ,Т), N(4 + Ъ) = IV (и'). (10)
3. Определение оптимального метода расчета
Характерные длины, на которых проявляются нелинейные и дисперсионные эффекты, т2 1
составляют Lp = -щ и Ln = ^рт- Логично выбирать шаг интегрирования много меньше каждой из этих длин. Можно, например, взять h = 0, 01 ■ min(Lp,Ln), однако по мере распространения импульса пиковая мощность импульса и его длительность могут существенно меняться. Так, в случае распространения солитона, например, третьего порядка, его пиковая мощность в разные моменты эволюции может различаться в 3 раза. Требуется
метод известен и используется, например, в [6].
В ходе разработки программы расчета были опробованы различные схемы численного решения:
h
step):
hh U (С + h,T)= exp( ^D) ■ exp(hN(0) ■ exp( ^D) ■ U (£, T), (11)
h = con st.
h
(8), (9), (10) (precise; constant step):
h h h U(Z + h,T)= exp(2D) ■ exp(^[N(0 + N(£, + h)]) ■ exp(^D) ■ U(£,T), (12)
N(£ + h) = N(U'),
hh U' (£ + h,T)= exp( ^D) ■ exp(hN(0) ■ exp( ^D) ■ U (Ц, T),
h = .
h
step):
hh U (i + h,T)= exp( ^D) ■ exp(hN(£)) ■ exp( ^D) ■ U (Ц, T), (13)
, ho
h =
(max |U (Ol)2' h
variable step):
h h h U(£ + h,T)= exp(^D) ■ exp(^[N(0 + N(£, + h)]) ■ exp(^D) ■ U(£,T), (14)
N(( + h) = N(U'),
)(M ho
(max|U(0\)2.
hh U' (£ + h,T)= exp( ^D) ■ exp(hN(0) ■ exp( ^D) ■ U (£, T),
h=
Были произведены расчеты по всем указанным схемам, а также по схеме, предложенной в [6], с целью определения оптимального по производительности метода. Для определения точности метода удобнее всего сравнивать численное решение с точным аналитическим. Так как уравнение распространения в общем виде (3) не имеет аналитических решений, для
оценки точности методов будем решать уравнение в случае, когда а = 0. Проводя аналогию, его иногда называют обобщенным нелинейным уравнением Шредингера (ОНУШ):
ВЛ -i rßjj
= --Signet) + iN 2| Л |2Л. (15)
Аналитическое решение этого уравнения при начальном условии Ло = у/Pq • sech(T/To) и при таких Ро и То, что 7PoTg/|ß21 = N2, является квадратом целого числа и солитоном порядка N. Оно периодично по z с периодом ад = или, в безразмерных координатах, с периодом ^s = В качестве функции ошибки будем использовать две метрики:
/ || Ло| — | Л(ад)||2dr
err"-" =1—тш*— (16)
т
егг|^| = |тах|Ло |—тах|Л (zs )||. (17)
Здесь Ло - комплексная амплитуда при г = 0. Как мы знаем из аналитического решения, она же - при г = ад. В качестве приемлемого уровня ошибки выбираем егг^ц = 10-3 и ориентировочное время расчета приводим именно для такого уровня. В статьях также часто используется данное предельное значение. Графики приводятся вплоть до 10-4.
Будем рассматривать импульсы длительностью То = 1 пс и центральной длиной волны Ао = 1550 нм, распространяющиеся в волокне с дисперсией групповых скоростей
пг2 1
р2 = -20 ^ и параметром 7 = 26км.Вт- Дисперсионная длина при таких параметрах составляет Ьр = 50 м, а период солитопа примерно ад = 78, 5 м.
Разбиение по Т должно быть таково, чтобы выполнялось йТ ^ Трщнм(%) V г € [0, ад]. В наших тестах будем брать йТ < 4 фс. Когда мы решаем уравнение (15), в таком мелком разбиении нет необходимости, однако при решении уравнения распространения (3) в
Т
функции к и (т).
и
4. Результаты вычислений и их анализ
Начнем с метода (11), с самого простого. Рассчитаем точность метода в зависимости от шага интегрирования для солитонов различных порядков. Шаг в методе (11) постоянный, будем выбирать его в виде к = е • тш(Ьд,Ь^), где е ^ 1. Начнем с солитона первого порядка, для него Ь^ = Ьр = 50 м.
•1.5 -1 -0.5 -1.5 -1 -0.5
log ю(^) toflio(^)
шаге для N =1. Метод (11) (coarse; constant step). Время счёта 0,1 с
На графике слева в логарифмических координатах изображены ошибки вычисления
грирования. Пунктиром обозначен уровень ошибки 10-3. На графике справа в логарифмических координатах приведено время расчета с соответствующим шагом. Треугольником и крестиком обозначено время расчета с точностью 10-3 (на основе интерполяции по левому графику).
Как видно из рис. 1, уже при шаге Ъ° = 0,1 ■ тш(£д) достигается точность 10-3,
при этом время расчета распространения на период солитона составляет всего ~ 0,1с. Из
Ъ°
мических координатах функции ошибки линейно зависят от шага, а тангенс угла наклона прямой равен 2. Приведем также аналогичный рисунок для солитона пятого порядка.
Рис. 2. Зависимость ошибки и времени вычисления от шага интегрирования по z при постоянном шаге для N = 5. Метод (11) (coarse; constant step). Время счёта 23 с
Здесь заданная точность достигается при шаге h0 ~ 0, 0073 ■ Ln, время расчета с заданной точностью уже около 23 секунд. Аналогичные расчеты были произведены для солито-нов 1 10 порядков, время расчета с заданной точностью приведено на рис. 3.
Рис. 3. Время расчета распространения солитопов различных порядков на расстояние zs с заданной точностью. Результаты метода (11) (coarse: constant step)
78, 5 м
Из рис. 3 видно, что время расчета, например, солитона 10 порядка составляет более 11 минут.
Метод (12) или precise; constant step, оказывается более медленным для солитонов порядка вплоть до 4, а начиная с 5 уже дает выигрыш во времени расчета с заданной точностью. На рис. 3 приведена точность метода в зависимости от шага интегрирования при расчете распространения солитона 5 порядка.
N = 5
Ъ
наклона кривых на левом графике составляет около 4, что соответствует методу 4 порядка точности. Однако для солитонов более высоких порядков средний наклон кривой уменьшается и составляет около 3. На рис. 5 приведена аналогичная зависимость для солитона 9 порядка.
шаге для N = 9. Метод (12) (precise; constant step). Время счёта 142,3 с
На следующем рисунке приведем зависимость времени вычисления с точностью 10-3 для солитонов различных порядков.
10-3
составляет 211 секунд, что более чем в 3 раза быстрее метода (11).
Метод (13) или coarse; variable step, дает очень хорошие результаты по производительности. Этим методом также было рассчитано распространение солитонов 1 10 порядков.
Начиная с солитона второго порядка, он быстрее метода (11) и для солитонов всех поряд-
10-3
для солитонов различных порядков приведена на рис. 7.
Рис. 6. Время расчета распространения солитонов различных порядков на расстояние zs ~ 78, 5 м с заданной точностью. Результаты метода (12) (precise: constant step)
Рис. 7. Время расчета распространения солитонов различных порядков на расстояние zs ~ 78, 5 м с заданной точностью. Результаты метода (13) (coarse: variable step)
Метод (14) или precise; variable step, является комбинацией методов (12) и (13). Можно было ожидать, что этот метод окажется наиболее эффективным.
Однако, как оказалось, это не так. Он проигрывает методам (12) и (13) и выигрывает лишь у метода (11). По крайней мере, при вычислении с точностью 10-3. На следующем рисунке приведена производительность метода (14) для солитонов различных порядков.
Вполне может оказаться, что этот метод эффективнее (вследствие теоретически более высокой точности) при расчете распространения более мощных импульсов. На практике практически не встречаются ситуации, когда энергия импульса в волоконном световоде
соответствует энергии солитона таких высоких порядков тального солитона).
(более 100 энергий фундамен-
Рис. 8. Время расчета распространения солитонов различных порядков на расстояние zs ~ 78, 5 м
с заданной точностью. Результаты метода (14) (precise; variable step)
Приведем сравнительный график производительности всех четырех методов.
123456789 10
N
Рис. 9. Результаты сравнитштьного теста всех четырех методов: о - метод (11) (coarse; constant step), А - метод (12) (precise; constant step), х - метод (13) (coarse; variable step), • - метод (14) (precise; variable step)
Время расчета с заданной точностью в упрощенной модели (4) и время в полной модели (12) существенно различаются, так как в полной модели каждый шаг содержит больше
данная точность, будет тем же самым.
Расчет задачи распространения в среде с положительной дисперсией требует гораздо меньше времени для достижения заданной точности, так как амплитуда резко снижается, также снижается роль нелинейности, а именно расчет нелинейного оператора является преимущественным местом накопления ошибок. Таким образом, при расчете задачи в среде с положительной дисперсией можно пользоваться любыми методами, хотя предпочтительными являются методы (12) и (13).
Таким образом, показано, что оптимальный выбор метода вычисления позволяет в несколько раз уменьшить время вычислений, что дает возможность получить рабочий приемлемый инструмент для теоретических исследований.
Литература
1. Agrawal G.P. Nonlinear Fiber Optics (Fifth Edition) // Academic Press. 2012. 629 p. P. 5358.
2. Blow K.J., Wood D. Theoretical description of transient stimulated Raman scattering in optical fibers // IEEE Journal of Quantum Electronics. 1989. V. 25, I. 12. P. 2665-2673.
3. Taha T.R., Ablowitz M.I. Analytical and numerical aspects of certain nonlinear evolution equations. II. Numerical, nonlinear Schrôdinger equation // Journal of Computational Physics. August 1984. V. 55, I. 2. P. 203-230.
4. Chang Q., Jia E., Sun, W. Difference Schemes for Solving the Generalized Nonlinear Schrôdinger Equation // Journal of Computational Physics. 20 January 1999. V. 148, I. 2. P. 397-415.
5. Liu X., Lee B. A fast method for nonlinear Schrôdinger equation // IEEE Photonics Technology Letters. November 2003. V. 15, I. 11. P. 1549-1551.
6. Sinkin О. V., Holzlôhner R., Zweck J., Menyuk C.R. Optimization of the split-step Fourier method in modeling optical-fiber communications systems // Journal of Lightwave Technology. January 2003. V. 21, I. 1. P. 61-68.
7. Shao J., Liang X., Kumar S. Comparison of Split-Step Fourier Schemes for Simulating Fiber Optic Communication Systems // IEEE Photonics JournalOpen Access. 1 August 2014. V. 6, I. 4. Article number 7200515.
8. Balac S., Fernandez A. Mathematical analysis of adaptive step-size techniques when solving the nonlinear Schrôdinger equation for simulating light-wave propagation in optical fibers // Optics Communications. 15 October 2014. V. 329. P. 1-9.
9. Mohammadi R. An exponential spline solution of nonlinear Schrôdinger equations with constant and variable coefficients // Computer Physics Communications. March 2014. V. 185, I. 3. P. 917-932.
10. Secondini M., Marsella D., Forestieri E. Enhanced split-step Fourier method for digital backpropagation // European Conference on Optical Communication, ECOC. 20 November 2014, Article number 6964122. September 2014, Code 109416.
11. Ba§han A., Uçar Y., Murat Yagmurlu N., Esen A. A new perspective for quintic B-spline based Crank-Nicolson-differential quadrature method algorithm for numerical solutions of the nonlinear Schrôdinger equation // European Physical Journal Plus. 1 January 2018. V. 133, I. 1. Article number 12.
12. Wang J., Liu X., Zhou, Y. A high-order accurate wavelet method for solving Schrôdinger equations with general nonlinearitv // Applied Mathematics and Mechanics (English Edition). 1 February 2018. V. 39, I. 2.
References
1. Agrawal G.P. Nonlinear Fiber Optics (Fifth Edition). Academic Press. 2012. 629 p. P. 53-58.
2. Blow K.J., Wood D. Theoretical description of transient stimulated Raman scattering in optical fibers. IEEE Journal of Quantum Electronics. 1989. V. 25, I. 12. P. 2665-2673.
3. Taha T.R., Ablowitz M.I. Analytical and numerical aspects of certain nonlinear evolution equations. II. Numerical, nonlinear Schrôdinger equation. Journal of Computational Physics. August 1984. V. 55, I. 2. P. 203-230.
4. Chang Q., Jia E., Sun, W. Difference Schemes for Solving the Generalized Nonlinear Schrôdinger Equation. Journal of Computational Physics. 20 January 1999. V. 148, I. 2. P. 397-415.
5. Liu X., Lee B. A fast method for nonlinear Schrôdinger equation. IEEE Photonics Technology Letters. November 2003. V. 15, I. 11. P. 1549-1551.
6. Sinkin O.V., Holzlôhner R., Zweck J., Menyuk C.R. Optimization of the split-step Fourier method in modeling optical-fiber communications systems. Journal of Lightwave Technology. January 2003. V. 21, I. 1. P. 61-68.
7. Shao J., Liang X., Kumar S. Comparison of Split-Step Fourier Schemes for Simulating Fiber Optic Communication Systems. IEEE Photonics JournalOpen Access. 1 August 2014. V. 6, I. 4. Article number 7200515.
8. Balac S., Fernandez A. Mathematical analysis of adaptive step-size techniques when solving the nonlinear Schrôdinger equation for simulating light-wave propagation in optical fibers. Optics Communications. 15 October 2014. V. 329. Pages 1-9.
9. Mohammadi R. An exponential spline solution of nonlinear Schrôdinger equations with constant and variable coefficients. Computer Physics Communications. March 2014. V. 185, 1.3. P. 917-932.
10. Secondini M., Marsella D., Forestieri E. Enhanced split-step Fourier method for digital backpropagation. European Conference on Optical Communication, ECOC. 20 November 2014, Article number 6964122. September 2014, Code 109416.
11. Ba§han A., Uçar Y., Murat Yagmurlu N., Esen A. A new perspective for quintic B-spline based Crank-Nicolson-differential quadrature method algorithm for numerical solutions of the nonlinear Schrôdinger equation. European Physical Journal Plus. 1 January 2018. V. 133, I. 1. Article number 12.
12. Wang J., Liu X., Zhou, Y. A high-order accurate wavelet method for solving Schrôdinger equations with general nonlinearitv. Applied Mathematics and Mechanics (English Edition). 1 February 2018. V. 39, I. 2.
Поступим в редакцию 31.05.2018
