fornberg(2020) weights based on hermite-based finite difference

src/derivative_operators/fornberg.jl

@@ -1,20 +1,27 @@
 #############################################################
-# Fornberg algorithm
+#= Fornberg algorithm
 
-# This implements the Fornberg (1988) algorithm (https://doi.org/10.1090/S0025-5718-1988-0935077-0)
-# to obtain Finite Difference weights over arbitrary points to arbitrary order.
+This implements the Fornberg (1988) algorithm (https://doi.org/10.1090/S0025-5718-1988-0935077-0)
+and hermite-based finite difference Fornberg(2020) algorithm (https://doi.org/10.1093/imanum/draa006)
+to obtain Finite Difference weights over arbitrary points to arbitrary order.
 
-function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real
-    #=
+Inputs:
         order: The derivative order for which we need the coefficients
         x0   : The point in the array 'x' for which we need the coefficients
         x    : A dummy array with relative coordinates, e.g., central differences
                need coordinates centred at 0 while those at boundaries need
                coordinates starting from 0 to the end point
+        dfdx : optional argument to consider weights of the first-derivative of function or not
+                if    dfdx == false (default kwarg), implies Fornberg(1988)
+                      dfdx == true,                  implies Fornberg(2020)
 
-        The approximation order of the stencil is automatically determined from
-        the number of requested stencil points.
-    =#
+    Outputs:
+        if dfdx == false (default kwarg),   _C : weights to approximate derivative of required order using function values only.
+                                 else,   _D,_E : weights to approximate derivative of required order using function and its first- 
+                                                 derivative values respectively.                                             
+
+=#
+function calculate_weights(order::Int, x0::T, x::AbstractVector; dfdx::Bool = false) where T<:Real
     N = length(x)
     @assert order < N "Not enough points for the requested order."
     M = order
@@ -56,5 +63,23 @@ function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real
     =#
     _C = C[:,end]
     _C[div(N,2)+1] -= sum(_C)
-    return _C
+     if dfdx == false
+        return _C
+     else
+        A = x .- x';
+        s = sum(1 ./ (A + I(N)), dims = 1) .- 1;
+        cp = factorial.(0:M);
+        cc = C./cp'
+        c̃ = zeros(N, M+2);
+        for k in 1:M+1
+           c̃[:,k+1] = sum(cc[:,1:k].*cc[:,k:-1:1], dims = 2);
+        end
+        E = c̃[:,1:M+1] - (x .- x0).*c̃[:,2:M+2];
+        D = c̃[:,2:M+2] + 2*E.*s';
+        D = D.*cp';
+        E = E.*cp';
+
+        _D = D[:,end];   _E = E[:,end]
+        return _D, _E
+    end
 end

src/docstrings.jl

@@ -8,15 +8,18 @@ CenteredDifference
 
 
 @doc """
-    calculate_weights(n::Int, x₀::Real, x::Vector)
+    calculate_weights(n::Int, x₀::Real, x::Vector; dfdx::Bool = false)
 
-Return a vector `c` such that `c⋅f.(x)` approximates ``f^{(n)}(x₀)`` for a smooth `f`.
+(Default kwarg `dfdx == false`)Return a vector `c` such that `c⋅f.(x)` approximates ``f^{(n)}(x₀)`` for a smooth `f`.
+(kwarg `dfdx == true`)Returnvectors `d` and `e` such that `d⋅f.(x) + e⋅f'.(x)` approximates ``f^{(n)}(x₀)`` for a smooth `f`.
 
 The points `x` need not be evenly spaced.
 
 The stencil `c` has the highest approximation order possible given values of `f` at `length(x)` points. More precisely, if `x` has length `m`, there is a function `g` such that ``g(y) = f(y) + O(y-x₀)^{m-n+?}`` and ``c⋅f.(x) = g'(x₀)``.
+The stencil `d` and `e` has the highest approximation order possible given values of `f` at `length(x)` points. More precisely, if `x` has length `m`, there is a function `g` such that ``g(y) = f(y) + O(y-x₀)^{m-n+?}`` and ``d⋅f.(x) + e⋅f'.(x) = g'(x₀)``.
 
-The algorithm is due to [Fornberg](https://doi.org/10.1090/S0025-5718-1988-0935077-0), with a [modification](http://epubs.siam.org/doi/pdf/10.1137/S0036144596322507) to improve stability.
+The algorithm is due to [Fornberg, 1988](https://doi.org/10.1090/S0025-5718-1988-0935077-0), with a [modification](http://epubs.siam.org/doi/pdf/10.1137/S0036144596322507) to improve stability,
+& [Fornberg, 2020](https://doi.org/10.1093/imanum/draa006)
 """
 calculate_weights
 

test/Misc/utils.jl

@@ -9,6 +9,34 @@ using ModelingToolkit
     @test DiffEqOperators.perpsize(zeros(2,2,3,2), 3) == (2, 2, 2)
 end
 
+@testset "finite-difference weights from fornberg(1988) & fornberg(2020)" begin
+    order = 2; z = 0.0; x = [-1, 0, 1.0];
+    @test DiffEqOperators.calculate_weights(order, z, x) == [1,-2,1]  # central difference of second-derivative with unit-step
+
+    order = 1; z = 0.0; x = [-1., 1.0];
+    @test DiffEqOperators.calculate_weights(order, z, x) == [-0.5,0.5] # central difference of first-derivative with unit step
+
+    order = 1; z = 0.0; x = [0, 1];
+    @test DiffEqOperators.calculate_weights(order, z, x) == [-1, 1] # forward difference
+
+    order = 1; z = 1.0; x = [0, 1];
+    @test DiffEqOperators.calculate_weights(order, z, x) == [-1, 1] # backward difference
+
+    # forward-diff of third derivative with order of accuracy == 3
+    order = 3; z = 0.0; x = [0,1,2,3,4,5]
+    @test DiffEqOperators.calculate_weights(order, z, x) == [-17/4,	71/4	,−59/2,	49/2,	−41/4,	7/4]
+
+    order = 3; z = 0.0; x = collect(-3:3)
+    d, e = DiffEqOperators.calculate_weights(order, z, x;dfdx = true)
+    @test d ≈ [-167/18000, -963/2000, -171/16,0,171/16,963/2000,167/18000]
+    @test e ≈ [-1/600,-27/200,-27/8,-49/3,-27/8,-27/200,-1/600]
+
+    order = 3; z = 0.0; x = collect(-4:4)
+    d, e = DiffEqOperators.calculate_weights(order, z, x;dfdx = true)
+    @test d ≈ [-2493/5488000, -12944/385875, -87/125 ,-1392/125,0,1392/125,87/125,12944/385875,2493/5488000]
+    @test e ≈ [-3/39200,-32/3675,-6/25,-96/25,-205/12, -96/25, -6/25,-32/3675,-3/39200]
+end
+
 @testset "count differentials 1D" begin
     @parameters t x
     @variables u(..)