Holoborodko2008: Smooth noise-robust differentiators

src/SummationByPartsOperators.jl

@@ -77,6 +77,7 @@ export periodic_central_derivative_operator, periodic_derivative_operator, deriv
        fourier_derivative_operator, spectral_viscosity_operator, super_spectral_viscosity_operator,
        legendre_derivative_operator
 
+export Fornberg1998, Holoborodko2008, BeljaddLeFlochMishraParés2017
 export MattssonNordström2004, MattssonSvärdNordström2004, MattssonSvärdShoeybi2008,
        Mattsson2012, Mattsson2014,
        MattssonAlmquistCarpenter2014Extended, MattssonAlmquistCarpenter2014Optimal

src/periodic_operators.jl

@@ -436,6 +436,146 @@ function periodic_derivative_coefficients(derivative_order, accuracy_order, left
 end
 
 
+
+"""
+    Holoborodko2008
+
+Coefficients of the periodic operators given in
+  Holoborodko (2008)
+  Smooth Noise Robust Differentiators.
+  http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
+"""
+struct Holoborodko2008 <: SourceOfCoefficients end
+
+function Base.show(io::IO, ::Holoborodko2008)
+    print(io,
+        "  Holoborodko (2008) \n",
+        "  Smooth Noise Robust Differentiators. \n",
+        "  http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/ \n")
+end
+
+"""
+    periodic_derivative_coefficients(source::Holoborodko2008, derivative_order, accuracy_order;
+                                     T=Float64, parallel=Val{:serial}(), 
+                                     stencil_width=accuracy_order+3)
+
+Create the `PeriodicDerivativeCoefficients` approximating the `derivative_order`-th
+derivative with an order of accuracy `accuracy_order` and scalar type `T` given
+by Holoborodko2008.
+The evaluation of the derivative can be parallised using threads by chosing
+parallel=Val{:threads}())`.
+"""
+function periodic_derivative_coefficients(source::Holoborodko2008, derivative_order, accuracy_order;
+                                          T=Float64, parallel=Val{:serial}(), 
+                                          stencil_width=accuracy_order+3)
+    method_exists = true
+    if derivative_order == 1
+        if accuracy_order == 2
+            if stencil_width == 5
+                lower_coef = SVector{2, T}([
+                    -2//8, -1//8
+                ])
+                upper_coef = -lower_coef
+                central_coef = T(0)
+            elseif stencil_width == 7
+                lower_coef = SVector{3, T}([
+                    -5//32, -4//32, -1//32
+                ])
+                upper_coef = -lower_coef
+                central_coef = T(0)
+            elseif stencil_width == 9
+                lower_coef = SVector{4, T}([
+                    -14//128, -14//128, -6//128, -1//128
+                ])
+                upper_coef = -lower_coef
+                central_coef = T(0)
+            elseif stencil_width == 11
+                lower_coef = SVector{5, T}([
+                    -42//512, -48//512, -27//512, -8//512, -1//512
+                ])
+                upper_coef = -lower_coef
+                central_coef = T(0)
+            else
+                method_exists = false
+            end
+        elseif accuracy_order == 4
+            if stencil_width == 7
+                lower_coef = SVector{3, T}([
+                    -39//96, -12//96, 5//96
+                ])
+                upper_coef = -lower_coef
+                central_coef = T(0)
+            elseif stencil_width == 9
+                lower_coef = SVector{4, T}([
+                    -27//96, -16//96, 1//96, 2//96
+                ])
+                upper_coef = -lower_coef
+                central_coef = T(0)
+            elseif stencil_width == 11
+                lower_coef = SVector{5, T}([
+                    -322//1536, -256//1536, -39//1536, 32//1536, 11//1536
+                ])
+                upper_coef = -lower_coef
+                central_coef = T(0)
+            else
+                method_exists = false
+            end
+        else
+            method_exists = false
+        end
+    elseif derivative_order == 2
+        if accuracy_order == 2
+            if stencil_width == 5
+                lower_coef = SVector{2, T}([
+                    0, 1//4
+                ])
+                upper_coef = lower_coef
+                central_coef = T(-2//4)
+            elseif stencil_width == 7
+                lower_coef = SVector{3, T}([
+                    -1//16, 2//16, 1//16
+                ])
+                upper_coef = lower_coef
+                central_coef = T(-4//16)
+            elseif stencil_width == 9
+                lower_coef = SVector{4, T}([
+                    -4//64, 4//64, 4//64, 1//64
+                ])
+                upper_coef = lower_coef
+                central_coef = T(-10//64)
+            else
+                method_exists = false
+            end
+        elseif accuracy_order == 4
+            if stencil_width == 7
+                lower_coef = SVector{3, T}([
+                    1//12, 5//12, -1//12
+                ])
+                upper_coef = lower_coef
+                central_coef = T(-10//12)
+            elseif stencil_width == 9
+                lower_coef = SVector{4, T}([
+                    -12//192, 52//192, 12//192, -7//192
+                ])
+                upper_coef = lower_coef
+                central_coef = T(-90//192)
+            else
+                method_exists = false
+            end
+        else
+            method_exists = false
+        end
+    else
+        method_exists = false
+    end
+    if method_exists == false
+        throw(ArgumentError("Method with derivative_order=$derivative_order, accuracy_order=$accuracy_order, stencil_width=$stencil_width not implemented/derived."))
+    end
+
+    PeriodicDerivativeCoefficients(lower_coef, central_coef, upper_coef, parallel, derivative_order, accuracy_order, source)
+end
+
+
 """
     PeriodicDerivativeOperator
 
@@ -575,6 +715,25 @@ end
     periodic_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N, left_offset, parallel)
 end
 
+"""
+    periodic_derivative_operator(source::Holoborodko2008, derivative_order, accuracy_order, 
+                                 xmin, xmax, N; parallel=Val{:serial}(), kwargs...)
+
+Create a `PeriodicDerivativeOperator` approximating the `derivative_order`-th
+derivative on a uniform grid between `xmin` and `xmax` with `N` grid points up
+to order of accuracy `accuracy_order` where the leftmost grid point used is
+determined by `left_offset`.
+The evaluation of the derivative can be parallised using threads by chosing
+`parallel=Val{:threads}())`.
+"""
+function periodic_derivative_operator(source::Holoborodko2008, derivative_order, accuracy_order, 
+                                      xmin, xmax, N; parallel=Val{:serial}(), kwargs...)
+    grid = linspace(xmin, xmax, N) # N includes two identical boundary nodes
+    coefficients = periodic_derivative_coefficients(source, derivative_order, accuracy_order;
+                                                    kwargs..., T=eltype(grid), parallel=parallel)
+    PeriodicDerivativeOperator(coefficients, grid)
+end
+
 """
     periodic_derivative_operator(derivative_order, accuracy_order, grid, left_offset=-(accuracy_order+1)÷2, parallel=Val{:serial}())
 

test/periodic_operators_test.jl

@@ -654,3 +654,137 @@ let T = Float64
     xplot, duplot = evaluate_coefficients(res, D)
     @test maximum(abs, duplot-dufunc.(xplot)) < tol
 end
+
+
+# Robust methods of Holoborodko
+let T = Float32
+    xmin = one(T)
+    xmax = 2*one(T)
+    N = 100
+    x1 = compute_coefficients(identity, periodic_derivative_operator(1, 2, xmin, xmax, N))
+
+    x0 = ones(x1)
+    x2 = x1 .* x1
+    x3 = x2 .* x1
+    x4 = x2 .* x2
+    x5 = x2 .* x3
+    x6 = x3 .* x3
+    x7 = x4 .* x3
+
+    res = zeros(x0)
+
+    # first derivative operators
+    der_order = 1
+    acc_order = 2
+    for stencil_width in 5:2:11
+        D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=stencil_width)
+        println(DevNull, D)
+        @test derivative_order(D) == der_order
+        @test accuracy_order(D) == acc_order
+        @test issymmetric(D) == false
+        A_mul_B!(res, D, x0)
+        @test all(i->abs(res[i]) < eps(T), stencil_width:length(res)-stencil_width)
+        A_mul_B!(res, D, x1)
+        @test all(i->res[i] ≈ x0[i], stencil_width:length(res)-stencil_width)
+        A_mul_B!(res, D, x2)
+        @test all(i->res[i] ≈ 2*x1[i], stencil_width:length(res)-stencil_width)
+    end
+    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=3)
+    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=13)
+
+    acc_order = 4
+    for stencil_width in 7:2:11
+        D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=stencil_width)
+        println(DevNull, D)
+        @test derivative_order(D) == der_order
+        @test accuracy_order(D) == acc_order
+        @test issymmetric(D) == false
+        A_mul_B!(res, D, x0)
+        @test all(i->abs(res[i]) < 50*eps(T),  stencil_width:length(res)-stencil_width)
+        A_mul_B!(res, D, x1)
+        @test all(i->res[i] ≈ x0[i],  stencil_width:length(res)-stencil_width)
+        A_mul_B!(res, D, x2)
+        @test all(i->res[i] ≈ 2*x1[i],  stencil_width:length(res)-stencil_width)
+        A_mul_B!(res, D, x3)
+        @test all(i->res[i] ≈ 3*x2[i],  stencil_width:length(res)-stencil_width)
+        A_mul_B!(res, D, x4)
+        @test all(i->res[i] ≈ 4*x3[i],  stencil_width:length(res)-stencil_width)
+    end
+    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=5)
+    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=13)
+
+    acc_order = 6
+    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=9)
+    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=13)
+
+    # second derivative operators
+    der_order = 2
+    acc_order = 2
+    for stencil_width in 5:2:9
+        D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=stencil_width)
+        println(DevNull, D)
+        @test derivative_order(D) == der_order
+        @test accuracy_order(D) == acc_order
+        @test issymmetric(D) == true
+        A_mul_B!(res, D, x0)
+        @test all(i->abs(res[i]) < eps(T),  stencil_width:length(res)-stencil_width)
+        A_mul_B!(res, D, x1)
+        @test all(i->abs(res[i]) < 100N*eps(T),  stencil_width:length(res)-stencil_width)
+        A_mul_B!(res, D, x2)
+        @test all(i->isapprox(res[i], 2*x0[i], atol=600N*eps(T)),  stencil_width:length(res)-stencil_width)
+        A_mul_B!(res, D, x3)
+        @test all(i->isapprox(res[i], 6*x1[i], atol=2000N*eps(T)),  stencil_width:length(res)-stencil_width)
+    end
+    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=3)
+    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=11)
+
+    acc_order = 4
+    for stencil_width in 7:2:9
+        D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=stencil_width)
+        println(DevNull, D)
+        @test derivative_order(D) == der_order
+        @test accuracy_order(D) == acc_order
+        @test issymmetric(D) == true
+        A_mul_B!(res, D, x0)
+        @test all(i->abs(res[i]) < 50N*eps(T),  stencil_width:length(res)-stencil_width)
+        A_mul_B!(res, D, x1)
+        @test all(i->abs(res[i]) < 1000N*eps(T),  stencil_width:length(res)-stencil_width)
+        A_mul_B!(res, D, x2)
+        @test all(i->isapprox(res[i], 2*x0[i], atol=5000N*eps(T)),  stencil_width:length(res)-stencil_width)
+        A_mul_B!(res, D, x3)
+        @test all(i->isapprox(res[i], 6*x1[i], atol=5000N*eps(T)),  stencil_width:length(res)-stencil_width)
+        A_mul_B!(res, D, x4)
+        @test all(i->isapprox(res[i], 12*x2[i], atol=7000N*eps(T)),  stencil_width:length(res)-stencil_width)
+        A_mul_B!(res, D, x5)
+        @test any(i->!(res[i] ≈ 30*x3[i]),  stencil_width:length(res)-stencil_width)
+    end
+    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=5)
+    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=11)
+
+    acc_order = 6
+    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=9)
+    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=13)
+
+    der_order = 3
+    acc_order = 2
+    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=5)
+    @test_throws ArgumentError D = periodic_derivative_operator(Holoborodko2008(), der_order, acc_order,
+                                         xmin, xmax, N, stencil_width=13)
+end