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

Step 2 — One bound state in a fixed well

Nuclear Physics

Before any self-consistency: one nucleon, one fixed Woods–Saxon potential, the s-wave radial Schrödinger equation for $u = r R$ :

- \frac{ℏ ^{2}}{2 m} u^{''} + U (r) u = ε u, u (0) = u (r_{m a x}) = 0.

This step carries the two pieces of machinery the whole tutorial runs on: the Hamiltonian is a tridiagonal matrix, and its ground state comes from shifted inverse iteration with a shift supplied by Gershgorin's theorem. Neither is exotic — but both deserve to be actually explained, not name-dropped.

Why the matrix is tridiagonal

Put $u$ on the grid: it becomes a vector of samples $u_{i} = u (r_{i})$ . The second derivative at point $i$ is approximated from the point and its two neighbors:

u^{''} (r_{i}) \approx \frac{u _{i - 1} - 2 u _{i} + u _{i + 1}}{h ^{2}} .

So when you write the Schrödinger equation at grid point $i$ , only three unknowns appear: $u_{i - 1}, u_{i}, u_{i + 1}$ . Row $i$ of the matrix $H$ therefore has exactly three nonzero entries — the diagonal $d_{i} = 2 B / h^{2} + U (r_{i})$ and one coupling $o = - B / h^{2}$ to each neighbor. For five grid points:

H = d_{1} o o d_{2} o o d_{3} o o d_{4} o o d_{5}

Every blank is a true zero. The reason is physical, not accidental: a derivative is a local operator — the curvature at $r_{i}$ doesn't care what the wavefunction does five femtometers away — and locality of the operator becomes sparsity of the matrix. That sparsity is the whole performance story. Solving a linear system with a full $N \times N$ matrix costs $O (N^{3})$ ; with a tridiagonal one it costs $O (N)$ — one forward sweep, one back-substitution (the Thomas algorithm). On our 140-point grid that's the difference between ~3 million operations and ~400.

The plan: inverse iteration

We want the lowest eigenvalue of $H$ and its eigenvector. Expand any starting vector $u$ in the (unknown) eigenvectors $v_{k}$ of $H$ : $u = \sum_{k} c_{k} v_{k}$ . Now apply $(H - σ I)^{- 1}$ once — that is, solve $(H - σ I) w = u$ , which is one O(N) Thomas solve. Each eigencomponent gets multiplied by its own factor:

w = k \sum \frac{c _{k}}{λ _{k} - σ} v_{k} .

The component whose eigenvalue $λ_{k}$ lies closest to the shift $σ$ has the smallest denominator, so it gets amplified the most. Repeat the solve a few dozen times, normalizing as you go, and everything else dies away: the vector converges to the eigenvector closest to $σ$ . So if we can find a $σ$ that is guaranteed to sit below the entire spectrum, the closest eigenvalue is the lowest one — the ground state, guaranteed, no luck involved.

Gershgorin: a lower bound you can read off the matrix

Gershgorin's circle theorem: every eigenvalue of a matrix lies within one of the discs centered at a diagonal entry $H_{ii}$ , with radius the sum of the magnitudes of that row's off-diagonal entries. The proof is three lines. Take $H v = λ v$ and look at the component $v_{i}$ with the largest magnitude. Row $i$ of the eigenvalue equation says

(λ - H_{ii}) v_{i} = j \neq = i \sum H_{ij} v_{j} ⟹ ∣ λ - H_{ii} ∣ \leq j \neq = i \sum ∣ H_{ij} ∣ \leq 1 \frac{∣ v _{j} ∣}{∣ v _{i} ∣} \leq j \neq = i \sum ∣ H_{ij} ∣.

Dividing by the biggest component makes every ratio at most one — that's the entire trick. Consequence: the number $min_{i} (H_{ii} - \sum_{j \neq = i} ∣ H_{ij} ∣)$ sits at or below every eigenvalue. A guaranteed lower bound on the whole spectrum, computed by just reading the matrix entries — no diagonalization, no iteration, O(N) arithmetic. That's our shift: the code takes this bound and subtracts 1 for safety.

Why this particular bound is tight — the cancellation

Here's the part that looks like luck and isn't. For our matrix the diagonal carries the kinetic term $+ 2 B / h^{2} \approx + 4147$ MeV — enormous. If the bound landed anywhere near $- 4000$ MeV, inverse iteration would crawl: the convergence rate per sweep is the ratio

contraction = \frac{λ _{1} - σ}{λ _{2} - σ},

and with $σ \approx - 4000$ that ratio would be $\approx (4000 - 31) / (4000 - 19) \approx 0.997$ — thousands of sweeps. But look at row $i$ : its two off-diagonals are each $- B / h^{2} \approx - 2074$ , so the Gershgorin radius is $\approx + 4147$ — and it cancels the kinetic diagonal almost exactly. What's left is just the potential: $H_{ii} - radius \approx U (r_{i})$ . The bound lands at the floor of the well. The run prints sigma: -51.43 — one MeV under the 51 MeV well depth, exactly as this argument predicts. With $λ_{1} = - 31.09$ and the next level near $- 19$ , the contraction is $(- 31.09 + 51.4) / (- 19 + 51.4) \approx 0.63$ per sweep — machine precision in a few dozen sweeps. The physics (kinetic scale cancels row-wise) hands the algorithm a tight shift for free.

The program

# Build your own nuclear DFT -- step 2: one bound state in a fixed well.
#
# Before any self-consistency, solve the simplest quantum problem in the file:
# a single nucleon in a FIXED Woods-Saxon potential. For a spherical potential
# the 3-D Schrodinger equation separates, and the s-wave (l = 0) radial part
# for u(r) = r R(r) is one-dimensional:
#
#     -(hbar^2/2m) u'' + U(r) u = eps u,      u(0) = u(rmax) = 0.
#
# On the mesh, u'' becomes the (u[i-1] - 2u[i] + u[i+1])/h^2 stencil, so H is
# a TRIDIAGONAL matrix. We want its lowest eigenpair, and there is an O(N)
# trick for that: pick a shift sigma below the whole spectrum (Gershgorin:
# diagonal minus the reach of the off-diagonals), then apply (H - sigma I)^-1
# over and over. Each application is a Thomas solve -- forward sweep, back
# substitution -- and the iteration converges to the lowest state because
# that's the eigenvalue closest to sigma.

pi = 3.141592653589793
hbar2_2m = 20.7355
h = 0.1
nr = 140
amass = 16.0

r = zeros(nr)
for i to nr { r[i] = (i + 1.0) * h }

# a standard fixed Woods-Saxon well: 51 MeV deep, R = 1.27 A^(1/3), a = 0.67
v0 = -51.0
rws = 1.27 * pow(amass, 1.0 / 3.0)
aws = 0.67
u_pot = zeros(nr)
for i to nr { u_pot[i] = v0 / (1.0 + exp((r[i] - rws) / aws)) }

# tridiagonal H: constant mass for now (the effective mass comes in step 5)
diag = zeros(nr)
off = zeros(nr)
for i to nr {
    diag[i] = 2.0 * hbar2_2m / (h * h) + u_pot[i]
    if i < nr - 1.0 { off[i] = 0.0 - hbar2_2m / (h * h) }
}

# Gershgorin lower bound: min over rows of (diagonal - |off-diagonals|).
# The kinetic +2/h^2 on the diagonal cancels against the two -1/h^2 couplings,
# so this lands near the bottom of the potential -- a tight shift.
sigma = diag[0]
for i to nr {
    g = diag[i]
    if i > 0 { g = g - abs(off[i - 1.0]) }
    if i < nr - 1.0 { g = g - abs(off[i]) }
    if g < sigma { sigma = g }
}
sigma = sigma - 1.0
show sigma           # ~ -52: just under the well floor, as promised

# inverse iteration with Thomas solves
u = zeros(nr)
w = zeros(nr)
tc = zeros(nr)
td = zeros(nr)
for i to nr { u[i] = 1.0 }
eps = 0.0
for sweep to 200 {
    # forward sweep of (H - sigma I) w = u
    tc[0] = off[0] / (diag[0] - sigma)
    td[0] = u[0] / (diag[0] - sigma)
    for i to nr {
        if i > 0 {
            m = diag[i] - sigma - off[i - 1.0] * tc[i - 1.0]
            if i < nr - 1.0 { tc[i] = off[i] / m }
            td[i] = (u[i] - off[i - 1.0] * td[i - 1.0]) / m
        }
    }
    # back substitution
    w[nr - 1.0] = td[nr - 1.0]
    k = nr - 2.0
    while k >= 0.0 {
        w[k] = td[k] - tc[k] * w[k + 1.0]
        k = k - 1.0
    }
    # normalize, measure the change, adopt
    nrm = 0.0
    for i to nr { nrm += w[i] * w[i] }
    nrm = sqrt(nrm)
    change = 0.0
    for i to nr {
        d = w[i] / nrm - u[i]
        change += d * d
        u[i] = w[i] / nrm
    }
    if change < 1.0e-24 { break }
}
# Rayleigh quotient: the eigenvalue of the converged vector
for i to nr {
    eps += diag[i] * u[i] * u[i]
    if i < nr - 1.0 { eps += 2.0 * off[i] * u[i] * u[i + 1.0] }
}
show eps             # the 1s level of the well, ~ -38 MeV: deeply bound

show u               # u(r) = r R(r): rises from 0, one bump, decays -- no nodes
show u_pot           # the well it lives in

Run it: sigma: -51.43 — just under the well floor, as the cancellation argument predicts — then eps: -31.09, the 1s level. The figure of

u (r)

rises from zero, makes one bump, and decays: a nodeless ground state.

Open in IDE ▸ · annotated source