# Generating finite-difference stencils

**URL:** <https://discourse.julialang.org/t/generating-finite-difference-stencils/85876>\
**Category:** Numerics\
**Tags:** differentiation\
**Created:** [August 17, 2022, 12:46pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876 "2022-08-17T12:46:13Z")\
**Posts on this page:** 1\
**Showing post:** 5

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 20, 2022, 4:19pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/5 "2022-08-20T16:19:17Z")

</div>

> [@dorn-gerhard](#):
>
> I am looking for a package that calculates the weights for any order of derivative (including the 0th derivative, corresponding to interpolation), approximated to any order of accuracy on an arbitrary grid in one dimension (as stated in that [paper: “Generation of Finite Difference Formulas on Arbitrarily Spaced Grids - Bengt Fornberg”](https://www.ams.org/journals/mcom/1988-51-184/S0025-5718-1988-0935077-0/S0025-5718-1988-0935077-0.pdf)

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.

```julia
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](https://github.com/mitmath/matrixcalc/tree/iap2024) — see [problem 2 of pset 2](https://github.com/mitmath/matrixcalc/blob/iap2024/psets/pset2sol.pdf).

---

_[View the full topic](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876)._