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

[arseniiv]

Я перепроверю вывод, но $i$ там, вроде, должно оставаться.

Спасибо, Mysterious Light! Кстати, в конце вашей анимации по мотивам моего кода на последних кадрах как раз видно, как модель разносит.

Пока поделюсь функцией раскраски — по-моему, она не очень ужасно отражает и фазу, и модуль:

AmplitudeCF[zoom_][z_] := Hue[Arg[z]/6.283185307179586`, Abs@Tanh[Abs[z]/zoom], .5 Tanh[.2 Abs[z]/zoom] + .5]

Этот код, наверно, тоже стоит поускорять, потому что он тоже выполняется $n^2$ раз…

Анимации будут! Сначала учту изменения.

[arseniiv]

Оказалось, что $\Delta x^2$ должны быть в знаменателях. Мне повезло, я никогда не менял $\Delta x = 1$.

Новый код:

MakeSchrodingerTimeStep[V_, m_, \[HBar]_, \[CapitalDelta]t_, \[CapitalDelta]x_] := Module[{a1, V1},
  a1 = I \[HBar] \[CapitalDelta]t/(2 m \[CapitalDelta]x^2);
  V1 = 1 - 2 I \[HBar] \[CapitalDelta]t/(m \[CapitalDelta]x^2) - (I \[CapitalDelta]t)/\[HBar] V;
  Function[{\[Psi]},
    ArrayPad[ListConvolve[{{0, a1, 0}, {a1, 0, a1}, {0, a1, 0}}, \[Psi]], 1] + V1 \[Psi]
  ]
]

Новый пример:

(* Улучшил цвет и немного ускорил *)
AmplitudeCF[zoom_] := Module[{f},
  f = Compile[{{z, _Complex}},
    Block[{absz = Abs[z] 1./zoom},
      {Arg[z]/6.283185307179586`, Tanh[absz], If[absz > 10., 0., 2. - Cosh[.6 absz]]}
    ]
  ];
  Hue @@ f[#] &
]
(* Пример *)
n = 100;
V = Array[0 &, {n, n}];
step = MakeSchrodingerTimeStep[V, 1., 1., .2, 1.];
\[Psi] = Array[Exp[-I 2 \[Pi] 6 #2/n] Exp[-((#1 - .5 n)^2 + (#2 - .5 n)^2)/(.3 n)]/Sqrt[2. \[Pi] (.3 n)] &, {n, n}];
First@Timing[steps = NestList[step, \[Psi], 180]] (* 8.469 *)
someSteps = steps[[;; ;; 5]]; (* отображение отнимает много времени *)
First@Timing[anim = ListAnimate[
  frames = ArrayPlot[#, ColorFunction -> AmplitudeCF[.1], ColorFunctionScaling -> False] & /@ someSteps,
  AnimationRunning -> False]
] (* 29.734 — и это пятая часть! *)
anim

Увы, разрисовывание получается дольше, чем вычисления. Без компиляции было ещё больше. Я и к вычислениям её собирался применить, но не понял, что было не так.

Пока это пример свободной частицы (хм, надо уже сделать преобразование от $[-1;1]^2$ в $1..n^2$)

$$\psi(y,x,t=0) = \frac1{\sqrt{2\pi\cdot0{,}3n}}\exp\left(-\frac{(x-n/2)^2 + (y-n/2)^2}{0{,}3n}\right)\exp(-12\pi iy/n),$$

и в конце можно наблюдать, как замечательно модель расходится:

Mysterious Light в сообщении #897942 писал(а):

P.S. Вам нулевые краевые условия принципиальны?

Даже не знаю. Можно вместо них поставить за краями бесконечно большой потенциал; надо будет посчитать, как изменить счёт (UPD: пока возился, подумал; не понял). Но на все случаи это не подойдёт.

-- Пт авг 22, 2014 19:54:24 --

(Белый — ноль, чёрный — значения модуля больше какого-то, посередине цветное в соответствии с аргументом.)

[arseniiv]

Текст примера сменим на более человекопонятный:

TableXYN[expr_, n_] := Table[expr, {y, 1, -1, -2/(n - 1)}, {x, -1, 1, 2/(n - 1)}] (* транспонирование и обратный порядок — с учётом отображения потом *)
n = 100;
V = TableXYN[0, n];
step = MakeSchrodingerTimeStep[V, 1., 1., .05, 1.];
\[Psi] = TableXYN[Exp[I 2 \[Pi] 4. x] Exp[-(x^2 + y^2)/.005]/Sqrt[2 \[Pi] .01], n];
First@Timing[steps = NestList[step, \[Psi], 600]]
someSteps = steps[[;; ;; 10]];
First@Timing[anim = ListAnimate[
  frames = ArrayPlot[#, ColorFunction -> AmplitudeCF[2.5], ColorFunctionScaling -> False] & /@ someSteps,
  AnimationRunning -> False]
]
anim

$$\psi(x,y,t=0) = \frac1{\sqrt{0{,}004\pi}}\exp\left(-(x^2 + y^2)/0{,}001 + 8\pi ix\right)$$

-- Пт авг 22, 2014 20:52:03 --

А вот при тех же параметрах и потенциале $10(x^2 + y^2)$ модель опять начинает расходиться.

-- Пт авг 22, 2014 20:52:58 --

На самом интересном месте. Что ж, надо будет ещё поиграть цифрами…

-- Пт авг 22, 2014 20:55:43 --

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

[Vince Diesel]

Для уравнения теплопроводности условие сходимости для явных схем, насколько я помню, это шаг по времени меньше чем шаг по $x$ в квадрате пополам. Не знаю, применимо ли это к уравнению Шредингера. Еще коэффициент при старших производных надо. А насчет колебаний, их не должно там быть? решение же нестационарное.

[arseniiv]

Те колебания вылезают с какой-то странной стороны. Например, при больших $\Delta t$ они запросто возникали в том месте, откуда волновой пакет практически ушёл. Явно видно как нарушается унитарность (можно вычислить Norm[Flatten[\[Psi]]] и сравнить с тем, что получается во время «зачернения» — отличия, по идее, должны быть не меньше какой-то зависящей от схемы и номера шага по времени — наверно, линейной по последнему).

[sithif]

arseniiv, ищите информацию по спектральному методу анализа устойчивости разностных схем. Например, эта ссылка.

[Xaositect]

У меня достаточно долго не расходится при $\Delta x = 10^{-2}$, $\Delta t = 10^{-7}$.

Расплывание колокола.
http://s9.postimg.org/to0bsj1h9/0_250_V0.gif
Пакет на торе
http://s29.postimg.org/a9qt2ln5x/1_250_V0.gif
Отражение от нулевого краевого условия (не знаю, насколько это правда, а на сколько артефакт численного метода)
http://s30.postimg.org/mxs90nskf/1_250_V0_Z.gif

И с потенциалом ($5000 r^2$)
Стоячая волна
http://s14.postimg.org/emqtnxbkf/0_250_V50.gif
Пакет
http://s30.postimg.org/755o35z6n/1_1000_V50_Z.gif
http://s1.postimg.org/m0vvtfi25/1_1000_V50_D.gif

[Munin]

Собственно, дальше картинки, которые хотелось бы увидеть:

Цитата:

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

[arseniiv]

Картинки классные, только раскраска фазы не совсем понятна.

[Xaositect]

arseniiv в сообщении #898771 писал(а):

Картинки классные, только раскраска фазы не совсем понятна.

Кажется, это кривая функция HSV->RGB, я не стал разбираться.

-- Сб авг 23, 2014 16:25:49 --

Да, это не было не в Mathematica, если что. И параметры начального условия немного отличаются от Вашего.

[Xaositect]

Вот более аккуратная картинка. Фаза кодируется оттенком (0 - красный, $\pi$ - циан), модуль - насыщенностью.
http://s1.postimg.org/ihzu0wk65/test.gif

Начальное условие $\phi_{0.1}(r) \exp(10x \pi i)$, масса $1$, потенциал $5000r^2$. ($\hbar = 1$). Квадрат $[-1, 1]^2$, Сетка $200\texttt{x}200$, $\Delta t = 10^{-7}$, время анимации $0.1$.

Вскоре после этого все-таки расходится, надо Рунге-Кутта попробовать написать.

[arseniiv]

(Оффтоп)

Если о функциях сопоставления $\mathbb C$ цветам, я бы ещё попытался цвета одного модуля сделать более одинаковой Y, надо будет поиграть со своей…

[arseniiv]

(Оффтоп)

Ай, ерунда получается.

[Утундрий]

Vince Diesel в сообщении #898436 писал(а):

Для уравнения теплопроводности условие сходимости для явных схем, насколько я помню, это шаг по времени меньше чем шаг по $x$ в квадрате пополам. Не знаю, применимо ли это к уравнению Шредингера.

Не применимо. Для УШ вообще любая явная схема неустойчива. Как считать, да чтобы не разваливалось, показано в Гугле и Самарском ( А.А. Самарский "Теория разностных схем", глава 5, §4 Нестационарное уравнение Шрёдингера).

[Munin]

(Оффтоп)

Хорошие соавторы Гугл и Самарский...