Generating finite-difference stencils

Specific Domains Numerics
5

I was thinking about it and it seems like a vastly simpler algorithm exists. I think the following code should return the desired stencil weights w_k such that:

f^{(m)}(x_0) ≈ \sum_{k=1}^n w_k f(x_k) + O(\max_k |x_k - x_0|^{n-m})

given any arbitrary set of n > m points x_k, the desired derivative order m, and the point x_0 where the derivative is desired.

function stencil(x::AbstractVector{<:Real}, x₀::Real, m::Integer)
    ℓ = 0:length(x)-1
    m in ℓ || throw(ArgumentError("order $m ∉ $ℓ"))
    A = @. (x' - x₀)^ℓ / factorial(ℓ)
    return A \ (ℓ .== m) # vector of weights w
end

I checked, and it returns the same results (up to roundoff error) as DiffEqOperators.calculate_weights for random inputs, so it seems fine?

Am I missing something? Probably it becomes badly conditioned for large n, but in practice it seems unlikely that anyone would need n > 10, and I’m not sure how well Fornberg does in that regime either…

PS. To derive the algorithm above, just Taylor expand \sum_{k=1}^n w_k f(x_k) around x_0 to order n-1, then set all of the f^{(\ell\ne m)} coefficients to zero while setting the f^{(m)} coefficient to 1. The result is an n\times n system of linear equations in the weights w_k.

PPS. We assigned the derivation of this formula as a homework exercise in our 2024 Matrix Calculus course — see problem 2 of pset 2.

8 Likes