From bd857c3ced129c69f6696d9b3a588ff0ed9be7ef Mon Sep 17 00:00:00 2001
From: Zander <alex.lemony@gmail.com>
Date: Wed, 17 Jul 2019 13:06:09 +0200
Subject: [PATCH 1/2] added error message

---
 src/derivative_operators/derivative_operator.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl
index 36ed4c6f4..34cd3f470 100644
--- a/src/derivative_operators/derivative_operator.jl
+++ b/src/derivative_operators/derivative_operator.jl
@@ -20,7 +20,7 @@ struct CenteredDifference{N} end
 function CenteredDifference{N}(derivative_order::Int,
                             approximation_order::Int, dx::T,
                             len::Int, coeff_func=nothing) where {T<:Real,N}
-
+    @assert approximation_order>1 "approximation_order must be greater than 1."
     stencil_length          = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2
     boundary_stencil_length = derivative_order + approximation_order
     dummy_x                 = -div(stencil_length,2) : div(stencil_length,2)
@@ -30,8 +30,8 @@ function CenteredDifference{N}(derivative_order::Int,
     deriv_spots             = (-div(stencil_length,2)+1) : -1
     boundary_deriv_spots    = boundary_x[2:div(stencil_length,2)]
 
-    stencil_coefs           = convert(SVector{stencil_length, T}, calculate_weights(derivative_order, zero(T), dummy_x))
-    _low_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, calculate_weights(derivative_order, oneunit(T)*x0, boundary_x)) for x0 in boundary_deriv_spots]
+    stencil_coefs           = convert(SVector{stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, zero(T), dummy_x))
+    _low_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, boundary_x)) for x0 in boundary_deriv_spots]
     low_boundary_coefs      = convert(SVector{boundary_point_count},_low_boundary_coefs)
     high_boundary_coefs     = convert(SVector{boundary_point_count},reverse(SVector{boundary_stencil_length, T}[reverse(low_boundary_coefs[i]) for i in 1:boundary_point_count]))
 

From 71a7cd448be2a9689a439af216914401cd11cc2a Mon Sep 17 00:00:00 2001
From: Zander <alex.lemony@gmail.com>
Date: Wed, 17 Jul 2019 13:06:47 +0200
Subject: [PATCH 2/2] added error benchmark

---
 test/error_benchmarking.jl | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)
 create mode 100644 test/error_benchmarking.jl

diff --git a/test/error_benchmarking.jl b/test/error_benchmarking.jl
new file mode 100644
index 000000000..5753b537a
--- /dev/null
+++ b/test/error_benchmarking.jl
@@ -0,0 +1,28 @@
+using DiffEqOperators, Plots
+
+n = 10000
+@show n
+x = range(0.0; length = n, stop = 2π)
+dx = x[2]-x[1]
+y = exp.(π*x)
+y_ = y[2:(end-1)]
+
+for dor in 1:4, aor in 2:6
+
+    D1 = CenteredDifference(dor,aor,dx,length(x))
+
+    #test result
+    @show dor
+    @show aor
+    #take derivatives
+    err_abs = abs.(D1*y .- (π^dor)*y_)
+    err_percent = 100*err_abs./abs.(((π^dor)*y_))
+    max_error = maximum(err_percent) # test operator with known derivative of exp(kx)
+    avg_error = sum(err_percent)/length(err_percent)
+
+    @show max_error
+    @show avg_error
+    plot(x[2:(end-1)], err_percent, title="Percentage Error, n=$n aor = $aor, dor = $dor")
+
+    #test result
+end