“Know how to solve every problem that has been solved.” “What I cannot create, I do not understand.” — Richard Feynman

Etude: 1D eigenvalue problems

Etudes

Don't worry, you don't need to play the piano or anything. An etude is a short piece of music designed to teach a technique. Eventually you play the pattern enough times to figure out how to play it. This page is all about using 1D finite differences to solve eigenvalue problems. Physics is filled with eigenvalue problems, and it's a shame that in linear algebra classes you don't get to solve more physics-y problems, since they often use quite a bit of linear algebra. Probably the physics background would take too much time in a math class. Still a loss though.

Everyone is using the word kernel these days. I don't know if it's a Linux tradition or what, but any time there's a serious routine you have to call it a kernel or people won't take you seriously. The other day I was adding 1 + 1 = 2 and I shipped it as a performant addition kernel. So this small program will also be a kernel. The kernel is about 12 lines of Python. You take ye-ole half-Laplacian ( $\frac{1}{2} \nabla^{2}$ ), discretize it on a grid, add a potential $V (x)$ , and hand the whole thing to np.linalg.eigh (a real kernel in my opinion (a true kernel if I do say so myself (a kernel of kernels))). The only thing that changes from problem to problem is the potential. The esteemed minds at notesoncomputing decided to choose the following systems:

Infinite well ( $V = 0$ ) — the calibration shot. We know the answer in closed form, so this is the test that the kernel works.
Harmonic oscillator ( $V = x^{2} /2$ ) — the textbook one. Gaussian wavefunctions, evenly spaced energies.
Anharmonic oscillator ( $V = x^{2} /2 + λ x^{4}$ ) — turn on the knob and watch perturbation theory miss.
Symmetric double well ( $V = - 2 x^{2} + 0.2 x^{4}$ ) — tunneling between two minima, with a splitting too small for any analytic series to handle cleanly.
Hydrogen radial ( $V = - 1/ r$ ) — the hydrogen atom, dropped onto a finite-difference grid.
And one bonus round — rerun the harmonic oscillator on a million-point grid, where the dense matrix won't fit in memory anymore. This is where Lanczos earns its keep.

The kernel

The Schrödinger equation on a finite interval, with the wavefunction vanishing at the ends, is

[- \frac{1}{2} \frac{d ^{2}}{d x ^{2}} + V (x)] ψ (x) = E ψ (x), ψ (x_{m i n}) = ψ (x_{m a x}) = 0.

Discretize the interval into $N$ evenly spaced points with spacing $Δ x$ . The second derivative becomes the three-point stencil

ψ^{''} (x_{i}) \approx \frac{ψ _{i - 1} - 2 ψ _{i} + ψ _{i + 1}}{( Δ x ) ^{2}} .

Plug that in and the differential equation collapses into a matrix eigenvalue equation $H ψ = E ψ$ , where $H$ is tridiagonal: the main diagonal carries the kinetic constant $+ 1/ (Δ x)^{2}$ plus the potential $V (x_{i})$ , and the two off-diagonals carry $- 1/ (2 (Δ x)^{2})$ . That's the whole derivation. Everything below is the same matrix with a different middle row.

import numpy as np

def solve_eigenvalue_1d(V_func, x_min, x_max, N, n_states=4):
    """Solve -1/2 d^2 psi/dx^2 + V(x) psi = E psi on [x_min, x_max].

    Dirichlet BCs (psi=0 at both ends — i.e., infinite well outside the
    interval, which is the standard "particle on a finite domain" setup).
    Uses a 3-point finite-difference stencil. Units where m = hbar = 1.
    """
    x = np.linspace(x_min, x_max, N)
    dx = x[1] - x[0]
    T_diag = np.full(N, 1.0 / dx**2)
    T_off  = np.full(N - 1, -0.5 / dx**2)
    V_diag = V_func(x)
    H = np.diag(T_diag + V_diag) + np.diag(T_off, 1) + np.diag(T_off, -1)
    eigvals, eigvecs = np.linalg.eigh(H)
    return x, eigvals[:n_states], eigvecs[:, :n_states]

Two things to notice now, before we use the kernel six times. First: the kernel knows nothing about quantum mechanics. It's a differential operator with Dirichlet boundary conditions. You could rename $E$ as the squared frequency of a standing wave on a string, or as the buckling load of a column, or as a mode of a heat kernel, and the same code would still solve the new problem. That's the breadth of the eigenvalue-problem template — it's about a self-adjoint operator on a 1D interval, not about quantum.

Second: every variation below changes one thing. The reader who knows the kernel cold can write any of these six in under a minute. That feeling — "I already know this, I just have to point the right function at the right potential" — is what the étude is for.

Variation 1 — infinite well (the calibration shot)

Set $V (x) = 0$ on $[0, π]$ . With the box walls at the boundary, this is the particle-in-a-box: a problem so embarrassingly tractable that anyone who has seen a quantum mechanics book in their life remembers the answer.

E_{n} = \frac{n ^{2}}{2}, n = 1, 2, 3, \dots

So we know exactly what the kernel should return: $0.5, 2, 4.5, 8$ . Anything else means a bug in the kernel.

# Variation 1: 1D infinite well on [0, pi], V = 0
x, E, _ = solve_eigenvalue_1d(lambda x: np.zeros_like(x), 0.0, np.pi, 800)
analytic = np.array([0.5 * n**2 for n in range(1, 5)])
print("finite difference:", E)
print("analytic n^2/2:   ", analytic)
print("relative error:   ", (E - analytic) / analytic)

finite difference: [0.4975056  1.99001475 4.47750449 7.95993655]
analytic n^2/2:    [0.5 2.  4.5 8. ]
relative error:    [-0.0049888  -0.00499263 -0.004999   -0.00500793]

The relative error is $- 0.005$ on every eigenvalue — the finite-difference scheme is biased low by a known amount. That bias scales like $- \frac{1}{12} (Δ x)^{2} \cdot E_{n}$ , the leading truncation error of the three-point stencil. You can read off the discretization error of your code from a problem with a closed form. That's what variation 1 is for. Never trust a numerical method on a hard problem until you've watched it solve an easy one.

Variation 2 — quantum harmonic oscillator

Change the potential to $V (x) = x^{2} /2$ . The walls are still there at $\pm 8$ , but they don't matter — the wavefunctions are Gaussian-localized near the origin and decay long before they reach the box edge. The closed-form answer is the textbook $E_{n} = n + \frac{1}{2}$ :

# Variation 2: 1D harmonic oscillator, V = x^2/2
x, E, _ = solve_eigenvalue_1d(lambda x: 0.5 * x**2, -8.0, 8.0, 800)
analytic = np.array([n + 0.5 for n in range(4)])
print("finite difference:", E)
print("analytic n + 1/2: ", analytic)
print("relative error:   ", (E - analytic) / analytic)

finite difference: [0.49998747 1.49993734 2.49983708 3.49968669]
analytic n + 1/2:  [0.5 1.5 2.5 3.5]
relative error:    [-2.50632459e-05 -4.17729136e-05 -6.51672028e-05 -8.95175154e-05]

Same kernel, different potential, errors down to $1 0^{- 5}$ . Why is this so much more accurate than variation 1? Because the wavefunctions of the oscillator are smooth Gaussians — they have far smaller second-derivative content at the boundaries than the abrupt sine waves of the box. The discretization error is largest where the wavefunction varies fastest. The same FD scheme gets two orders of magnitude more accurate just because the problem it's being pointed at is gentler.

The lesson: the kernel's accuracy isn't a property of the kernel. It's a property of the kernel-times-the-problem. The kernel is fine; the box problem is jagged.

Variation 3 — anharmonic oscillator (perturbation theory check)

Take the same oscillator and bend it: $V (x) = x^{2} /2 + λ x^{4}$ . For small $λ$ first-order perturbation theory says the ground-state energy moves up by the matrix element of $λ x^{4}$ in the unperturbed ground state:

E_{0} (λ) \approx \frac{1}{2} + λ ⟨ 0∣ x^{4} ∣0 ⟩ = \frac{1}{2} + \frac{3}{4} λ + O (λ^{2}) .

The kernel doesn't care that it's a perturbation problem — it just diagonalizes the new matrix. Read out the true $E_{0} (λ)$ for a series of $λ$ values and compare to the PT prediction:

# Variation 3: anharmonic oscillator, V = x^2/2 + lam * x^4
for lam in [0.0, 0.05, 0.1, 0.2]:
    x, E, _ = solve_eigenvalue_1d(
        lambda x, l=lam: 0.5 * x**2 + l * x**4, -8.0, 8.0, 800
    )
    # First-order perturbation theory: <0|lam x^4|0> = 3 lam / 4
    pt_E0 = 0.5 + 0.75 * lam
    print(f"lambda = {lam:.2f}:  E_0 (FD) = {E[0]:.6f},  E_0 (PT, 1st) = {pt_E0:.6f}")

lambda = 0.00:  E_0 (FD) = 0.499987,  E_0 (PT, 1st) = 0.500000
lambda = 0.05:  E_0 (FD) = 0.532627,  E_0 (PT, 1st) = 0.537500
lambda = 0.10:  E_0 (FD) = 0.559129,  E_0 (PT, 1st) = 0.575000
lambda = 0.20:  E_0 (FD) = 0.602383,  E_0 (PT, 1st) = 0.650000

At $λ = 0.05$ the PT estimate is off by 1%. By $λ = 0.2$ it's 8% high. The asymptotic series overestimates the energy because the higher-order terms have alternating signs and the series is divergent — Bender and Wu showed that the radius of convergence is exactly zero. (For more on resumming divergent series, the perturbation-methods section is where to look.) The numerical answer doesn't care about any of this; it's just an eigenvalue of a matrix.

What variation 3 buys you: a way to see perturbation theory breaking. You can dial $λ$ from "small" to "not small" and watch the analytic estimate drift away from the truth. That discrepancy is a meter for how big "small" actually is.

Variation 4 — symmetric double well (tunneling)

Flip the curvature of the harmonic oscillator and stabilize it with a quartic: $V (x) = - 2 x^{2} + 0.2 x^{4}$ . Two minima at $x = \pm 5$ , with a barrier between them. The wavefunctions of the lowest states are almost-degenerate symmetric/antisymmetric pairs: one even, one odd, with a tiny gap set by tunneling through the barrier.

# Variation 4: symmetric double well, V = -2 x^2 + 0.2 x^4
x, E, _ = solve_eigenvalue_1d(
    lambda x: -2.0 * x**2 + 0.2 * x**4, -8.0, 8.0, 1200, n_states=6
)
print("lowest 6 eigenvalues:", E)
print(f"tunneling splitting E_1 - E_0 = {E[1] - E[0]:.6e}")
print(f"next pair splitting E_3 - E_2 = {E[3] - E[2]:.6e}")

lowest 6 eigenvalues: [-3.64161741 -3.63990568 -1.25321343 -1.15131969  0.36885582  1.17695177]
tunneling splitting E_1 - E_0 = 1.711729e-03
next pair splitting E_3 - E_2 = 1.018937e-01

The lowest two eigenvalues differ by $1.7 \times 1 0^{- 3}$ while the next pair already differ by $0.1$ . That ratio of sixty is tunneling. The instanton-style estimate gives the splitting as $Δ E \sim e^{- S /ℏ}$ where $S$ is the Euclidean action along the bounce trajectory through the barrier. Going from a single-well to a double-well changes nothing about the numerical method — it's the same kernel — but produces qualitatively new physics: degeneracy lifting, parity-paired states, exponentially small splittings.

Worth feeling: in variations 1 and 2 the question was "how close is the answer to the analytic one." In variation 4 there is no clean analytic answer for the splitting — instanton theory is itself an asymptotic, semiclassical estimate. The numerical eigenvalue is the ground truth. The method has gone from "verifiable against a closed form" to "the closed form is verified against it." Same code, flipped relationship.

Variation 5 — hydrogen radial

The 3D hydrogen atom, after separating out the angular part, leaves a 1D radial equation for $u (r) = r R (r)$ :

[- \frac{1}{2} \frac{d ^{2}}{d r ^{2}} + V_{eff} (r)] u (r) = E u (r), V_{eff} (r) = - \frac{1}{r} + \frac{ℓ ( ℓ + 1 )}{2 r ^{2}} .

Strip the centrifugal term (set $ℓ = 0$ ). What's left is exactly the kernel's template, just with $V (r) = - 1/ r$ and the domain shifted to $r > 0$ . The boundary $u (0) = 0$ is automatic — set the grid to start a hair off the origin (e.g. $r_{m i n} = 0.01$ ) to keep the potential finite. The closed-form spectrum is the Rydberg formula $E_{n} = - \frac{1}{2 n ^{2}}$ :

# Variation 5: hydrogen radial (l = 0). V(r) = -1/r.
# In u-coordinates, psi(r) = u(r) / r. The boundary u(0) = 0 handles the
# singularity, provided we start the grid slightly off the origin.
def V_hydrogen(r):
    return -1.0 / r

r, E, U = solve_eigenvalue_1d(V_hydrogen, 0.01, 50.0, 4000, n_states=4)
analytic = np.array([-0.5 / n**2 for n in range(1, 5)])
print("finite difference 1s..4s:", E)
print("analytic -1/(2n^2):      ", analytic)
print("relative error:          ", (E - analytic) / analytic)

finite difference 1s..4s: [-0.50508463 -0.12563423 -0.05574334 -0.03128502]
analytic -1/(2n^2):       [-0.5        -0.125      -0.05555556 -0.03125  ]
relative error:           [0.01016926 0.00507388 0.00338018 0.00112048]

Same code as the harmonic oscillator. Different physical problem. The reader who's spent any time on atomic physics will have seen the Rydberg series derived from the analytic Coulomb-bound-state algebra — factor the radial equation, find associated Laguerre polynomials, quantize the energy through a regularity condition at infinity. None of that is what happens here. Here, the bound states are matrix eigenvalues, full stop. The discrete Rydberg spectrum is what falls out of a finite-difference grid. The structural form of the answer — a bound-state ladder converging to $E = 0$ — is the same answer the algebra produces, but the route is unrecognizable.

Two artifacts worth flagging. First: the FD error here is much bigger than for the oscillator (around 1% on the 1s, dropping to 0.1% on higher states), because the Coulomb wavefunctions have a sharp cusp at the origin that the smooth FD stencil can't resolve. If you want chemistry-grade accuracy from this method, you need a logarithmic grid or a different basis. Second: pushing $r_{m i n}$ very close to zero introduces a spurious "fall to the origin" state — the discrete potential at $r = r_{m i n}$ goes to $- 1/ r_{m i n}$ , artificially deep, and the FD scheme finds a fake eigenstate trapped against that wall. The fix is the offset we used. Both of these are discretization stories. The continuum problem doesn't have either issue.

Variation 6 — scaling: dense vs sparse Lanczos

The kernel above stores $H$ as a full $N \times N$ matrix. That's fine for $N = 800$ (5 MB) and even $N = 10, 000$ (800 MB if you're patient). For a 3D problem with a 100³ grid — a routine ask in chemistry or materials simulation — you'd be storing $(1 0^{6})^{2} = 1 0^{12}$ floating-point numbers, which is eight terabytes of RAM you don't have. The kernel doesn't scale.

But the matrix is tridiagonal. It has $3 N$ nonzero entries, not $N^{2}$ . Store only the nonzeros, give an iterative method a way to compute $H v$ for any vector $v$ , and you never materialize $H$ at all. The Lanczos algorithm uses that operator and produces the extreme eigenvalues in $O (N)$ memory and $O (k N)$ work per eigenpair, where $k$ is the number of Lanczos iterations (typically dozens, not thousands).

     N |  dense memory (full H) | sparse memory (tridiag) | dense vs sparse flop ratio
------ | ---------------------- | ----------------------- | --------------------------
   800 |              0.005 GB |              0.019 MB  |                       2e+04
10,000 |              0.800 GB |              0.240 MB  |                       3e+06
100,000|             80.000 GB |              2.400 MB  |                       3e+08
1e6    |          8,000.000 GB |             24.000 MB  |                       3e+10

Same eigenvalue problem. Same answer. But variation 6 is the moment when "switch to sparse storage and an iterative solver" goes from a style preference to a physical necessity. Dense at $N = 1 0^{6}$ is impossible on Earth's largest computer. Sparse-Lanczos at the same $N$ fits on a phone.

This is the move that distinguishes the model in undergraduate textbooks from how production quantum-mechanical codes are written. Casida's equation in TDDFT, eigenstates of an actual ab-initio Hamiltonian, density-of-states calculations in band theory, the eigenstates of a discretized Schrödinger operator on a real 3D molecular grid — none of these can afford dense diagonalization. They all live or die by iterative subspace methods (Lanczos, Davidson, LOBPCG) acting on a never-materialized Hamiltonian. Same skeleton, different basis for the matrix.

What just happened

Six problems. Same twelve lines. The variations swapped, in order:

nothing (calibration against a closed form),
the potential (smooth wavefunctions instead of jagged ones),
the potential with a parameter (perturbation theory's domain of validity, made visible),
the curvature sign (qualitatively new physics: tunneling),
the meaning of the coordinate (radial instead of Cartesian — a 3D problem reduced to 1D),
the matrix representation (dense to sparse — the move that makes the method tractable).

The kernel never changed. The physics changed every time. That asymmetry — small input change, large output change — is the affordance the eigenvalue-problem template gives you. Once you've internalized it, you stop seeing six problems. You see one problem with a parameter slot, and the next time a Schrödinger-shaped differential equation arrives wearing a costume, you reach for the same dozen lines.

A few patterns to take with you:

Discretize, then diagonalize. Any self-adjoint linear operator on an interval, with sensible boundary conditions, becomes a real symmetric matrix. Eigenvalues of the operator are (approximately) eigenvalues of the matrix. The recipe doesn't care whether the operator came from quantum mechanics, elasticity, acoustics, or graph Laplacians.
Calibrate against a closed form first. Variation 1 isn't there because anyone needs the box energies — they're three keystrokes by hand. It's there to verify that the code finds them. The first thing every numerical method should produce is a known answer.
Smooth problems are easier than jagged ones. Wavefunctions with sharp features (cusps, kinks, hard walls) require finer grids or different bases. The accuracy you get isn't a property of the method — it's a property of how well the method's assumptions match the problem.
The same code can produce qualitatively different physics. Variation 4's tunneling splitting and variation 5's Rydberg ladder come from the same routine; the potential decides whether the spectrum is a near-degenerate pair or a hydrogenic series. If you find yourself writing two separate programs for two problems that look this similar, you're probably missing the abstraction.
The matrix's structure decides whether the method is possible. Tridiagonality isn't a stylistic choice — it's what makes Lanczos work and what lets the method scale to a million points. Recognizing the sparsity pattern of $H$ is more useful than knowing six special functions.

Graduation problem

Try this without looking back. 2D harmonic oscillator: $V (x, y) = (x^{2} + y^{2}) /2$ , with the wavefunction confined to a square $[- 8, 8]^{2}$ . The analytic spectrum is $E_{n_{x}, n_{y}} = n_{x} + n_{y} + 1$ , so the lowest eigenvalues you should see are 1, 2, 2, 3, 3, 3, 4, 4, 4, 4 — with the degeneracy pattern $g_{n} = n + 1$ .

The mechanical question is "extend the 1D stencil to 2D." The pedagogically interesting question is "what changes in the kernel, and what stays the same?" Hint: the 2D Laplacian on an $N \times N$ grid produces an $N^{2} \times N^{2}$ matrix with five nonzero diagonals — the matrix Kronecker-sums of two 1D Laplacians. That extension is the entry point to every higher-dimensional finite-difference code on the site.

If you build the 2D solver and the degeneracies come out wrong, almost certainly the bug is in how you indexed the 2D grid into a 1D eigenvector. The math is fine; the bookkeeping is brittle. That's also the lesson.

Related on this site

The Lanczos algorithm gets its own page at Lanczos iteration; the closely-related Arnoldi iteration handles the non-symmetric case. For the linear-response generalization — same eigenvalue-problem skeleton, but non-Hermitian and twice as big — see Casida's equation. The single-particle template that turns interacting many-body problems into eigenvalue problems is on single-particle models. The full path from "discretize an operator" to "Schrödinger ground state via brute-force FD" is in the from-scratch H2 series. For why perturbation series around variation 3 diverge, the Borel-Padé and related pages are where to look.