made first order actually first order and added fourth order

pynumdiff/finite_difference/__init__.py

@@ -1,5 +1,5 @@
 """This module implements some common finite difference schemes
 """
-from ._finite_difference import first_order, second_order
+from ._finite_difference import first_order, second_order, fourth_order
 
-__all__ = ['first_order', 'second_order'] # So these get treated as direct members of the module by sphinx
+__all__ = ['first_order', 'second_order', 'fourth_order'] # So these get treated as direct members of the module by sphinx

pynumdiff/finite_difference/_finite_difference.py

@@ -1,3 +1,4 @@
+"""This is handy for this module https://web.media.mit.edu/~crtaylor/calculator.html"""
 import numpy as np
 from pynumdiff.utils import utility
 from warnings import warn
@@ -25,29 +26,10 @@ def first_order(x, dt, params=None, options={}, num_iterations=None):
         return _iterate_first_order(x, dt, num_iterations)
 
     dxdt_hat = np.diff(x) / dt # Calculate the finite difference
-    dxdt_hat = np.hstack((dxdt_hat[0], dxdt_hat, dxdt_hat[-1])) # Pad the data
-    dxdt_hat = np.mean((dxdt_hat[0:-1], dxdt_hat[1:]), axis=0) # Re-finite dxdt_hat using linear interpolation
+    dxdt_hat = np.hstack((dxdt_hat, dxdt_hat[-1])) # using stencil -1,0, you get expression for previous value
 
     return x, dxdt_hat
 
-
-def second_order(x, dt):
-    """Second-order centered difference method
-
-    :param np.array[float] x: data to differentiate
-    :param float dt: step size
-
-    :return: tuple[np.array, np.array] of\n
-             - **x_hat** -- estimated (smoothed) x
-             - **dxdt_hat** -- estimated derivative of x
-    """
-    dxdt_hat = (x[2:] - x[0:-2]) / (2 * dt)
-    first_dxdt_hat = (-3 * x[0] + 4 * x[1] - x[2]) / (2 * dt)
-    last_dxdt_hat = (3 * x[-1] - 4 * x[-2] + x[-3]) / (2 * dt)
-    dxdt_hat = np.hstack((first_dxdt_hat, dxdt_hat, last_dxdt_hat))
-    return x, dxdt_hat
-
-
 def _x_hat_using_finite_difference(x, dt):
     """Find a smoothed estimate of the true function by taking FD and then integrating with trapezoids
     """
@@ -56,7 +38,6 @@ def _x_hat_using_finite_difference(x, dt):
     x0 = utility.estimate_initial_condition(x, x_hat)
     return x_hat + x0
 
-
 def _iterate_first_order(x, dt, num_iterations):
     """Iterative first order centered finite difference.
 
@@ -79,3 +60,43 @@ def _iterate_first_order(x, dt, num_iterations):
     x_hat, dxdt_hat = first_order(x, dt)
 
     return x_hat, dxdt_hat
+
+
+def second_order(x, dt):
+    """Second-order centered difference method, with special endpoint formulas.
+
+    :param np.array[float] x: data to differentiate
+    :param float dt: step size
+
+    :return: tuple[np.array, np.array] of\n
+             - **x_hat** -- estimated (smoothed) x
+             - **dxdt_hat** -- estimated derivative of x
+    """
+    dxdt_hat = np.zeros(x.shape)
+    dxdt_hat[1:-1] = x[2:] - x[:-2]
+    dxdt_hat[0] = -3 * x[0] + 4 * x[1] - x[2]
+    dxdt_hat[-1] = 3 * x[-1] - 4 * x[-2] + x[-3]
+    dxdt_hat /= 2*dt
+
+    return x, dxdt_hat
+
+
+def fourth_order(x, dt):
+    """Fourth-order centered difference method, with special endpoint formulas.
+
+    :param np.array[float] x: data to differentiate
+    :param float dt: step size
+
+    :return: tuple[np.array, np.array] of\n
+             - **x_hat** -- estimated (smoothed) x
+             - **dxdt_hat** -- estimated derivative of x
+    """
+    dxdt_hat = np.zeros(x.shape)
+    dxdt_hat[2:-2] = (8*(x[3:-1] - x[1:-3]) - x[4:] + x[:-4])
+    dxdt_hat[0] = -25*x[0] + 48*x[1] - 36*x[2] + 16*x[3] - 3*x[4]
+    dxdt_hat[1] = -3*x[0] - 10*x[1] + 18*x[2] - 6*x[3] + x[4]
+    dxdt_hat[-2] = 3*x[-1] + 10*x[-2] - 18*x[-3] + 6*x[-4] - x[-5]
+    dxdt_hat[-1] = 25*x[-1] - 48*x[-2] + 36*x[-3] - 16*x[-4] + 3*x[-5]
+    dxdt_hat /= 12*dt
+
+    return x, dxdt_hat

pynumdiff/tests/test_diff_methods.py

@@ -3,12 +3,12 @@
 from pytest import mark
 from warnings import warn
 
+from ..finite_difference import first_order, second_order, fourth_order
 from ..linear_model import lineardiff, polydiff, savgoldiff, spectraldiff
 from ..total_variation_regularization import velocity, acceleration, jerk, iterative_velocity, smooth_acceleration, jerk_sliding
 from ..kalman_smooth import constant_velocity, constant_acceleration, constant_jerk
 from ..smooth_finite_difference import mediandiff, meandiff, gaussiandiff, friedrichsdiff, butterdiff, splinediff
-from ..finite_difference import first_order, second_order
-# Function aliases for testing cases where parameters change the behavior in a big way
+# Function alias for testing a case where parameters change the behavior in a big way
 def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs)
 
 dt = 0.1
@@ -33,6 +33,7 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs)
     (first_order, {}), # empty dictionary for the case of no parameters. no params -> no diff in new vs old
     (iterated_first_order, {'num_iterations':5}), (iterated_first_order, [5], {'iterate':True}),
     (second_order, {}),
+    (fourth_order, {}),
     (lineardiff, {'order':3, 'gamma':5, 'window_size':11, 'solver':'CLARABEL'}), (lineardiff, [3, 5, 11], {'solver':'CLARABEL'}),
     (polydiff, {'polynomial_order':2, 'window_size':3}), (polydiff, [2, 3]),
     (savgoldiff, {'polynomial_order':2, 'window_size':4, 'smoothing_win':4}), (savgoldiff, [2, 4, 4]),
@@ -60,23 +61,29 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs)
 # big ol' table by the method, then the test function, then the pair of quantities we're comparing.
 error_bounds = {
     first_order: [[(-25, -25), (-25, -25), (0, 0), (1, 1)],
-                  [(-25, -25), (-13, -14), (0, 0), (1, 1)],
-                  [(-25, -25), (0, 0), (0, 0), (1, 0)],
+                  [(-25, -25), (-13, -13), (0, 0), (1, 1)],
                   [(-25, -25), (0, 0), (0, 0), (1, 1)],
+                  [(-25, -25), (1, 0), (0, 0), (1, 1)],
                   [(-25, -25), (2, 2), (0, 0), (2, 2)],
                   [(-25, -25), (3, 3), (0, 0), (3, 3)]],
-    iterated_first_order: [[(-8, -9), (-11, -11), (0, -1), (0, 0)],
-                           [(-6, -6), (-6, -7), (0, -1), (0, 0)],
-                           [(-1, -1), (0, 0), (0, -1), (0, 0)],
-                           [(0, 0), (1, 0), (0, 0), (1, 0)],
-                           [(1, 1), (2, 2), (1, 1), (2, 2)],
-                           [(1, 1), (3, 3), (1, 1), (3, 3)]],
+    iterated_first_order: [[(-8, -9), (-25, -25), (0, 0), (1, 1)],
+                           [(-6, -6), (-6, -7), (0, 0), (1, 1)],
+                           [(1, 0), (1, 0), (1, 1), (1, 1)],
+                           [(1, 0), (1, 1), (1, 0), (1, 1)],
+                           [(2, 2), (3, 2), (2, 2), (3, 2)],
+                           [(2, 2), (3, 3), (2, 2), (3, 3)]],
     second_order: [[(-25, -25), (-25, -25), (0, 0), (1, 1)],
                    [(-25, -25), (-13, -13), (0, 0), (1, 1)],
                    [(-25, -25), (-13, -13), (0, 0), (1, 1)],
                    [(-25, -25), (0, -1), (0, 0), (1, 1)],
                    [(-25, -25), (1, 1), (0, 0), (1, 1)],
                    [(-25, -25), (3, 3), (0, 0), (3, 3)]],
+    fourth_order: [[(-25, -25), (-25, -25), (0, 0), (1, 1)],
+                   [(-25, -25), (-13, -13), (0, 0), (1, 1)],
+                   [(-25, -25), (-13, -13), (0, 0), (1, 1)],
+                   [(-25, -25), (-2, -2), (0, 0), (1, 1)],
+                   [(-25, -25), (1, 0), (0, 0), (1, 1)],
+                   [(-25, -25), (2, 2), (0, 0), (2, 2)]],
     lineardiff: [[(-6, -6), (-5, -6), (0, -1), (0, 0)],
                  [(0, 0), (1, 1), (0, 0), (1, 1)],
                  [(1, 0), (2, 2), (1, 0), (2, 2)],