diff --git a/src/SummationByPartsOperators.jl b/src/SummationByPartsOperators.jl index d8fb4f7e..aa63b0ec 100644 --- a/src/SummationByPartsOperators.jl +++ b/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 diff --git a/src/periodic_operators.jl b/src/periodic_operators.jl index ee85baf6..e8de3c21 100644 --- a/src/periodic_operators.jl +++ b/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}()) diff --git a/test/periodic_operators_test.jl b/test/periodic_operators_test.jl index ddb74e47..2843b126 100644 --- a/test/periodic_operators_test.jl +++ b/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