QM Wavepackets

Crank-Nicholson Calculation

by Kevin J. McCann
September 15, 1999

Mathematica Preliminaries

Clear["Global`*"]

Off[General::"spell"];
Off[General::"spell1"];

Needs["Utilities`Notation`"];
Needs["Graphics`Colors`"];
Needs["LinearAlgebra`Tridiagonal`"];
Needs["Graphics`Graphics3D`"];

\(\(SetOptions[Plot, 
	Frame -> True, 
	
        FrameStyle -> AbsoluteThickness[2], 
	
        GridLines -> Automatic, 
	
        PlotStyle -> {{Red, AbsoluteThickness[3]}, {Black, 
              AbsoluteThickness[3]}, {Blue, 
              AbsoluteThickness[3]}}
];\)\)

\(\(SetOptions[ListPlot, 
	Frame -> True, 
	
        FrameStyle -> AbsoluteThickness[2], 
	
        GridLines -> Automatic, 
	
        PlotStyle -> {Red, 
            AbsoluteThickness[
              3]}, 
\ \ \ \ \ \ \ \ PlotJoined -> 
          True
];\)\)

$DefaultFont = {"Arial-Bold", 12};

{$LicenseID, $Version, $ReleaseNumber}

My system: {L2221-7779,4.0 for Microsoft Windows (April 21, 1999),0} PentiumII 450 MHz, running WNT.

Symbols

Symbolize[(Subscript[k, o])];Symbolize[(Subscript[x, ∞])];Symbolize[(Subscript[x, o])];Symbolize[(Subscript[N, x])];Symbolize[(Subscript[N, t])];

Crank-Nicholson

Parameters and initial condition

ℏ = 1.0;
m = 1.0;
σ = 5.0;
τ = 2  (σ^2);
(Subscript[x, o]) = -50.0;
(Subscript[k, o]) = 1.0;
L = 200.0;
T = 100.0;
(Subscript[N, x]) = 1001;
(Subscript[N, t]) = 401;

Δx  =  
        (L/((Subscript[N, x]) - 1));

Δt  =  
        (T/((Subscript[N, t]) - 1));

ρ  =  ((I ℏ Δt)/(2  m (((Δx))^2)));

Remove[Ψ];
Ψ[x_, 
        t_] := ((((1/(2  π (σ^2)))))^(1/
              4)) (1/(Sqrt[1  +  I (t/τ)])) (E^(I ((Subscript[k, o]) (x - 
                        (Subscript[x, o])) - ((ℏ ((Subscript[k, o])^2))/(2  m)) 
                    t))) (E^(-((((x - ((Subscript[x, o]) + ((ℏ (Subscript[k, o]))/m) t)))^2)/(4  (σ^2) (1 + I t/τ)))))

Set up the x-grid

x = 
        Table[-(L/2) + (j - 1) Δx, {j, 1, 
            (Subscript[N, x])}];

ListPlot[Transpose[{x, Abs[Ψ[x, 0]]}], 
        PlotRange -> All];

U is a matrix which will hold the wavepacket for (x,t).

Initialize U:

U = {Ψ[x, 0]};

Set up the Crank-Nicholson coefficients for the free-particle wavepacket:

offDiag  =  Table[-ρ, {j, 2, (Subscript[N, x]) - 2}];

diag  =  Table[2 + 2  ρ, {j, 2, (Subscript[N, x]) - 1}];

{(Subscript[N, x]), Δx, (Subscript[N, t]), Δt, ρ}

Solution

Now run the thing for (Subscript[N, t]) time steps

Timing[
      For[j = 2, 
        j ≤ 
          (Subscript[N, t]), 
χ  =  ρ (Take[
                    U[[j - 1]], {3, (Subscript[N, x])}] + 
                  Take[U[[j - 1]], {1, (Subscript[N, x]) - 2}]) + (2 - 
                  2  ρ) 
              Take[U[[j - 1]], {2, (Subscript[N, x]) - 1}]; 
u = 
          TridiagonalSolve[offDiag, diag, 
            offDiag, χ]; 
AppendTo[u, 
          0.0]; 
PrependTo[u, 
          0.0]; 
AppendTo[U, 
          u]; 
j++;
]
]

For my system, the above gave {438.26 Second,Null}.

Length[U]

Generate an InterpolatingFunction for \(TraditionalForm\`\(\(|\)\(ψ \((x, t)\)\)\( | \^2\)\)\)

f  =  
      ListInterpolation[
        Abs[U]^2 // Transpose, {{-L/2, L/2}, {0, T}}]

Plot3D[f[x, t], {x, -L/2, L/2}, {t, 0, T}, 
        ViewPoint -> {4, 10, 5}, PlotPoints -> 50, 
        PlotRange -> {0, 0.1}, Lighting -> False, 
        ColorFunction -> Hue];

Check with the analytic answer

Plot[{f[x, 0], f[x, 100], 
          Abs[Ψ[x, 100]]^2}, {x, -L/2, L/2}, 
        PlotRange -> All];

Reflected Wavepacket

Parameters and initial conditions

x =. ;

Remove[V];
V[x_] := If[x ≥ 0, 100.0, 0.0]

SetAttributes[V, Listable]

{V[-0.001], V[0.001]}

ℏ = 1.0;
m = 1.0;
σ = 5.0;
τ = 2  (σ^2);
(Subscript[x, o]) = -50.0;
(Subscript[k, o]) = 1.0;
L = 200.0;
T = 100.0;
(Subscript[N, x]) = 1001;
(Subscript[N, t]) = 401;

Δx  =  
        (L/((Subscript[N, x]) - 1));

Δt  =  
        (T/((Subscript[N, t]) - 1));

ρ  =  ((I ℏ Δt)/(2  m (((Δx))^2)));

Set up the x-grid

x = 
        Table[-(L/2) + (j - 1) Δx, {j, 1, 
            (Subscript[N, x])}];

Initialize U:

U = {Ψ[x, 0]};

Set up the Crank-Nicholson coefficients for the free-particle wavepacket:

offDiag  =  Table[-ρ, {j, 2, (Subscript[N, x]) - 2}];

diag  =  
        Table[2 + 2  ρ + (I/ℏ) V[x[[j]]], {j, 
            2, (Subscript[N, x]) - 1}];

Solution

Timing[
      For[j = 2, 
        j ≤ 
          (Subscript[N, t]), 
χ  =  ρ (Take[
                    U[[j - 1]], {3, (Subscript[N, x])}] + 
                  Take[U[[j - 1]], {1, (Subscript[N, x]) - 2}]) + 
            Table[(2 - 
                    2  ρ - (I/ℏ) 
                      V[x[[n]]]) U[[j - 1]][[n]], {n, 2, 
                (Subscript[N, x]) - 1}]; 
u = 
          TridiagonalSolve[offDiag, diag, 
            offDiag, χ]; 
AppendTo[u, 
          0.0]; 
PrependTo[u, 
          0.0]; 
AppendTo[U, 
          u]; 
j++;
]
]

For my system, the above gave {505.156 Second,Null}.

Length[U]

g = 
      ListInterpolation[
        Abs[U]^2 // Transpose, {{-L/2, L/2}, {0, T}}]

x =.

Plot3D[g[x, t], {x, -L/2, 0}, {t, 0, T}, 
        ViewPoint -> {-0.6, 3.7, 4.3}, PlotPoints -> 100, 
        PlotRange -> {0, 0.1}, Lighting -> False, 
        ColorFunction -> Hue, Mesh -> False];

ContourPlot[g[-x, t], {x, 0, L/2}, {t, 0, T}, 
        PlotPoints -> 100, ColorFunction -> Hue];

Plot[{g[x, 0], g[x, 100], 
          Abs[Ψ[-x, 100]]^2}, {x, -L/2, 0}, 
        PlotRange -> All];

Probability Conservation

After all, this is why we are using Crank-Nicholson

NIntegrate[g[x, 0], {x, -L/2, L/2}]

NIntegrate[g[x, 50], {x, -L/2, L/2}]

NIntegrate[g[x, 100], {x, -L/2, L/2}]