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

[arseniiv]

Взяв двумерное уравнение Шрёдингера с потенциалом

$$i\hbar\psi'_t = V\psi - \frac{\hbar^2}{2m}(\psi''_{xx} + \psi''_{yy}})$$

и поиграв в конечные разности, написал такое:

MakeSchrodingerTimeStep[V_, m_, \[HBar]_, \[CapitalDelta]t_, \[CapitalDelta]x_] := Module[
  {a1, a2, V1},
  a1 = I \[HBar] \[CapitalDelta]x^2 \[CapitalDelta]t/(2 m);
  a2 = 2 I \[HBar] \[CapitalDelta]x^2 \[CapitalDelta]t/m;
  V1 = (-((I \[CapitalDelta]t)/\[HBar]) # &) /@ V;
  Function[\[Psi],
   Block[{xnm1, ynm1, \[Psi]1},
    \[Psi]1 = Array[0 &, Dimensions[\[Psi]]];
    {xnm1, ynm1} = Dimensions[\[Psi]] - 1;
    Do[
     \[Psi]1[[xi, yi]] = a1 (\[Psi][[xi - 1, yi]] + \[Psi][[xi + 1, yi]] + \[Psi][[xi, yi - 1]] + \[Psi][[xi, yi + 1]]) + \[Psi][[xi, yi]] (1 - a2 + V1[[xi, yi]]),
     {xi, 2, xnm1}, {yi, 2, ynm1}];
    \[Psi]1
    ]
   ]
  ]

И используется оно примерно так:

n = 100;
V = Array[5 ((#1/n - .5)^2 + (#2/n - .5)^2) &, {n, n}];
step = MakeSchrodingerTimeStep[V, 1., 1., .3, 1.]; (* V, m, \[HBar], \[CapitalDelta]t, \[CapitalDelta]x *)
\[Psi] = Array[Sin[6.28 4 #2/n] Exp[-((#1 - .3 n)^2 + (#2 - .3 n)^2)/(.3 n)]/Sqrt[2. \[Pi] (.3 n)] &, {n, n}];
grs = {};
Do[
 AppendTo[grs,
  ArrayPlot[\[Psi], ColorFunctionScaling -> False, 
   ColorFunction -> AmplitudeCF[.05] (* моя *)]];
 \[Psi] = step[\[Psi]],
 {80}]
ListAnimate[grs, AnimationRunning -> False]

Как бы оптимизировать? Долго слишком (несколько минут приведённый пример работает). Да ещё и, боюсь, я написал что-то неправильно. Ну и в общем код может быть где-то неестественным — сам себя не поправишь, приму все замечания.

-- Ср авг 20, 2014 03:50:36 --

Кстати, разбирающиеся в численных методах, может быть, подскажут, какие соотношения должны тут выполняться между параметрами, чтобы не было того, что тут появляется на последних шагах по времени? Что-то я как будто соблюдал, но, видимо, или не то, или не так, или того недостаточно.

[Aritaborian]

В уравнении Шрёдингера и численных методах я разбираюсь чуть лучше, чем в разведении страусов и древней истории Мадагаскара, но позволю себе наивный вопрос (раз уж «чуть лучше» ;-) Зачем вы используете процедурное программирование (Do)? Здесь же наверняка можно применить что-то нативное для Wolfram Language. Возможно, конструкцию Sow...Reap.

[arseniiv]

С нею как раз беда, этот коан я так ещё и не раскусил. Да и выполнится ли она быстрее? (одно ядро с Hyper-Treading). Хотя ещё есть MapIndexed, который тут вписывается… Проверю его. Вру.

Видимо, придётся писать потом на другом языке.

[Aritaborian]

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

Видимо, придётся писать потом на другом языке.

Не смейте предавать Математику! ;-) Мы же с вами прекрасно знаем, что «Нет ничего невозможного» — лозунг, лучше всего подходящий именно этой системе. Что-нибудь придумается, но нужно подумать.

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

С нею как раз беда, этот коан я так ещё и не раскусил.

Честно признаться, я тоже, хотя были задачи, где эта конструкция подошла бы лучше всего.

[arseniiv]

Aritaborian в сообщении #897860 писал(а):

Что-нибудь придумается, но нужно подумать.

Что я ещё вижу — компиляцию, но она, вроде, для функции, преобразующей матрицу, хоть и чисел, не допускается.

Aritaborian в сообщении #897860 писал(а):

Не смейте предавать Математику!

У меня есть C# и MathNet.Numerics — это серьёзное оправдание!

Нестабильность после 60-80-й итерации весьма удручает. Как там всё же надо анализировать разностные схемы?.. Помогает ли повышение порядка?

[Aritaborian]

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

Что я ещё вижу — компиляцию, но она, вроде, для функции, преобразующей матрицу, хоть и чисел, не допускается.

Никогда серьёзно не вникал в Compile (не нужно было), но теперь прочёл в справке:

Цитата:

Compile handles numerical functions, matrix operations, procedural programming constructs, list manipulation functions, and functional programming constructs, etc.

Или это всё не о том?

[Mysterious Light]

arseniiv, я вижу несколько способов ускорить Ваш код:
1. Использовать Compile для частичной компиляции кода.
2. Переписать в другом стиле (с императивного на функциональное или др.), другими методами.
3. Использовать мощный инструмент символьных вычислений.
С последним я ещё не разобрался, в первый — неинтересно. Оставляю второй вариант, который у меня работает в 10-15 раз шустрее Вашего:

Step[\[Psi]_] := ArrayPad[a1 ListConvolve[{
       {0, 1, 0},
       {1, 0, 1},
       {0, 1, 0}
      } , \[Psi]], 1] + \[Psi] (1 - a2 + V1);

(полностью)

m = 1.;
\[HBar] = 1.;
\[CapitalDelta]t = .3;
\[CapitalDelta]x = 1.;

n = 100;
V = Array[5 ((#1/n - .5)^2 + (#2/n - .5)^2) &, {n, n}];
\[Psi]0 = 
  Array[Sin[6.28 4 #2/n] Exp[-((#1 - .3 n)^2 + (#2 - .3 n)^2)/(.3 n)]/
      Sqrt[2. \[Pi] (.3 n)] &, {n, n}];

a1 = (I \[HBar] \[CapitalDelta]x^2 \[CapitalDelta]t)/(2 m);
a2 = 2 I \[HBar] \[CapitalDelta]x^2 \[CapitalDelta]t/m;
V1 = -((I \[CapitalDelta]t)/\[HBar]) # & /@ V;

Step[\[Psi]_] := ArrayPad[a1 ListConvolve[{
       {0, 1, 0},
       {1, 0, 1},
       {0, 1, 0}
      } , \[Psi]], 1] + \[Psi] (1 - a2 + V1);

First[Timing[grp = NestList[Step, \[Psi]0, 80]]]
ListAnimate[ArrayPlot[#, ColorFunctionScaling -> False] & /@ grp]

Все параметры всегда можно замкнуть в Module-Function, на скорости не сказывается.
Обратите внимание, Times автоматически отображается на списки: Times[list1, list2] == MapThread[Times, {list1, list2}].
Также можно менять ядро свёртки, тем самым меняя разностную схему.

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

[Munin]

А как насчёт результаты повыкладывать? :-)

[Mysterious Light]

Munin, если Вы мне замечание сделали, то я подкреплю циферками:

(преамбула)

m = 1.;
\[HBar] = 1.;
\[CapitalDelta]t = .3;
\[CapitalDelta]x = 1.;
n = 100;
V = Array[5 ((#1/n - .5)^2 + (#2/n - .5)^2) &, {n, n}];
\[Psi]0 = 
  Array[Sin[6.28 4 #2/n] Exp[-((#1 - .3 n)^2 + (#2 - .3 n)^2)/(.3 n)]/
      Sqrt[2. \[Pi] (.3 n)] &, {n, n}];
a1 = I \[HBar] \[CapitalDelta]x^2 \[CapitalDelta]t/(2 m);
a2 = 2 I \[HBar] \[CapitalDelta]x^2 \[CapitalDelta]t/m;
V1 = (-((I \[CapitalDelta]t)/\[HBar]) # &) /@ V;
s = 80; (* число шагов *)

Исходный код

(Оффтоп)

Step = Function[\[Psi],
   Block[{xnm1, ynm1, \[Psi]1},
    \[Psi]1 = Array[0 &, Dimensions[\[Psi]]];
    {xnm1, ynm1} = Dimensions[\[Psi]] - 1;
    Do[
     \[Psi]1[[xi, yi]] = 
      a1 (\[Psi][[xi - 1, yi]] + \[Psi][[xi + 1, yi]] + \[Psi][[xi, 
            yi - 1]] + \[Psi][[xi, yi + 1]]) + \[Psi][[xi, yi]] (1 - 
          a2 + V1[[xi, yi]]),
     {xi, 2, xnm1}, {yi, 2, ynm1}];
    \[Psi]1]
   ];

\[Psi] = \[Psi]0;
grs = {};
First[Timing[Do[AppendTo[grs, \[Psi]]; \[Psi] = Step[\[Psi]], {s}]]]

Out= 11.6287

Заменим императивный цикл Do-Append-Set на Nest

(Оффтоп)

First[Timing[grs = NestList[Step, \[Psi]0, s]]]

Out= 11.6647

Заменим в исходном коде Function на правило SetDelay

(Оффтоп)

Step[\[Psi]_] :=
  Block[{xnm1, ynm1, \[Psi]1},
   \[Psi]1 = Array[0 &, Dimensions[\[Psi]]];
   {xnm1, ynm1} = Dimensions[\[Psi]] - 1;
   Do[
    \[Psi]1[[xi, yi]] = 
     a1 (\[Psi][[xi - 1, yi]] + \[Psi][[xi + 1, yi]] + \[Psi][[xi, 
           yi - 1]] + \[Psi][[xi, yi + 1]]) + \[Psi][[xi, yi]] (1 - 
         a2 + V1[[xi, yi]]),
    {xi, 2, xnm1}, {yi, 2, ynm1}];
   \[Psi]1];

\[Psi] = \[Psi]0;
grs = {};
First[Timing[Do[AppendTo[grs, \[Psi]]; \[Psi] = Step[\[Psi]], {s}]]]

Out= 11.4567

Вывод: от этих двух факторов ничего не зависит!

А вот теперь попробуем переписать тело функции.

(Оффтоп)

Step[\[Psi]_] := 
  ArrayPad[a1 ListConvolve[{{0, 1, 0}, {1, 0, 1}, {0, 1, 0}}, \[Psi]],
     1] + \[Psi] (1 - a2 + V1);

\[Psi] = \[Psi]0;
grs = {};
First[Timing[Do[AppendTo[grs, \[Psi]]; \[Psi] = Step[\[Psi]], {s}]]]

Out= 1.80411

Профит налицо, на порядок.

Хотя как по мне, основное достоинство последнего кода в том, что вся схема компактно выражена без излишеств. В частности, хорошо обозревается структура конечно-разностной схемы.

[Munin]

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

[Mysterious Light]

Позволю себе вставить картиночку по мотивам кода ТС. Единственное, у ТС разукраска немножко другая, собственная.

Размер: 100 х 100 точек, 80 кадров
У меня считалось полторы минуты.

[Aritaborian]

Mysterious Light, GIF-анимация на форуме не отображается, имейте в виду.
(Но можно кликнуть ;-)

[Munin]

Mysterious Light
Я так понимаю, это свободная частица и неподвижный пакет?

[Mysterious Light]

Aritaborian в сообщении #898001 писал(а):

Mysterious Light, GIF-анимация на форуме не отображается, имейте в виду.
(Но можно кликнуть

Спасибо, я уж думал было, что это мой браузер глючит: по ссылке показывает гиф-анимацию, а на форуме — нет.

Munin, за физической интерпретацией обращаться к ТС.
Насколько я вижу, решается вышеуказанное уравнение Ш. с начальным условием

$$\psi_0(x,y) = \sin\left(2\pi \frac{4y}n\right) \exp\left(-\frac{(x-0.3n)^2 + (y-0.3n)^2}{0.3 n}\right)$$

для потенциала

$$V(x,y) = 5 \left(\frac xn-\frac12\right)^2 + 5 \left(\frac yn-\frac12\right)^2$$

[Vince Diesel]

Чисто экспериментально, если убрать $i$ в поределении V1, и уменьшить шаг по $t$ в $4$ раза, то облако движется в точку минимума потенциала, где и застывает.