Full Text
Алгоритмы решения обратной задачи рассеяния для модели Манакова [1]
1. ВВЕДЕНИЕ
Векторный вариант нелинейного уравнения Шрёдингера, учитывающий одновременно эффекты дисперсии, нелинейности и поляризации оптического излучения, известен как модель Манакова [1]. Эта модель состоит из пары связанных нелинейных уравнений Шрёдингера с дисперсией второго порядка и кубической (керровской) нелинейностью для двух оптических поляризаций. Исследование поляризационных эффектов, осложненных дисперсией и нелинейностью, имеет большое значение для современной нелинейной оптики и оптических технологий. Модель Манакова находит важные применения в нелинейной физической оптике для описания явлений оптического просветления и самофокусировки. В последние годы эта модель особенно востребована для изучения нелинейно-дисперсионных и поляризационных эффектов, возникающих при распространении оптического излучения по волоконным линиям связи [2].
Векторное нелинейное уравнение Шрёдингера (НУШ), называемое моделью Манакова, имеет следующий вид:
(1)
где – мнимая единица, – двухкомпонентный вектор решения, содержащий две нормированные поляризационные компоненты поля и , – пространственная координата, – время. Параметр принимает значения +1 или 1, для дефокусирующего и фокусирующего случаев, соответственно. Векторное уравнение Шрёдингера (1) принадлежит к нетривиальному классу интегрируемых нелинейных уравнений, изучаемых методом обратной задачи рассеяния (ОЗР) [3, 4]. В рамках метода ОЗР решение нелинейного уравнения сводится к исследованию линейных спектральных задач рассеяния (прямых и обратных), ассоциированных с исходным нелинейным уравнением.
Метод ОЗР – выдающееся достижение современной математической физики, позволившее аналитически исследовать целый ряд нелинейных эволюционных волновых уравнений (см. [4]), оказался мощным инструментом решения задачи Коши для этих нелинейных уравнений. Численная реализация метода ОЗР на основе алгоритмов решения обратных задач рассеяния позволяет развивать эффективные алгоритмы решения нелинейной задачи Коши, без каких-либо итераций.
Сущность метода ОЗР состоит в том, что нелинейному эволюционному уравнению ставится в соответствие система линейных уравнений, где в качестве коэффициентов используются решения нелинейного уравнения. Задачи, формулируемые как обратные задачи рассеяния для такой системы, позволяют находить решения нелинейного уравнения, при произвольном значении эволюционного параметра . В случае модели Манакова эта линейная система называется системой Манакова и имеет следующий вид:
(2)
где – комплексное собственное значение, – трехмерный собственный вектор, а символ * означает комплексное сопряжение. В дефокусирующем случае – спектр собственных значений системы Манакова содержит только непрерывную компоненту, соответствующую диспергирующим волновым пакетам. В более сложном фокусирующем случае к непрерывному спектру добавляется дискретный, соответствующий солитонным решениям модели Манакова. Определение (спектральных) данных рассеяния по заданному решению представляет собой прямую задачу рассеяния, а нахождение решения по данным рассеяния, при фиксированном эволюционном параметре , является обратной задачей рассеяния. Для достаточно быстро убывающих на бесконечности решений модели Манакова, эволюция непрерывной и дискретной компонент данных рассеяния сводится к фазовым сдвигам [1, 3, 4].
Таким образом, решение задачи Коши дается тремя процедурами.
- По заданному начальному условию, путем решения прямой задачи рассеяния, находятся начальные спектральные данные задачи.
- Умножение спектральных компонент на их фазовые эволюционные множители дает конечный спектр задачи.
- Решение обратной задачи рассеяния для полученных спектральных данных дает решение задачи Коши для интегрируемого нелинейного эволюционного уравнения.
Описанная схема решения напоминает применение метода Фурье для решения линейных дифференциальных уравнений и часто называется методом нелинейного преобразования Фурье.
Данная работа посвящена наиболее сложному этапу схемы решения – обратной задаче рассеяния, алгоритмам ее решения для модели Манакова.
Метод ОЗР для нелинейных уравнений с локализованными (убывающими на бесконечности) решениями выводит из системы Манакова систему из девяти линейных интегральных уравнений Гельфанда – Левитана – Марченко (ГЛМУ) [5, 6]. Для решения обратной задачи рассеяния достаточно трех интегральных уравнений системы ГЛМУ, которые (для левой задачи рассеяния и случая конечного носителя – ) для системы Манакова запишем в следующем виде [7]:
(3)
где , – ядра интегральных уравнений, . Решение ГЛМУ дает синтезирующее соотношение: .
Для случая более простого скалярного нелинейного уравнения Шрёдингера, не учитывающего поляризацию, ассоциированного с системой Захарова – Шабата [3], в нашей лаборатории на основе ГЛМУ были разработаны эффективные вычислительные алгоритмы тёплицева внутреннего окаймления (Toeplitz Inner Bordering – TIB) [8, 9]. Эти алгоритмы основаны на численном решении сеточного дискретного аналога ГЛМУ с помощью алгоритма окаймления типа Левинсона [10]. Их численная эффективность обусловлена использованием тёплицевой симметрии дискретизованных ГЛМУ. В алгоритме решения обратной задачи второго порядка точности используется оригинальный прием перенос ряда слагаемых, содержащих неизвестные, в правую часть, восстанавливающий тёплицеву структуру матрицы. Алгоритмы TIB были успешно применены для расчета волоконных брэгговских решеток [11, 12] и решения обратных задач для уравнения Гельмгольца [13]. Эти алгоритмы открывают перспективы для разработки новых подходов к передаче информации по оптическим сетям, основанных на интегрируемости НУШ [14, 15].
Мы не касаемся здесь более простых прямых задач рассеяния. Основной целью работы является сравнение двух новых блочных алгоритмов, один из которых является дальнейшим развитием скалярного алгоритма TIB, а другой основан на тёплицевом разложении матрицы и алгоритме окаймления Тыртышникова [16], для решения обратных задач рассеяния, ассоциированных с моделью Манакова.
Сначала мы рассмотрим алгоритм первого порядка точности аппроксимации интегралов ГЛМ, затем будут детально представлены два алгоритма второго порядка точности, где нарушена тёплицева структура дискретизованной системы уравнений. В заключительной части работы на примере точного аналитического решения – векторного солитона Манакова, представлено сравнение точности и скорости этих двух алгоритмов.
2. АЛГОРИТМ ПЕРВОГО ПОРЯДКА ТОЧНОСТИ АППРОКСИМАЦИИ
Опишем сначала общий подход первого порядка точности, основанный на блочной форме дискретизованных ГЛМУ. Введем дискретную вычислительную сетку путем разбиения рабочего интервала на равных частей:
(4)
Здесь – шаг сетки. Сеточные функции вводятся посредством аппроксимации интегралов в ГЛМУ методом прямоугольников:
Дискретизация интегральных ГЛМУ (3) с первым порядком точности приводит к набору вложенных систем линейных алгебраических уравнений (СЛАУ). Так как ядра ГЛМУ (3) зависят от разности аргументов, матрицы СЛАУ состоят из девяти квадратных тёплицевых блоков. Несмотря на то, что все эти блоки тёплицевы, матрицы системы, в отличие от скалярного случая, не обладают тёплицевой симметрией. Для фиксированной координаты дискретизованные ГЛМУ (3) представляют собой систему линейных уравнений размерности . При изменении размера системы от 1 до получаются системы, вложенные одна в другую, т.е. каждая система является окаймлением системы на единицу меньшей размерности. Для численного решения ОЗР необходимо решить все полученные вложенные системы и взять последние элементы решения: . Метод исключения Гаусса для вложенных систем ГЛМУ требует арифметических операций, что может потребовать применения суперкомпьютеров для решения задачи, при реальных размерах задачи в несколько тысяч отсчетов. Однако в случае первого порядка точности матрица системы дискретизованных ГЛМУ может быть представлена в виде блочной матрицы, имеющей тёплицеву симметрию, содержащую блоки матрицы размером . Блочным алгоритмам окаймления типа Левинсона [10] требуется всего операций для решения всех вложенных систем с блочно-тёплицевой матрицей. Для решения тёплицевых систем разработаны эффективные “супербыстрые” алгоритмы (см. [17]), требующие всего арифметических операций, не использующие идею окаймления. Однако для решения вложенных друг в друга блочно-тёплицевых СЛАУ алгоритмы окаймления оказываются в итоге более быстрыми, поскольку решают каждую из вложенных систем всего за арифметических операций.
Запишем систему уравнений 1-го порядка точности в следующем матричном виде:
(5)
где использованы следующие обозначения:
(6)
(7)
Здесь и далее символ обозначает блочное транспонирование: трехкомпонентные вектор-столбцы , переставляются, оставаясь столбцами, квадратные блоки матрицы также переставляются без внутренних изменений.
Блочно-тёплицева матрица содержит блоки и , размером , и их эрмитовы сопряжения:
(8)
где символ † означает эрмитово сопряжение, а также введены следующие обозначения:
(9)
(10)
Далее в этом и следующем разделах будем опускать верхний индекс для всех матриц и векторов. Оставим его только для матрицы системы , вектора решения , верхней и нижней строк обратной матрицы.
Для блочно-тёплицевой матрицы существуют известные алгоритмы обращения, основанные на блочных вариантах алгоритма Левинсона, см. например, [18]. При этом, для обращения тёплицевой матрицы достаточно рекуррентно найти верхнюю и нижнюю блочные строки матрицы, обратной тёплицевой. Эти строки являются основными вспомогательными массивами в алгоритмах окаймления типа Левинсона. Обозначим элементы этих строк на шаге блочного алгоритма следующим образом:
, . (11)
Удобно ввести следующие матрицы, размером , для более компактной записи алгоритма: которые определяются следующими суммами:
(12)
Эти суммы, по существу, составляют основной цикл алгоритма, выполняемый для каждого от 1 до . Основным шагом алгоритма окаймления является вычисление первой и последней блочных строк обратной матрицы на основе их значений , на предыдущем шаге. Алгоритмы типа Левинсона находят новые строки как линейную комбинацию строк предыдущего шага:
(13)
где 0 блок из нулей, размером . Умножив первую блочную строку обратной матрицы (13) на первый и последний блочный столбцы матрицы , получим систему уравнений на коэффициенты
Здесь единичный блок размером . Решение этой системы имеет вид:
(14)
Проделав то же самое с последней блочной строкой обратной матрицы, получим выражения для :
(15)
3. АЛГОРИТМ А ВТОРОГО ПОРЯДКА ТОЧНОСТИ
Второй порядок точности нарушает тёплицеву структуру блочной матрицы системы ГЛМУ, так как формула трапеций приводит к изменению первого и последнего блочных столбцов матрицы (а также блоков ее главной диагонали). Нарушение тёплицевой структуры матрицы формально не позволяет использовать для решения блочные алгоритмы окаймления типа Левинсона. Идея решения “испорченной” системы уравнений была предложена в работе авторов [9]. Она опирается на перенос части слагаемых системы, содержащих неизвестные, в правую часть СЛАУ. Убирая из системы уравнение для , вновь приходим к системе вложенных матриц размером , . После процедуры переноса ряда слагаемых систем в правую часть СЛАУ, матрицы систем снова становятся тёплицевыми, как в (8), однако в правой части появляются неизвестные. Сделав замену , можно исключить из системы первое блочное уравнение и вновь прийти, как в (5), к последовательности вложенных систем с блочно-тёплицевыми матрицами размером , . В результате преобразования диагональные блоки матрицы системы и вектора-блоки правой части принимают вид:
В правой части системы появились первый и два последних неизвестных элемента решения , , которые, собственно, и требуется найти. Умножив первую и последнюю строки обратной матрицы на правую часть системы ГЛМУ, мы получаем пару уравнений для пары неизвестных – для первого и последнего элементов блочного вектора неизвестных:
Из шести уравнений этой системы мы выделили более простую подсистему всего из трех уравнений, содержащую компоненты , , вектора решения:
(16)
Решение системы (16) находится по формулам Крамера, что практически не влияет на оценку числа арифметических операций, которая составляет .
3.1. Схема блочного алгоритма а
Для вычисляем начальное значение для вспомогательного массива: . Затем вычисляем элементы матриц и , размером : , а также находим матрицы и : (находим только те элементы матриц, которые используются в алгоритме).
Находим m-ю компоненту вектора решения , , используя решение системы (16) (это результат работы алгоритма на каждом шаге).
Вычисляем блочные коэффициенты , используя формулы (14) и (15).
Определяем следующие значения вспомогательных массивов, т.е. строк и , используя (13).
Находим необходимые элементы новых матриц , вычислив суммы (12).
Увеличиваем на единицу и переходим к шагу 2, пока .
4. АЛГОРИТМ B ВТОРОГО ПОРЯДКА ТОЧНОСТИ
Алгоритм B также аппроксимирует систему уравнений ГЛМ (3) системой линейных алгебраических уравнений, при помощи введения сеточных функций и замены интегралов конечными суммами. При этом используется та же численная сетка (4). В этом алгоритме второго порядка сеточные функции вводятся посредством аппроксимации интегралов в ГЛМУ методом трапеций, но немного по иному, чем в алгоритме А. Учитывая, что подынтегральные выражения состоят из двух множителей: известных ядер , и вектора неизвестных функций , множитель формулы трапеций можно приписывать не элементам матриц, а соответствующим элементам искомых сеточных векторов, которые определены следующим образом:
Такой выбор сеточных функций приводит к системе вложенных друг в друга СЛАУ. Каждая из этих систем имеет сравнительно простую структуру матрицы. Перестановкой строк и столбцов эту матрицу можно преобразовать к матрице, близкой к блочно-тёплицевой, состоящей из ( блоков размерности :
где по-прежнему обозначено
Тёплицеву структуру блочной матрицы нарушают лишь левый верхний и правый нижний блоки, которые имеют следующий вид:
Система уравнений по-прежнему имеет вид (5), но теперь включает и уравнение для индекса 0. Искомый вектор и вектор правой части имеют вид
При решение записывается в виде:
(17)
Как было предложено в работе [19], матрицу СЛАУ размерности можно представить в виде
(18)
4.1. Алгоритм Тыртышникова обращения матриц, близких к тёплицевым
Применительно к матрицам, близким к тёплицевым, Тыртышников [16] разработал общий алгоритм обращения, основанный на тёплицевом разложении – представлении такой матрицы в виде суммы парных произведений левых и правых треугольных тёплицевых матриц [20]. Напомним, что тёплицево разложение не единственно и может иметь разное число слагаемых (длину разложения). Наименьшая возможная для данной матрицы длина называется тёплицевым рангом матрицы. Этот ранг не превышает двух для тёплицевой матрицы, а для произвольной матрицы может достигать ее размерности. Нахождение тёплицева разложения является в общем случае нетривиальной задачей, но если оно найдено, решение СЛАУ может быть получено по рекуррентным формулам, основанным на методе окаймления. Алгоритм не требует коммутативности операции умножения, что позволяет обобщить его на блочные матрицы. Алгоритм B, предложенный в работе [16], содержит арифметических операций, где – длина тёплицева разложения, – размер блока (в нашем случае ).
В работе [19] алгоритм Тыртышникова был применен для решения дискретизованной задачи рассеяния с семейством матриц (18). Тёплицево разложение для первых трех слагаемых очевидно, так как они являются тёплицевыми треугольными. Четвертое слагаемое – матрица, содержащая всего два ненулевых элемента, порождает еще три члена тёплицева разложения:
(19)
Здесь обозначает левую треугольную блочно-тёплицеву матрицу размерности , составленную из элементов вектора по правилу: первый столбец матрицы совпадает с второй столбец – с , сдвинутым вниз на один элемент, третий столбец – на два элемента и т.д. Освободившиеся места сверху заполняются нулями.
Алгоритм обращения матрицы [16] основан на методе окаймления. Он последовательно решает усеченные подсистемы, в которых матрицы являются ведущими подматрицами исходной матрицы СЛАУ, а вектор правой части состоит из соответствующего числа первых элементов исходного. При решении обратной задачи рассеяния с первым порядком аппроксимации последовательность систем уравнений при различных значениях параметра совпадает с последовательностью вложенных подсистем алгоритма обращения тёплицевой матрицы. Таким образом, обратная задача рассеяния сводится к однократному обращению матрицы, при , методом окаймления. В случае второго порядка аппроксимации, как уже отмечалось, теплицева структура матриц нарушается. Как видно из (19), все векторы , , за исключением могут быть получены окаймлением , , начиная с . Алгоритм окаймления второго порядка точности (алгоритм B) для решения этой проблемы рассматривает вспомогательную задачу, которая отличается от исходной тем, что в тёплицевом разложении (19) . Тогда вспомогательная задача на шаге будет окаймлением вспомогательной задачи на шаге : , . Исходная задача на шаге также будет окаймлением вспомогательной задачи на шаге , но другого вида: , , , . Вспомогательная и исходная задачи должны решаться совместно, что приводит к небольшому объему дополнительных вычислений, не меняя асимптотической оценки числа операций. Основные шаги (схема) алгоритма приведены ниже.
4.2. Cхема блочного алгоритма B
1. Решение при дается формулами (17).
2. При требуется инициализация процесса, поскольку матрица для не является окаймлением матриц для :
3. Завершение вычислений при , а также все последующие шаги при определяются общими формулами:
Здесь все величины представляют собой блоки , или вектор столбцы, состоящие из таких блоков, за исключением , которые состоят из блоков размером . Стрелка вниз означает сдвиг элементов вектора вниз с добавлением нулевого блока в качестве первого элемента, стрелка вверх – аналогичный сдвиг в противоположном направлении, длина вектора при этом увеличивается на единицу. Вектор дает решение обратной задачи рассеяния для точки .
5. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ
Для проверки алгоритма использовалось точное аналитическое решение модели Манакова – векторный солитон Манакова. Солитон, как решение уравнений ГЛМ, для комплексного собственного значения дискретного спектра , имеет следующий вид [1, 7]:
где задает положение центра солитона, – фазы компонент солитона, – угол поляризации. Cолитон Манакова является решением ГЛМУ для следующего ядра:
(20)
Обратная задача восстановления солитона Манакова по ядру ГЛМ решалась на интервале при следующем наборе параметров солитона: Численное моделирование для алгоритмов A и B показало, что увеличение вдвое размера сетки уменьшает ошибку расчета в 4 раза, что подтверждает второй порядок точности расчетов, при этом время расчета увеличивается в 4 раза, что означает квадратичную зависимость числа арифметических операций алгоритмов от размера задачи .
Результаты численного моделирования представлены на фиг.1.
Фиг. 1.а график зависимости ошибки восстановления векторного солитона; б график зависимости времени расчета восстановления векторного солитона.Квадратами показаны расчетные значения для алгоритма A, кружками для алгоритма B.
На фиг.1а показана зависимость натурального логарифма ошибки восстановления векторного солитона Манакова от количества точек расчета в логарифмическом (по основанию 2) масштабе. Прямые линии соответствуют подгонке функцией Совпадение расчета с подгонкой подтверждает второй порядок точности обоих алгоритмов. При этом алгоритм B оказался примерно 1.8 раза точнее алгоритма A, при прочих равных условиях. Однако алгоритм A решает обратную задачу рассеяния почти в 2.2 раза быстрее алгоритма B. В подтверждение этого на фиг. 1б представлена зависимость натурального логарифма времени расчета (в секундах) обратной задачи рассеяния в зависимости от двоичного логарифма числа расчетных точек. Сплошными линиями показана подгонка функцией . Расчеты проводились на компьютере под управлением ОС Windows 7, 64 бит, с процессором Intel Core i7 Q720.
Для трансляции программы использовался Фортран: GNU Fortran (Rev1, Built by MSYS2 project) 12.2.0, с системными библиотеками MINGW, 64 бит.
6. ЗАКЛЮЧЕНИЕ
Рассмотрены два численных алгоритма решения обратной задачи рассеяния для модели Манакова векторного нелинейного уравнения Шрёдингера. Оба алгоритма опираются на блочные варианты дискретизованных интегральных уравнений Гельфанда-Левитана-Марченко. Решение задачи для первого порядка точности аппроксимации сводится к обращению ряда вложенных друг в друга блочных тёплицевых матриц блочным вариантом метода окаймления Левинсона. Повышение точности аппроксимации нарушает тёплицеву структуру блочных матриц. Проблема построения алгоритма повышенного (второго) порядка точности по-разному решена в представленных алгоритмах. Алгоритм А использует восстановление тёплицевой структуры матриц системы, путем переноса ряда слагаемых СЛАУ в правую часть. Затем в алгоритме А используется блочный вариант алгоритма окаймления Левинсона. В алгоритме B применяется тёплицево разложение для блочной матрицы, близкой к тёплицевой и алгоритм окаймления Тыртышникова.
Представленные алгоритмы соответствуют расчетному интервалу координат . Однако оба алгоритма инвариантны относительно выбора рабочего интервала и начала координат, как и большинство уравнений и алгоритмов математической физики. При изменении интервала меняется только расчетная сетка, но сам алгоритм не меняется. Численные тесты, проведенные путем сравнения результатов расчетов с известными точными аналитическими решениями (векторный солитон Манакова, секанс-потенциал), подтвердили устойчивость и второй порядок точности предложенных алгоритмов.
На примере векторного солитона Манакова приводятся результаты сравнения скорости и точности расчетов представленных алгоритмов. Алгоритм B оказался 1.8 раза точнее алгоритма A, при этом алгоритм A решает обратную задачу рассеяния почти в 2.2 раза быстрее алгоритма B. Алгоритм A значительно проще для кодирования, чем алгоритм B, однако алгоритм B имеет перспективы дальнейшего обобщения на еще более высокие порядки точности аппроксимации интегральных уравнений ГЛМ.
[1] Работа выполнена при финансовой поддержке РНФ (грант 22-22-00653).