X = 2.0
H = 0.1
rung "The slope of a secant" {
#| The derivative of $f$ at $x$ is the limit, as $h \to 0$, of the slope of the line
#| through $(x, f(x))$ and $(x + h, f(x + h))$. Before the limit that slope is an ordinary
#| number, and for $f(x) = x^3$ the algebra gives it exactly:
#| $\frac{(x+h)^3 - x^3}{h} = 3x^2 + 3xh + h^2$.
#| The derivative is the $3x^2$; the rest, $3xh + h^2$, is the error, and it shrinks like $h$:
#| halve $h$ and the error roughly halves.
#|
#| `forward` takes the function as an argument, so each estimate on this ladder is
#| written once and works for any function.
#| $f(x) = x^3$, the function every rung on this ladder differentiates.
def cube(t) {
return t * t * t
}
#| $\frac{f(x+h) - f(x)}{h}$: the slope of the secant from $x$ to $x + h$, wrong by $O(h)$.
def forward(f, x, h) {
"The forward difference: the slope of the secant from x to x + h."
return (f(x + h) - f(x)) / h
}
exact = 3.0 * X * X
show forward(cube, X, H)
test "the forward difference of x^3 is 3x^2 + 3xh + h^2, exactly" {
assert_near(forward(cube, X, H), exact + 3.0 * X * H + H * H, 1.0e-12)
}
test "halving h halves the error, to within 1%" {
e1 = forward(cube, X, H) - exact
e2 = forward(cube, X, H / 2.0) - exact
assert_near(e1 / e2, 2.0, 0.02)
}
}
rung "A secant centered on the point" {
#| Take the secant from $x - h$ to $x + h$ instead. In $(x+h)^3 - (x-h)^3$ every even
#| power of $h$ cancels, and $\frac{(x+h)^3 - (x-h)^3}{2h} = 3x^2 + h^2$.
#| The $3xh$ term is gone. What is left is $h^2$, so halving $h$ quarters the error.
#| Rung 1's two tests still run here, unchanged, and still pass; a ladder requires
#| every earlier test to keep passing. The new estimate gets a new name, `centered`,
#| so `forward` and its tests stay as they were.
#| $\frac{f(x+h) - f(x-h)}{2h}$: the odd powers of $h$ survive, the even ones cancel, so the error is $O(h^2)$.
def centered(f, x, h) {
"The centered difference: the slope of the secant from x - h to x + h."
return (f(x + h) - f(x - h)) / (2.0 * h)
}
show centered(cube, X, H)
test "the centered difference of x^3 is 3x^2 + h^2, exactly" {
assert_near(centered(cube, X, H), exact + H * H, 1.0e-12)
}
test "halving h quarters the error" {
e1 = centered(cube, X, H) - exact
e2 = centered(cube, X, H / 2.0) - exact
assert_near(e1 / e2, 4.0, 1.0e-9)
}
}
rung "Combining two centered slopes to cancel the h^2 term" {
#| The centered error is $h^2$ at step $h$ and $h^2/4$ at step $h/2$. So four times the
#| finer estimate, minus the coarser one, over three, has no $h^2$ term at all:
#|
#| $$\frac{4\left(3x^2 + \tfrac{h^2}{4}\right) - \left(3x^2 + h^2\right)}{3} = 3x^2.$$
#|
#| For a cubic that is the derivative
#| exactly, to rounding. This is Richardson extrapolation, and it uses the same idea as the
#| previous rung: combine two estimates so that their leading error terms cancel.
#| $\frac{4 D_{h/2} - D_h}{3}$, where $D_h$ is the centered slope at step $h$: the $h^2$ terms cancel and the error is $O(h^4)$.
def richardson(f, x, h) {
"Combine the centered slopes at h and h/2 so the h^2 error cancels."
return (4.0 * centered(f, x, h / 2.0) - centered(f, x, h)) / 3.0
}
show richardson(cube, X, H)
test "for x^3 the extrapolated slope is 3x^2 to rounding" {
assert_near(richardson(cube, X, H), exact, 1.0e-12)
}
}
rung "On a function where nothing is exact" {
#| For $e^x$ the Taylor series never ends, so none of the three is exact and what remains
#| is the order of each: halve $h$ and the forward error falls by about 2, the centered by
#| about 4, the extrapolated by about 16. At $x = 1$, $h = 0.1$ the three errors are
#| $0.14$, $0.0045$ and $5.7 \times 10^{-7}$. On a log-log plot they are lines of slope 1, 2
#| and 4, and the figure draws them from $h = 0.4$ down to $h = 0.0125$.
ex = fn(t) -> exp(t)
e = exp(1.0)
f1 = forward(ex, 1.0, H) - e
f2 = forward(ex, 1.0, H / 2.0) - e
c1 = centered(ex, 1.0, H) - e
c2 = centered(ex, 1.0, H / 2.0) - e
r1 = richardson(ex, 1.0, H) - e
r2 = richardson(ex, 1.0, H / 2.0) - e
show f1
show c1
show r1
test "halving h shrinks the three errors by factors of 2, 4 and 16, each to within 2%" {
assert_near(f1 / f2, 2.0, 0.04)
assert_near(c1 / c2, 4.0, 0.01)
assert_near(r1 / r2, 16.0, 0.02)
}
NH = 6
hs = zeros(NH)
err_f = zeros(NH)
err_c = zeros(NH)
err_r = zeros(NH)
for i to NH {
h = 0.4 / pow(2.0, i)
hs[i] = h
err_f[i] = abs(forward(ex, 1.0, h) - e)
err_c[i] = abs(centered(ex, 1.0, h) - e)
err_r[i] = abs(richardson(ex, 1.0, h) - e)
}
figure("Error of three derivative estimates of $e^x$ at $x = 1$", "h", "|error|")
xscale("log")
yscale("log")
plot(hs, err_f, "forward: slope 1")
plot(hs, err_c, "centered: slope 2")
plot(hs, err_r, "extrapolated: slope 4")
test "the figure's steps run from h = 0.4 down to h = 0.0125, five halvings" {
assert_near(hs[0], 0.4, 1.0e-15)
assert_near(hs[NH - 1], 0.0125, 1.0e-15)
}
test "the log-log slopes are 1, 2 and 4 between the two finest steps" {
dh = log(hs[NH - 1]) - log(hs[NH - 2])
assert_near((log(err_f[NH - 1]) - log(err_f[NH - 2])) / dh, 1.0, 0.03)
assert_near((log(err_c[NH - 1]) - log(err_c[NH - 2])) / dh, 2.0, 0.01)
assert_near((log(err_r[NH - 1]) - log(err_r[NH - 2])) / dh, 4.0, 0.01)
}
}
# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors