From 332d628744a0670234585053dbe32a3e82e0c4db Mon Sep 17 00:00:00 2001 From: danieljfarrell Date: Mon, 12 Aug 2013 14:33:11 +0900 Subject: ENH: Improve accuracy of numpy.gradient at edges * numpy.gradient has been enhanced to use a second order accurate one-sided finite difference stencil at boundary elements of the array. Second order accurate central difference are still used for the interior elements. The result is a fully second order accurate approximation of the gradient over the full domain. * The one-sided stencil uses 3 elements each with a different weight. A forward difference is used for the first element, dy/dx ~ -(3.0*y[0] - 4.0*y[1] + y[2]) / (2.0*dx) and backwards difference is used for the last element, dy/dx ~ (3.0*y[-1] - 4.0*y[-2] + y[-3]) / (2.0*dx) * Because the datetime64 datatype cannot be multiplied a view is taken of datetime64 arrays and cast to int64. The gradient algorithm is then applied to the view rather than the input array. * Previously no dimension checks were performed on the input array. Now if the array size along the differentiation axis is less than 2, a ValueError is raised which explains that more elements are needed. If the size is exactly two the function falls back to using a 2 point stencil (the old behaviour). If the size is 3 and above then the higher accuracy methods are used. * A new test has been added which validates the higher accuracy. Old tests have been updated to pass. Note, this should be expected because the boundary elements now return different (more accurate) values. --- numpy/lib/tests/test_function_base.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) (limited to 'numpy/lib/tests/test_function_base.py') diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py index 61113ae7a..f52eb5fbe 100644 --- a/numpy/lib/tests/test_function_base.py +++ b/numpy/lib/tests/test_function_base.py @@ -466,7 +466,7 @@ class TestGradient(TestCase): '1910-10-12', '1910-12-12', '1912-12-12'], dtype='datetime64[D]') dx = np.array( - [-5, -3, 0, 31, 61, 396, 731], + [-7, -3, 0, 31, 61, 396, 1066], dtype='timedelta64[D]') assert_array_equal(gradient(x), dx) assert_(dx.dtype == np.dtype('timedelta64[D]')) @@ -477,11 +477,23 @@ class TestGradient(TestCase): [-5, -3, 10, 12, 61, 321, 300], dtype='timedelta64[D]') dx = np.array( - [2, 7, 7, 25, 154, 119, -21], + [-3, 7, 7, 25, 154, 119, -161], dtype='timedelta64[D]') assert_array_equal(gradient(x), dx) assert_(dx.dtype == np.dtype('timedelta64[D]')) + def test_second_order_accurate(self): + # Testing that the relative numerical error is less that 3% for + # this example problem. This corresponds to second order + # accurate finite differences for all interior and boundary + # points. + x = np.linspace(0, 1, 10) + dx = x[1] - x[0] + y = 2 * x ** 3 + 4 * x ** 2 + 2 * x + analytical = 6 * x ** 2 + 8 * x + 2 + num_error = np.abs((np.gradient(y, dx) / analytical) - 1) + assert_(np.all(num_error < 0.03) == True) + class TestAngle(TestCase): def test_basic(self): -- cgit v1.2.1