[Mathematica 8] Оптимизация модели ур-я Шрёдингера

[Утундрий]

Трёхслойка как-то идеологически сложнее двухслойки. Я бы порекомендовал играющим играться примерно в следующем направлении: одномерная без потенциальной энергии --> одномерная с потенциальной энергией --> двумерная без потенциальной энергии -->… и т.д.

Для любителей острых ощущений - посчитать уравнение Дирака)

[Munin]

Утундрий в сообщении #902194 писал(а):

посчитать уравнение Дирака)

2+1. Ага. Его ещё написать надо...

[arseniiv]

А что, 3+1 так неочевидно снижается до 2+1? (Не пробовал.) (Я всё ещё оффлайн.)

[Munin]

С векторами-тензорами - очевидно. Но уравнение Дирака - это спиноры :-)

[arseniiv]

А, я почему-то перепутал его с Клейном—Гордоном…

[kw_artem]

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

(Оффтоп)

TMA[V_, m_, \[HBar]_, \[CapitalDelta]t_, \[CapitalDelta]x_] := 
 Module[{a, b, P, Q, d, ans}, a = \[HBar]/(4 m \[CapitalDelta]x^2);
  b = V/\[HBar] + 2 a - I/\[CapitalDelta]t;
  Function[{\[Psi] }, d = ListConvolve[{1, -2, 1}, \[Psi]]; 
   PrependTo[d, d[[1]]]; AppendTo[d, d[[-1]]];
   d = I \[Psi]/\[CapitalDelta]t - a d; P = Array[0 &, n]; Q = P; 
   ans = Array[0 &, n + 1];
   P[[1]] = a/b[[1]]; Q[[1]] = -d[[1]]/b[[1]];
   Do[P[[i]] = a/(b[[i]] - a P[[i - 1]]), {i, 2, n - 1}];
   Do[Q[[i]] = (a Q[[i - 1]] - d[[i]])/(b[[i]] - a P[[i - 1]]), {i, 2,
      n}];
   Do[ans[[i]] = P[[i]] ans[[i + 1]] + Q[[i]], {i, n, 1, -1}]; 
   Most[ans]]]

(*Пример*)
n = 100;
V = Array[0&, n];

step = TMA[V, 100., 1., .1, .1];
\[Psi]0 = 
 Array[Exp[I 2 \[Pi] 6 #/40] Exp[-(# - 100)^2/(40)]/
     Sqrt[2. \[Pi] (40)] &, 
  n];
First@Timing[steps = NestList[step, \[Psi]0, p = 1000]] (*8.469*)
someSteps = Abs /@ steps[[;; ;; Floor[p/200]]];
First@Timing[
  anim = ListAnimate[
    frames = 
     ListPlot[{#, k ArrayPad[V, -1]}, Joined -> True, 
        ColorFunction -> Automatic, PlotRange -> 0.1] & /@ someSteps, 
    AnimationRunning -> False, DefaultDuration -> 13]] 
anim

,а на основе этой одномерной -- двухмерную (с "переменными направлениями"):

(Оффтоп)

MakeSchrodingerTimeStep[V_, 
  m_, \[CapitalDelta]t_, \[CapitalDelta]x_] := 
 Module[{a, b, tmp}, a = \[HBar]/(4 m \[CapitalDelta]x^2); 
  b = V/\[HBar] + 2 a - I/\[CapitalDelta]t; tmp = Array[0 &, {n, n}];
  Function[{\[Psi]}, 
   Do[tmp[[i]] = 
     LineTMA[\[Psi][[i]], a, b[[i]], \[CapitalDelta]t], {i, 1, n}]; 
   Do[tmp[[All, i]] = 
     LineTMA[tmp[[All, i]], a, b[[All, i]], \[CapitalDelta]t], {i, 1, 
     n}]; tmp]]

LineTMA := Module[{P, Q, d, ans},
  Function[{\[Psi], a, b, \[CapitalDelta]t}, 
   d = ListConvolve[{1, -2, 1}, \[Psi]]; PrependTo[d, d[[1]]]; 
   AppendTo[d, d[[-1]]];
   d = I \[Psi]/\[CapitalDelta]t - a d; P = Array[0 &, n]; Q = P; 
   ans = Array[0 &, n + 1];
   P[[1]] = a/b[[1]]; Q[[1]] = -d[[1]]/b[[1]];
   Do[P[[i]] = a/(b[[i]] - a P[[i - 1]]), {i, 2, n - 1}];
   Do[Q[[i]] = (a Q[[i - 1]] - d[[i]])/(b[[i]] - a P[[i - 1]]), {i, 2,
      n}];
   Do[ans[[i]] = P[[i]] ans[[i + 1]] + Q[[i]], {i, n, 1, -1}]; 
   Most[ans]]]

\[HBar] = 1.; n = 100;
V = Array[0 &, {n, n}];
step = MakeSchrodingerTimeStep[V, 100., 1., .1];
\[Sigma] = 1.5; p0 = .75;
\[Psi]0 = 
 Array[0.5 Exp[
      I p0 #2/\[HBar] - ((#1 - 50)^2 + (#2 - 25)^2)/
         4/\[Sigma]^2]/(2 \[Pi] \[Sigma]^2)^(1/4) &, {n, n}];
First@Timing[steps = NestList[step, \[Psi]0, p = 100]] (*8.469*)
someSteps = steps[[;; ;; Ceiling[p/100]]];
Norm[Flatten[#]] & /@ {\[Psi]0, Last@steps}

AmplitudeCF[zoom_] := 
 Module[{f}, 
  f = Compile[{{z, _Complex}}, 
    Block[{absz = Abs[z] 1./zoom}, {Arg[z]/6.283185307179586`, 
      1 - Tanh[0.5 absz], Tanh[2 absz]}]];
  Hue @@ f[#] &]

First@Timing[
  anim = ListAnimate[
    frames = 
     ArrayPlot[#, ColorFunction -> AmplitudeCF[.1], 
        ColorFunctionScaling -> False] & /@ someSteps, 
    AnimationRunning -> False, DefaultDuration -> 7]]
anim 

Разные начальные условия и вид потенциала:

V = Array[0 &, {n, n}];(*Array[0&,{n,n}];*)
(*Array[If[((#1-50)^2+(#2-50)^2)<25^2,0,5]&,{n,n}];*)(*Array[0.002((#\
1-50)^2+(#2-50)^2)&,{n,n}];*)(*5(Array[If[49<#2<52,1,0]&,{n,n}]-\
ArrayPad[{{1,1},{1,1},{1,1},{0,0},{0,0},{0,0},{0,0},{1,1},{1,1},{1,1}}\
,{{n/2-5},{n/2-1}}]);*)
(*Array[If[n/2<#<n/2+5,1.4,0]&,n+2];*)(*5(Array[If[49<#2<52,1,0]&,{n,\
n}]-ArrayPad[{{1,1},{1,1}(*,{1,1},{1,1},{1,1},{1,1}*)},{{n/2-1},{n/2-\
1}}]);*)

\[Psi]0 = 
  Array[Exp[I 1 #2] + Exp[-I 1 #2] &, {n, 
    n}];(*Array[wk/.{x\[Rule]#2,y\[Rule]#1}&,{n,n}];*)(*Array[Exp[I 2 \
\[Pi] 3 #2/40] Exp[-((#1-50)^2+(#2-50)^2)/(2 20)]/Sqrt[2. \[Pi] \
(20)]&,{n,n}];*)
(*Array[0.5Exp[I p0 #2/\[HBar]-((#1-50)^2+(#2-25)^2)/4/\[Sigma]^2]/(2\
\[Pi] \[Sigma]^2)^(1/4)&,{n,n}];*)

Последняя, как предполагалось, считает дольше, но может можно код как-то оптимизировать, например, через функцию Compile, с ней я не очень разобрался.

Картинки выложу чуть позже.

[kw_artem]

Забыл сказать, модель сделана по схеме:

$$i\left( {\psi - \psi ^{ - t} } \right) + \sigma \left( {\psi _{ + x} - 2\psi + \psi _{ - x} } \right) = 0$$

[Vince Diesel]

А команда NDSolve краевые задачи для уравнения Шредингера разве не считает?

[kw_artem]

Обещанные анимации (swf):
Свободная частица: анимация;
Странное взаимодействие ВФ с границей области сетки, приводящее к ее отражению: анимация;
Колебания ВФ внутри квадратичного потенциала: анимация;
Столкновение с потенциальной стенкой: анимация;
То же, но с повышенной чувствительностью к амплитуде ВФ: анимация;
Столкновение с отн. невысокой потенциальной стеной: анимация;
Падение с потенциальной стенки: анимация;
Имитация двухщелевого эксперимента. Картинка получилась очень даже похожей на ту, что приводил Munin. анимация;
Внутри круговой ямы с вертикальными стенками: анимация;
Рассеяние на круглом барьере: маленький импульс налетающей частицы -- анимация, большой -- анимация.

Хоть модель построена так, чтобы максимально сохранялась норма, но в некоторых примерах, она все равно сильно изменяется (до 25%), например, в квадратичном потенциале.

-- 30.04.2015, 22:54 --

Vince Diesel, надо будет попробовать! Но самому интереснее!

[Munin]

А как бы их скачать? :-)

[Утундрий]

kw_artem в сообщении #1009716 писал(а):

Хоть модель построена так, чтобы максимально сохранялась норма, но в некоторых примерах, она все равно сильно изменяется (до 25%), например, в квадратичном потенциале.

Вообще-то, нет. У вас простая неявная схема, а она унитарность не сохраняет. То что вам нужно идёт у Кунина начиная с формулы $(7.30)$.

[kw_artem]

Утундрий в сообщении #1009728 писал(а):

kw_artem в сообщении #1009716 писал(а):

Хоть модель построена так, чтобы максимально сохранялась норма, но в некоторых примерах, она все равно сильно изменяется (до 25%), например, в квадратичном потенциале.

Вообще-то, нет. У вас простая неявная схема, а она унитарность не сохраняет. То что вам нужно идёт у Кунина начиная с формулы $(7.30)$.

Подумал, что не так, а оказалось, я неправильно сообщил какая у меня схема. А схема вот какая (точно!):

$$i\left( {\psi - \psi ^{ - t} } \right) + \frac 1 2 \sigma \left(\left( {\psi _{ + x} - 2\psi + \psi _{ - x} } \right) + \left({\psi^{-t} _{ + x} - 2\psi^{-t} + \psi^{-t} _{ - x} \right)}\right) = 0$$

-- 30.04.2015, 23:55 --

Вот https://yadi.sk/d/-Nwxp_1SgMorB. Прога для просмотра тоже в архиве.

[Munin]

Чего осталось:

Munin в сообщении #898662 писал(а):

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

[arseniiv]

kw_artem
Как хорошо, что вы продолжили — и сделали нормально! По крайней мере, много нормальнее, чем был мой вклад в этой теме. Да ещё и кодом поделились. (Потом поизучаю!)

[ArtDen]

kw_artem, вообще здорово, но у меня пара вопросов.
Во первых, почему свободная частица расходится? В той теме мне удалось добиться для свободной частицы, чтобы она пролетала от одной границы поля до другой без искажений. Правда, я там забыл похвастаться анимацией ))

Во вторых, как долго всё это считается? Мне удалось сделать расчёт практически в реальном времени для поля 1024x1024. И кстати, методы из книжки Кунина в моём случае не работали на больших сетках из-за очень большой погрешности. Правда, максимум, что я осилил - это неявная схема, решаемая при помощи метода редукции.