Numerical Schrödinger Equation in 1D: a Crank-Nicolson Solver

1. Discretising Space

Divide the domain [x_min, x_max] into N grid points spaced by dx = (x_max − x_min)/(N−1). The complex wave function ψⱼ = ψ(xⱼ, t) is stored as two parallel Float64Array buffers, real and imaginary — never as a single complex type, since JavaScript has none built in:

const N  = 400;
const xMin = -20, xMax = 20;
const dx = (xMax - xMin) / (N - 1);
const x  = Float64Array.from({ length: N }, (_, j) => xMin + j * dx);
const psiRe = new Float64Array(N);
const psiIm = new Float64Array(N);

// Initial Gaussian wave packet: psi(x,0) = exp(-(x-x0)^2/(4*sigma^2)) * exp(i*k0*x)
function initGaussian(x0, sigma, k0) {
  let norm = 0;
  for (let j = 0; j < N; j++) {
    const env = Math.exp(-(x[j] - x0) ** 2 / (4 * sigma * sigma));
    psiRe[j] = env * Math.cos(k0 * x[j]);
    psiIm[j] = env * Math.sin(k0 * x[j]);
    norm += (psiRe[j] ** 2 + psiIm[j] ** 2) * dx;
  }
  const invSqrtNorm = 1 / Math.sqrt(norm);
  for (let j = 0; j < N; j++) { psiRe[j] *= invSqrtNorm; psiIm[j] *= invSqrtNorm; }
}

2. The Finite-Difference Hamiltonian

The second spatial derivative in Ĥ = −ℏ²/(2m)∂²/∂x² + V(x) is approximated by the standard three-point central difference:

Written as a matrix, Ĥ is tridiagonal: −α on the two off-diagonals, (2α + Vⱼ) on the diagonal. This sparsity is the entire reason the whole scheme stays fast — a direct dense solve would cost O(N³) per time step; the tridiagonal structure brings it down to O(N).

3. Why Explicit Euler Fails

The obvious first attempt — ψ(t+dt) = ψ(t) − i(dt/ℏ)Ĥψ(t) — is the explicit (forward) Euler method applied to iℏ∂ψ/∂t = Ĥψ. It is unconditionally unstable: because Ĥ is Hermitian with real eigenvalues E, the exact time evolution factor e^−iEdt/ℏ has magnitude exactly 1, but its first-order Taylor approximation 1 − iEdt/ℏ has magnitude √(1+(Edt/ℏ)²) > 1 for any non-zero dt. Every mode grows exponentially, however small the time step.

4. The Crank-Nicolson Scheme

Crank-Nicolson averages the Hamiltonian's action at the old and new time steps, giving a scheme that is second-order accurate in dt and exactly unitary (norm-preserving) for any time step:

The right-hand side is a simple tridiagonal matrix-vector product (explicit, cheap). The left-hand side requires solving a tridiagonal linear system each step — still O(N), but a genuine solve rather than a direct evaluation.

5. Solving the Tridiagonal System: the Thomas Algorithm

A tridiagonal system Ax = d (sub-diagonal a, diagonal b, super-diagonal c) is solved in a single forward-then-back sweep, the Thomas algorithm — Gaussian elimination specialised to bandwidth 1:

// Complex tridiagonal solve via Thomas algorithm.
// a, b, c, d are arrays of {re, im} complex numbers; b/d are overwritten.
function thomasSolve(a, b, c, d, N) {
  // Forward sweep: eliminate sub-diagonal
  for (let j = 1; j < N; j++) {
    const w = cdiv(a[j], b[j - 1]);         // w = a[j] / b[j-1]
    b[j] = csub(b[j], cmul(w, c[j - 1]));  // b[j] -= w * c[j-1]
    d[j] = csub(d[j], cmul(w, d[j - 1]));  // d[j] -= w * d[j-1]
  }
  // Back substitution
  const x = new Array(N);
  x[N - 1] = cdiv(d[N - 1], b[N - 1]);
  for (let j = N - 2; j >= 0; j--) {
    x[j] = cdiv(csub(d[j], cmul(c[j], x[j + 1])), b[j]);
  }
  return x;
}
// cmul/cdiv/csub: trivial complex arithmetic helpers on {re, im} objects
function cmul(p, q) { return { re: p.re*q.re - p.im*q.im, im: p.re*q.im + p.im*q.re }; }
function csub(p, q) { return { re: p.re - q.re, im: p.im - q.im }; }
function cdiv(p, q) { const d2 = q.re*q.re + q.im*q.im; return { re: (p.re*q.re + p.im*q.im)/d2, im: (p.im*q.re - p.re*q.im)/d2 }; }

6. Boundary Conditions and Potentials

Two boundary choices are common. Hard walls fix ψ(x_min) = ψ(x_max) = 0 for all t (equivalent to an infinite square well containing the whole grid) — simplest to implement, but wave packets reflect unphysically off the domain edges. Absorbing boundaries add a smooth imaginary potential −iΓ(x) near the edges (a "complex absorbing potential", CAP) that damps ψ before it reaches the boundary, mimicking an open, infinite domain:

// Potential: finite square well/barrier + absorbing edges
function buildPotential(x, N, wellDepth, wellWidth, gamma0 = 2.0, absorbWidth = 3) {
  const V = new Float64Array(N);
  const Vi = new Float64Array(N);  // imaginary part: absorbing layer
  for (let j = 0; j < N; j++) {
    V[j] = Math.abs(x[j]) < wellWidth / 2 ? -wellDepth : 0;
    const distFromEdge = Math.min(x[j] - x[0], x[N-1] - x[j]);
    if (distFromEdge < absorbWidth) {
      Vi[j] = -gamma0 * (1 - distFromEdge / absorbWidth) ** 2;
    }
  }
  return { V, Vi };
}

7. Full Working Code

Putting it all together: build the tridiagonal Ĥ, assemble the (I ± iĤdt/2ℏ) matrices, and step forward each frame.

function createSolver({ N, dx, hbar, mass, V, dt }) {
  const alpha = hbar * hbar / (2 * mass * dx * dx);
  const r = { re: 0, im: dt / (2 * hbar) };   // i*dt/(2*hbar) split into re/im

  // Diagonal of H (tridiagonal): 2*alpha + V[j]; off-diagonals: -alpha
  const diagH = Array.from({ length: N }, (_, j) => 2 * alpha + V[j]);

  function step(psiRe, psiIm) {
    // Right-hand side: (I - i*H*dt/2hbar) psi, tridiagonal matvec
    const dRe = new Float64Array(N), dIm = new Float64Array(N);
    for (let j = 0; j < N; j++) {
      let hRe = diagH[j] * psiRe[j], hIm = diagH[j] * psiIm[j];
      if (j > 0)     { hRe -= alpha * psiRe[j-1]; hIm -= alpha * psiIm[j-1]; }
      if (j < N - 1) { hRe -= alpha * psiRe[j+1]; hIm -= alpha * psiIm[j+1]; }
      // psi - i*(dt/2hbar)*H*psi = psi_re + H_im*(dt/2hbar), psi_im - H_re*(dt/2hbar)
      dRe[j] = psiRe[j] + r.im * hIm;
      dIm[j] = psiIm[j] - r.im * hRe;
    }
    // Left-hand side matrix (I + i*H*dt/2hbar): tridiagonal, solved with Thomas algorithm
    const a = Array.from({ length: N }, () => ({ re: 0, im: -r.im * -alpha }));
    const b = Array.from({ length: N }, (_, j) => ({ re: 1, im: r.im * diagH[j] }));
    const c = a.slice();
    const d = Array.from({ length: N }, (_, j) => ({ re: dRe[j], im: dIm[j] }));
    const x = thomasSolve(a, b, c, d, N);
    for (let j = 0; j < N; j++) { psiRe[j] = x[j].re; psiIm[j] = x[j].im; }
  }
  return { step };
}

8. Verifying Accuracy: Normalisation and Energy

Two cheap diagnostics catch nearly every discretisation bug before it becomes a visual artefact:

• Normalisation: Σ|ψⱼ|²dx should stay extremely close to 1 every frame (drift beyond ~10⁻⁶ per step without an absorbing boundary means a sign error, an asymmetric matrix, or dt too large relative to the finest energy scale in Ĥ).
• Energy: ⟨Ĥ⟩ = Σψⱼ*(Ĥψ)ⱼdx should be constant in time for a closed system (no absorbing layer, no explicit time-dependence in V) — Crank-Nicolson conserves this exactly in exact arithmetic, so drift here signals a coding bug rather than a scheme limitation.