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']}
48 The returned :class:`Data` block shares its x vector with the
49 input data. The ith df/dx value in the returned data is
52 (df/dx)[i] = (SUM_j w[j] f[i+j]) / h
54 where ``h = x[i+1]-x[i]`` is the x coordinate spacing (assumed
55 constant) and ``j`` ranges over the keys of `weights`.
57 There standard schemes translate as follows:
59 ======== ====================== ===================
60 scheme formula weights
61 ======== ====================== ===================
62 forward ``(f[i+1]-f[i])/h`` ``{0:-1,1:1}``
63 backward ``(f[i]-f[i-1])/h`` ``{0:1,-1:-1}``
64 central ``(f[i+1]-f[i-1])/2h`` ``{-1:-0.5,1:0.5}``
65 ======== ====================== ===================
67 The default scheme is central differencing.
71 The boundary conditions could be configurable in principle. The
72 current scheme just extrapolates virtual points out to negative
75 f[i<0] = 2*f[0] - f[-i]
77 With analogous treatment for `i > data.shape[0]`. This ensures that
78 `f[i]-f[0]` is odd about `i=0`, which keeps derivatives smooth.::
80 f[i] - f[0] = f[0] - f[-i] == -(f[-i] - f[0])
82 output = Data((data.shape[0],2), dtype=data.dtype)
83 output.info = copy.copy(data.info)
84 output.info['columns'] = [
85 data.info['columns'][x_col],
86 'deriv %s with respect to %s' \
87 % (data.info['columns'][f_col], data.info['columns'][x_col]),
89 h = data[1,x_col] - data[0,x_col]
91 for i,w in weights.items():
92 chunk = numpy.roll(w*data[:,f_col], -i)
93 if i > 0: # chunk shifted down, replace the high `i`s
94 zero = len(chunk) - 1 - i
95 for j in range(1,i+1):
96 chunk[zero+j] = 2*chunk[zero] - chunk[zero-j]
97 elif i < 0: # chunk shifted up, replace the low `i`s
99 for j in range(1,zero+1):
100 chunk[zero-j] = 2*chunk[zero] - chunk[zero+j]
103 output[:,0] = data[:,x_col]
104 output[:,1] = sum(chunks)