The Crank-Nicolson technique below is obtained from:


Options[CrankNicolson] = {MaxIterations -> 5, 
   Tolerance -> Automatic};
CrankNicolson /: 
 NDSolve`InitializeMethod[CrankNicolson, stepmode_, sd_, rhs_, state_,
   OptionsPattern[CrankNicolson]] := 
 Module[{prec, rtol, maxit}, maxit = OptionValue[MaxIterations];
  prec = state@"WorkingPrecision";
  rtol = OptionValue[Tolerance];
  If[rtol === Automatic, rtol = 10^(-prec*3/4)];
  CrankNicolson[maxit, rtol]]
CrankNicolson[maxit_, rtol_]["Step"[f_, h_, t0_, x0_, f0_]] := 
  Module[{J, LU, t1 = t0 + h, x1, f1, residual, err, done = False, 
    tol = rtol, count = 0}, x1 = x0 + h f0;
   f1 = f[t1, x1];
   x1 = x0 + (h/2) (f0 + f1);
   J = f["JacobianMatrix"[t1, x1]];
   LU = IdentityMatrix[Length[x1], SparseArray] - (h/2) J;
   LU = LinearSolve[LU];
   While[(count <= maxit) && ! done, f1 = f[t1, x1];
    residual = x1 - x0 - (h/2)*(f0 + f1);
    err = Norm[residual, Infinity];
    If[err < tol, done = True
     (*else*), x1 = x1 - LU[residual];
   If[count > maxit, Message[CrankNicolson::cvmit, maxit];
    x1 = $Failed];
   {x1, f1}];
CrankNicolson[___]["StepInput"] = {"F"["T", "X"], "H", "T", "X", "XP"};
CrankNicolson[___]["StepOutput"] = {"X", "XP"};
CrankNicolson[___]["DifferenceOrder"] := 2;
CrankNicolson[___]["StepMode"] := "Fixed";

Using it to solve and graph the Klein-Gordon equation:

usol = First[
   u /. NDSolve[{D[u[t, x], t, t] - D[u[t, x], x, x] + u[t, x] == 0, 
      u[0, x] == Exp[-((x^2)/8)]*Exp[-I*1.5*x]/(4*Pi)^0.25, 
      Derivative[1, 0][u][0, x] == 0, u[t, -100] == u[t, 100]}, 
     u, {t, 0, 100}, {x, -100, 100}, 
     Method -> {"DoubleStep", Method -> CrankNicolson}]];
p1 = Plot[Re[usol[0, x]], {x, -100, 100}, PlotRange -> All, 
   PlotStyle -> Red];
p2 = Plot[Im[usol[0, x]], {x, -100, 100}, PlotRange -> All];
Show[p1, p2]
p3 = Plot[Re[usol[50, x]], {x, -100, 100}, PlotRange -> All, 
   PlotStyle -> Red];
p4 = Plot[Im[usol[50, x]], {x, -100, 100}, PlotRange -> All];
Show[p3, p4]
p5 = Plot[Re[usol[100, x]], {x, -100, 100}, PlotRange -> All, 
   PlotStyle -> Red];
p6 = Plot[Im[usol[100, x]], {x, -100, 100}, PlotRange -> All];
Show[p5, p6]

Do the opposite-moving wave packets correspond to different signs of momentum or different signs of energy? The Klein-Gordon equation has both positive energy/frequency and negative energy/frequency solutions, and I need to get each type of solution by itself. How can this be done?

  • 1
    $\begingroup$ Your initial condition sets the time derivative to zero. This means that the energy expectation value of the wave packet is zero - so it's neither positive nor negative. Since plane waves solve the KG equation and are also energy eigenfunctions, you could construct a wave packet solution by superimposing only those plane waves that have the desired sign of the frequency (energy). Then the energy expectation value of the wave packet would have the same sign as its Fourier components. I don't see this as a Mathematica question, though. It seems to be all about choosing the initial conditions. $\endgroup$
    – Jens
    Commented Jun 30, 2018 at 5:18

1 Answer 1


To answer the question with minimal modifications to the original code, I would just change the initial condition for the time derivative to one of these:

Derivative[1, 0][u][0, x] == 2 I Exp[-((x^2)/8)]*Exp[-I*1.5*x]


Derivative[1, 0][u][0, x] == -2 I Exp[-((x^2)/8)]*Exp[-I*1.5*x]

The initial condition for u[0, x] contains a factor that provides a momentum expectation value of roughly $-5.3$. More precisely, if I define the initial wave packet as

f[x_] := Exp[-((x^2)/8)]*Exp[-I*1.5*x]

Then the momentum expectation value is

-I NIntegrate[Conjugate[#] D[#, x] &[f[x]], {x, -∞, ∞}]

-5.31736 - 5.22818*10^-9 I

The imaginary part can be neglected. In my tests I also dropped the factor 1/(4*Pi)^0.25 from the initial conditions because I care only about the signs and not the magnitudes of the expectation values.

We can distinguish positive from negative energy by checking whether the wave packet predominantly propagates in the direction of the momentum (for positive energy) or opposite to it (for negative energy). These two cases can be approximately generated by the two initial conditions for the time derivative above. With the first choice I get a wave packet that moves predominantly toward positive x, i.e., opposite to the momentum direction. This corresponds to a negative energy:

negative energy

The second choice of initial condition gives a similar plot but with the larger wave-packet component moving toward negative x, corresponding to positive energy.

My goal here was to make only a small modification to the given code, not to produce exact results. The main point is that you have to give the wave packet a nonzero time derivative. It should also be chosen such that the energy expectation value is real-valued, which is why I added the imaginary prefactor on the first two lines above.

  • $\begingroup$ Thank you Jens, that is very helpful. $\endgroup$ Commented Jun 30, 2018 at 21:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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