This site is version 0.0.1. Most of the code on this site is drafted by an agent and needs review. I apologize
for the AI-speak, we are working on rewriting it!
Adam Zettel · September 25, 2026 · written with an agent
Numerov's method integrates y′′=y outward from y(0)=1,y(h)=e−h, aiming at the decaying solution e−x. By x=9 it is off by a factor of 30, and by x=18 by 2×109. Rewrite decaying so that it follows e−x to Numerov's own accuracy, within 2.7×10−5, on both intervals.
y′′=y,y(x)=e−x:outward, the error grows as e2x
Numerov's method is the standard way to integrate y′′=f(x)y, which is
the form of the radial Schrödinger equation. For y′′=k2y on a grid of
spacing h it is the three-term recurrence
accurate to fourth order in h. The
equation has two solutions, e−x and ex, and the recurrence carries
both. The two starting values are never exactly on the decaying one (they
differ from the recurrence's own decaying solution by its truncation error),
so a little of ex is always mixed in. Going outward, that part grows by
e2x relative to the part you want. In the region where a bound state's
wavefunction decays, every outward integration has this problem.
The function must return y0,…,ym at xi=ih, normalized to
y0=1. The tests compare it with e−x at every point. Their
tolerances come from posts/refs/make-numerov-find-the-decaying-solution.py,
which runs the working method in numpy: its largest relative error is
8.87×10−6 for h=0.15 on [0,9],5.6×10−7 for
h=0.075, and 1.84×10−5 for h=0.15 on [0,18]. The
tests allow about 1.5 times those. Halving h divides the error by 15.9,
the factor 16 of a fourth-order method.
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.
This is a problem: its tests fail until the functions in it are written. Edit a cell and run it to try; the cells below it show whether the tests pass.
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 36–88
#| The decaying solution of $y'' = y$ at $x_i = ih$, $i = 0, \ldots, m$, with $y_0 = 1$.#| As written it runs Numerov outward from $y_0 = 1$, $y_1 = e^{-h}$.defdecaying(m, h) {
c =1.0- h * h /12.0
d =1.0+5.0* h * h /12.0
y =zeros(m +1)
y[0] =1.0
y[1] =exp(-h)
for i from1to m {
y[i +1] = (2.0* d * y[i] - c * y[i -1]) / c
}
return y
}
#| The largest relative difference between y and $e^{-x}$ on the grid of spacing h.defmax_rel_error(y:Arr, h) {
e =0.0for i tolen(y) {
e =max(e, abs(y[i] /exp(-i * h) -1.0))
}
return e
}
M_SHORT =60
M_LONG =120
H =0.15
H_HALF =0.075
err_short =max_rel_error(decaying(M_SHORT, H), H)
show err_short
# each test's bound over the reference's error on the same grid
tolerance_margin = [1.3e-5/8.873863e-6, 8.0e-7/5.576447e-7, 2.7e-5/1.835749e-5]
show tolerance_margin
test"y(0) = 1" {
assert_near(decaying(M_SHORT, H)[0], 1.0, 1.0e-15)
}
test"with h = 0.15 on [0, 9] every point is within 1.3e-5 of e^(-x), where numpy's inward run reaches 8.87e-6" {
assert_lt(max_rel_error(decaying(M_SHORT, H), H), 1.3e-5)
}
test"with h = 0.075 on [0, 9] every point is within 8e-7, where numpy reaches 5.58e-7" {
assert_lt(max_rel_error(decaying(M_LONG, H_HALF), H_HALF), 8.0e-7)
}
test"with h = 0.15 on [0, 18], where outward is off by 2e9, every point is within 2.7e-5 (numpy reaches 1.84e-5)" {
assert_lt(max_rel_error(decaying(M_LONG, H), H), 2.7e-5)
}
test"halving h divides the error by 12 to 20, as a fourth-order method's should (numpy: 15.9)" {
ratio =max_rel_error(decaying(M_SHORT, H), H) /max_rel_error(decaying(M_LONG, H_HALF), H_HALF)
assert_gt(ratio, 12.0)
assert_lt(ratio, 20.0)
}
# SPDX-License-Identifier: MIT OR Apache-2.0# SPDX-FileCopyrightText: 2026 The Viv authors