Численное решение нелинейного уравнения Шредингера в декартовой системе координат Текст научной статьи по специальности «Физика»

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

А.А. Дегтярев, А.Е. Деркач Самарский государственный аэрокосмический университет

Аннотация

Для решения уравнения Шредингера, описывающего распространение электромагнитной волны в нелинейной среде, строится разностная схема типа Писмена-Рекфорда с итерационным уточнением, имеющая квадратичный порядок точности по всем трем пространственным переменным. Для исследования схемы рассмотрены частные случаи задачи, допускающие аналитическое решение. Проведено сопоставление аналитического и численного решений. Приведены результаты численного моделирования процесса распространения электромагнитной волны в нелинейной среде, позволяющие наблюдать как явление самофокусировки светового пучка, так и его самодефокусировки, что полностью согласуется с теорией электромагнитных волн [1].

Отмечается возможность применения разностной схемы для исследования процессов распространения электромагнитной волны в средах с пространственно-зависимым показателем преломления.

Введение

Как известно [1], нелинейное уравнение Шредингера является частным случаем волнового уравнения в параболическом приближении, записанном с учетом эффекта самовоздействия. Эффект самовоздействия проявляется при распространении оптического излучения в средах с кубичной нелинейностью (поляризация пропорциональна напряженности электрического поля в третьей степени).

С учетом этого эффекта уравнение Шредингера в декартовой системе координат запишем как [1]

dU

dz 2kn

-A DU + ■

о

iknH 2n0

-|U|2 U = 0, (1)

-b. < x <bL, -b. < y < ^ о < z < L, 2 2 2 2

где Л d = Л x +Л y =

^d 2

d

2 Л

dx2 dy2

оператор Лапла-

са в декартовой системе координат,

и - напряженность электрического поля, к - волновой вектор, По - показатель преломления среды, пнл - изменение показателя преломления под действием поля распространяющейся волны.

Отметим, что ось 2 совпадает с направлением распространения волны, а оси X и У лежат в плоскости, перпендикулярной оси 2.

Для получения замкнутой краевой задачи дополним уравнение (1) следующими краевыми условиями

= ф{х, y )

bx/2 = 0 /2 = 0 . = о

(2)

У=-Ly /2

y=Ly /2

= о

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

Без учета нелинейных эффектов уравнение Шредингера принимает вид

dU + _J_

dz 2kn

.ЛDU=0

(3)

0

В дальнейшем уравнение (1) будем называть нелинейным уравнением Шредингера, а уравнение (3) - линейным уравнением.

1. Расчетная конечно-разностная схема

с итерационным уточнением Для численного решения системы (1) - (2) построим разностную схему Писмена-Рэкфорда [2]. Определим сетку

Zk = k • hz ,

k = 0, K,

hz = L/K

xi = i • hx, i = -1/2,1/2,

hx = Lx /I

у= 1 • Ну, ] = -3/2,3/2, Ну = Ьу /3 , Кх, Ку - шаги дискретизации по переменным 2, х и у.

Обозначим через ик сеточный аналог численного решения и (х1, у 1, ).

Оператор Лапласа аппроксимируем следующим разностным оператором [2]:

, . ик,, - 2ик + и+ ,

л К ик = 1-11_у_1+11

ЛхиЧ~ К 2 .

"х

Задаче (1)-(2) поставим в соответствие следующую двухслойную нелинейную разностную схему:

J-rk + 1/2 Tjk Uij - Uij +_ l

0,5hz

2kn,

^+1/2 )+

2n0

rrk+1 jrk+1/2

UiJ - UV + l

+ ikПнл Uk+l/lf Uk+l/2 = 0

0,5hz

2kn

in- (

hj + ^U>+)+ (4)

+ iknrn Uk+1/2 I2 uk+1 /2 = 0 2n0 I 11 I 11

U0 = ф(г)

U-I/2 j = UI/2 j = Ui-J/2 = UiJ/2 = 0

+

z=0

0

Первое уравнение системы разностных уравнений (4) является нелинейным. Для нахождения неиз-

« л. ттк+1/2

вестной функции иа можно использовать итерационный метод последовательных приближений в сочетании с разностной прогонкой по х и у:

^к+1/2

и- и к I

0,5кх

2кп0

1кп„

( +к+1/2 . лкхх и у +лИууик

(5),

„п0

к+1/2

и„

*+1 и,

к+1/2

= 0

, к+1/2

= ик

^ у ■

причем и у

Таким образом, система (4) примет вид:

*+1к+1/2 и,

- и.

0,5кх

у , — +-

2кп0

,кп2

( *+1к+1/2 лкх и а

+л^уик

у у

2п0

, к+1/2

и„

*+1 и,

к+1/2

= 0

ттк+1 ттк+1/2

иа - иа + _ 1

0,5кх

2кп,

к<ик+112 +лк;ик+1)+

(6)

+ 1кп2л |и<к+1/2|2 ик+1/2 = 0

2п0 \ 11 I 11

, к+1/2

и,0 = ф(г), и ,1 = и,

ик/2 а = и1 /2 а = и,к-г/2 = ик

= 0

I/21 ~ /21 - /2 _ /2

В ходе теоретических исследований удалось доказать, что схема (6) имеет порядок аппроксимации О (¿х2, ¿у2, к] ).

В результате проведения вычислительных экспериментов получено подтверждение сходимости разностной схемы (6). Процесс итерационного уточнения обнаружил быструю сходимость: для уточнения решения до величин порядка 10-4 требуются 2-3 итерации.

2. Вычислительные эксперименты для среды с

постоянным показателем преломления

Рассмотрим сначала аналитическое решение линейного уравнения.

Будем решать линейное уравнение методом разделения переменных. В этом случае одно из частных решений задачи (1) -(2) запишется в виде

и (х, у, г ) = ехр 1

(( 2

71 71 +

\ \

V V

Ь2 Ь

012

/ У

Б1П

(х ^

V Ьх у

Б1П

( \

пу

V Ьу У

где а = ■

2кщ

а начальное условие

(пх А ( ^

ТЕС

(7)

ф(х, у) = БШ

V Ьх У

Б1П

лу

V Ь У

2.1. Результаты экспериментальных исследований линейного уравнения

Пучок, описываемый формулой (7), является стационарным, так как (х, у, 2) = !(х, у,0) , то есть распределение интенсивности в плоскости ОХУ не зависит от расстояния 2. Это распределение изображено на рис. 1.

Рис. 1. Распределение интенсивности волны в линейной среде. На рис. 2-4 приведены результаты численного решения линейного уравнения Шредингера на различных расстояниях от входа в среду.1

Для проведения численных расчетов были использованы следующие значения параметров: Ьх=Ьу =50 мкм; Х=0,63 мкм, причем Х=2л/к; п0 =1.

Рис. 2. Распределение интенсивности волны на расстоянии Ь от входа в среду (Ь = 100мкм, пх = 1000, пх = 70, пу =70).

Рис. 3. Распределение интенсивности волны на расстоянии Ь от входа в среду (Ь = 500мкм, пх = 5000, пх =70, пу =70).

1 В дальнейшем использованы следующие обозначения: пх, пу, пх- количество интервалов разбиений по осям х, у и

2 соответственно.

+

2

+

+

2

+

Рис. 4. Распределение интенсивности волны на расстоянии Ь от входа в среду (Ь = 1 мм, пх =5000, пх = 70, пу = 70).

Из приведенных рисунков видно, что на расстояниях порядка миллиметра распределения интенсивности волны, полученные с помощью численного решения уравнения Шредингера практически совпадают с аналитически рассчитанным распределением. Среднеквадратическая ошибка при этом не превосходит 5%.

3.2. Численное решение нелинейного уравнения Шредингера

Самофокусировка пучка

На рис. 5-7 приведены результаты численного решения нелинейного уравнения Шредингера сосле-дующими параметрами:

Ьх= Ьу = 50 мкм;

X = 0,63 мкм, причем Х=2л/к;

п0 = 1;

пнл =0,001, причем п

Рис. 5. Распределение интенсивности волны на расстоянии Ь от входа в среду (Ь = 100 мкм, пх = 1000, пх =50, пу = 50).

Рис. 6. Распределение интенсивности волны на расстоянии Ь от входа в среду (Ь =300 мкм, пх = 5000, пх = 50, пу = 50).

Рис. 7. Распределение интенсивности волны на расстоянии Ь от входа в среду

(Ь = 500 мкм, пх = 10000, пх =50, пу = 50).

В работе [1] показано, что при пнл >0 происходит фокусировка пучка, то есть смещение энергии пучка к центру. Это явление можно наблюдать и на рисунках 5-7.

Самодефокусировка пучка

На рис. 8-11 приведены результаты численного решения нелинейного уравнения Шредингера со следующими параметрами:

Ьх= Ьу = 50 мкм;

X = 0,63 мкм, причем Х=2л/к;

п = 1;

„ 2

=0,001, причем пнл<0.

Рис. 8. Распределение интенсивности волны на расстоянии Ь от входа в среду (Ь =200 мкм, пх = 1000, пх = 50, пу = 50).

□

Рис. 9. Распределение интенсивности волны на расстоянии Ь от входа в среду (Ь =500 мкм, пх = 1000, пх = 50, пу = 50).

нл

Рис. 10. Распределение интенсивности волны на расстоянии Ь от входа в среду (Ь = 800 мкм, П2 = 1000, пх = 50, пу = 50).

Рис. 11. Распределение интенсивности волны на расстоянии Ь от входа в среду

(Ь = 1000 мкм, т = 1000, пх = 50, пу = 50).

В работе [1] показано, что при пнл <0 происходит самодефокусировка пучка в нелинейной среде, то есть смещение энергии пучка в периферийную зону. Приведенные выше графики численного моделирования демонстрируют этот факт.

3. Вычислительные эксперименты для среды с показателем преломления, распределенным по параболическому закону

Пусть показатель преломления зависит от координаты х следующим образом

(

п0 = п

1 - 2Д

2 Л

(8)

причем пнл = 0.

Заметим [3], что при х<а

1 -Д|-

2 Л

(9)

Будем искать решение уравнения в виде и (х, у, 2) = А( ))(х )С (у). (10)

Подставляя (10) в уравнение (1), получаем фак-торизованное уравнение

21к дА 1 (1 1

---=—I —Л хВ +--Л уС | = -/,

А д2 п0 VВ х С у 1

(11)

где / - произвольная константа.

Таким образом, А(2) = К ехр| 1 -2^2 |. (12)

Преобразуя (11), получаем

1 1 2

- Л хВ + Х =- - Л уС = Я 2;

В С

2

где g - произвольная константа. С учетом (13)

С = ).

(13)

(14)

Представим В(х) в виде произведения двух функций

В(х )= X (х)ехр

( х2 Л

(15)

После подстановки (15) и (9) в (13), и осуществив замену переменных, получаем

X (х) = НА

(Л

V 0 у

где Нм - многочлен Эрмита, причем

2 N = ( - g 2 )) -1,

т0 =-

2а

(16)

(17)

(18)

Таким образом, с учетом (15), (14), (12) и (9) получаем искомое решение

(х>/21

и (х, у, 2) = К ехр( 121 siп(yg )Н„

т0

V 0 у

х ехр

( х2 Л

(19)

3.1.Результаты численного решения линейного уравнения Шредингера Для вычислительных экспериментов выберем, например, моду Гаусса-Эрмита с номером N=11.

Для выполнения граничного условия необхо-

(± ьЛ ^

димо, чтобы Ны

± Ьх42

т0

= 0 . Пусть, например,

= ±12,

(20)

так как Ны (± 12)«±10-16.

Для выполнения граничного условия при у=Ьу получаем

П

g =

Ьу

(21)

Из (17) и (18) можно найти / и р = -Д-

На рис. 12 изображено распределение интенсивности в пучке, описываемом формулой (19).

х

а

т

0

пп

0

а

а

На рис. 13 - 14 приведены результаты численного моделирования линейного уравнения со следующими параметрами: Я = 50 мкм;

X = 0,63 мкм, причем Х=2я/&; п1 = 1.

Как видно из приведенных графиков, результат численного решения уравнения Шредингера практически полностью совпадает с результатом аналитического решения.

Рис. 12. Распределение интенсивности волны в линейной среде.

Рис.13. Распределение интенсивности волны на расстоянии Ь от входа в среду (Ь=1 мм, пх =1000, пх = 80, пу = 80).

Рис.14. Распределение интенсивности волны на расстоянии Ь от входа в среду (Ь = 10 мм, пх =50000, пх = 80, пу = 80).

Заключение

Проведенные исследования подтвердили высокую точность численного решения задачи на расстояниях порядка миллиметра. Так для модели среды, не учитывающей эффект самовоздействия, численное решение практически совпадает с аналитическим, рассчитанным по формуле (7).

Для среды с "параболическим профилем" показателя преломления моды Гаусса-Эрмита, рассчитанные с помощью схемы (6), практически не претерпевают изменений в процессе распространения волны на расстояния порядка десятков миллиметров, что соответствует теории [3].

Для модели среды, учитывающей эффект самовоздействия, численный расчет позволяет наблюдать как явление самофокусировки излучения, так и его самодефокусировки, что спрогнозировано теоретическими оценками, например, в работе [1].

Литература

1. Виноградова М.Б., Руденко О.В., Сухору-ков А.П. Теория волн // М., Наука. 1979.

2. Самарский А. А. Теория разностных схем // М., Наука. 1977.

3. Адамс М. Введение в теорию оптических волноводов // М., Мир. 1984.