7
$\begingroup$

The book "Advanced Visual Quantum Mechanics" by Thaller includes Mathematica software packages for the numerical solution of the Klein-Gordon equation and the Dirac equation (subject to user-defined boundary conditions), and the plotting of these solutions. However, this software is now out-of-date. Does anyone know of more current software packages for doing the same things?

8 July 2018 update: I have found a way to reproduce the results of Thaller's software packages using simpler and faster Mathematica code:

Options[CrankNicolson] = {MaxIterations -> 10, 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];
     count++;]];
   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";
(*psi at t=0 centered at x=0:*)

psi0[x_] = 0.3943 E^(1.5 I x - x^2/8);
(*dpsi/dt at t=0:*)

dpsidtfit0[x_] = 
  0.7199 E^(-0.122 x^2) (-I Cos[1.611 x] + Sin[1.611 x]);
xl = -20;
xr = +240;
tf = 240;
usol = First[
  u /. NDSolve[{D[u[t, x], t, t] - D[u[t, x], x, x] + u[t, x] == 0, 
     u[0, x] == psi0[x], Derivative[1, 0][u][0, x] == dpsidtfit0[x], 
     u[t, xl] == 0, u[t, xr] == 0}, u, {t, 0, tf}, {x, xl, xr}, 
    Method -> {"DoubleStep", Method -> CrankNicolson}]]
pr0 = Plot[Re[usol[0, x]], {x, -10, 10}, PlotRange -> All, 
   PlotStyle -> Red];
pi0 = Plot[Im[usol[0, x]], {x, -10, 10}, PlotRange -> All, 
   PlotStyle -> Blue];
Show[pr0, pi0]
rho[t, x] = -(I/
    2) (usol[t, x] Conjugate[Derivative[1, 0][usol][t, x]] - 
     Conjugate[usol[t, x]] Derivative[1, 0][usol][t, x]);
Plot3D[rho[t, x], {t, 0, 240}, {x, -20, 240}, PlotRange -> All, 
 PlotPoints -> 500, AxesLabel -> Automatic]

The Crank-Nicolson algorithm was copied from: http://reference.wolfram.com/language/tutorial/NDSolvePlugIns.html

$\endgroup$
6
  • 6
    $\begingroup$ Could you please include some sample code or hyperlink to the reference you mention? Currently, it reads as a general question on "how to solve..." That would probably help lots of users provide you with feedback. Does this help with the Klein-Gordon equation and perhaps this? As for other software that could be used for solution - have you tried Mathematica's NDSolve option first? Perhaps, the python community already has code on this? $\endgroup$ – dearN Jun 25 '17 at 21:58
  • 1
    $\begingroup$ I was unable to find any code in the book. Could you please tell us the page number? $\endgroup$ – zhk Jun 26 '17 at 0:06
  • $\begingroup$ The original Mathematica code is on a CD attached to the back cover of the book. Additional code can be downloaded from vqm.uni-graz.at/index.html. $\endgroup$ – Michael B. Heaney Jun 26 '17 at 1:11
  • 1
    $\begingroup$ Edit your question and include the relevant equation with initial/boundary conditions. It seems you are looking for the solution of two PDE's, so ask about it in two questions. $\endgroup$ – zhk Jun 26 '17 at 1:43
  • 1
    $\begingroup$ @MichaelB.Heaney I am voting to reopen your question. If it does get reopened, then it would be best for you to transfer your new code into an answer. $\endgroup$ – MarcoB Jul 9 '18 at 12:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.