Free Fall

Maria Daghofer,  Wolfgang von der Linden
Institut fuer Theoretische Physik
Technische Universitaet Graz
A - 8010 Graz, 
    Austria
Projekt Multimediale Lehre des
Zentralen Informatikdienstes der TU Graz

2001 - 12 - 13
Null

Parameters

Begin["fallQM`"]; 
    Clear["fallQM`*"]; 
Needs["Graphics`MultipleListPlot`"];

Approximate Solution:

Size of discrete Hamilton-matrix nd and maximal height to be included in the discrete Matrix xmaxd:

Lowering nd increases discretization error and speed. Too small a value for 
xmaxd gives false eigenvalues and eigenfunctions; very large xmaxd cause 
uninteresting regions to be included in calculation.

Dynamics

Center x0 and width σ0 of initial gaussian wave packet:

x0 should be well below xmaxd!

Number of periods Nper for which the evaluation is to be calculated,  Number of timesteps Nt0 per period, Number of eigenstates Mmaxe included in spectral theorem:

Explanation

Approximate solution

Definition

\(M[n_, xmax_] := 
      Block[{H, Δ, a}, Δ = xmax/n; 
        a = 1/Δ\^2; 
				H = 
          Table[If[i == j, 2  a + \ i\ Δ, 
              If[Abs[i - j] == 1, \(-a\), 0]], {i, 1, n}, {j, 1, 
              n}]
			]\)

Meaning of Parameters:
n:		dimension of Hemline-matrix
xmax:	maximal x-coordinate included in discrete Hamiltonian

Display Matrix

M[10, 1. ] // MatrixForm

Definition:

Ψ[n , xmax , m] plots the (m^th)eigenvector of the discrete matrix (with dimension n and maximal  x-coordinate xmax), Eigen[n ,xmax ] displays a few eigenvalues and expectation values for the x-coordinate and  compares them to the exact solution. (Caution: To execute Eigen, the exact solution must first be defined, see below.)

\(\(Eigen[n_, xmax_] := 
        Block[{x, 
            ew, Δ, Φ}, 
					Δ = xmax/n; 
					eigsy = 
            Eigensystem[N[M[n, xmax], 16]]; 
					ew = 
            eigsy[\([1]\)]; 
					ListPlot[ew, 
            PlotStyle -> {Hue[1]}, 
            PlotLabel -> "Eigenvalues of Discrete Hamiltonmatrix"]; 
					Print[
            TableForm[{{\*""(Subscript[λ, 1]) ="", λ\_1 = \ 
                    N[ew[\([n]\)], 8], "  exakt:", \(-ExactEigenvalue[
                      1]\)}, 
					{\*""(Subscript[λ, 2]) ="", λ\_2 = \ 
                    N[ew[\([n - 1]\)], 8], "  exakt:", \(-ExactEigenvalue[
                      2]\)}, 
					{\*""(Subscript[λ, 3]) ="", λ\_3 = \ 
                    N[ew[\([n - 2]\)], 8], "  exakt:", \(-ExactEigenvalue[
                      3]\)}}]]; 
					x\  = \ 
            Table[i\ Δ, {i, 1, 
                n}]; 
					Φ\  = \ 
            eigsy[\([2, n]\)]; \ \ \ \ \ \ ξ\_1\  = \ 
            x . Abs[Φ]^2/\((Abs[Φ]\  . 
                  Abs[Φ])\); 
					Φ\  = \ eigsy[\([2, n - 1]\)]; ξ\_2\  = \ 
            x . Abs[Φ]^2/\((Abs[Φ]\  . 
                  Abs[Φ])\); \ \ 
					Φ\  = \ eigsy[\([2, n - 2]\)]; ξ\_3\  = \ 
            x . Abs[Φ]^2/\((Abs[Φ]\  . 
                  Abs[Φ])\); \ \ 
					Print[TableForm[{{\*""(Subscript[ξ, 1]) ="", ξ\_1, "  exakt:", 
                  ExactMean[
                    1]}, 
					{\*""(Subscript[ξ, 2])  ="", ξ\_2, "  exakt:", 
                  ExactMean[
                    2]}, 
					{\*""(Subscript[ξ, 3])  ="", ξ\_3, "  exakt:", 
                  ExactMean[
                    3]}}]];
					];\)\)
Ψ[n_, xmax_, m_] := 
        Block[{x, Δ}, Δ = xmax/n; 
          x = Table[i Δ, {i, 1, n}]; 
          t = Eigensystem[
                N[M[n, xmax]]][[2]]; 
p1ymax = 
            Max[t[[m]]^2/Δ]; 
p1xmax = 
            Max[xmax, 8]; 
p1 = 
            ListPlot[Transpose[{x, t[[m]]^2/Δ}], 
              PlotRange -> {{0, p1xmax}, {0, p1ymax}}, 
              AxesLabel -> {"", "ϵ"}, 
              PlotStyle -> {PointSize[0.015], Hue[1]}, 
              DisplayFunction -> 
                Identity];
];

Meaning of parameters:
n:		dimension of Hemline-matrix
Δ:		step size  of finite differences in spacial derivatives
m:		quantum number

Exact solution

Block[{z, y}, 
      DSolve[-y''[z] + z y[z] == 0, y[z], z]]

Block[{z, airygraph}, 
      airygraph = 
        Plot[{AiryAi[z], AiryBi[z]}, {z, -3, 2}, 
          DisplayFunction -> Identity]; 
      Show[airygraph , 
        Graphics[{Text["Ai(z)", {-2.1, 0.4}, 
              TextStyle -> {FontColor -> Hue[0.6]}], 
            Text["Bi(z)", {0.8, 1.4}, 
              TextStyle -> {FontColor -> Hue[1]}]}], 
        PlotRange -> {{-3, 2}, {-0.5, 2}}, 
        DisplayFunction -> $DisplayFunction, 
        AxesLabel -> {"z", ""}];
]

Block[{z, airygraph},airygraph  = 
              Plot[AiryAi[z], {z, -13, 10}, 
                PlotStyle -> {Thickness[0.01], Hue[0.6]}, 
                DisplayFunction -> Identity];
airygraph  = 
              Show[airygraph , 
                Plot[(Sqrt[(1/Pi)]) 
                    Cos[2/3  (-z)^(3/2) - Pi/4]/
                      Abs[z]^0.25, {z, -13, 0}, 
                  PlotStyle -> {Hue[1]}, 
                  DisplayFunction -> Identity]];
Show[airygraph ,Graphics[{Text["Ai(z)", {2.5, 0.2}, 
                        TextStyle -> {FontColor -> Hue[0.6]}],Text["(Sqrt[(1/π)]) cos((2/3) (((-z))^((3/2))) - (π/4))/|z\!\(\( | \^\(1/4\)\)\)",{6, 0.75},TextStyle -> {FontSize -> 8, 
                              FontColor -> Hue[1]}]}],PlotRange -> {{-13, 10}, {-0.5, 
                      1}},DisplayFunction -> $DisplayFunction,AxesLabel -> {"z", ""}];
]

Definition:

\(\(ExactEigenvalue[m_] := 
        Block[{x0, xy}, 
							x0\  = \ 
            N[\((\((2  m - 1\/2)\) \(π\ 3\)\/4)\)\^\(2/3\)]; 
							xy /. 
            FindRoot[
              AiryAi[xy] == 
                0, {xy, \(-x0\)}]
						];\)\)
\(\(ExactMean[m_] := \ 
        Block[{λ, 
            x}, 
					λ = \(-ExactEigenvalue[
                m]\); 
					N[∫\_\(0 \(\(\ \)\(\ \)\)\)\%\(\(∞\)\(\ \ \)\)x\ \ \ \(AiryAi[x - λ]\^2\) \[DifferentialD]x]\/N[∫\_0\%∞\ \ \ AiryAi[x - λ]\^2\ \ \ \[DifferentialD]x]
			\ \ \ \ 	\ \ \ ];\)\)
\(\(Compare[n_, xmax_, nEig_] := 
        Block[{x, λ}, 
								Ψ[n, xmax, 
            n + 1 - nEig]; 
								λ\ \ \ \ \ \ \ \  = \(-ExactEigenvalue[
                nEig]\); 
								norm\  = 
            N[∫\_0\%∞ AiryAi[x - λ]\^2\ \ \ \[DifferentialD]x]; 
								p2\ \ \ \ \  = 
            Plot[AiryAi[x - λ]\^2/norm, {x, 0, p1xmax}, 
              PlotStyle -> {Hue[0.6]}, 
              DisplayFunction -> 
                Identity]; 
								Show[{p1, p2}, 
            AxesLabel -> {"x", \*""|Ψ(x)\!\(\( | \^2\)\)""}, PlotRange -> {{0, p1xmax}, {0, p1ymax}}, 
            DisplayFunction -> $DisplayFunction];
							];\)\)

Meaning of parameters:
m:		quantum number
n:		dimension of Hemline-matrix
xmax:	maximal x-coordinate included in d1scretized region

Groundstate and low-lying excited states

This Section gives a comparison of the results obtained with a discrete  Matrix and the exact results for the lowest three eigenstates.
First a plot  of all eigenvalues plus numerical values of lowest  3 eigenvalues and  corresponding averages < (Subscript[Ψ, n]) | 
        x | (Subscript[Ψ, n])>:

Eigen[nd, xmaxd];

Comparison of wave-functions Finite Differences  <-> Exact solution for quantum numbers n=1,2,3

1) Groundstate

Compare[nd, xmaxd, 1];

2) First excited state

Compare[nd, xmaxd, 2];

3) Second excited state

Compare[nd, xmaxd, 3];

Dynamics

Definition:

dyn[n, xmax, Nt, Nper, x0, σ] computes the animation of a Gaussian wave packet centered at x0 with a  width of σ using a discrete matrix with dimension n and maximal  x-coordinate xmax. Nt timesteps are checkmated (with OneStep) for a total of Nper "periods." A "period" corresponds to the motion  of a classical particle to the ground and back up to its starting point.
DynMean calculates the dynamics of the expectation value ⟨x ⟩.
ExactDyn[n, xmax, Nt, Nper, x0, σ, Mmax] simulates the dynamics using the first Mmax exact eigenfunctions and  eigenvalues, n and xmax are only needed for plot resolution.

\(\(\(OneStep[DM_, x_, xmax_, Φ_, n_] := 
        Block[{p1, p2, xmean, 
            x2mean}, 
								\ \ xmean = \ 
            x\  . Abs[Φ]^2/\((Abs[Φ]\  . 
                  Abs[Φ])\); 
								\ \ x2mean = \ \@\(\((x\ ^2)\) . Abs[Φ]^2/\((Abs[Φ]\  . Abs[Φ])\) - xmean^2\); 
								\ \ app\ \ \ \ \  = \ {{xmean, 0}, {xmean, 1}, {xmean - x2mean, 1}, {xmean + x2mean, 
                1}}; 
								\ \ yyy\ \ \ \ \ \  = \ Join[Transpose[{x, Abs[Φ]\^2}], 
              app]; 
							\ \ \ \ \ \ \ \ \ \ \ ListPlot[yyy, PlotRange -> {{0, xmax}, {0, 1}}, PlotStyle -> {Hue[1]}, 
            PlotJoined -> True, 
            AxesLabel -> {"x", \*""|Φ\!\(\( | \^2\)\) ""}]; 
							\ \ \ \ \ \ \ \ \ \ \ \ DM . Φ];\)\(
\)
    \)\)
\(\(dyn[n_, xmax_, Nt_, Nper_, x0_, σ_] := 
        Block[{Δ = xmax/n, dt = Nper*2 \@ x0/Nt}, 
          H = N[M[n, xmax]]; 
									x = 
            Table[i\ Δ, {i, 1, 
                n}]; 
									Φ = E\^\(-\(\((x - x0)\)\^2\/\(2\ σ\^2\)\)\)\/\((π\ σ)\)\^\(1/4\)\ ; 
									\ DM = 
            MatrixExp[I\ dt\ H]; 
									\ dynpl = 
            Table[\(Φ = 
                  OneStep[DM, x, xmax, Φ, n]\ ;\)\ , {i, 1, 
                Nt}]; 
SelectionMove[EvaluationNotebook[], 
            After, EvaluationCell]; 
SelectionMove[
            EvaluationNotebook[], Next, 
            CellGroup]; 
FrontEndTokenExecute["OpenCloseGroup"];
							];\)\)
\(\(DynMean[n_, xmax_, Nt_, Nper_, x0_, σ_] := 
        Block[{tmax = N[Nper*2 \@ x0], dt = tmax/Nt, 
            x, Δ = N[xmax/n], Φ, DM, tavg, tarr, 
            tp, p1}, 
									H\ \ \ \ \ \ \ \ \ \ \ \ \  = 
            M[n, xmax]; 
									x\ \ \ \ \ \ \ \ \ \ \ \ \  = 
            Table[i\ Δ, {i, 1, 
                n}]; 
									Φ\ \ \ \ \ \ \ \ \ \ \ \ \  = \ E\^\(-\(\((x - x0)\)\^2\/\(2\ σ\^2\)\)\)\/\((π\ σ)\)\^\(1/4\)\ ; 
									DM\ \ \ \ \ \ \ \ \ \  = 
            MatrixExp[I\ dt\ H]; 
									tavg\ \ \ \ \  = 
            Table[Φ = DM . Φ; 
              x\  . Abs[Φ]^2/\((Abs[Φ]\  . 
                    Abs[Φ])\), {i, 1, 
                Nt}]; 
									tarr\ \ \ \ \  = \ 
            Table[it\ dt, {it, 1, 
                Nt}]; 
									tp\ \ \ \ \ \ \ \ \ \  = \ \@x0; 
									p1\ \ \ \ \ \ \ \ \ \  = 
            Plot[x0 - \((Mod[t - tp, 2\ tp] - tp)\)\^2, {t, 0, 
                tmax}, 
								\ \ \ \ \ PlotStyle -> {Hue[0.3, 0.8, 0.5]}, 
              DisplayFunction -> 
                Identity]; 
									p2\ \ \ \ \ \ \ \ \ \ \  = 
            Show[p1, 
              ListPlot[
                Transpose[{tarr, 
                    tavg}], 
							\ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {{0, tmax}, {0, x0*1.1}}, 
                PlotStyle -> {PointSize[0.01], 
                    Hue[1]}, 
								\ \ DisplayFunction -> Identity], \ 
              DisplayFunction -> 
                Identity]; 
									Show[p2, 
            Graphics[
								\ \ \ \ {Text["class.", {tmax*1.05, \((x0 - \((Mod[tmax - tp, 2\ tp] - tp)\)\^2)\)*1.05}, 
								TextStyle -> {FontColor -> Hue[0.3, 0.8, 
                          0.5]}], 
								\ \ Text["QM", {tmax*1.05, 0.95*tavg[\([Nt]\)]}, 
                  TextStyle -> {FontColor -> 
                        Hue[
                          1]}]}], 
								\ \ \ \ \ \ \ PlotRange -> {{0, tmax*1.1}, {0, x0*1.1}}, 
            AxesLabel -> {"t", "<x>(t)"}, 
            DisplayFunction -> $DisplayFunction];
							];\)

\)
\(\(\(ExactDyn[n_, xmax_, Nt_, Nper_, x0_, σ_, 
          Mmax_] := 
							Block[{Δ = 
              xmax/n, tp, tavg, tarr, x, zeroes, normq, coeff, y, 
            sq = 2\ σ\^2, tmax = N[Nper*2 \@ x0], 
            dt = tmax/Nt}, 
									Off[
            NIntegrate::ncvb]; 
									Off[
            General::unfl]; 
									x\ \ \ \ \ \ \ \ \ \ \ \  = 
            Table[i\ Δ, {i, 
                n}]; 
									zeroes = 
            Table[\(-ExactEigenvalue[m]\), {m, 
                Mmax}]; 
									normq\ \ \ \  = \ NIntegrate[
              AiryAi[xx - zeroes[\([1]\)]]\^2, {xx, 
                0, ∞}]\ ;  (*norm\^2\ of\ 1\^st\ eigenfct*) 
									coeff\ \  = \ \ NIntegrate[AiryAi[xx - zeroes[\([1]\)]]\ E\^\(-\(\((xx - x0)\)\^2\/sq\)\), {xx, 0, ∞}]\/\(\(normq\ \((π\ σ\^2)\)\^\(1/4\)\)\(
\)
              \); 
									Do[
										y[m - 1]\  = \ \ AiryAi[x - zeroes[\([m - 1]\)]]\ *\ 
                coeff; 
										normq\ \ \ \ \ \ \ \  = normq\  + \ \ NIntegrate[
                  AiryAi[xx]\^2, {xx, \(-zeroes[\([m]\)]\), \(-zeroes[\([m - 
                            1]\)]\)}]\ ; 
										coeff\ \ \ \ \ \  = \ \ NIntegrate[AiryAi[xx - zeroes[\([m]\)]]\ E\^\(-\(\((xx - x0)\)\^2\/sq\)\), {xx, 0, ∞}]\/\(normq\ \((π\ σ\^2)\)\^\(1/4\)\);
									\ \ , {m, 2, 
              Mmax}]; 
									y[
            
            Mmax]\  = \ \ AiryAi[x - zeroes[\([Mmax]\)]]\ *\ 
              coeff; 
									tavg\  = \ 
            Table[0, {it, 0, 
                Nt}]; 
									yy\ \ \  = \ 
            x*0; 
									Do[
            yy\  = \ yy\  + \ \ \ y[m], {m, 1, 
              Mmax}]; 
										t0\ \ \  = 
            Transpose[{x, 
                Abs[yy]\^2}]; 
										xa\ \ \  = \ x\  . Abs[yy]\^2/\((Abs[yy]\  . 
                  Abs[yy])\); 
										va0\  = \ \@\(\(\(x\)\(\ \)\)\^2 . Abs[yy]\^2/\((Abs[yy]\  . Abs[yy])\) - xa\^2\); 
										ta0\  = \ {{xa, 
                0}
									, {xa, 
                1}}; 
										tv0\  = \ {{xa - 
                  va0, 1}, {xa + va0, 
                1}}; 
									exdynpl = {Null}; 
									Do[
										yy\ \ \  = \ 
              x*0; 
										Do[
              yy\  = \ yy\  + \ \ Exp[\(-I\)\ dt*\ it*
                        zeroes[\([m]\)]]\ *\ \ y[m], {m, 1, 
                Mmax}]; 
										t1\ \ \  = \ 
              Transpose[{x, 
                  Abs[yy]\^2}]; 
										xa\ \ \  = \ x\  . Abs[yy]\^2/\((Abs[yy]\  . 
                    Abs[yy])\); 
										va\ \ \  = \ \@\(\(\(x\)\(\ \)\)\^2 . Abs[yy]\^2/\((Abs[yy]\  . Abs[yy])\) - xa\^2\); 
										ta\ \ \  = \ {{xa, 0}, {xa, 
                  1}}; 
										tv\ \ \  = \ {{xa - va,  .99}, {xa + 
                    va,  .99}}; 
										exdynpl = Append[
                exdynpl, 
										MultipleListPlot[t1, ta, tv, t0, ta0, tv0, 
                  PlotJoined -> 
                    True, 
											PlotStyle -> {Hue[1], Hue[1], Hue[1], Hue[ .6], Hue[ .6], 
                      Hue[ .6]}, 
											SymbolShape -> None, PlotRange -> {{0, xmax}, {0, 1}}, 
                  
                  AxesLabel -> {"x", \*""|Φ\!\(\( | \^2\)\)""}]]; 
										tavg[\([it + 1]\)] = 
              x\  . Abs[yy]\^2/\((Abs[yy]\  . 
                    Abs[yy])\);
									, {it, 0, 
              Nt}]; 
SelectionMove[EvaluationNotebook[], 
            After, EvaluationCell]; 
SelectionMove[
            EvaluationNotebook[], Next, 
            CellGroup]; 
FrontEndTokenExecute["OpenCloseGroup"]; 
									On[NIntegrate::ncvb]; 
          On[General::unfl]; 
									tarr\ \ \ \ \ \  = \ Table[
              it\ dt, {it, 0, 
                Nt}]; 
									tp\ \ \ \ \ \ \ \ \ \  = \ N[\@x0]; 							\ \ 
\ \ \ \ \ \ \ p1\ \ \ \ \ \ \ \ \ \  = 
            Plot[x0 - \((Mod[t - tp, 2\ tp] - tp)\)\^2, {t, 0, 
                tmax}, 
								\ \ \ \ \ PlotStyle -> {Hue[0.3, 0.8, 0.5]}, 
              DisplayFunction -> 
                Identity]; 
									p2\ \ \ \ \ \ \ \ \ \ \  = 
            Show[p1, 
              ListPlot[
                Transpose[{tarr, 
                    tavg}], 
							\ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {{0, tmax}, {0, x0*1.1}}, 
                PlotStyle -> {PointSize[0.01], 
                    Hue[1]}, 
								\ \ DisplayFunction -> Identity], \ 
              DisplayFunction -> 
                Identity]; 
									p2 = 
            Show[{p1, p2}, 
              Graphics[
								\ \ \ \ {Text["class.", {tmax*1.05, \((x0 - \((Mod[tmax - tp, 2\ tp] - tp)\)\^2)\)*1.05}, 
								TextStyle -> {FontColor -> Hue[0.3, 0.8, 
                            0.5]}], 
								\ \ Text["QM", {tmax*1.05, 0.95*tavg[\([Nt]\)]}, 
                    TextStyle -> {FontColor -> 
                          Hue[1]}]}], 
								\ \ \ \ \ \ \ PlotRange -> {{0, tmax*1.1}, {0, x0*1.1}}, 
              AxesLabel -> {"t", "<x>(t)"}, 
              DisplayFunction -> 
                Identity];
];\)\(
\)
    \)\ 
\(start = 
        ButtonBox[
          StyleBox["Run Animation", "Title", 
            FontFamily -> "Helvetica"], 
          ButtonFunction :> 
            FrontEndExecute[{SelectionMove[EvaluationNotebook[], Before, 
                    ButtonCell]
SelectionMove[
                    EvaluationNotebook[], Previous, 
                    CellGroup]
SelectionAnimate[
                    SelectedNotebook[], 
                    AnimationDisplayTime -> simtime]}], 
          ButtonEvaluator -> Automatic, 
          Active -> True];\)\)

Meaning of parameter:
n:		dimension of Hemline-matrix
		in DynExact only needed  for plot-resolution
xmax:	maximal x-coordinate included in discretized region
Nt:		# of time steps
Nper:	Number of Periods 
x0:		 center of initial Gaussian wave-package
σ:		width of initial  Gaussian wave-package
Mmax:		max # of states included in spectral theorem

Time-evolution of Gaussian wave-package

dyn[nd, xmaxd, Nt0, Nper0, 
        x0, σ0];
\(start /. simtime -> \ 0.1 // DisplayForm	\)

If you want to save this animation as an animated GIF file, you have to 
uncomment the following cell. Just delete the (* and the *), saving may take 
fairly long.

Time-evolution of <x>

\(\(\(DynMean[nd, xmaxd, Nt0, Nper0, x0, σ0];\)\(	\)\)\)

We see that the quantum mechanical particle never reaches the ground. The 
reason for this behavior is the uncertainty principle: Both x and momentum 
would be zero and thus certain.

The second difference is, that the quantum mechanical particle never again 
reaches its original position.

Save plot to ps-file

SetDirectory["~/mathem"];
Export["meanfinite.gif", p2, "GIF"];

Continuum dynamics (no finite differences!)

ExactDyn[nd, xmaxd, Nt0, Nper0, x0, σ0, Mmaxe]; 
    start /. simtime ->  0.1 // DisplayForm

If you want to save this animation as an animated GIF file, you have to 
uncomment the following cell.

The red curves show the development of the wave function, the  expectation value  and the variance of the particle's position., the blue  curves give the initial wave packet and expectation value. If the blue wave  function does not look Gaussian, Mmaxe is probably too small!
Comparison to classical motion:

p2 = 
        Show[{p2}, DisplayFunction -> $DisplayFunction];

Save plot to ps-file

SetDirectory["~/mathem"];
Export["meanfall.gif", p2, "GIF"];