diff --git a/src/derivative_operators/fornberg.jl b/src/derivative_operators/fornberg.jl index 0b24ae706..eba4201dc 100644 --- a/src/derivative_operators/fornberg.jl +++ b/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 diff --git a/src/docstrings.jl b/src/docstrings.jl index d2ef73efe..b0887e51d 100644 --- a/src/docstrings.jl +++ b/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 diff --git a/test/Misc/utils.jl b/test/Misc/utils.jl index f3dfac3f5..9f16c852b 100644 --- a/test/Misc/utils.jl +++ b/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(..)