О разностном методе нахождения соб ственных мод нелинейного уравнения Шредингера Текст научной статьи по специальности «Математика»

3. Banks H. Т., К а р р е 1 F., Wang С. Weak solutions and differentiability for size structured population models // Int. Ser. Numer. Math. 1991. 100. P. 35-50.

4. Ackleh A. S., De n g K. Monotone method for first order nonlocal hyperbolic initial-boundary value problems // Applic. Analys. 1997. 67. P. 283-293.

5. Sinko J. W., Streifer W. A new model for age-sized structure for a population // Ecology. 1967. 48. P. 910-918.

6. Денисов A. M.,Макеев A. С. Итерационные методы решения обратной задачи для одной модели популяции // ЖВМиМФ. 2004. 44. № 8. С. 1480-1489.

7. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М: Наука, 1974.

Поступила в редакцию 01.03.05

УДК 517.958:535.14

С. А. Варенцова, В. А. Трофимов

О РАЗНОСТНОМ МЕТОДЕ НАХОЖДЕНИЯ СОБСТВЕННЫХ МОД НЕЛИНЕЙНОГО УРАВНЕНИЯ ШРЕДИНГЕРА1

(кафедра вычислительных методов факультета ВМиК)

1. Введение. Нахождение солитонных режимов распространения световых импульсов и волно-водного распространения пучков представляет собой проблему, важную для многих приложений, что объясняет постоянный интерес к данному классу задач [1-11]. Как известно, солитонное решение (решение в виде солитона) соответствующего нелинейного уравнения оптики обладает тем свойством, что оно не изменяется либо во времени (сохраняется пространственное распределение интенсивности световой волны), либо вдоль координаты, совпадающей с направлением распространения оптического излучения в нелинейной среде. Если эта среда — оптическое волокно, то солитон представляет собой импульс неизменной во времени формы. В этом случае наиболее существенным приложением солитонов являются задачи передачи информации оптическими методами.

Следует подчеркнуть, что построение эффективных методов нахождения солитонных решений соответствующих уравнений до сих пор остается актуальной задачей, так как общих методов решения не существует. Можно лишь упомянуть метод обратной задачи рассеяния [1], который позволяет вычислить спектр нелинейного уравнения Шредингера (НУШ) для случая слабой нелинейности, однако построить форму солитона на его основе в общем случае очень сложно. Еще один его существенный недостаток заключается в невозможности обобщения на многомерные задачи. Других конструктивных и достаточно универсальных методов построения финитных солитонов по пространству (времени) с ограниченной энергией, на наш взгляд, нет. Именно по этой причине наиболее известным в литературе является солитон, имеющий форму ch-2(i), или "темный" солитон. Финитные солитоны более высокого порядка (более сложной формы) практически неизвестны. Естественно, тривиальный случай, представляющий собой композицию достаточно удаленных друг от друга солитонов вида ch~ (i), не обсуждается.

Отметим, что в отличие от большинства упомянутых выше работ нами рассматривается задача нахождения солитонных решений НУШ на конечном отрезке. Между тем с математической точки зрения солитоны представляют собой собственные функции соответствующего НУШ с заданными краевыми условиями. Поэтому очевидной является попытка обобщить имеющиеся методы нахождения собственных функций и собственных значений на данный класс задач. Однако оказалось, что непосредственный перенос рекомендаций по нахождению собственных значений и функций линейных

1 Работа выполнена при поддержке РФФИ (проект № 02-1-727).

задач с помощью компьютерного моделирования неприемлем. В частности, для линейных задач на собственные значения, заданных на бесконечном (полубесконечном) интервале, известна рекомендация об увеличении длины отрезка, на котором решается задача, для более точного нахождения собственного значения. Важно подчеркнуть, что для НУШ уже при малых значениях коэффициента нелинейности этот способ не дает положительного результата. Поэтому цель настоящей работы заключается в выявлении закономерностей построения собственных мод любого порядка НУШ в зависимости от коэффициента, характеризующего нелинейность.

Следует также подчеркнуть имеющееся различие между собственными модами НУШ и солитона-ми. Так, найденные собственные функции НУШ в общем случае могут и не быть соответствующими его солитонными решениями. Причина этого заключается в выбранной нами постановке задачи: ищутся финитные решения на конечном отрезке времени или пространственной координаты с нулевыми краевыми условиями. При этом не требуется выполнения равенства нулю первой производной по этой координате (или по времени) в граничных точках. Как следствие этого при подстановке собственной моды в качестве начального условия в НУШ она может изменяться вдоль координаты распространения, если увеличить (или уменьшить) область по поперечной координате по сравнению с ее длиной, на которой найден солитон. Следовательно, чтобы собственная функция НУШ была бы его солитонным решением для любой области по поперечной координате (или по времени), необходимо выполнение условия равенства нулю ее производной совместно с обращением в нуль самой функции на границах области, в пределах которой находится собственная функция. Однако в настоящей работе такая постановка задачи не обсуждается, так как основное внимание уделено закономерностям нахождения собственных функций НУШ. Тем не менее полученные результаты имеют физическое приложение к задачам нахождения собственных мод нелинейных плоских волноводов, которые широко используются при передаче информации оптическими методами.

2. Постановка задачи. Как известно, процесс распространения оптического излучения в среде с кубичной нелинейностью описывается следующим безразмерным нелинейным уравнением Шредин-гера:

ЗА ,д2А . ,„|2

иг ох1

с начальными и граничными условиями

- + гО-^-т + га\А\гА = й, 2>0, 0 < ж < Ь, (1)

г=0

= А0(ж), А =0. (2)

х=0,Ь

Здесь А(х, г) — комплексная амплитуда импульса, нормированная на максимальное значение, х — безразмерная поперечная координата, г — координата, вдоль которой распространяется импульс, Ь — характерный поперечный размер области, коэффициент И пропорционален коэффициенту дифракции светового пучка. Ниже для определенности будем считать его равным единице: И=1. Коэффициент а характеризует отношение начальной мощности импульса к характерной мощности самовоздействия. Отметим, что положительное значение а соответствует самофокусировке пучка (компрессии, или сжатию импульса), а его отрицательное значение — дефокусировке пучка (декомпрессиии, или увеличению длительности импульса).

В соответствии с целью данной работы будем искать решение задачи (1), (2) в следующем виде:

А(х,г) = ф(х)е~гХ*. (3)

Подставляя функцию (3) в уравнение (1), получаем нелинейную задачу на собственные значения:

в2ф(х йх2

+ а\ф(х)\2ф(х) = Хф(х), 0 < ж < Ь, ф(0) = ф(Ь) = 0. (4)

Покажем действительность параметра А и собственных функций ф(х). Для этого умножим уравнение (4) на ф*(х), сопряженное к нему — на ф(х), сложим их и проинтегрируем по ж от 0 до Ь. Получим следующее интегральное тождество:

ь ь

В силу действительности левой части тождества (5) и коэффициента при А в его правой части параметр А — также действительное число. Следовательно, (4) — дифференциальное уравнение с действительными коэффициентами, а ф(х) — действительная функция. Поэтому далее знак модуля у функции ф будет опущен.

Нетрудно видеть, что при а = 0 задача (4) является линейной и имеет известное решение:

у = , ¿ = 1,2,.... (6)

Отметим также, что из непрерывной зависимости собственных значений А& от коэффициентов уравнения (4) А^ = Afc(a) (см., например, [12]) следует, что существует интервал (0, а*), в котором все собственные значения задачи (4) будут отрицательны и близки к решению (6) линейной задачи.

В заключение данного пункта заметим, что во многих работах (см., например, [10]), посвященных нахождению солитонных решений уравнения (1), вместо (3) предлагается искать солитоны в виде

A(x,z) = = ф{х)е~1Х*,

а

что делает возможным фиксировать значение коэффициента перед нелинейным слагаемым в уравнении (4). На наш взгляд, такая замена не всегда целесообразна, поскольку она не позволяет изучать и использовать свойства собственных значений и векторов задачи (4), например при малых значениях параметра а.

3. Разностная схема. Задачу (4) будем решать численно. Для этого на отрезке [О, Ь] введем равномерную сетку

иь = = г = 0,..., N + 1, 1г = Ь/Щ и рассмотрим разностную схему [13] с итерационным процессом вида

^й + «№?)Ч?+1 = А-+г = 1,..., ЛГ, (7)

Фо+1 = Ф^+г = 0, * = 0,1,....

Нетрудно показать, что схема (7) имеет второй порядок аппроксимации. На каждой итерации в разностная задача (7) является линейной задачей на собственные значения

Ь 8ф8+1 = А 8+1 ф'+1 (8)

с трехдиагональной симметричной матрицей оператора 1/, который определяется следующим образом:

I/ = А + а1(ф8)2.

Здесь Ау1 = ухх,г-, I — единичный оператор, (ф(, ..., ф%-) = фв■ Решение задачи (8) ищется методом бисекции, а соответствующие собственные векторы вычисляются методом прогонки.

Начальное приближение, соответствующее в = 0, строится специальным образом, в зависимости от значения параметра а. В случае слабой нелинейности а ^ 1 в качестве начального приближения для вычисления А&, фь целесообразно выбирать соответствующую собственную функцию линейной (а = 0) разностной задачи (7):

С = вт . (9)

Однако этот случай малоинтересен как для практики, так и с точки зрения вычислительной математики из-за того, что при малых значениях а задача (7) является квазилинейной: нелинейность оказывает слабое влияние на ее решение, и его можно получить аналитически, например методом обратной задачи рассеяния. Нас будет интересовать случай а 1.

Окончание итерационного процесса происходит при выполнении условия:

Ш+1 - Фк\\с < ЖНс+Ъ, £1>£2 > Существенно, что при компьютерном моделировании контролировалась также абсолютная погрешность решения в норме С:

Ф* =

d Фк(х) + ^\фк{х)\2 фк{х) - Акфк{х)

dx2

Задача (4) нелинейна по собственной функции, поэтому ее спектр зависит от выбранной нормировки собственной функции. Так как в исходном уравнении Шредингера (1) амплитуда А(х,г) уже нормирована на максимальное значение, поэтому будем решать разностную задачу (7), фиксируя

тах|^(ж)| = 1.

4. Результаты расчетов. Проведенные расчеты показали, что для всех рассмотренных значений длины отрезка Ь и параметра нелинейности а задача (4) обладает дискретным спектром. При малых значениях коэффициента нелинейности а ^ 0,1 для значений длины отрезка Ь = 0,5 -т- 10 вычисленные собственные числа нелинейной разностной задачи (7) являются отрицательными и близки к собственным числам линейной (а = 0) разностной задачи. Построенный итерационный процесс сходится для всех рассмотренных номеров собственных чисел: к ^ 50. Сравнение нескольких первых собственных чисел на вложенных сетках показало, что они (А^) вычисляются со вторым порядком точности. Так как первую собственную функцию можно записать аналитически (см., например, [1-3, 9]), то она использовалась нами для тестирования разностной схемы.

Для практики расчетов собственных функций и собственных чисел представляет большой интерес изучение влияния длины отрезка Ь на характер сходимости итерационного процесса и вид спектра. Компьютерные эксперименты показывают, что при Ь = 0,5 -т- 2,0 итерационный процесс (7) сходится для всех а ^ 100 и номеров к = 1 ^ 50 собственных значений и функций, например при к = 0,01. Важно подчеркнуть, что собственные значения А^ с ростом параметра а также растут, при этом значение А1 стремится к величине а/2. Так как собственные числа зависят от длины выбранного отрезка, то в таблице для различных Ь приведены значения акр, при превышении которых для первого собственного значения выполнено условие |А1 — а/2\ ^ Ю-2. Таблица наглядно иллюстрирует, что с ростом длины отрезка Ь значение акр, при котором первое собственное значение "склеивается" с а/2, быстро уменьшается.

Зависимость акр от длины отрезка Ь при условии А1 ~ а/2

ь 2 5 10 20

«кр 36 6 0,9 0,2

Для отрезков с длиной Ь > 2 собственные значения А& с номерами к = 2 10 при увеличении коэффициента нелинейности а, так же как и Х\, стремятся к величине а/2, образуя группу близко расположенных (возможно, кратных) собственных чисел. Это обстоятельство вызывает существенные трудности при вычислении собственных значений и функций с номерами к > 1. Так, в условиях близости нескольких А& использование упомянутого выше точного к-го решения (9) разностной задачи (7) при а = 0 в качестве начального приближения уже не дает сходимости к нужному собственному значению и собственной функции. Поэтому предлагается использовать другой подход, при котором начальное приближение, обеспечивающее сходимость, конструируется специальным образом. Так, на отрезке меньшей длины, где не имеет место "склеивание" собственных значений с разными номерами (например, при Ь = 1), вычисляется первая собственная функция ф\ для достаточно большого значения а, что обеспечивает нужную "узкую" форму полученной функции. Затем она несколько раз воспроизводится на отрезке большей длины, имитируя соответствующую собственную функцию с нужным числом осцилляций. Сконструированная подобным образом функция использовалась в качестве начального приближения для вычисления А&, ф^, на отрезках длины Ь = 2, 5, 10, 20 для тех значений а, при которых А& близки к а/2.

Отметим, что для сконструированного начального приближения итерационный процесс сходился, если выполнялось условие "близости" начального приближения к решению:

||Ф1 ~ Фк\\с ^ г,

где величина е = е(Ь,к). Так, вычислительные эксперименты показали, что, например, для второго собственного значения к = 2 при Ь = 5 величина е должна удовлетворять условию е ^ 0,05, в то время как для Ь ^ 10 требуется уже выполнение условия е ^ 0,0001. Таким образом, если первые к собственных значений образуют группу близко расположенных чисел, при наличии хорошего начального приближения возможно их разделение. К сожалению, нам не во всех случаях удавалось построить такое приближение.

В качестве примера на рис. 1 представлено построенное для отрезка L = 5 и а = 25 начальное приближение ф®, которое сначала вычислялось при L = 1, а = 125, и вторая собственная функция Ф2, полученная с помощью данного начального приближения. Важно подчеркнуть, что выбор а (при фиксированном размере области начального приближения L = 1) обусловлен длиной рассматриваемой области и номером собственной функции. В этом случае можно оценить ширину "отдельного" начального приближения и соответственно выбрать а. Для сравнения на рис. 1 пунктиром изображена функция ф(х) = sin (^ip), при использовании которой в качестве начального приближения сходимость итерационного процесса отсутствовала.

Рис. 1. Собственная функция ф2(х) (сплошная кривая), вычисленная на отрезке Ь = 5 для а = 25; построенное на основе предлагаемого подхода начальное приближение (штрих-пунктир); точное решение разностной задачи

(7) при а = 0 (пунктир)

Отметим также, что для отрезков Ь = 10, 20 имеют место трудности уже при вычислении первого собственного значения А1 в силу того, что первые собственные значения "склеиваются" в группу близко расположенных собственных чисел при значительно меньших значениях а. Однако сходимость итерационного процесса можно улучшить, увеличивая шаг сетки к. Так, для Ь = 10 при к = 0,001 итерационный процесс не сходится при а > 50. Увеличение шага сетки до к = 0,01 позволяет вычислить Х\, ф\ вплоть до а = 70, а при дальнейшем увеличении к до 0,05 сходимость имеет место даже при а = 84. В случае больших значений отрезка Ь, например при Ь = 20, к = 0,002, сходимость отсутствует при а > 10. При этом увеличение шага сетки до значения к = 0,02 позволяет вычислять Х\, ф\ до значений а = 20. Такая зависимость связана с известным фактом: с ростом шага сетки уменьшается число воспроизводимых на сетке гармоник. Это в свою очередь уменьшает максимальный номер собственной функции разностной задачи, что упрощает проблему отделения собственных чисел друг от друга.

Что касается собственных чисел и векторов с номерами к ^ 10, то при их нахождении вычислительных трудностей не наблюдается для достаточно больших интервалов по параметру а и при специальном выборе начального приближения, а также шага сетки. Так, для Ь = 10 значения А&, фь с номерами к = 20, 30, 40, 50 устойчиво вычислялись при а ^ 100. Десятое собственное значение и собственная функция Аю, ^ю вычислялись вплоть до а = 40. При увеличении Ь до 20 значения А&, фк с номерами к = 30, 40, 50 вычислялись вплоть до а = 40. В качестве примера на рис. 2,а приведена функция -0Ю, вычисленная на отрезке Ь = 10 при а = 20, а на рис. 2,6 — Ф15, вычисленная соответственно при Ь = 20, а = 10. Данный рисунок наглядно показывает, что собственные функции существенно отличаются от набора последовательных функций вида сЬ_2(ж).

Отметим также, что во всех проведенных численных экспериментах абсолютная погрешность в норме С при подстановке в уравнение (4) составляла в зависимости от шага к величину, не превосходящую Ю-4 для далеко отстоящих друг от друга собственных значений и 0,05 — для собственных чисел, близких к а/2.

5. Заключение. Приведенные в работе результаты позволяют сделать следующие выводы и дать некоторые рекомендации по численному нахождению собственных функций и собственных значений НУШ для практически важного случая а > 1 (в том числе а 1).

Во-первых, в отличие от линейного уравнения Шредингера, рассматриваемого на бесконечном (полубесконечном) интервале, при решении которого разностным методом необходимо увеличивать отрезок Ь, чтобы отделить собственные значения уравнения с кубичной нелинейностью, отрезок Ь необходимо уменьшать. Причина этого заключается в стремлении А& при а > акр(Ь) к одному значению, равному а/2. При этом акр(Ь) уменьшается с ростом отрезка Ь.

Во-вторых, с ростом номера собственного значения для обеспечения сходимости итерационного метода необходимо выбирать начальное приближение достаточно близким к соответствующей собственной функции. Для этого в работе предложен один из возможных способов его построения, основанный на свойстве первой собственной функции: уменьшении области ее локализации с ростом параметра нелинейности.

В-третьих, целесообразно учитывать известное свойство разностных схем: номер максимальной рассчитываемой гармоники определяется шагом сетки. Поэтому для ограничения количества собственных функций с собственными числами, близкими к а/2, необходимо использовать "грубые" сетки. В этом случае можно определить собственные функции с номерами к 6 Затем, уменьшая шаг сетки, определить следующую группу собственных функций с номерами к 6 [к* + 1 ,&**]. Потом при необходимости повторить указанную процедуру.

Таким образом, предложенный в работе метод расчета собственных функций и собственных значений уравнения Шредингера с кубичной нелинейностью показал свою высокую эффективность для практически интересного случая а £ [1,100]. С его помощью найдены собственные функции с номерами к ^ 50. Важно подчеркнуть, что найденные собственные функции использовались в качестве начального условия А(х, 0) = фк{х) Для НУШ (1) для проверки достоверности полученных результатов. Нами рассмотрена их эволюция вдоль координаты z на отрезках 1 ^ Lz ^ 20. Расчеты подтвердили, что для всех рассмотренных значений параметра нелинейности а найденные собственные функции действительно сохраняли свою форму (|А(ж,0)|2 = |А(ж,,г)|2) на трассах распространения, существенно превышающих дифракционную длину отдельного субпучка при условии неизменности размера области по поперечной координате, на которой они были найдены. Однако при увеличении области по поперечной координате для собственных функций фк{х) с А& < а/2 профиль пучка в процессе распространения не сохранялся. Причина этого обсуждалась во введении.

Заметим также, что предложенный метод решения нелинейной задачи на собственные значения может быть обобщен на случай двух координат.

СПИСОК ЛИТЕРАТУРЫ

1. Захаров В.Е.,Манаков С.В., Новиков С.П.,Питаевский Л.П. Теория солитонов: Метод обратной задачи. М.: Наука, 1980.

2. Ablowitz М. Y., Segur Н. Solitons and the inverse scattering transformation. Phyladelphia: SIAM, 1981.

3. Kivshar Yu.S., Pelinovsky D.E. Self-focusing and transverse instabilities of solitary waves // Physics Reports. 2000. 33. P. 117-195.

4. Маймистов А. И. Распространение предельно коротких электромагнитных импульсов в нелинейной среде. Некоторые модели // Квант, электр. 2000. 30. № 4. С. 287-304.

5. Montesinos G.D., Perez-Garsia V.M., Torres P.J. Stabilization of solitons of the multidimensional nonlinear Shrodinger equation: Matter-wave breathers // Physica D. 2004. 191. P. 193-210.

6. Лыс а к T.M., Трофимов В. А. О возможности солитоноподобного режима двухволнового распространения фемтосекундных импульсов в оптическом волокне в условиях ГВГ // Оптика и спектроскопия. 2003. 94. № 4. С. 808-814.

7. Liu X., Beckwitt К., Wise F.W. Two-dimensional optical spatiotemporal solitons in quadratic media // Phys. Rev. E. 2000. N 1. P. 1328-1340.

8. Liu X., Beckwitt K., Wise F.W. Transverse instability of optical spatiotemporal solitons in quadratic media // Phys. Rev. Lett. 2000. 85. N 9. P. 1871-1874.

9. Akhmediev N.N., Mitzkevich N.V. Extremly high degree of Ar-soliton pulse compression in an optical fiber // IEEE J. Quant. Electron. 1991. 37. N 3. P. 849-857.

10. A n k i e w i с z A., Krolikowski W., Akhmediev N.N. Partially coherent solitons of variable shape in a slow Kerr-like medium: Exact solutions // Phys. Rev. E. 1999. 59. N 5. P. 6079-6087.

11. Громов E.M., Таланов В. И. Высшие приближения теории дисперсии нелинейных волн в однородных и неоднородных средах // Изв. РАН. Сер. физическая. 1996. 60. № 12. С. 16-28.

12. Коллатц Л. Задачи на собственные значения. М.: Наука, 1968.

13. Варенцова С. А., Пономарева Е.В., Трофимов В. А. О расчете собственных значений и собственных функций одномерного уравнения Шредингера на адаптивных сетках // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2000. № 3. С. 23-28.

Поступила в редакцию 12.05.04