1
$\begingroup$

I don't know Mathematica at all, but I'm interesting in replicating some code to reproduce this animated gif.

enter image description here

Mathematica code is this:

(* Source code written in Mathematica 6.0 by Steve Byrnes, Feb. 2011. This source code is public domain. *)
(* Shows classical and quantum trajectory animations for a harmonic potential. Assume m=w=hbar=1. *)
ClearAll["Global`*"]
(*** Wavefunctions of the energy eigenstates ***)
psi[n_, x_] := (2^n*n!)^(-1/2)*Pi^(-1/4)*Exp[-x^2/2]*HermiteH[n, x];
energy[n_] := n + 1/2;
psit[n_, x_, t_] := psi[n, x] Exp[-I*energy[n]*t];
(*** A random time-dependent state ***)
SeedRandom[1];
CoefList = Table[Random[]*Exp[2 Pi I Random[]], {n, 0, 4}];
CoefList = CoefList/Norm[CoefList];
Randpsi[x_, t_] := Sum[CoefList[[n + 1]]*psit[n, x, t], {n, 0, 4}];
(*** A coherent state (or "Glauber state") ***)
CoherentState[b_, x_, t_] := Exp[-Abs[b]^2/2] Sum[b^n*(n!)^(-1/2)*psit[n, x, t], {n, 0, 15}];
(*** Make the classical plots...a red ball anchored to the origin by a gray spring. ***)
classical1[t_, max_] := ListPlot[{{max Cos[t], 0}}, PlotStyle -> Directive[Red, AbsolutePointSize[15]]];
zigzag[x_] := Abs[(x + 0.25) - Round[x + 0.25]] - .25;
spring[x_, left_, right_] := (.9 zigzag[3 (x - left)/(right - left)])/(1 + Abs[right - left]);
classical2[t_, max_] := Plot[spring[x, -5, max Cos[t]], {x, -5, max Cos[t]}, PlotStyle -> Directive[Gray, Thick]];
classical3 = ListPlot[{{-5, 0}}, PlotStyle -> Directive[Black, AbsolutePointSize[7]]];
classical[t_, max_, label_] := Show[classical2[t, max], classical1[t, max], classical3, 
   PlotRange -> {{-5, 5}, {-1, 1}}, Ticks -> None, Axes -> {False, True}, PlotLabel -> label, AxesOrigin -> {0, 0}];
(*** Put all the plots together ***)
SetOptions[Plot, {PlotRange -> {-1, 1}, Ticks -> None, PlotStyle -> {Directive[Thick, Blue], Directive[Thick, Pink]}}];
MakeFrame[t_] := GraphicsGrid[
   {{classical[t + 2, 1.5, "A"], classical[t, 3, "B"]},
    {Plot[{Re[psit[0, x, t]], Im[psit[0, x, t]]}, {x, -5, 5}, PlotLabel -> "C"], 
     Plot[{Re[psit[1, x, t]], Im[psit[1, x, t]]}, {x, -5, 5}, PlotLabel -> "D"]},
    {Plot[{Re[psit[2, x, t]], Im[psit[2, x, t]]}, {x, -5, 5}, PlotLabel -> "E"], 
     Plot[{Re[psit[3, x, t]], Im[psit[3, x, t]]}, {x, -5, 5}, PlotLabel -> "F"]},
    {Plot[{Re[Randpsi[x, t]], Im[Randpsi[x, t]]}, {x, -5, 5}, PlotLabel -> "G"], 
     Plot[{Re[CoherentState[1, x, t]], Im[CoherentState[1, x, t]]}, {x, -5, 5}, PlotLabel -> "H"]}
    }, Frame -> All, ImageSize -> 300];
output = Table[MakeFrame[t], {t, 0, 4 Pi*96/97, 4 Pi/97}];
SetDirectory["C:\\Users\\Steve\\Desktop"]
Export["test.gif", output]

In particular I'm stucked in this line:

psit[n_, x_, t_] := psi[n, x] Exp[-I*energy[n]*t];

First question: I first thought to translate it with a multiplication where the multiplicand is psi[n, x] and the multiplier is Exp[-I*energy[n]*t], because I interpret both the space and the star * with the multiplication but.. is it true or does they have a different meaning?

EDIT

Answer to the first question is (due to comment): space and star have the same meaning, multiplication.

Second question: with f(x, t) as Im[psit[1, x, t]] in my JavaScript translation I have those results:

x:0,  t: 0,   f(x, t):  0
x:-2, t: 0,   f(x, t):  0
x:2,  t: 0,   f(x, t):  0
x:0,  t: 0.5, f(x, t):  -485738.84954226436
x:-2, t: 0.5, f(x, t):  -4.717559812725077e+136
x:2,  t: 0.5, f(x, t):  3.917168783701828e+49
x:0,  t: 1,   f(x, t):  -710819.418779301
x:-2, t: 1,   f(x, t):  -6.903572006434811e+136
x:2,  t: 1,   f(x, t):  5.732297593069263e+49

Can someone plot me results in Mathematica please so I have value to test again?

Just for sake of curiosity this is my JavaScript code

$\endgroup$
  • 1
    $\begingroup$ both the space and the asterisk are multiplication in that example $\endgroup$ – Jason B. May 15 at 22:19
  • $\begingroup$ @JasonB. thanks! can you see the edits? $\endgroup$ – nkint May 15 at 22:28
3
$\begingroup$

Surprisingly, code written for version 6 still works in version 12 (albeit very slowly). Test case:

psi[n_, x_] := (2^n*n!)^(-1/2)*Pi^(-1/4)*Exp[-x^2/2]*HermiteH[n, x];
energy[n_] := n + 1/2;
psit[n_, x_, t_] := psi[n, x] Exp[-I*energy[n]*t];
 Flatten[
 Table[{x, t, Im[psit[1, x, t]]}, {x, -2, 2, 2}, {t, 0, 1, .5}], 1]

(*Out[]= {{-2, 0., 0.}, {-2, 0.5, 0.195985}, {-2, 1., 0.2868}, {0, 0.,
   0.}, {0, 0.5, 0.}, {0, 1., 0.}, {2, 0., 0.}, {2, 
  0.5, -0.195985}, {2, 1., -0.2868}}*)

 % // TableForm

fig1

I changed the color and some parameters

psi[n_, x_] := (2^n*n!)^(-1/2)*Pi^(-1/4)*Exp[-x^2/2]*HermiteH[n, x];
energy[n_] := n + 1/2;
psit[n_, x_, t_] := psi[n, x] Exp[-I*energy[n]*t];

SeedRandom[1];
CoefList = Table[Random[]*Exp[2 Pi I Random[]], {n, 0, 4}];
CoefList = CoefList/Norm[CoefList];
Randpsi[x_, t_] := Sum[CoefList[[n + 1]]*psit[n, x, t], {n, 0, 4}];

CoherentState[b_, x_, t_] := 
  Exp[-Abs[b]^2/2] Sum[b^n*(n!)^(-1/2)*psit[n, x, t], {n, 0, 5}];

classical1[t_, max_] := 
  ListPlot[{{max Cos[t], 0}}, 
   PlotStyle -> 
    Directive[RGBColor[1, 0, 1, .5], AbsolutePointSize[15]]];
zigzag[x_] := Abs[(x + 0.25) - Round[x + 0.25]] - .25;
spring[x_, left_, 
   right_] := (.9 zigzag[3 (x - left)/(right - left)])/(1 + 
     Abs[right - left]);
classical2[t_, max_] := 
  Plot[spring[x, -5, max Cos[t]], {x, -5, max Cos[t]}, 
   PlotStyle -> Directive[Gray, Thick]];
classical3 = 
  ListPlot[{{-5, 0}}, 
   PlotStyle -> Directive[Black, AbsolutePointSize[7]]];
classical[t_, max_, label_] := 
  Show[classical2[t, max], classical1[t, max], classical3, 
   PlotRange -> {{-5, 5}, {-1, 1}}, Ticks -> None, 
   Axes -> {False, True}, PlotLabel -> label, AxesOrigin -> {0, 0}];

SetOptions[
  Plot, {PlotRange -> {-1, 1}, Ticks -> None, 
   PlotStyle -> {Directive[Thick, RGBColor[1, 0, 1, .5]], 
     Directive[Thick, RGBColor[0, 1, 0, .5]]}}];
MakeFrame[t_] := 
  GraphicsGrid[{{classical[t + 2, 1.5, "A"], 
     classical[t, 3, 
      "B"]}, {Plot[{Re[psit[0, x, t]], Im[psit[0, x, t]]}, {x, -5, 5},
       PlotLabel -> "C"], 
     Plot[{Re[psit[1, x, t]], Im[psit[1, x, t]]}, {x, -5, 5}, 
      PlotLabel -> "D"]}, {Plot[{Re[psit[2, x, t]], 
       Im[psit[2, x, t]]}, {x, -5, 5}, PlotLabel -> "E"], 
     Plot[{Re[psit[3, x, t]], Im[psit[3, x, t]]}, {x, -5, 5}, 
      PlotLabel -> "F"]}, {Plot[{Re[Randpsi[x, t]], 
       Im[Randpsi[x, t]]}, {x, -5, 5}, PlotLabel -> "G"], 
     Plot[{Re[CoherentState[1, x, t]], 
       Im[CoherentState[1, x, t]]}, {x, -5, 5}, PlotLabel -> "H"]}}, 
   Frame -> All, ImageSize -> 300];
output = Table[MakeFrame[t], {t, 0, 4 Pi*96/97, 4 Pi/97}];
SetDirectory["C:\\Users\\username\\Desktop"]
Export["test.gif", output, AnimationRepetitions -> Infinity]

fig1

$\endgroup$
  • 1
    $\begingroup$ @b3m2a1 there was a question about a test case. But, then I added the animation by slightly changing the parameters. $\endgroup$ – Alex Trounev May 16 at 13:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.