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

src/discretization/fornberg.jl

@@ -1,20 +1,23 @@
-#############################################################
-# 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.
-
-function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real
-    #=
+#= Fornberg algorithm
+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.
+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)
+    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.                                             
+=#
 
-        The approximation order of the stencil is automatically determined from
-        the number of requested stencil points.
-    =#
+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
@@ -58,5 +61,23 @@ function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real
     if order != 0
         _C[div(N,2)+1] -= sum(_C)
     end
-    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

test/pde_systems/MOLfornberg_weights.jl

@@ -0,0 +1,30 @@
+using Test, LinearAlgebra
+using MethodOfLines
+
+@testset "finite-difference weights from fornberg(1988) & fornberg(2020)" begin
+    order = 2; z = 0.0; x = [-1, 0, 1.0];
+    @test MethodOfLines.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 MethodOfLines.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 MethodOfLines.calculate_weights(order, z, x) == [-1, 1] # forward difference
+
+    order = 1; z = 1.0; x = [0, 1];
+    @test MethodOfLines.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 MethodOfLines.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 = MethodOfLines.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 = MethodOfLines.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

test/runtests.jl

@@ -17,6 +17,7 @@ const is_TRAVIS = haskey(ENV,"TRAVIS")
 
     if GROUP == "All" || GROUP == "Integration Tests"
         @time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLtest.jl") end
+        @time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLfornberg_weights.jl") end
         #@time @safetestset "MOLFiniteDifference Interface: Linear Convection" begin include("pde_systems/MOL_1D_Linear_Convection.jl") end
         @time @safetestset "MOLFiniteDifference Interface: 1D Linear Diffusion" begin include("pde_systems/MOL_1D_Linear_Diffusion.jl") end
         @time @safetestset "MOLFiniteDifference Interface: 1D Non-Linear Diffusion" begin include("pde_systems/MOL_1D_NonLinear_Diffusion.jl") end