Crank-Nicolson scheme for the Schrödinger equation

Question

I am trying to solve the 1-D Schrödinger equation with the Crank-Nicolson method. I use the atomic units:

• time < $10^{-14}$ => $413$ a.u.t.
• length $30 \times 10^{-9}$ m => $567$ Bohr radii

---

a = 50; 
σ = a/10;
p[x_, 0] := Exp[-(x^2)/(2 σ^2)]; (* граничные условия *)
p[±a, t_] := 0; (* н.у. *)
NN = 100;
h = 1/NN;
repl = x -> ((i*6*a)/NN - 3*a); 

eq = D[Subscript[p, i][t], t] == I/2(Subscript[p, i - 1][t] + 2 Subscript[p, i][t] - 
  Subscript[p, i + 1][t])/h^2(*+D[u[x],x]*Subscript[p,i][t]*)/. repl; 

Table[eq, {i, 0, NN}];
boundary = {Subscript[p, 0][t] == 0, Subscript[p, NN][t] == 0}; 
p0[x_] := Exp[-(x^2)/(2 σ^2)]; 
Cauchy = Table[Subscript[p, i][0] == p0[x] /. repl, {i, 1, NN - 1}]; (*задача Коши*)

eqns = Join[Table[eq, {i, 1, NN - 1}], boundary, Cauchy];
sol = NDSolve[N[eqns], 
N[Table[Subscript[p, i][t], {i, 0, NN}]], {t, 0, 413}];

and get this error:

NDSolve::mconly: For the method IDA, only machine real code is available. 
Unable to continue with complex values or beyond floating-point exceptions. >>

I want to get something like this:

sol = NDSolve[{I D[u[t, x], t] == (-1/2)D[u[t, x], {x, 2}], 
     u[0., x] == Exp[-(x^2./(2*σ^2))], u[t, a] == 0, 
     u[t, -a] == 0}, u, {t, 0, 4130}, {x, -a, a}, (*MaxStepSize->0.01,*)
     AccuracyGoal -> 3, PrecisionGoal -> 3];

Animate[Plot[Evaluate[Abs[u[t, x] /. First[sol]]^2], {x, -a, a}, 
     PlotRange -> {0, 1}], {t, 0, 413}]

What's wrong with my code and how to fix it?

Answer 1

Several issues here.

1. Your code isn't an implementation of Crank–Nicolson method, but a implementation of method of lines. The way for setting Crank–Nicolson method inside NDSolve has been included in this tutorial, in the rest part of this answer I'll simply fix your code.

2. There're several simple mistakes in your code:

(1) The step size is wrong, h = 1/NN should be h = (2 a)/NN.

(2) The transformation rule is wrong. repl = x -> ((i*6*a)/NN - 3*a) should be repl = x -> (i 2 a)/NN - a. BTW, though not quite necessary here, one can make use of RescalingTransform to obtain the transform relation like this:

    RescalingTransform[{{0, NN}}, {{-a, a}}][{i}]
    (* {-a + (2 a i)/NN} *)

(3) The difference formula is wrong. Subscript[p, i - 1][t] + 2 Subscript[p, i][t] - Subscript[p, i + 1][t] should be Subscript[p, i - 1][t] - 2 Subscript[p, i][t] + Subscript[p, i + 1][t].

3. The DAE solver of NDSolve seems not to be strong enough. The boundary conditions needs to be modified to ODEs, or NDSolve will spit out mconly warning and fails. This is curtly discussed in this tutorial, particularly the part about Boundary Conditions.

The following is the fixed code:

a = 50;
σ = a/10;
p[x_, 0] := Exp[-(x^2/(2 σ^2))];
NN = 100;
h = (2 a)/NN;
repl = x -> (i 2 a)/NN - a;
eq = D[Subscript[p, i][t], t] == 
       (I (Subscript[p, i - 1][t] - 2 Subscript[p, i][t] + 
               Subscript[p, i + 1][t]))/(2 h^2);
boundary = {Subscript[p, 0][t] == 0, Subscript[p, NN][t] == 0};
(* Modify the original b.c.s to almost equivalent ODEs,
   this is also the processing strategy inside NDSolve. *)
diffboundary = With[{sf = 1}, Map[D[#, t] + sf # &, boundary, {2}]];
p0[x_] := Exp[-(x^2/(2 σ^2))];
Cauchy = Table[
  Subscript[p, i][0] == p0[x] /. repl // Evaluate, {i, 0, NN}(*Notice the modification*)]; 
eqns = {Table[eq, {i, 1, NN - 1}], diffboundary, Cauchy};
sollst = NDSolveValue[eqns, Table[Subscript[p, i], {i, 0, NN}], {t, 0, 413}, 
    MaxSteps -> Infinity]; // AbsoluteTiming
rebuild = ListInterpolation[
    Developer`ToPackedArray@#["ValuesOnGrid"] & /@ # // 
     Transpose, {#[[1]]["Coordinates"][[1]], Array[# &, NN + 1, {-a, a}]}] &@sollst
Manipulate[Plot[Abs[rebuild[t, x]]^2, {x, -a, a}, PlotRange -> {0, 1}], {t, 0, 413}]