Notebooks · posts/oscillator-levels

How far a finite-difference Hamiltonian's levels sit below n + 1/2

Adam Zettel · September 21, 2026 · written with an agent

Put $H = - \frac{1}{2} \partial_{x}^{2} + \frac{1}{2} x^{2}$ on a grid and diagonalize: the levels come out at $n + \frac{1}{2}$ , a little low. How low is predictable. First-order perturbation theory gives $- \frac{h ^{2}}{32} (2 n^{2} + 2 n + 1)$ , and for the six lowest levels the measured shift matches it to within half a percent (the test allows 4%).

H^= - \frac{1}{2} \partial_{x}^{2} + \frac{1}{2} x^{2}, E_{n} = n + \frac{1}{2}

The three-point second difference is $\partial_{x}^{2} + \frac{h ^{2}}{12} \partial_{x}^{4} + \dots$ , so the grid Hamiltonian is the true one plus $- \frac{h ^{2}}{24} p^{4}$ . To first order the level moves by the expectation of that, and for the oscillator $⟨ n ∣ p^{4} ∣ n ⟩ = \frac{3}{4} (2 n^{2} + 2 n + 1)$ . The grid's error is a small perturbation of the Hamiltonian, so the perturbation theory used for the Stark shift predicts it.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 18–76

N = 160
XMAX = 8.0
H = 2.0 * XMAX / (N + 1.0)
X = linspace(-XMAX + H, XMAX - H, N)

A = zeros(N, N)
for i to N {
    A[i, i] = 1.0 / (H * H) + 0.5 * X[i] * X[i]
    if i + 1.0 < N {
        A[i, i + 1.0] = -0.5 / (H * H)
        A[i + 1.0, i] = -0.5 / (H * H)
    }
}
levels = zeros(N)
states = zeros(N, N)
eig_sym(A, levels, states)

LEVELS = 6
measured = zeros(LEVELS)
predicted = zeros(LEVELS)
for n to LEVELS {
    measured[n] = levels[n] - (n + 0.5)
    predicted[n] = -(H * H / 32.0) * (2.0 * n * n + 2.0 * n + 1.0)
}
show H
show measured
show predicted

#| Column n of the eigenvector matrix, as a wavefunction normalized on the grid.
def wavefunction(n) {
    psi = zeros(N)
    for i to N {
        psi[i] = states[i, n] / sqrt(H)
    }
    if psi[N / 2.0 + 6.0] < 0.0 {
        psi = -psi
    }
    return psi
}

figure("The three lowest eigenfunctions", "$x$", "$\\psi$", -6.0, 6.0, -0.8, 0.8)
plot(X, wavefunction(0.0), "$n = 0$")
plot(X, wavefunction(1.0), "$n = 1$")
plot(X, wavefunction(2.0), "$n = 2$")

test "the ground state is at 1/2, to a part in a thousand" {
    assert_near(levels[0], 0.5, 5.0e-4)
}
test "each of the first six levels is low, because the grid's p^4 term is negative" {
    assert_lt(max(measured), 0.0)
}
test "each shift is within 4% of what perturbation theory predicts" {
    for n to LEVELS {
        assert_near(measured[n], predicted[n], 0.04 * abs(predicted[n]))
    }
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

N = 160
XMAX = 8.0
H = 2.0 * XMAX / (N + 1.0)
X = linspace(-XMAX + H, XMAX - H, N)

A = zeros(N, N)
for i to N {
    A[i, i] = 1.0 / (H * H) + 0.5 * X[i] * X[i]
    if i + 1.0 < N {
        A[i, i + 1.0] = -0.5 / (H * H)
        A[i + 1.0, i] = -0.5 / (H * H)
    }
}
levels = zeros(N)
states = zeros(N, N)
eig_sym(A, levels, states)

LEVELS = 6
measured = zeros(LEVELS)
predicted = zeros(LEVELS)
for n to LEVELS {
    measured[n] = levels[n] - (n + 0.5)
    predicted[n] = -(H * H / 32.0) * (2.0 * n * n + 2.0 * n + 1.0)
}
show H
show measured
show predicted

#| Column n of the eigenvector matrix, as a wavefunction normalized on the grid.
def wavefunction(n) {
    psi = zeros(N)
    for i to N {
        psi[i] = states[i, n] / sqrt(H)
    }
    if psi[N / 2.0 + 6.0] < 0.0 {
        psi = -psi
    }
    return psi
}

figure("The three lowest eigenfunctions", "$x$", "$\\psi$", -6.0, 6.0, -0.8, 0.8)
plot(X, wavefunction(0.0), "$n = 0$")
plot(X, wavefunction(1.0), "$n = 1$")
plot(X, wavefunction(2.0), "$n = 2$")

test "the ground state is at 1/2, to a part in a thousand" {
    assert_near(levels[0], 0.5, 5.0e-4)
}
test "each of the first six levels is low, because the grid's p^4 term is negative" {
    assert_lt(max(measured), 0.0)
}
test "each shift is within 4% of what perturbation theory predicts" {
    for n to LEVELS {
        assert_near(measured[n], predicted[n], 0.04 * abs(predicted[n]))
    }
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors