1626435587-55f52a4de97976f3c6215fa7c103f544 (844241), страница 18
Текст из файла (страница 18)
Альтернативный вывод нелинейного уравнения Шрёдингера без использования приближения медленно меняющейся амплитуды приведён в статье[A15].4.5. Численные методыПрежде чем переходить к рассмотрению различных методов численного интегрирования нелинейного уравнения Шрёдингера (112), обратим внимание на сходство его линейной части = 12 2 с уравнением теплопроводности. Как мы помним из гл. 3, основная сложность построения численного решения уравнения теплопроводностибыла связана с построением устойчивых численных схем. В этой связиначнём наше рассмотрение именно с вопросов устойчивости.Несложно заметить, что, в отличие от уравнения теплопроводности,ˆ = − 2 2 в правой части уравненияспектр линейного оператора 2(112) целиком лежит на мнимой оси, что соответствует набегу фазыпри распространении волн вместо диссипации энергии в задаче о теплопроводности. Как следствие, разностный аналог уравнения Шрёдингера напрямую не попадает под определение жёстких систем [1, с.
103].Тем не менее, вполне очевидно, что появление множителя не влияетна ключевые выводы, сделанные в гл. 3 для уравнения теплопроводности. А именно, при использовании простейших явных разностных схемдля уравнения (112) шаг численного интегрирования по должен выбираться достаточно малым для предотвращения взрывного роста амплитуды самой высокочастотной гармоники на сетке exp(max ). Приэтом проявляется характерная черта жёстких систем: малость шага ℎчисленного интегрирования определяется не масштабом, на которомэволюционирует физическое решение, но свойствами уравнения и сетки, используемой для построения численного решения.Как и в случае уравнения теплопроводности, гармоника с частотойНайквиста характеризуется самым большим по модулю собственнымзначением:⃒⃒⃒ 2 ⃒ 2 |2 |ˆ⃒ = , max || = ⃒ max ⃒⃒ =,(125)22 2где — шаг сетки по .
Максимальное собственное значение (125) быстро (квадратично) возрастает с уменьшением шага сетки . При использовании явных схем это приводит к необходимости соответствующегоуменьшения шага ℎ сетки по , ℎ ∝ 2 . Ниже мы рассмотрим альтернативное решение — использование безусловно устойчивых численных92схем. Хотя для уравнения (112) можно использовать неявные схемы,подобно тому как это было сделано в гл. 3 для уравнения теплопроводности, гораздо интереснее рассмотреть два принципиально иныхподхода, позволяющих построить эффективные безусловно устойчивыечисленные схемы.4.5.1. Метод расщепления по физическим процессамЗапишем уравнение (112) в видеˆ +ˆ ), = ((126)ˆ иˆ — дисперсионный и нелинейный операторы соответственно:где ˆ = − 2 2 ,ˆ = ||2 .(127)2Будем обозначать посредством и Ψ точное и численное решение уравнения (126) соответственно.Основная идея метода расщепления по физическим процессам33 заключается в разделении шага ℎ численного интегрирования по координате на две или более части, на каждой из которых учитываетсятолько дисперсия либо только нелинейность.
Похожий приём использовался нами в локально одномерном методе для решения уравнения теплопроводности (см. п. 3.7 на с. 82), хотя, как мы увидим далее, данныеметоды основаны на различных механизмах обеспечения устойчивостичисленного решения.Самый простой способ продолжения численного решения из точки до +ℎ методом расщепления может быть записан следующим образом:(︀ )︀(︀ )︀ˆ exp ℎˆ Ψ(),Ψ( + ℎ) = exp ℎ(128)тогда как точное решение эволюционирует на шаге ℎ по закону(︁ (︀)︀)︁ˆ +ˆ ().( + ℎ) = exp ℎ (129)точныйКлючевым моментом являетсяучёт действия дисперсии вФурье-представлении, см.
п. 4.1 и формулу (113) на с. 86. Использование точной аналитической формулы вместо разностной схемы даётабсолютную (не зависящую от шага сетки ℎ) устойчивость численнойсхемы. Подставляя (113) в (128), имеем:(︂)︂(︀ )︀ˆ Ψ(),Ψ( + ℎ) = ˆ − exp2 2 ℎ ˆ + exp ℎ233Step-split Fourier method в англоязычной литературе.93где ˆ + и ˆ − — прямое и обратное дискретное преобразование Фурье.Применение эффективных программных реализаций алгоритмов быстрого преобразования Фурье позволяет выполнять вычисления быстрее,чем при использовании разностных схем.В отличие от дисперсионного оператора, нелинейность не приводитк развитию неустойчивости, поэтому действие нелинейного оператораˆ ) может быть реализовано любым способом: как точно, в соотexp(ℎветствии с формулой (116), так и с помощью какой-либо разностнойсхемы для уравнения = ||2 .Таким образом, для выполнения шага численного интегрированияв соответствии с (128) необходимо выполнить следующие шаги (символ ← в выражениях ниже соответствует оператору присваивания впрограмме):ˆ )Ψ();1) нелинейный шаг: Ψ() ← exp(ℎ2) прямое преобразования Фурье: Ψ̃() ← ˆ + Ψ();3) линейный шаг в Фурье-представлении: Ψ̃() ← Ψ̃() exp(︀ 2 2 2)︀ℎ ;4) обратное преобразования Фурье: Ψ() ← ˆ − Ψ̃().Исследуем погрешность схемы (128).
В случае, если действие опеˆ ) в (128) вычисляется точно (Ψ() ← Ψ() exp(|Ψ|2 ℎ)),ратора exp(ℎединственным источником ошибки (не считая неизбежных погрешноˆ +ˆ )) на суперпостей округления) является замена оператора exp(ℎ(ˆ exp(ℎˆ ). Данная замена приводит к ошибке, связаннойзицию exp(ℎ)ˆ и ˆ . Вычислим погрешностьс некоммутативностью34 операторов схемы (128) на одном шаге, вычитая из точного решения ( + ℎ) численное, построенное в соответствии c (128):[︁]︁^^^ ^^^1 = ( + ℎ) − ℎ ℎ () = ℎ(+ ) − ℎ ℎ ().Расписывая разность в квадратных скобках по определению через рядТейлора для экспоненты, имеем:^ ^ ] легко понять из физических соображений:Отличие от нуля коммутатора [,как мы видели в п.
4.1, действие дисперсии приводит к увеличению длительностиимпульсов и пропорциональному уменьшению их пиковой мощности . Посколькунелинейный набег фазы = пропорционален , перестановка ^ и ^ будетприводить к изменению ⇒ ^ ^ ̸= ^ ^ , что и требовалось показать.3494[︂2 (︀(︀)︀)︀ˆ +ˆ +ℎ ˆ +ˆ 2 + ...−1+ℎ 2(︂)︂ (︂)︂]︂22ˆ+ℎ ˆ2 + . . .ˆ+ℎ ˆ 2 + . . .
() =− 1 + ℎ1 + ℎ22)︁2 (︁ℎˆˆ −ˆˆ () + (ℎ3 ).=21 =Таким образом, погрешность численной схемы (128) на одном шагеесть (ℎ2 ). При интегрировании по отрезку длины = const с использованием /ℎ шагов погрешность возрастает до (ℎ), так что схема(128) обеспечивает первый порядок точности по шагу -сетки.Точность метода расщепления можно повысить за счёт использования симметричной схемы:(︂)︂(︂)︂(︀ )︀ℎˆℎˆˆΨ( + ℎ) = exp exp ℎ exp Ψ().(130)22Расписывая аналогичным образом экспоненты от операторов черезсумму рядов Тейлора, несложно показать, что погрешность симметричной схемы (130) на одном шаге есть:[︁]︁1^ ^^^ 1 ^1 = ℎ(+ ) − 2 ℎ ℎ 2 ℎ () = (ℎ3 ).Следовательно, схема (130) имеет второй порядок точности по шагусетки ℎ.
Абсолютная устойчивость схемы (130) достигается за счёт использования точного (аналитического) ответа для эволюции решенияˆ.под действием дисперсии групповых скоростей Формально сравнивая формулы (128) и (130), можно прийти к выводу о том, что повышение порядка точности достигается ценой дополнительного действия дисперсионного оператора на каждом шаге.Однако в действительности объём вычислений при переходе от (128)к (130) почти не возрастает. В самом деле, для выполнения шаговпо формуле (130) необходимо подействовать на численное решение Ψоператором)︁(︁)︁(︁ 11^^ 1 ^ ^^^ −1 ℎ^ 1 ^(2)= 2 ℎ ℎ ℎ 2 ℎ .ˆ×ℎ = 2 ℎ ℎ 2 ℎ(131)ˆ в конце и началеКак видно из (131), группировка операторов exp( 21 ℎ)промежуточных шагов позволяет записать оператор эволюции ˆ×ℎ в95виде суперпозиции нелинейных и + 1 дисперсионных операторов.По сравнению со схемой первого порядка (128), -кратное применеˆ exp(ℎˆ )) , произведение в форние которой даёт оператор (exp(ℎ)муле (131) содержит всего на один оператор больше и, соответственно,требует одного дополнительного быстрого преобразования Фурье, чтопочти незаметно при большом числе шагов .
Данная ситуация вполнеаналогична переходу от квадратурной формулы правых (левых) прямоугольников к формуле трапеций, также позволяющей повысить наединицу порядок точности ценой всего одного дополнительного вычисления подынтегральной функции, см. [1, с. 59].Реализовав в программном коде преобразования численного решеˆ и exp(ℎˆ ), совмещённые с вызония под действием операторов exp(ℎ)вом функций для быстрого преобразования Фурье, достаточно простозапрограммировать метод расщепления четвёртого порядка точностипо ℎ [A15]. Для этого нужно, используя симметричную схему (130),сделать 4 шага ℎ вперёд, один шаг 2ℎ в обратном направлении и ещё 4шага ℎ вперёд. Соответствующее преобразование численного решенияΨ может быть записано в виде:(2)(2)(2)Ψ( + 6ℎ) = ˆ4×ℎ ˆ1×(−2ℎ) ˆ4×ℎ Ψ(),(132)(2)где оператор ˆ×ℎ определён в соответствии с выражением (131).
Группировка дисперсионных операторов в (132) позволяет сократить объёмвычислений до 19 преобразований на шаге 6ℎ:(︁)︁(︁)︁ℎ ^ℎ ^ℎ ^^ ℎ ^ 3^^^ ℎ ^ 3^ ℎ ^(4)ˆ1×6ℎ = 2 ℎ 2 ℎ − 2 −2ℎ − 2 ℎ 2 ℎ 2 .Это приблизительно втрое увеличивает объём вычислений в расчёте наодин шаг сетки ℎ, однако позволяет повысить общую скорость моделирования благодаря увеличению ℎ (сокращению числа шагов).4.5.2. Переход к представлению взаимодействияРассмотрим ещё один подход к построению абсолютно устойчивойчисленной схемы для интегрирования нелинейного уравнения Шрёдингера (112). Данный метод [A16] аналогичен переходу к, или представлению Дирака, используемому в квантовой механике наряду с представлениями Шрёдингера и Гейзенберга.Основная идея, лежащая в основе данного метода: выполнить замену переменной и перейти к новой искомой функции, которая бы непредставлениювзаимодействия96ˆ:эволюционировала под действием дисперсионного оператора35 (︀)︀ˆ , = exp − ( − 0 )(133)где — решение нелинейного уравнения Шрёдингера в представлениивзаимодействия, 0 = const — точка, в которой представления совпадают.