3 """The `calculus` module provides functions for calculating
4 derivatives and integrals of discrete data.
11 from ..curve import Data
14 def derivative(data, x_col=0, f_col=1, weights={-1:-0.5, 1:0.5}):
15 """Calculate the discrete derivative (finite difference) of
16 data[:,f_col] with respect to data[:,x_col].
22 >>> d = Data((5,2), dtype=numpy.float,
23 ... info={'columns':['x', 'x**2']})
24 >>> for i in range(5):
33 >>> dd = derivative(d)
40 >>> pprint.pprint(dd.info)
41 {'columns': ['x', 'deriv x**2 with respect to x']}
49 The returned :class:`Data` block shares its x vector with the
50 input data. The ith df/dx value in the returned data is
53 (df/dx)[i] = (SUM_j w[j] f[i+j]) / h
55 where ``h = x[i+1]-x[i]`` is the x coordinate spacing (assumed
56 constant) and ``j`` ranges over the keys of `weights`.
58 There standard schemes translate as follows:
60 ======== ====================== ===================
61 scheme formula weights
62 ======== ====================== ===================
63 forward ``(f[i+1]-f[i])/h`` ``{0:-1,1:1}``
64 backward ``(f[i]-f[i-1])/h`` ``{0:1,-1:-1}``
65 central ``(f[i+1]-f[i-1])/2h`` ``{-1:-0.5,1:0.5}``
66 ======== ====================== ===================
68 The default scheme is central differencing.
73 These could be configurable in principle. The current scheme just
74 extrapolates virtual points out to negative `i` following::
76 f[i<0] = 2*f[0] - f[-i]
78 With analogous treatment for `i > data.shape[0]`. This ensures that
79 `f[i]-f[0]` is odd about `i=0`, which keeps derivatives smooth.::
81 f[i] - f[0] = f[0] - f[-i] == -(f[-i] - f[0])
83 output = Data((data.shape[0],2), dtype=data.dtype)
84 output.info = copy.copy(data.info)
85 output.info['columns'] = [
86 data.info['columns'][x_col],
87 'deriv %s with respect to %s' \
88 % (data.info['columns'][f_col], data.info['columns'][x_col]),
90 h = data[1,x_col] - data[0,x_col]
92 for i,w in weights.items():
93 chunk = numpy.roll(w*data[:,f_col], -i)
94 if i > 0: # chunk shifted down, replace the high `i`s
95 zero = len(chunk) - 1 - i
96 for j in range(1,i+1):
97 chunk[zero+j] = 2*chunk[zero] - chunk[zero-j]
98 elif i < 0: # chunk shifted up, replace the low `i`s
100 for j in range(1,zero+1):
101 chunk[zero-j] = 2*chunk[zero] - chunk[zero+j]
104 output[:,0] = data[:,x_col]
105 output[:,1] = sum(chunks)