# Crank-Nicolson scheme for the Schrödinger equation

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?

• Back in the late 90s I recall Rob Knapp writing an example for CN for Schrodinger. Please search the Wolfram library and I think you should be able to locate it. Commented Dec 13, 2013 at 2:18
• @MikeHoneychurch I tried to search but can't find Ron Knapp's example, could you search if you have time? I'm very interested in this. But I find this example by Terry Robb, which solves the TDSE using a external Fortran subroutine SCHROED. And SCHROED uses a Crank-Nicholson scheme. Commented Mar 6, 2014 at 4:51

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].

1. 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[
DeveloperToPackedArray@#["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}]
`