# 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:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![dorn-gerhard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dorn-gerhard/32/27153_2.png) [@dorn-gerhard](https://discourse.julialang.org/u/dorn-gerhard)\
**Post date:** [August 17, 2022, 12:46pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/1 "2022-08-17T12:46:13Z")

</div>

Regarding the second suggested approach in the topic title: **discrete derivatives** (without interpolating data)  
Is there a Julia package that performs discrete derivatives (e.g. using some [higher order stencils](https://de.wikipedia.org/wiki/Differenzenquotient#Gew%C3%B6hnliche_Differenzenquotienten_h%C3%B6herer_Ableitungs-_und_Fehlerordnung_2)) also for non-uniform grids?

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)

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [August 20, 2022, 8:33am UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/2 "2022-08-20T08:33:20Z")

</div>

FWIW, here are some links to codes based on Fornberg’s (1998 and 1988) articles: [Julia code](https://github.com/bekaiser/Fornberg/blob/master/Lagrange_poly_derivatives.jl) and [Python code](https://github.com/bjodah/finitediff)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 20, 2022, 10:31am UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/3 "2022-08-20T10:31:49Z")

</div>

Here’s a maintained version of Fornberg:

> <https://github.com/SciML/DiffEqOperators.jl/blob/master/src/derivative_operators/fornberg.jl>

With prebuilt operators:

> **[GitHub - SciML/DiffEqOperators.jl: Linear operators for discretizations of...](https://github.com/SciML/DiffEqOperators.jl)**
>
> Linear operators for discretizations of differential equations and scientific machine learning (SciML) - GitHub - SciML/DiffEqOperators.jl: Linear operators for discretizations of differential equa...

Of course, doing operators is a slow way to implement it so I wouldn’t recommend it, but the package exists.

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [August 20, 2022, 2:28pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/4 "2022-08-20T14:28:36Z")

</div>

Isn’t that package in the middle of being deprecated?

---

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

---

<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:26pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/6 "2022-08-20T16:26:36Z")

</div>

> [@stevengj](#):
>
> I was thinking about it and it seems like vastly simpler algorithm exists

In fact, simply by passing `Integer` and `Rational` arguments (so that `/` becomes exact), the _same_ code can give _exact_ rational results. For example:

```julia
julia> stencil(-4:4, 0//1, 4)
9-element Vector{Rational{Int64}}:
    7//240
   -2//5
  169//60
 -122//15
   91//8
 -122//15
  169//60
   -2//5
    7//240

```

matching the last line of Table 1 in the Fornberg paper:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/8/2850737d9b37838b243c8aed1fdcc3765843f80e.png)  
(And I can reproduce the other tables in Fornberg as well.)

Gotta love generic programming!

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [August 20, 2022, 5:28pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/7 "2022-08-20T17:28:14Z")

</div>

This is the kind of example that should be used to “sell” Julia! The claims of raw speed rarely rise to the point where doubters can be persuaded, but this is quite another matter.

Well done!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 20, 2022, 5:43pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/8 "2022-08-20T17:43:41Z")

</div>

> [@stillyslalom](#):
>
> Isn’t that package in the middle of being deprecated?

Yeah, just like ODE.jl. It’ll stay alive for at least another decade as something I tell people not to do.

It’s just for most applications (i.e. anything nonlinear) using a “stencil” calculation is simply worse than fusing the nonlinearity in and doing a straight flat loop. By large margin. We’ve demonstrated this over and over again, so now we’re at the point where we basically say, we’ll keep the stencil thing because it’s what people learn in school, but please don’t do that under the vague notion that it’s somehow high-performance. Yes, you can get a “fast” implementation of a stencil, but it’s still worse than naively fusing, unless your equations are actually linear.

> [@stevengj](#):
>
> I checked, and it returns the same results (up to roundoff error) as `DiffEqOperators.calculate_weights` for random inputs, so it seems fine?

Huh, cool. I would add that the last term should just impose the summation condition. There’s some things out there that say it matters over a large integration, but if you’re seeing the same results up to round-off error then 🤷.

---

<div class="post-metadata">

**Author:** ![tobydriscoll](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tobydriscoll/32/1843_2.png) [@tobydriscoll](https://discourse.julialang.org/u/tobydriscoll)\
**Post date:** [August 26, 2022, 2:30pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/9 "2022-08-26T14:30:40Z")

</div>

> [@stevengj](#):
>
> Am I missing something? Probably it becomes badly conditioned for large n nn , but in practice it seems unlikely that anyone would need n \> 10 n\>10n \> 10 , and I’m not sure how well Fornberg does in that regime either…

What a fun little code!

I’d be cautious about this version, though. When I try for `x=1:5` to evaluate the first derivative at an endpoint, it’s lost 2 digits of accuracy relative to Chris’ code (and exact rational results). It’s not terrible, but given how sensitive FD methods are to roundoff already, it’s not great, either. It also loses exact antisymmetry at the ends, which I suppose is just unaesthetic.

I’d argue that it’s safe only in the `Rational` case. @ChrisRackauckas’s code can turn out rational results, too, however.

Fornberg’s 1988 intro says that “the main anticipated application” of his method is to “dynamically changing grids.” I’m not sure how that ever played out…

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [September 14, 2022, 7:06pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/10 "2022-09-14T19:06:23Z")

</div>

@stevengj’s stencil code turned out to be an incredibly painless way to calculate wave speeds from wave arrival times at unevenly-spaced pressure transducers with measurement uncertainties:

```julia
function wavespeed(xt, N=1)
    map(axes(xt, 1)) do i
        window = max(1, i - N):min(size(xt, 1), i + N)
        mapreduce(*, +, stencil(xt.t[window], xt.t[i], 1), xt.x[window])
    end
end

```

```julia
julia> xt
10×3 DataFrame
 Row │ name x t
     │ String7 Measurem… Measurem…        
─────┼────────────────────────────────────────
   1 │ PT1 0.63±0.004 -0.002128±1.0e-5
   2 │ PT2 1.22±0.004 -0.001209±1.0e-5
   3 │ PT3 2.072±0.004 0.0±1.0e-5
   4 │ PT4 3.952±0.004 0.002638±1.0e-5
   5 │ PT5 4.606±0.004 0.003565±1.0e-5
   6 │ PT6 5.207±0.004 0.004423±1.0e-5
   7 │ PT7 5.536±0.004 0.004879±1.0e-5
   8 │ PT8 6.087±0.004 0.00567±1.0e-5
   9 │ PT9 6.607±0.004 0.006409±1.0e-5
  10 │ PT10 6.862±0.004 0.006772±1.0e-5

julia> wavespeed(xt)
10-element Vector{Measurement{Float64}}:
 642.0 ± 12.0
 669.1 ± 5.7
 707.2 ± 5.9
 707.4 ± 8.6
 702.9 ± 6.5
 714.0 ± 15.0
 712.0 ± 15.0
 700.2 ± 7.5
 703.0 ± 19.0
 702.0 ± 31.0

```

Thanks!

---

<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:** [September 14, 2022, 8:44pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/11 "2022-09-14T20:44:35Z")

</div>

> [@stillyslalom](#):
>
> @stevengj’s stencil code turned out to be an incredibly painless way to calculate wave speeds from wave arrival times at unevenly-spaced pressure transducers with measurement uncertainties:

Cool! But as @tobydriscoll pointed out above, my code is susceptible to roundoff errors for floating-point calculations (probably because my matrix `A` is badly conditioned), especially as the number of points becomes larger, whereas `DiffEqOperators.calculate_weights` is better-behaved and otherwise equivalent, so I would tend to recommend the latter in practice (though I don’t know if it handles `Measurement` values?), especially for window sizes larger than a few points.

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [September 14, 2022, 8:51pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/12 "2022-09-14T20:51:53Z")

</div>

For measurements where the inherent uncertainty (we’re talking about an eyeball and a tape measure here) is many orders of magnitude larger than floating-point roundoff errors, I’ll take your simple, `Measurement`-compatible approach over adding a dependency with an active deprecation warning.

---

<div class="post-metadata">

**Author:** ![JeffreySarnoff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeffreysarnoff/32/1980_2.png) [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)\
**Post date:** [September 14, 2022, 8:54pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/13 "2022-09-14T20:54:22Z")

</div>

Is there Julia code that does this for nonuniform spacing, one-sided differences?

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [September 14, 2022, 8:58pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/14 "2022-09-14T20:58:57Z")

</div>

> [@stevengj](#):
>
> though I don’t know if it handles `Measurement` values?

It seems to work:

```julia
using DiffEqOperators: calculate_weights

function wavespeed2(xt, N=1)
    map(axes(xt, 1)) do i
        window = max(1, i - N):min(size(xt, 1), i + N)
        mapreduce(*, +, calculate_weights(1, xt.t[i], xt.t[window]), xt.x[window])
    end
end

V2 = wavespeed2(xt)
10-element Vector{Measurement{Float64}}:
 642.0 ± 12.0
 669.1 ± 5.7
 707.2 ± 5.9
 707.4 ± 8.6
 702.9 ± 6.5
 714.0 ± 15.0
 712.0 ± 15.0
 700.2 ± 7.5
 703.0 ± 19.0
 702.0 ± 31.0

julia> wavespeed2(xt) - wavespeed(xt)
10-element Vector{Measurement{Float64}}:
        0.0 ± 0.0
        0.0 ± 2.7e-15
 -1.137e-13 ± 1.4e-15
  4.547e-13 ± 9.5e-15
        0.0 ± 1.3e-14
        0.0 ± 3.8e-14
  -4.55e-13 ± 3.3e-14
   4.55e-13 ± 3.2e-14
   1.82e-12 ± 1.1e-13
        0.0 ± 0.0

```

---

<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:** [September 14, 2022, 9:31pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/15 "2022-09-14T21:31:34Z")

</div>

> [@JeffreySarnoff](#):
>
> Is there Julia code that does this for nonuniform spacing, one-sided differences?

Both my code and `DiffEqOperators.calculate_weights` handle one-sided differences as a special case.

---

<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:** [September 14, 2022, 9:40pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/16 "2022-09-14T21:40:58Z")

</div>

> [@stillyslalom](#):
>
> For measurements where the inherent uncertainty (we’re talking about an eyeball and a tape measure here) is many orders of magnitude larger than floating-point roundoff errors,

Right. Probably it should be fine in that context up to about 20 points, for which the condition number of my A is about 10^8 (for centered differences) and so it should still be accurate to half of double precision. Beyond that the `factorial` function will throw an integer-overflow error (at least, as the code is currently written), so at least it will fail loudly.

---

<div class="post-metadata">

**Author:** ![mastrof](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mastrof/32/17443_2.png) [@mastrof](https://discourse.julialang.org/u/mastrof)\
**Post date:** [September 22, 2022, 1:39pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/17 "2022-09-22T13:39:19Z")

</div>

> [@ChrisRackauckas](#):
>
> It’s just for most applications (i.e. anything nonlinear) using a “stencil” calculation is simply worse than fusing the nonlinearity in and doing a straight flat loop. By large margin. We’ve demonstrated this over and over again, so now we’re at the point where we basically say, we’ll keep the stencil thing because it’s what people learn in school, but please don’t do that under the vague notion that it’s somehow high-performance. Yes, you can get a “fast” implementation of a stencil, but it’s still worse than naively fusing, unless your equations are actually linear.

Can I ask what is meant here by “fusing the nonlinearity in and doing a straight flat loop”?  
I’m genuinely interested in this but my knowledge of numerical integration is limited to the finite difference stuff you learn in introductory uni courses

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 22, 2022, 10:08pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/18 "2022-09-22T22:08:41Z")

</div>

See [Code Optimization for Differential Equations · DifferentialEquations.jl](https://diffeq.sciml.ai/stable/tutorials/faster_ode_example/#Example-Accelerating-Linear-Algebra-PDE-Semi-Discretization)

---

<div class="post-metadata">

**Author:** ![ohno](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohno/32/38937_2.png) [@ohno](https://discourse.julialang.org/u/ohno)\
**Post date:** [July 23, 2024, 9:56pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/19 "2024-07-23T21:56:17Z")

</div>

I developed a new package based on your post. Thank you.

> **[GitHub - ohno/FiniteDifferenceMatrices.jl: A Julia package for discrete...](https://github.com/ohno/FiniteDifferenceMatrices.jl)**
>
> A Julia package for discrete approximations of differential operators - ohno/FiniteDifferenceMatrices.jl

---

<div class="post-metadata">

**Author:** ![CosmoProf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cosmoprof/32/15516_2.png) [@CosmoProf](https://discourse.julialang.org/u/CosmoProf)\
**Post date:** [July 23, 2024, 10:39pm UTC](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/20 "2024-07-23T22:39:44Z")

</div>

Thank you,  
I already created a stencil like that by hand, but the physics is actually more complicated and I had to make a more complex stencil (leading edge different than trailing edge, anti-reflection factors, etc).