CC BY
56
Вычислительные технологии
Том 20, № 3, 2015
Использование аппроксимации Паде для решения систем нелинейных уравнений Шрёдингера с помощью метода расщепления по физическим процессам
И. С. Чеховской
Новосибирский государственный университет, Россия
Институт вычислительных технологий СО РАН, Новосибирск, Россия
Контактный e-mail: [email protected]
Для метода Фурье расщепления по физическим процессам разработана модификация на основе аппроксимации Паде матричной экспоненты. Эта модификация расширяет возможность применения метода Фурье на системы связанных нелинейных уравнений Шредингера, описывающих распространение световых импульсов в многоядерных и многомодовых оптических волокнах. Предложена параллельная реализация данного численного алгоритма на вычислительных системах с общей памятью.
Ключевые слова: метод Фурье расщепления по физическим процессам, аппроксимация Паде, нелинейное уравнение Шредингера, многоядерные оптические волокна, многомодовые оптические волокна, нелинейная волоконная оптика.
Введение
Система связанных нелинейных уравнений Шредингера (НУШ) широко используется при математическом моделировании в различных областях нелинейной волоконной оптики [1,2]. С ее помощью может быть описана эволюция комплексных огибающих импульсов, распространяющихся вдоль оптического волокна с сильной связью между соседними модами света. Примером таких волокон, где нужно учитывать взаимодействие оптических импульсов друг с другом, являются многоядерные оптические волокна (multicore fibers — MCF), представляющие собой набор близко расположенных волноводов. Аналогичное взаимодействие возникает при распространении нескольких мод света в одном волноводе (многомодовые оптические волокна) [3,4].
Для решения систем связанных НУШ (через нелинейную часть) разработаны конечно-разностные схемы [5,6], в том числе повышенного порядка аппроксимации [7], псевдоспектральные методы [8], мультисимплектические методы расщепления [9]. Также был предложен параллельный алгоритм SSFM для решения связанных через нелинейность систем НУШ [10] и набор спектральных методов разных порядков аппроксимации [11]. Однако в силу отсутствия линейных связей между комплексными переменными при численном решении рассматриваемых уравнений не возникало проблем с реализацией алгоритмов для решения линейной части уравнений, рассмотренных в работах [5-11].
© ИВТ СО РАН, 2015
Существует много численных методов решения нелинейного уравнения Шрёдингера с одной переменной. В классе спектральных алгоритмов выделяется высокоточный метод Фурье расщепления по физическим процессам (split-step Fourier method — SSFM), в котором используется предположение о том, что при распространении оптического импульса по волоконному световоду линейные и нелинейные эффекты действуют независимо. Этот метод обладает большей точностью по сравнению с конечно-разностными методами благодаря использованию преобразования Фурье для вычисления второй временной производной, поскольку это обеспечивает экспоненциальную скорость сходимости по шагу временной переменной.
В работе рассматривается периодическая начально-краевая задача для системы связанных НУШ
на комплексные функции Ак = Ак(г^), к =1, ...,М, где г Е [0,Ь] — пространственная переменная, а Ь Е [-Т, Т] — временная. Вещественные коэффициенты и 7 постоянны. Матрица С = (Ск,т) задает линейные связи между функциями Ак. При обобщении ББЕМ на системы НУШ вида (1) возникают определенные сложности, вызванные необходимостью учета линейных связей между неизвестными функциями Ак.
В настоящей работе предлагается использовать аппроксимацию Паде с 13-м порядком точности для вычисления линейного экспоненциального оператора системы связанных НУШ. Данный оператор возникает при попытке расширения применения стандартного ББЕМ на системы нелинейных уравнений Шредингера с линейными связями между переменными. Такая модификация численного алгоритма ББЕМ оказывается удобной для распараллеливания. Ее использование позволяет эффективно сократить время расчета, так как основное время при решении системы НУШ затрачивается не на алгоритм быстрого преобразования Фурье, который обладает низкой эффективностью распараллеливания, а на матрично-векторные операции, которые, как известно, хорошо распараллеливаются.
Численная аппроксимация Паде различных функциональных зависимостей широко применяется благодаря простоте вычисления и более высокому порядку сходимости по сравнению с аппроксимацией, например, рядом Тейлора [12]. Так, в работе [13] для численного решения уравнения Гельмгольца с помощью ББЕМ использовалось приближение нелинейного оператора с помощью аппроксимации Паде. В ряде разностных схем для НУШ высокого порядка применяется аппроксимация Паде для вычисления линейного оператора [14,15]. Аппроксимация Паде применялась и для расчета линейного оператора при численном решении системы связанных НУШ с помощью метода Рунге — Кутты [16].
1. Метод Фурье расщепления по физическим процессам
(1)
Ак (0,t) = А0>к (t), Ак(z,T) = Ак(z, -Т)
(3)
Данный численный алгоритм является одним из представителей спектральных численных методов, обладающих высокой точностью расчета по временной переменной [17,18].
Рассматриваемый метод Фурье расщепления по физическим процессам был адаптирован для системы связанных НУШ специальным образом. Перепишем сначала систему (1) в операторном виде:
8 А
= (И + М)А, (4)
где А = (Ах,..., А^)т — вектор из комплексных переменных, И — линейный дифференциальный оператор, задающий дисперсию и линейные связи между различными переменными Ак, к = 1...М, а N — оператор, отвечающий за нелинейность. Для рассматриваемой системы уравнений эти операторы имеют следующую структуру:
" = ' (С-, (5)
N4 = ё1аё{г т1 Ах |2, г т|А2|2 ,...,г^у |АМ |2} , (6)
где I — тождественный оператор, С — матрица линейных связей. В ББЕМ приближенное решение получают, предполагая, что при распространении оптического поля на один шаг к по эволюционной переменной г линейные и нелинейные эффекты действуют независимо. В частности, распространение от г до г + к осуществляется в два этапа. На первом этапе действует только нелинейность и оператор И = 0. На втором этапе учитываются только линейные эффекты и поэтому нелинейный оператор N4 = 0. На практике для достижения второго порядка аппроксимации используется следующая симметризованная форма алгоритма:
А(г + к, г) и ехр(^1 к^ ехр(кЁ) ехр кЖ^ А(г, г), (7)
которая при многократном применении для получения итогового решения на расстоянии Ь = кМх может быть представлена в виде
А(Ь, ^ и ехр ^-1 кN^ exp(кN) ехр(кГ)^ ехр к^ А(0, Ь). (8)
При обобщении ББЕМ на систему НУШ (4) вычисление нелинейного оператора ехр(ак IV), а = {0.5, -0.5}, в силу его диагональности тривиальным образом распадается на расчет отдельных действий нелинейного оператора для каждой переменной Ак, к = 1,...,N. Поскольку линейный оператор И представляет собой недиагональную матрицу, возникает проблема реализации линейного шага, а именно вычисления матричной экспоненты ехр(кИ). Была использована следующая схема ее вычисления. Для упрощения нахождения второй производной по Ь оператор ехр(к И) вычислялся в фурье-пространстве с помощью представления
ехр(кИ) А(г, Ь) = ехр[кЬ(-1 ш)ЩА(г, Ь), (9)
где Р^, обозначает операцию применения преобразования Фурье, А(х, ¿) — промежуточное значение вектора из N переменных, которое может быть записано в виде
, *) = (П ехр(к1)) ехр(кГ))^ ехр к^ А(0, Ь), (10)
а оператор О(-гш) получен из дисперсионного оператора (5) заменой производной д/дЬ на — гш, где ш — круговая частота в фурье-пространстве. Полученная матрица ехр[кИ(—гш)} является постоянной, поэтому в случае использования равномерной сетки по г ее вычисление достаточно осуществить один раз перед началом счета.
Расчет матричной экспоненты может быть проведен различными способами [19]. Хорошо зарекомендовал себя на практике метод численной аппроксимации Паде для экспоненты (например, по формуле 13-го порядка точности) [19,20],
Re(x)
(2 ■ 6 - j)!6! Xj
j=0 (2 ■ 6)!(6 - j)! "J
L
L
(2 ■ 6 - j)!6! (-X))
= (2 ■ 6)!(6 - j)! j!
exp[X], (11)
в которой используется приближение экспоненты рациональной функцией с многочленами 6-й степени в числителе и знаменателе, чем и обеспечивается порядок точности 0(Х13).
При вычислении матричной экспоненты удобно использовать соотношение ех = (ех/а)а для перенормировки линейного оператора (scaling and squaring), где X — матрица размера N х N, а а — некоторая ненулевая константа. Такой подход позволяет
и \
значительно увеличить точность вычисления матричной экспоненты е^ по сравнению с прямым применением формулы численного приближения (11) [19]. Вначале выбирается некоторое число а, такое, чтобы норма матрицы Х/а была порядка 1. Затем вычисляется аппроксимация Паде (11) для нормированной матрицы, которая возводится в степень а.
Предлагаемый численный алгоритм универсален и может быть применен к системам НУШ с матрицами любого вида. При моделировании динамики оптических импульсов в многоядерных световодах возможен другой подход, использующий для диа-гонализации матрицы линейных связей С разложение в ряд Фурье по дискретному пространству ядер. Стоит заметить, что такое разложение в ряд Фурье по к может быть использовано лишь в случае периодических граничных условий по ядрам, что не выполняется в общем случае. Также применение дискретного разложения в ряд Фурье для диагонализации матрицы С возможно в случае хорошей локализации решения по пространству, например в задачах об образовании пространственно-временного оптического солитона. Однако предложенный алгоритм с использованием аппроксимации Паде был разработан для применения в задачах расчета динамики набора временных солитонов, введенных в структуру, состоящую из относительно небольшого числа ядер, и их последующего сжатия в один мощный оптический импульс. Кроме того, представленный алгоритм эффективно распараллеливается по ядрам, что будет представлено далее, тогда как ускорение работы вычислительного алгоритма с использованием преобразования Фурье по пространству затруднительно.
2. Параллельная реализация метода и использованные технологии
В программной реализации представленного алгоритма применялась специализированная математическая библиотека Intel MKL. В частности, использовались функции из компонентов DFT, BLAS, LAPACK и VML. Преобразование Фурье выполнялось с помощью параллельной версии быстрого преобразования Фурье (FFT) из компонента DFT,
имеющего, как известно, асимптотическую сложность 0(Mt logMt). Для этого число узлов по временной переменной Mt всегда равнялось числу 2 в некоторой степени. Действие оператора, соответствующего кубической нелинейности в системе НУШ, рассчитывалось с помощью функций компонента VML (операции над векторами), в частности поэлементного нахождения экспоненты векторов и поэлементного возведения в квадрат.
Матричная экспонента, вычисляемая с помощью аппроксимации Паде (11), находилась с использованием функции zgemm (умножение матрицы на матрицу) из компонента BLAS, а также функций из компонента LAPACK, выполняющих обращение матриц. Сама матричная экспонента хранилась в оперативной памяти в виде двумерного комплексного массива размера N2 х Mt, где N — число уравнений в системе (1). Умножение вектора неизвестных величин на матричную экспоненту осуществлялось не прямым способом по определению (т. е. не как умножение Mt раз матрицы размера N х N на вектор длины N), а с помощью N2 поэлементных умножений и сложений одной из N2 строк длины Mt матричной экспоненты на временное распределение одной из N компонент вектора FtA(z) длиной Mt (рис. 1). Такой подход существенно сокращает количество вызовов функций из библиотеки MKL, а также улучшает доступ к памяти. Операции поэлементного сложения и умножения векторов выполнялись высокооптимизированными функциями из компонента VML.
Матрица линейных связей С обычно является разреженной. Однако матричная экспонента уже представляет собой плотную матрицу, а значит, умножение на нее вектора неизвестных проводится за 0(N2) операций. Таким образом, вычисление линейного шага — самая трудозатратная операция при достаточно большом количестве уравнений в исходной системе (7). Стоит заменить, что для расчета систем НУШ с малым числом уравнений не существует более быстрых алгоритмов.
Распараллеливание алгоритма с использованием технологии OpenMP позволяет существенно сократить время расчета на вычислительных системах с общей памятью. Для эффективной загрузки всех вычислительных потоков удобно применять разделение по обрабатываемым данным при расчете умножения матричной экспоненты на вектор неизвестных, а также при вычислении действия нелинейного оператора. Так как
Рис. 1. Схема умножения матричной экспоненты ехр[Ь1)(—гш)] на промежуточное значение временного распределения вектора переменных ^А(^) для расчета действия линейных эффектов. А!(г) — результат умножения
Рис. 2. Схема разделения обрабатываемых алгоритмом данных по т потокам OpenMP
Mt = 2га, при числе потоков ОрепМР, также являющемся числом 2 в некоторой степени т, эффективно разделение обрабатываемых данных, при котором каждому потоку достаются М4/2т узлов по временной переменной (рис. 2). Ускорение алгоритма, а также его эффективность на некоторых тестах рассмотрены в следующем разделе.
На основе этого алгоритма был реализован и запатентован вычислительный программный комплекс для решения задач распространения оптических импульсов в многоядерных волноводах с расположенными по окружности ядрами "МиШСогеР1Ьег-1" (И.С. Чеховской, О.В. Штырина, М.П. Федорук, свидетельство о государственной регистрации программы для ЭВМ № 2015610111). Результаты численных экспериментов по нелинейному сжатию и сложению импульсов в многоядерных оптических волокнах с использованием данного программного комплекса представлены в работе [21].
3. Тестовые расчеты, анализ алгоритма
Вычисления проводились на сервере с общей памятью HP DL980 G7 кластера НГУ, содержащем восемь 8-ядерных процессоров Intel Xeon X7560 с тактовой частотой 2266 МГц и оперативной памятью 2 097 152 МБ.
Метод Фурье расщепления по физическим процессам с аппроксимацией Паде матричной экспоненты тестировался с параметрами 7 =1 и = —1 на решении скалярного НУШ в виде фундаментального солитона
1
2 ) cosh(i)'
которое также является решением для системы НУШ
A(z,t)=exp( , (12)
В А 1 г) 2 А
г~д.V = — 2— (Ак+1 — 2^ + А—) — ' к = 1>->м> (13)
если для каждой переменной Ак распределение (12) задать в качестве начального. Данная система НУШ рассматривалась в работе [21] в качестве модели для описания распространения оптических импульсов в многоядерных оптических волокнах с круговым расположением ядер. Матрица линейных связей такой системы уравнений имеет вид
/-2 1 0 0 1-210
С
0 0 0 0 \ 1 0 0 0
0 0 0 0 0 0
1 0
0 1 -2 1 0 0 1 -2)
:14)
Расчеты проводились в области (0 < г < 10) х (-20 < £ < 20) для системы из шести уравнений (М = 6). Размер сетки по временной переменной М4 равнялся 29. С -норма ошибки аппроксимации
5 = тах тах |Ж!,- — (гпЛ4) 1 <к<Ы о<з<мг '
:15)
и коэффициент ее убывания при различном числе узлов по пространственной переменной г приведены в табл. 1.
На системе НУШ (13) также проанализировано ускорение параллельной реализации представленного алгоритма при различном числе потоков OpenMP на сетке Мх х М4 =
Таблица 1: Значение ошибки аппроксимации 8 на стационарном решении и коэффициент ее убывания в зависимости от размера сетки Мх
М2 200 400 800 1600 3200 6400
5 2.2939e-03 5.742^-04 1.4360e-04 3.5903e-05 8.9759e-06 2.2440e-06
К — 3.99 4.00 4.00 4.00 4.00
Таблица 2: Время расчетов в зависимости от количества потоков OpenMP, с
Число Количество потоков OpenMP
уравнений N 1 2 4 8 16
6 13.88 7.83 4.74 3.21 2.60
19 111.92 55.04 30.16 16.59 10.21
16
о в
В
£
' 10.96
/ 6.74 ^Ъ.ЪЪ
1.00 / г<
2 4
Число ПОТОКОВ
16
Линейное ускорение
N=6 19
Рис. 3. Ускорение выполнения вычислений при увеличении количества потоков OpenMP
2000 х 4096. Зависимость ускорения от числа потоков при N = 6 и N =19 приведена на рис. 3. Время расчетов представлено в табл. 2.
4. Заключение
С развитием телекоммуникационных технологий, таких, как, например, многоядерные оптические волокна, возникла необходимость численного решения нового класса систем связанных нелинейных уравнений Шрёдингера. Предложенная модификация широко известного метода Фурье расщепления по физическим процессам позволяет эффективно применять данный метод для численного решения систем НУШ с линейными связями между переменными. Матричная экспонента, естественным образом возникающая при обобщении метода Фурье на такие системы уравнений, может быть найдена с помощью численной аппроксимации Паде. В работе продемонстрирована высокая точность предложенного численного алгоритма и показана возможность его эффективного распараллеливания на вычислительных системах с общей памятью.
Благодарности. Работа выполнена при финансовой поддержке Российского научного фонда (Соглашение № 14-21-00110).
Список литературы / References
[1] Kivshar, Y.S., Agrawal, G.P. Optical solitons: from fibers to photonic crystals. 5th ed. San Diego: Academic Press, 2003. 540 p.
[2] Agrawal, G. Applications of Nonlinear Fiber Optics. Burlington, Ma: Academic Press, 2008. 508 p.
[3] Richardson, D.J., Fini, J.M., Nelson, L.E. Space-division multiplexing in optical fibres // Nature Photonics. 2013. Vol. 7, No. 5. P. 354-362.
[4] Mumtaz, S., Essiambre, R., Agrawal, G.P. Nonlinear propagation in multimode and multicore fibers: Generalization of the manakov equations // Journal of Lightwave Technology.
2013. Vol. 31, No. 3. P. 398-406.
[5] Wang, T. Maximum norm error bound of a linearized difference scheme for a coupled nonlinear schrodinger equations // Journal of Computational and Applied Mathematics. 2011. Vol. 235, No. 14. P. 4237-4250.
[6] Wang, D., Xiao, A., Yang, W. A linearly implicit conservative difference scheme for the space fractional coupled nonlinear schrodinger equations // Journal of Computational Physics.
2014. Vol. 272. P. 644-655.
[7] Ma, Y., Kong, L., Hong, J., Cao, Y. High-order compact splitting multisymplectic method for the coupled nonlinear schrodinger equations // Computers & Mathematics with Applications. 2011. Vol. 61(2). P. 319-333.
[8] Dehghan, M., Taleei, A. A Chebyshev pseudospectral multidomain method for the soliton solution of coupled nonlinear Schrodinger equations // Computer Physics Communications. 2011. Vol. 182(12). P. 2519-2529.
[9] Chen, Y., Zhu, H., Song, S. Multi-symplectic splitting method for the coupled nonlinear schrodinger equation // Computer Physics Communications. 2010. Vol. 181(7). P. 1231-1241.
[10] Taha, T.R., Xu, X. Parallel split-step fourier methods for the coupled nonlinear Schrodinger type equations // Journal of Supercomputing. 2005. Vol. 32(1). P. 5-23.
[11] Wang, S., Wang, T., Zhang, L. Numerical computations for n-coupled nonlinear Schrodinger equations by split step spectral methods // Applied Mathematics and Computation. 2013. Vol. 222. P. 438-452.
[12] Бейкер Д., Грейвс-Моррис П., Рахманов Е.А., Суетин С.П. Аппроксимации Паде: Основы теории. Обобщения и приложения. М.: Мир, 1986.
Baker, G., Graves-Morris, P. Pade approximants. Part I: Basic theory. Part II: Extensions and applications. Encyclopedia of Mathematics and Its Applications. 1981. Vol. 13-14.
[13] Zhang, L., Rector, J.W., Hoversten, G.M., Fomel, S. Split-step complex pade-fourier depth migration // Geophysical Journal International. 2007. Vol. 171(3). P. 1308-1313.
[14] Smadi, M., Bahloul, D. A compact split step pade scheme for higher-order nonlinear schrodinger equation (hnls) with power law nonlinearity and fourth order dispersion // Computer Physics Communications. 2011. Vol. 182(2). P. 366-371.
[15] Xu, Z., He, J., Han, H. Semi-implicit operator splitting pade method for higher-order nonlinear schrodinger equations // Applied Mathematics and Computation. 2006. Vol. 179(2). P. 596-605.
[16] Bhatt, H.P., Khaliq, A.Q. Higher order exponential time differencing scheme for system of coupled nonlinear schrodinger equations // Applied Mathematics and Computation. 2014. Vol. 228. P. 271-291.
[17] Agrawal, G. Nonlinear Fibre Optics. 4th ed. San Diego: Academic Press, 2013.
[18] Turitsyn, S.K., Bale, B.G., Fedoruk, M.P. Dispersion-managed solitons in fibre systems and lasers // Physics Reports. 2012. Vol. 521(4). P. 135-203.
[19] Moler, C., Van Loan, C. Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later // SIAM Review. 2003. Vol. 45(1). P. 3-49.
[20] Higham, N.J. The scaling and squaring method for the matrix exponential revisited // SIAM Journal on Matrix Analysis and Applications. 2005. Vol. 26(4). P. 1179-1193.
[21] Rubenchik, A.M., Chekhovskoy, I.S., Fedoruk, M.P., Shtyrina, O.V., Turitsyn, S.K.
Nonlinear pulse combining and pulse compression in multi-core fibers // Optics Letters. 2015. Vol. 40(5). P. 721-724.
Поступила в 'редакцию 6 апреля 2015 г., с доработки —15 апреля 2015 г.
Using Pade approximation for solving systems of nonlinear Schrodinger equations by the split-step Fourier method
Chekhovskoy, Igor S.
Novosibirsk State University, Novosibirsk, 630090, Russia
Institute of Computational Technologies SB RAS, Novosibirsk, 630090, Russia
Corresponding author: Chekhovskoy, Igor S., e-mail: [email protected]
Purpose. This paper presents a new modification of the well-known split-step Fourier method. The development of telecommunication technologies caused the necessity to solve numerically new types of systems of coupled nonlinear Schrodinger
© ICT SB RAS, 2015
108
H. C. HexoBCKoA
equations. For devices with strong coupling between the modes of light, such as multi-core optical fibers, there is a need to consider the linear coupling between the various modes. Presented numerical algorithm can solve the NSE system, taking into account such linear coupling.
Methodology. Pade approximation for the numerical calculation of the matrix exponential was used. The matrix exponential appears in the generalization of the split-step Fourier method to the system of coupled NSE. Pade approximation formula uses the approximation of exponential by rational functions with polynomials of degree 6 in the numerator and the denominator and having the 13th-order accuracy.
Findings. The new modification of split-step Fourier method and its parallel implementation based on the Intel MKL library are represented. In this paper we demonstrate the high accuracy of the proposed numerical algorithm on the solution of scalar NSE in the form of fundamental soliton. The possibility of efficient parallelization of the algorithm for computing systems with shared memory was demonstrated by the figures with the acceleration of the algorithm on different numbers of threads.
Originality/value. The proposed new numerical algorithm can effectively solve the problem of propagation of light pulses in multi-core and multi-mode optical fibers. The possibility of creating its parallel implementation allows using this algorithm on a large spectrum of computing systems with shared memory.
Keywords: split-step Fourier method, Pade approximation, nonlinear Schrodinger equation, multi-core fibers, multi-mode fibers, nonlinear fiber optics.
Acknowledgements. This research was supported by Russian Science Foundation (Agreement No. 14-21-00110).
Received 6 April 2015 Received in revised form 15 April 2015
