1 # Copyright (C) 2010 W. Trevor King <wking@drexel.edu>
3 # This file is part of Hooke.
5 # Hooke is free software: you can redistribute it and/or modify it
6 # under the terms of the GNU Lesser General Public License as
7 # published by the Free Software Foundation, either version 3 of the
8 # License, or (at your option) any later version.
10 # Hooke is distributed in the hope that it will be useful, but WITHOUT
11 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
13 # Public License for more details.
15 # You should have received a copy of the GNU Lesser General Public
16 # License along with Hooke. If not, see
17 # <http://www.gnu.org/licenses/>.
19 """The `calculus` module provides functions for calculating
20 derivatives and integrals of discrete data.
27 from ..curve import Data
30 def derivative(x_data, f_data, weights={-1:-0.5, 1:0.5}):
31 """Calculate the discrete derivative (finite difference) of
32 f_data with respect to x_data.
37 >>> d = Data((5,2), dtype=numpy.float,
38 ... info={'columns':['x', 'x**2']})
39 >>> for i in range(5):
48 >>> dd = derivative(x_data=d[:,0], f_data=d[:,1])
50 Data([ 1., 2., 4., 6., 7.])
57 The returned :class:`Data` block shares its x vector with the
58 input data. The `i`\th df/dx value in the returned data is
61 (df/dx)[i] = (SUM_j w[j] f[i+j]) / h
63 where ``h = x[i+1]-x[i]`` is the x coordinate spacing (assumed
64 constant) and ``j`` ranges over the keys of `weights`.
66 There standard schemes translate as follows:
68 ======== ====================== ===================
69 scheme formula weights
70 ======== ====================== ===================
71 forward ``(f[i+1]-f[i])/h`` ``{0:-1,1:1}``
72 backward ``(f[i]-f[i-1])/h`` ``{0:1,-1:-1}``
73 central ``(f[i+1]-f[i-1])/2h`` ``{-1:-0.5,1:0.5}``
74 ======== ====================== ===================
76 The default scheme is central differencing.
80 The boundary conditions could be configurable in principle. The
81 current scheme just extrapolates virtual points out to negative
84 f[i<0] = 2*f[0] - f[-i]
86 With analogous treatment for `i > data.shape[0]`. This ensures that
87 `f[i]-f[0]` is odd about `i=0`, which keeps derivatives smooth.::
89 f[i] - f[0] = f[0] - f[-i] == -(f[-i] - f[0])
91 output = Data(f_data.shape, dtype=f_data.dtype)
92 h = x_data[1] - x_data[0]
94 for i,w in weights.items():
95 chunk = numpy.roll(w*f_data, -i)
96 if i > 0: # chunk shifted down, replace the high `i`s
97 zero = len(chunk) - 1 - i
98 for j in range(1,i+1):
99 chunk[zero+j] = 2*chunk[zero] - chunk[zero-j]
100 elif i < 0: # chunk shifted up, replace the low `i`s
102 for j in range(1,zero+1):
103 chunk[zero-j] = 2*chunk[zero] - chunk[zero+j]