Notebooks · posts/finite-differences-and-spectral-methods

Finite differences and spectral methods

Adam Zettel · October 14, 2026 · written with an agent

Two ways to differentiate a function you only know at grid points. A central difference uses the two neighbors; a spectral derivative uses every point. On a smooth periodic function the spectral one reaches rounding error with a few dozen points. Put a jump in the function and, away from the jump, the central difference is the accurate one.

f^{'} (x_{j}) \approx \frac{f _{j + 1} - f _{j - 1}}{2 h} vs f^{'} (x_{j}) \approx k \sum ik f^_{k} e^{ik x_{j}}

Take a periodic function sampled at $N$ equally spaced points, $x_{j} = j h$ with $h = 2 π / N$ , and ask for its derivative at those points. The central difference $(f_{j + 1} - f_{j - 1}) /2 h$ looks at two neighbors. Its error is $\frac{h ^{2}}{6} f^{'''}$ , so doubling $N$ divides it by 4. The spectral derivative fits the trigonometric polynomial through all $N$ samples and differentiates that exactly: transform, multiply mode $k$ by $ik$ , transform back. Every sample contributes to every derivative value.

A smooth function

For $e^{s i n x}$ , which is smooth and periodic, the Fourier coefficients fall off faster than any power of $k$ , and so does the spectral derivative's error. The two lines for $e^{s i n x}$ in the figure show it: the spectral error reaches rounding by 32 points, while the central difference is still losing a factor of 4 per doubling.

[claim: spectral-derivative]

A jump

Now add a unit step to $sin x$ , so that the function jumps at $x = 0$ and $x = π$ . Next to a jump neither method is right, and the largest error over all points grows with $N$ for both. Away from the jumps they behave differently. The central difference at a point $π /4$ or farther from a jump only sees $sin x$ , and its error falls like $h^{2}$ as before. The spectral derivative at that point still includes the jump: in the spectral differentiation matrix the weight of a sample at distance $d$ is $\frac{1}{2} cot (d /2)$ in size, and that does not shrink as $N$ grows. The two lines marked far points show it: the central difference's error keeps falling, the spectral one stays near $0.7$ .

[claim: central-difference-near-a-jump]

What the two claims cover

The two claims, [claim: spectral-derivative] and [claim: central-difference-near-a-jump], cover one smooth function and one function with jumps, both periodic, on equally spaced points. Functions on an interval that are not periodic, which spectral methods handle with Chebyshev points, and the cost of each method per time step, are not tested here.

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 52–190

#: pi
PI_K = 4.0 * atan(1.0)
SIZES = [8, 16, 32, 64, 128]
#: the points the jump comparison uses: at least FAR from either jump
FAR = PI_K / 4.0

def smooth(x: Num) -> Num {
    return exp(sin(x))
}

def smooth_slope(x: Num) -> Num {
    return cos(x) * exp(sin(x))
}

#| sin x, plus 1 on [0, pi): jumps at x = 0 and x = pi.
def jumpy(x: Num) -> Num {
    step = 0.0
    if x < PI_K {
        step = 1.0
    }
    return sin(x) + step
}

#| Distance from x in [0, 2 pi) to the nearer jump of jumpy.
def jump_distance(x: Num) -> Num {
    return min(abs(x - PI_K), min(x, 2.0 * PI_K - x))
}

#| Samples of the smooth function (which = 0) or the jumpy one (which = 1) on n points.
def samples(which: Int, n: Int) {
    h = 2.0 * PI_K / n
    v = zeros(n)
    for j to n {
        if which == 0 {
            v[j] = smooth(j * h)
        }
        else {
            v[j] = jumpy(j * h)
        }
    }
    return v
}

#| The spectral derivative of the samples v: a direct discrete Fourier transform, mode k times i k, back.
def spectral_slope(v) {
    n: Int = len(v)
    h = 2.0 * PI_K / n
    re = zeros(n)
    im = zeros(n)
    for k to n {
        for j to n {
            re[k] = re[k] + v[j] * cos(k * j * h)
            im[k] = im[k] - v[j] * sin(k * j * h)
        }
    }
    d = zeros(n)
    for j to n {
        for k to n {
            #: the wavenumber of mode k, the Nyquist mode dropped
            w = 0.0
            if 2 * k < n {
                w = k * 1.0
            }
            if 2 * k > n {
                w = k - n * 1.0
            }
            d[j] = d[j] + w * (-im[k] * cos(k * j * h) - re[k] * sin(k * j * h)) / n
        }
    }
    return d
}

#| The central difference of the periodic samples v.
def central_slope(v) {
    n: Int = len(v)
    h = 2.0 * PI_K / n
    return [(v[(j + 1) % n] - v[(j + n - 1) % n]) / (2.0 * h) for j to n]
}

#| The largest error of a derivative d on n points: all points for the smooth function, the far ones for the jumpy one.
def worst_error(which: Int, d) -> Num {
    n: Int = len(d)
    h = 2.0 * PI_K / n
    worst = 0.0
    for j to n {
        if which == 0 {
            worst = max(worst, abs(d[j] - smooth_slope(j * h)))
        }
        else if jump_distance(j * h) >= FAR {
            worst = max(worst, abs(d[j] - cos(j * h)))
        }
    }
    return worst
}

n_sizes: Int = len(SIZES)
ns = zeros(n_sizes)
smooth_spectral = zeros(n_sizes)
smooth_central = zeros(n_sizes)
jump_spectral = zeros(n_sizes)
jump_central = zeros(n_sizes)
for k to n_sizes {
    n: Int = int(SIZES[k])
    ns[k] = n * 1.0
    smooth_spectral[k] = max(worst_error(0, spectral_slope(samples(0, n))), 1.0e-16)
    smooth_central[k] = worst_error(0, central_slope(samples(0, n)))
    jump_spectral[k] = worst_error(1, spectral_slope(samples(1, n)))
    jump_central[k] = worst_error(1, central_slope(samples(1, n)))
}
show smooth_spectral
show smooth_central
show jump_spectral
show jump_central

figure("Largest derivative error against the number of points", "N", "largest error")
xscale("log")
yscale("log")
plot(ns, smooth_spectral, "spectral, exp(sin x)")
plot(ns, smooth_central, "central difference, exp(sin x)")
plot(ns, jump_spectral, "spectral, sin x + step, far points")
plot(ns, jump_central, "central difference, sin x + step, far points")

test "the spectral error on exp(sin x) is below 1e-12 by N = 32" {
    assert_lt(smooth_spectral[2], 1.0e-12)
}
test "the central difference error falls by 3.9 to 4.1 per doubling from N = 32 on both functions" {
    for k from 2 to n_sizes - 1 {
        assert_gt(smooth_central[k] / smooth_central[k + 1], 3.9)
        assert_lt(smooth_central[k] / smooth_central[k + 1], 4.1)
        assert_gt(jump_central[k] / jump_central[k + 1], 3.9)
        assert_lt(jump_central[k] / jump_central[k + 1], 4.1)
    }
}
test "the spectral error at the points far from the jumps stays above 0.7" {
    assert_gt(min(jump_spectral[1:]), 0.7)
}

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

#: pi
PI_K = 4.0 * atan(1.0)
SIZES = [8, 16, 32, 64, 128]
#: the points the jump comparison uses: at least FAR from either jump
FAR = PI_K / 4.0

def smooth(x: Num) -> Num {
    return exp(sin(x))
}

def smooth_slope(x: Num) -> Num {
    return cos(x) * exp(sin(x))
}

#| sin x, plus 1 on [0, pi): jumps at x = 0 and x = pi.
def jumpy(x: Num) -> Num {
    step = 0.0
    if x < PI_K {
        step = 1.0
    }
    return sin(x) + step
}

#| Distance from x in [0, 2 pi) to the nearer jump of jumpy.
def jump_distance(x: Num) -> Num {
    return min(abs(x - PI_K), min(x, 2.0 * PI_K - x))
}

#| Samples of the smooth function (which = 0) or the jumpy one (which = 1) on n points.
def samples(which: Int, n: Int) {
    h = 2.0 * PI_K / n
    v = zeros(n)
    for j to n {
        if which == 0 {
            v[j] = smooth(j * h)
        }
        else {
            v[j] = jumpy(j * h)
        }
    }
    return v
}

#| The spectral derivative of the samples v: a direct discrete Fourier transform, mode k times i k, back.
def spectral_slope(v) {
    n: Int = len(v)
    h = 2.0 * PI_K / n
    re = zeros(n)
    im = zeros(n)
    for k to n {
        for j to n {
            re[k] = re[k] + v[j] * cos(k * j * h)
            im[k] = im[k] - v[j] * sin(k * j * h)
        }
    }
    d = zeros(n)
    for j to n {
        for k to n {
            #: the wavenumber of mode k, the Nyquist mode dropped
            w = 0.0
            if 2 * k < n {
                w = k * 1.0
            }
            if 2 * k > n {
                w = k - n * 1.0
            }
            d[j] = d[j] + w * (-im[k] * cos(k * j * h) - re[k] * sin(k * j * h)) / n
        }
    }
    return d
}

#| The central difference of the periodic samples v.
def central_slope(v) {
    n: Int = len(v)
    h = 2.0 * PI_K / n
    return [(v[(j + 1) % n] - v[(j + n - 1) % n]) / (2.0 * h) for j to n]
}

#| The largest error of a derivative d on n points: all points for the smooth function, the far ones for the jumpy one.
def worst_error(which: Int, d) -> Num {
    n: Int = len(d)
    h = 2.0 * PI_K / n
    worst = 0.0
    for j to n {
        if which == 0 {
            worst = max(worst, abs(d[j] - smooth_slope(j * h)))
        }
        else if jump_distance(j * h) >= FAR {
            worst = max(worst, abs(d[j] - cos(j * h)))
        }
    }
    return worst
}

n_sizes: Int = len(SIZES)
ns = zeros(n_sizes)
smooth_spectral = zeros(n_sizes)
smooth_central = zeros(n_sizes)
jump_spectral = zeros(n_sizes)
jump_central = zeros(n_sizes)
for k to n_sizes {
    n: Int = int(SIZES[k])
    ns[k] = n * 1.0
    smooth_spectral[k] = max(worst_error(0, spectral_slope(samples(0, n))), 1.0e-16)
    smooth_central[k] = worst_error(0, central_slope(samples(0, n)))
    jump_spectral[k] = worst_error(1, spectral_slope(samples(1, n)))
    jump_central[k] = worst_error(1, central_slope(samples(1, n)))
}
show smooth_spectral
show smooth_central
show jump_spectral
show jump_central

figure("Largest derivative error against the number of points", "N", "largest error")
xscale("log")
yscale("log")
plot(ns, smooth_spectral, "spectral, exp(sin x)")
plot(ns, smooth_central, "central difference, exp(sin x)")
plot(ns, jump_spectral, "spectral, sin x + step, far points")
plot(ns, jump_central, "central difference, sin x + step, far points")

test "the spectral error on exp(sin x) is below 1e-12 by N = 32" {
    assert_lt(smooth_spectral[2], 1.0e-12)
}
test "the central difference error falls by 3.9 to 4.1 per doubling from N = 32 on both functions" {
    for k from 2 to n_sizes - 1 {
        assert_gt(smooth_central[k] / smooth_central[k + 1], 3.9)
        assert_lt(smooth_central[k] / smooth_central[k + 1], 4.1)
        assert_gt(jump_central[k] / jump_central[k + 1], 3.9)
        assert_lt(jump_central[k] / jump_central[k + 1], 4.1)
    }
}
test "the spectral error at the points far from the jumps stays above 0.7" {
    assert_gt(min(jump_spectral[1:]), 0.7)
}

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