CC BY
12
8. Andreev V. B. Numerical solution of singularly perturbed boundary value problems with power and logarithmic boundary layers // Analytical and Computational Methods for Convection-Dominated and Singularly Perturbed Problems / Eds.: J.J.H. Miller, G.I. Shishkin, L. Vulkov. N.Y.: Nova Science Publishers, 2000. P. 1-5.
Поступила в редакцию 05.05.06
УДК 517.958:535.14
П. В. Догадушкин, В. А. Трофимов
КОНСЕРВАТИВНАЯ РАЗНОСТНАЯ СХЕМА ДЛЯ ЗАДАЧИ РАСПРОСТРАНЕНИЯ ФЕМТОСЕКУНДНОГО ИМПУЛЬСА В СЛОИСТОЙ ПОГЛОЩАЮЩЕЙ СРЕДЕ С НЕОТРАЖАЮЩИМИ КРАЕВЫМИ УСЛОВИЯМИ1
(кафедра вычислительных методов факультета ВМиК)
1. Введение. Компьютерное моделирование распространения световых импульсов в слоистой (периодически) поглощающей среде широко анализируется в задачах взаимодействия лазерного излучения с веществом. Это прежде всего распространение оптического излучения в фотонных кристаллах [1]. Подобные же задачи имеют место при рассмотрении конденсата Бозе-Эйнштейна [2-4] и при анализе конверсии оптических частот в средах с периодическими свойствами [5, 6]. Как правило, рассмотрение этих проблем базируется на нелинейном уравнении Шредингера с периодическими параметрами, характерной чертой эволюции решения которого является многократное отражение световых импульсов на границах раздела сред. В результате этого временной интервал, на котором анализируется решение, становится достаточно большим. Как следствие этого нулевые краевые условия для комплексной амплитуды необходимо задавать на границе, которая значительно удалена от нелинейной среды. Таким образом, эффективность компьютерного моделирования существенно снижается, так как решение вне области, занятой периодической средой, не представляет большого интереса. Поэтому необходимо уменьшить расчетную область, используя неотражающие краевые условия. Они применяются для многих задач нелинейного взаимодействия лазерного излучения с веществом [7-9]. Однако неотражающие краевые условия для нелинейного уравнения Шредингера практически не обсуждались в литературе, за исключением [10], где предложен один из возможных типов данных условий для прозрачной среды. В настоящей же работе рассматривается задача распространения светового импульса в периодически поглощающей среде, для которой предложены несколько видов неотражающих краевых условий и проведено их сравнение.
Следует подчеркнуть, что при построении неотражающих краевых условий для данного круга задач важным требованием является непревышение максимума модуля комплексной амплитуды отраженной волны некоторого значения, чтобы отсутствовало влияние отраженной волны на динамику взаимодействия падающей волны с нелинейным (и даже линейным) кристаллом. Эмпирически установлено, что максимальная амплитуда отраженной волны не должна превышать 0,1% от падающей амплитуды. Известные же в литературе неотражающие краевые условия, применяемые, например, для линейных уравнений Максвелла, дают коэффициент отражения по амплитуде, равный 3-5%. Таким образом, строящиеся нами неотражающие условия позволяют в 30-50 раз уменьшить коэффициент отражения.
1 Работа выполнена при частичной финансовой поддержке РФФИ (грант № 05-01-507).
2. Постановка задачи. Процесс распространения оптического излучения в периодической, линейно (или нелинейно) поглощающей среде в одномерном случае описывается следующим волновым уравнением:
д2Е 1 д2Е 4тг д2Р 1 д(<тЕ) „
= + Р = (1)
В (1) Е{х,£) — напряженность электрического поля, с — скорость света, х — координата, вдоль которой распространяется световой импульс, £ — время, а характеризует поглощение среды, и — частота, х'1' — линейная восприимчивость, енл описывает нелинейную добавку к поляризации среды.
Для записи уравнения относительно медленно изменяющейся во времени амплитуды А представим напряженность электрического поля в виде [11]
Е(х,Ь) = ^[А(х,Ь)е^ + к.с.], (2)
где символы к.с. означают комплексное сопряжение. После подстановки (2) в (1) и проведения стандартных упрощений, соответствующих методу медленно изменяющихся амплитуд, пренебрегая вторыми производными по времени, получим следующее уравнение:
*£.ы1М + {ек2_.к(г)А = 0} (3)
где к = ^ — волновое число, А — длина волны, е = п2 = 1+47ГХ'1' — диэлектрическая проницаемость линейной среды.
Введем следующие безразмерные переменные х = = —, А = -г-, параметры г и Ат определим
Л Т А-тп
ниже. В этих переменных уравнение (3) для случая линейного распространения (енл = 0, в дальнейшем рассматриваем именно этот случай) запишется в виде
дА ,^д2А (. ¡кст\ — ^ 1 ст с аст . .
Если в качестве нормирующего множителя г по времени выбрать период волны Т, то коэффициенты
1 г _ <та
4тге ' 0 2е
уравнения (4) будут соответственно равны И = 5о = и оно примет вид (черту над перемен-
ными для простоты опустим)
дА д2 А
— + {0— + 605(х)А + тА = 0. (5)
В (5) для общности коэффициент поглощения записан как функция продольной координаты х. Нетрудно видеть, что при этом предыдущий вывод уравнения не изменяется.
Введя еще одну новую функцию А = Ае~гж*, подставив ее в (5) и опустив в дальнейшем волну над комплексной амплитудой, получим следующее безразмерное уравнение:
дА д2 А
— + + 805(х)А = 0, ¿>0, 0 < ж < (6)
где Ь — граница области, в пределах которой анализируется взаимодействие оптического импульса со средой. При этом в периодически поглощающей среде функция 5(х) имеет вид
+ (7)
I 0, если х £
где С к — область (пит), в которой имеет место поглощение световой энергии, Ир — число питов. Вне питов физической ситуации соответствует равенство нулю функции 8{х).
Граничные и начальные условия для комплексной амплитуды с учетом ее нормировки имеют вид
А(0, = А(Ь, = 0, (8)
А(х, 0) = А0(ж) = ехр < - ( -—— ) - ¡х(х - х0) > , (9)
где хо — центр распределения комплексной амплитуды светового пучка, а — продольный размер пучка (он также характеризует длительность импульса), выраженный в единицах длины волны, ^ = кХ = 2тг — безразмерное волновое число.
3. Неотражающие краевые условия. Недостатком постановки задачи (6)-(9) являются нулевые граничные условия (8), которые ограничивают эффективность компьютерного моделирования из-за необходимости задания большой области по координате ж, чтобы оптическое излучение не дошло до ее границ на протяжении всего времени расчетов. Поэтому для ее повышения целесообразно поставить неотражающие краевые условия для исходной задачи. Для их построения рассмотрим эквивалентную (6) систему уравнений относительно действительной и = Ее А и мнимой V = 1т А частей комплексной функции А при отсутствии поглощающих питов:
ди ду ^д2 и г . .
Ъ = *>0' 0<ж<х- (10)
Решение этой системы представляет собой, в частности, волну и = Щ соз(о^—хж), V = —С/о Xх)-,
распространяющуюся вдоль оси ж со скоростью Их (или —Ох) и частотой и = Их2 ■ Следовательно, волне, распространяющейся в положительном направлении ж, соответствует граничное условие
ди ^ ди ду ^ дю . .
или
ди „ 9 дю „ 9 .
— = £>хЧ ^ = -Вх2Щ ® = 0, ¿>0. (12)
Условие (12) означает, что на границе световая волна представляет собой гармоническую волну с частотой .Ох2 - Возвращаясь снова к комплексной амплитуде, условия (11) и (12) запишем соответственно в виде
дА дА
дА
— - ШХ2А = 0, ж = 0. (14)
Необходимо отметить, что условие (14) не содержит производной по пространственной координате, т.е. не представляет собой распространяющуюся волну. Поэтому непосредственно для построения неотражающих краевых условий оно непригодно. Однако сочетание условий (13) и (14) позволяет получить хороший результат. Заметим также, что еще одно краевое условие получается при вычитании (14) из (13):
дА
— -гХА = 0, ж = 0, (15)
которое означает, что на границе волна имеет волновое число х-
Для правой границы ж = Ь запишем аналогичные (13) и (15) условия
дА дА
Ж + = ж = £, (16)
дА
— + ixA = 0, х = Ь. (17)
Следует подчеркнуть, что записать неотражающие граничные условия для исходных уравнений можно и на основе введения новой функции А = Аег%х, А — медленно изменяющаяся во времени и по пространственной координате комплексная амплитуда, распространяющаяся в отрицательном направлении оси ж. Перейдя в (6) к А и отбрасывая вторую производную от нее, запишем уравнение на левой границе среды:
д~А д~А _
— -2£>х^-¿Х2£>А = 0, Ж = 0, ¿>0. (18) Возвращаясь снова к исходной функции А, получим
дА дА
— -2£>Х^ + ^Х2£>А = 0, Ж = 0, ¿>0. (19)
Аналогичное неотражающее условие на границе х = Ь имеет вид
дА дА
— + 2£>х^ + *Х2£>А = 0, ¿>0, х = Ь. (20)
Нетрудно видеть, что (19) можно получить вычитанием из (13) выражения (15), предварительно умноженного на Их- Следовательно, полученные краевые условия (19), (20) являются линейной комбинацией условий (13), (15) и (16), (17) соответственно. Поэтому они описывают как распространение волны, так и то, что она имеет частоту Ох2 и волновое число Их- Далее среди записанных краевых условий будем рассматривать только условия (13), (15) и (19) (и соответствующие им условия на правой границе) и проводить компьютерное моделирование только для них. Заметим, что во введенных выше безразмерных переменных параметры Их = ^ и Ох2 = 7 определяются только свойствами вещества.
4. Инварианты распространения. На примере инварианта энергии продемонстрируем различия между рассматриваемыми условиями (13), (15) и (19). Для записи инварианта уравнение (6) умножим на А*, а сопряженное к нему — на А, сложим их, а затем полученную сумму проинтегрируем по координате х от 0 до Ь. В результате получим следующее интегральное соотношение, описывающее процесс изменения энергии среды:
ь
а г „ /ял ял*
о
(21)
которое показывает, что изменение энергии оптического излучения в среде равно потоку его энергии через ее границы.
Выразив ^ и из условий (13), (15), (19) и (16), (17), (20), для \¥\ получим соответственно три выражения:
i (дА _ дА* А
= ^ (■т-4 -1г-4
х— Ь
+1 (а4л- - 94-А
х=0
Цг1 = 20Х(\А\2 +|А|2 ),
V х=Ь ж=0/
х— Ь
= 0. (22)
х=0
Последнее выражение, таким образом, подтверждает то, что (19) и (20) являются линейной комбинацией выражений (13), (15) и (16), (17).
Для записи второго инварианта умножим (6) на ^г, а сопряженное к нему — на вычтем второе из первого, а затем проинтегрируем полученное уравнение по продольной координате от 0 до Ь. В результате получим
т =
д_ т
дА* л дА . -—А - —А* ) йх + дх дх
дА Л* дЛ л
-—А--—А
дь дь
ь дА 2 Ь
- 1У2 = 0, 1У2 = -2Ш
0 дх 0
(23)
Выразив член 2 через граничные условия (13), (15), (19) и (16), (17), (20), для него соответственно запишем следующие соотношения:
1¥2 =
2г дА 2
вх 2 т
1¥2 = -2ШХ2 |А|2
ь
1 1У2
0
на да*
9(
Ох2
дА
т
- ъВх2 |А|2
(24)
да 9(
из первого уравнения второе, проинтегрируем полученное выражение по продольной координате от 0 до Ь. В результате получим
дА
дх
-о 1У~ I—— дЛ* 9А
3 ' 3 дх дЬ дх дЬ
С учетом граничных условий имеем следующие выражения:
2 дА 2 2 дА 2
Ох т т х=0
тж, • , лдА* л*дА
^ = ** И-зГ - л а"
х— Ь
. . дА* л^дА
х=0
дА
т
1х ( лдА* л„дА + — А—--А* —
2 V дЬ дЬ
и.
дА
т
+н №
2 V дг
л*дА
(26)
х=0
соответственно для (13), (15), (19) и (16), (17), (20).
5. Разностная схема. Для решения задачи (6) с соответствующими начальными и граничными условиями введем в области С = {0 ^ ж ^ 0 ^ £ ^ равномерную сетку
ш = шх х^, их = {х1 = 11ъ, I = 0,1,..., Л^, кНх = Ь}, = = ]т, 1 = 0,1,..., Л^, тNt =
Определим следующие сеточные функции и введем безындексные обозначения (для простоты сохраним прежние обозначения):
0,5 0,5 I
А = А1 = А(х1,г1), А = А1 = А(Х1,г1+1), А = А\ = - (А(х1,г1) + А(х1,г1+1)) ,
а,ха=А1+^-2А^ + А1-^, ах=А1+^~А1'\ а-х=А1-~А1-^, 5 = 51 = 5(х1).
На сетке и запишем следующую разностную схему:
^—^ + 10К-ХХА+ 508>А=0, 0<1<НХ, 0<j^Nt, (27)
г
с начальными условиями вида
А(х1, 0) = ехр | — (Х'~ХЛ _ ¿х(Жг _ Жо) 1 о^/^ЛГ*. (28)
Для каждого из рассматриваемых граничных условий (13), (15), (19) и (16), (17), (20) запишем соответственно разностные уравнения
А - А °>5 А - А °>5
--БХАх = 0, 1 = 0, -+ 0ХА-х = О, 1 = МХ, (29)
т т
0,5 0,5 0,5
Ах-1ХА= 0, 1 = 0, А£ + 1ХА=0, 1 = МХ, (30)
А — А °'5 °'5 А — А °'5 °'5
--2ВХАх + 1х2ОА =0, 1 = 0, -+ 2ВХА-х + 1Х2ВА =0, 1 = МХ. (31)
т т
Нетрудно заметить, что (29)—(31) аппроксимируют исходные краевые условия с порядком О (г2 + К) относительно точки (ж/, ¿^+0,5)- Если же условия (13), (15), (19) и (16), (17), (20) аппроксимировать
симметричными относительно точки (ж/+о,5 5 ¿7+0,5) выражениями
0,5 0,5 0,5 0,5
Ао - Ар + ¿! - Ах _ ^Аг-Ар _ АМт - АМх + АМх_г - АМх_г „ К К I _ .. г,.))
2т Х к _и' 2т + Х к _и'
0,5 0,5 . . 0,5 0,5 . .
А1 - А0 гу /0,5 0,5 \ Адг - Адг _1 ¿у / о,5 0,5 \
——т[А• + *)=<>. + Т <33)
Л°-А°2+г-4'--4' - + ^ (X + Ж) = О,
0,5 0,5 . .
Адг - Адг + Адг -1 - Адг -1 „ Адг Адг —^ Ц)Х - ( " Г. 0,5 \
- 2г -^^ • '21)\- ь ■ I >\ • >\ I I <)• (34)
то тогда получим второй порядок О (г2 + Ь2) относительно этой точки.
Приведем также разностные краевые условия, которые аппроксимируют исходные краевые условия со вторым порядком относительно точки (x^,tj+otrJ). Например, для (19) и (20) они имеют вид
А — А °'5 °'5 А — А °>5 °>5
(1 -гХЬ)--2ВХАх + 1Х2ВА = 0, / = 0, (1 + гХй)-+ 2£>уА^у2£> А = О, 1 = НХ.{Щ
т т
6. Сравнение эффективности граничных условий. Эффективность неотражающих краевых условий (29)-(35) оценивалась на основе компьютерного моделирования. В проведенных экспериментах шаги сетки и по времени и пространству выбирались согласованно так, чтобы выполнялось равенство ^ = В процессе расчетов вычислялись коэффициенты отражения световой энергии
и максимальной интенсивности от границ среды по формуле
= \А(хьЬ)\2, (36)
где Аотр — амплитуда отраженной волны, tj для Яе выбирается таким образом, чтобы импульс полностью прошел через границу. Заметим, что необходимость определения последней величины обусловлена тем, что нелинейное взаимодействие светового импульса со средой зависит именно от интенсивности волны.
В качестве иллюстрации полученные результаты представлены в табл. 1 в виде зависимости коэффициентов отражения от шагов сетки для различных видов краевых условий. Как видим, при уменьшении шагов в 2 раза коэффициент отражения изменяется практически в 4 раза для условий (29)—(31). Коэффициенты отражения для условий (32)-(35) слабо зависят от шагов сетки. Из таблицы следует, что среди условий с первым порядком аппроксимации по пространственной координате на границе наилучшими являются условия (31), а со вторым порядком — (34).
Необходимо подчеркнуть еще одну важную черту. В построенных краевых условиях используется волновое число у, которое присутствует в начальном распределении комплексной амплитуды. Однако при х < Ю у светового импульса, дошедшего до границы, волновое число существенно изменяется по сравнению с первоначальным значением, что приводит к значительному росту коэффициентов отражения. Для его уменьшения в граничных условиях целесообразно использовать локальное значение волнового числа у/, определеное вблизи границ. Для примера соответствующие значения приведены в табл. 2 для схемы (27)-(29). Так, при начальном распределении с волновым числом у = 1,0 дошедшая до границы волна имеет локальное волновое число, равное у/ = 3,074. Поэтому использование локального волнового числа приводит к уменьшению на порядок (и более) энергии и максимальной амплитуды отраженной от границы волны.
]%е —
У, I А>тр(ж/, £7)1 1=0
ЛГ
Е |А>тр(ж/,о)| 1=0
Таблица 1
Зависимость доли отраженной энергии и максимальной интенсивности от шагов сетки ш при х = 80, О = 0,001, а = 0,2
Условие h т Re Ri
(29) 0,005 0,03125 0,0102 0,0101
0,0025 0,015625 0,002662 0,002680
0,00125 0,0078125 0,0007819 0,0008576
(30) 0,005 0,03125 0,0102 0,009858
0,0025 0,015625 0,002662 0,002656
0,00125 0,0078125 0,0007819 0,0008447
(31) 0,005 0,03125 0,01007 0,009763
0,0025 0,015625 0,002505 0,002404
0,00125 0,0078125 0,00062573 0,00060107
(32) 0,005 0,03125 0,00032407 0,00023843
0,0025 0,015625 0,00016451 0,00012471
0,00125 0,0078125 0,00015638 0,00011451
(33) 0,005 0,03125 0,00020828 0,00015624
0,0025 0,015625 0,00016084 0,00011806
0,00125 0,0078125 0,00015703 0,00011155
(34) 0,005 0,03125 0,00009859 0,00009292
0,0025 0,015625 5,59 х 10"b 4,76 x 10"b
0,00125 0,0078125 2,7 х 10"' 1,4 x 10"'
(35) 0,005 0,03125 0,0987 0,009531
0,0025 0,015625 0,002415 0,002309
0,00125 0,0078125 0,00062034 0,0005997
Таблица 2
Сравнение характеристик отраженной волны для условий (29) с XI или х на сетке с шагами h = 0,005, т = 0,1 и параметрами D = 0,001, а = 0,2
X Re Ri
xi X xi X
0 0,0146 0,984 0,0135 0,95
1 0,0076 0,452 0,0071 0,447
2 0,0024 0,013 0,0022 0,012
6 0,0004 0,006 0,0004 0,006
10 0,0001 0,0015 0,0001 0,0014
СПИСОК ЛИТЕРАТУРЫ
1. Joannopoulos J.D.,Meadl R. D.,Winn J.N. Photonic crystals: molding the flow of light. N. Y.: Princeton, 1995.
2. Pitaevskii L.P., Stringari S. Theory of Bose-Einstein condensation in trapped gase // Reviews of Modern Physics. 1999. 71. N 3. P. 463-512.
3. Abdullaev F.K., Gammal A., Kamchatnov A.M., Tomio L. Dynamics of bright matter wave solitons in a Bose-Einstein condensate // Int J. of Modern Physics. B. 2005. 19. N 22. P. 3414-3473.
4. Gubeskys A., Malomed B.A., Merhasin I.M. Two-component gap solitons in two- and one-dimensional Bose-Einstein condensates // Phys. Review. A. 2006. 73. 023607.
5. Junqiang Sun, Deming Lin, Dexiu Jeuang. Difference-frequency generation among ultrashort optical pulses in quasi-phase—matching waveguides // Optical and Quantum Electronics. 2004. 36. P. 577-587.
6. Junqiang Sun. Wavelength conversion between picosecond pulses in periodically Poled NiNbOs waveguides using pulsed pumping // Microwave and Optical Technology Letters. 2005. 46. N 5. P. 464466.
7. Alpert B.,Greengard L., H agst ram J. Nonreflecting boundary conditions for the time-dependent wave equation // J. Comput. Physics. 2002. 180. P. 270-296.
8. Turkel E., Yefet A. Absorbing PML boundary layers for wave-like equations // Applied Numerical Mathematics. 1998. 27. P. 533-557.
9. Fibich G., Tsynkov S. High order two-way artificial boundary conditions to nonlinear wave propagation with back-scattering // J. Comput. Physics. 2001. 171. N 2. P. 632-677.
10. Те реш и h Е.Б., Трофимов В. А., Федотов M. В. Консервативная разностная схема для задачи распространения фемтосекундного импульса в нелинейном фотонном кристалле с неотражающими краевыми условиями // ЖВМиМФ. 2006. 46. № 1. С. 161-171.
11. Трофимов В. А. Инварианты распространения фемтосекундных импульсов в фотонных кристаллах // ЖВМиМФ. 2001. 41. № 9. С. 1429-1433.
УДК 519.2
О. В. Шестаков
ОЦЕНКА ТОЧНОСТИ ВОССТАНОВЛЕНИЯ ФУНКЦИИ
ПО ЕЕ ЭКСПОНЕНЦИАЛЬНОМУ ПРЕОБРАЗОВАНИЮ РАДОНА
ПРИ ИСПОЛЬЗОВАНИИ КОНЕЧНОГО ЧИСЛА ПРОЕКЦИЙ1
(кафедра математической статистики факультета ВМиК)
1. Введение. Проблема обращения экспоненциального преобразования Радона, которое отличается от классического преобразования Радона наличием экспоненциального множителя у интегрируемой функции, возникает в задачах однофотонной эмиссионной томографии [1]. По полному набору экспоненциальных проекций можно единственным образом восстановить функцию, описывающую объект. Однако на практике известно лишь конечное число проекций. В этом случае единственность не гарантирована даже когда коэффициент поглощения в экспоненциальном множителе равен нулю (т.е. при обращении классического преобразования Радона). Для этого случая доказана теорема, утверждающая, что если известны проекции объекта по конечному числу направлений, то существует другой объект, имеющий такой же носитель и такие же проекции по данным направлениям, что и первый объект, и совершенно произвольную структуру на любом компактном множестве, содержащемся внутри носителя первого объекта. Эта теорема приводит к так называемому парадоксу компьютерной томографии, поскольку томографы позволяют получать пригодные для практической цели изображения. Этот парадокс решается в работах [2] и [3] путем сглаживания искомой функции плотностью стандартного нормального распределения, которая выступает в роли регуляризатора. Также получены оценки близости функций, имеющих одинаковые или близкие проекции по конечному набору заданных направлений. В данной работе такие оценки получены для экспоненциального преобразования Радона.
2. Экспоненциальное преобразование Радона. Пусть функция /(2:1,2:2), описывающая объект, имеет компактный носитель. В данной работе мы будем полагать, что этим носителем является круг II с единичным радиусом и центром в начале координат. Кроме того, будем считать, что функция /(2:1, х2) неотрицательна и нормирована, т.е. интеграл от нее по всей области, где она отлична от нуля, равен единице. Множество всех таких функций обозначим за Г¡7. Экспоненциальным преобразованием Радона для коэффициента поглощения д > 0 и функции /(2:1,2:2) называется преобразование вида
Поступила в редакцию
22.03.06
ОС
RiifiWi = / efit f(s cos Ф + £ sin s sin cp — t cos cp) dt, s G R, G [0, 2îr).
— ОС
1 Работа выполнена при финансовой поддержке РФФИ (гранты 05-01-00535, 04-01-00671).
