This series starts from a numerical question: when a continuous Schrodinger equation becomes a finite matrix problem, what has actually been approximated? DVR and finite difference can both converge to the same continuum physics, but their finite-dimensional pictures are very different.

The common target

All notes in this series are organized around the time-dependent Schrodinger equation

\[ i\hbar\frac{d}{dt}\psi(t)=\hat H\psi(t),\qquad \hat H=-\frac{\hbar^2}{2m}\partial_x^2+V(x), \]

with initial condition \(\psi(0)=\psi_0\). If \(\hat H\) is time independent, the exact formal solution is

\[ \psi(t)=\exp\!\left(-\frac{i}{\hbar}\hat Ht\right)\psi_0. \]

The numerical task has two layers: approximate the continuous coordinate Hilbert space, then apply the exponential time-evolution operator accurately enough for the observable of interest.

Continuous operator viewpoint

In the continuum, the Hamiltonian is an unbounded self-adjoint operator on a domain inside \(L^2(\mathbb{R})\). That statement sounds abstract, but it carries a practical warning: the true object is infinite dimensional and may contain both bound and scattering parts of the spectrum.

A numerical method replaces this operator by a finite one. The relevant question is not whether the matrix is "the Hamiltonian"; it is which limiting process makes the matrix approach the continuum Hamiltonian.

DVR as projection

A discrete variable representation can be understood as a projection of the continuous Hamiltonian onto a finite basis:

\[ P_N=\sum_{k=1}^N|\phi_k\rangle\langle\phi_k|,\qquad H_N=P_N\hat H P_N. \]

When \(P_N\) approaches the identity strongly and the chosen finite subspace becomes dense enough in the Hamiltonian domain, the projected dynamics approximates the continuum dynamics. In computation this translates to a more concrete rule: the grid must cover the physically visited region and resolve the shortest wavelength that matters.

Finite difference as an auxiliary lattice

Finite difference starts from a local derivative stencil. For the second derivative, the standard central stencil is

\[ (D_{\Delta x}u)_j=\frac{u_{j+1}-2u_j+u_{j-1}}{\Delta x^2}. \]

This produces a tridiagonal kinetic matrix. In a sinc DVR, by contrast, the kinetic matrix is dense but structured:

\[ T_{ij}=\frac{\hbar^2}{2m\Delta x^2}(-1)^{i-j} \begin{cases} \pi^2/3,& i=j,\\ 2/(i-j)^2,& i\ne j. \end{cases} \]

The finite-difference stencil is local and easy to scale. The DVR kinetic operator is nonlocal, but it retains a spectral character that is valuable for small quantum benchmark problems.

A small benchmark

As a first check, compare both discretizations on the harmonic oscillator with \(\hbar=m=\omega=1\). The exact energies are \(E_n=n+1/2\), so the numerical error can be read directly from the eigenvalues of the finite Hamiltonian.

Harmonic oscillator energy benchmark comparing sinc DVR and finite difference
Harmonic-oscillator benchmark on the same coordinate interval. In the left panel, blue circles are the sinc DVR absolute energy errors and orange squares are the second-order finite-difference errors for the first fourteen states on a fixed grid. In the right panel, the blue and orange curves show the ground-state error as the grid spacing \(\Delta x\) is reduced. The DVR errors sit near floating-point precision for this smooth bound-state problem, while the finite-difference stencil converges more slowly because its derivative approximation is only local and second order.

This benchmark is not a claim that DVR always wins. It says something narrower and more useful: for smooth low-dimensional reference problems, the dense sinc kinetic matrix can reach spectral accuracy with modest grid sizes. That is why it is a good basis for the later operator and correlation-function notes.

Exact finite-dimensional evolution

After a DVR Hamiltonian matrix has been assembled, real-time evolution can be performed by diagonalization, matrix exponentiation, or a stable Krylov method. A simple explicit finite-difference update, for comparison, has the recurrence form

\[ i\hbar\frac{\psi_j^{n+1}-\psi_j^n}{\Delta t} = -\frac{\hbar^2}{2m} \frac{\psi_{j+1}^n-2\psi_j^n+\psi_{j-1}^n}{(\Delta x)^2} +V_j\psi_j^n. \]

Its accuracy is tied simultaneously to \(\Delta x\) and \(\Delta t\). DVR is not automatically cheaper, but it gives a clean finite Hamiltonian whose spectral structure can be inspected directly.

Practical takeaway

I use DVR in the later notes because the target systems are small enough for dense matrices, and the desired outputs are reference-quality quantum wavepacket dynamics and Kubo correlation functions. The cost is larger memory and matrix algebra; the gain is a transparent operator representation.

The next note zooms in on the matrix objects themselves: kinetic blocks, local potentials, side projectors, and flux commutators.

Code used in this note

  • dvr_fd_benchmark.py builds the harmonic-oscillator benchmark and regenerates the figure above.
  • dvr_kubo_minimal.py collects the shared sinc DVR, finite-difference, wavepacket, Kubo, side, and flux routines used across the series.

DVR series map

  1. Part I: Wavepacket Dynamics Methods - discretization choices and the harmonic-oscillator benchmark.
  2. Part II: Operator Matrices in the DVR Basis - matrix structure behind kinetic, side, and flux operators.
  3. Part III: Wavepacket Dynamics with DVR - exact finite-dimensional propagation of a prepared wavepacket.
  4. Part IV: Density-Matrix Ensemble Evolution - thermal ensembles, filtering, and absorber bookkeeping.
  5. Part V: Kubo-Transformed Position Correlation - spectral thermal correlation functions.
  6. Part VI: Flux-Side Kubo Correlation - rate-like observables from side and flux operators.
  7. Part VII: Excited-Ground Kubo Population Correlation - local adiabatic projectors and population-resolved transfer diagnostics.