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:
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.