[SciPy-Dev] Code review for trapz update

[fizyxnrd]

Alan G Isaac <alan.isaac <at> gmail.com> writes:

> 
> On 3/10/2015 12:48 AM, K. N. wrote:
> > I've proposed a change to trapz, so that bin widths (`x` array 
deltas) are always considered positive. 
> This allows monotonically decreasing `x`
> > sequences to be used, in addition to monotonically increasing 
sequences.
> 
> Can you elaborate please.  Monontically decreasing sequences can 
already be used.
> Aren't you losing the sign of the area?
> 
> Alan Isaac
> 


Trapezoid rule is an approximation of $\int y(x) dx$, replacing integral 
with sum by $\sum_{i=1}^{N-1} (y_{i+1} + y{i}) * (x_{i+1} - x{i}) / 2$.  
Written in this form, it might seem that we want to keep track of the 
sign of $x_{i+1} - x{i}$.  However, this construction always assumes 
that the $x_i$ are an ordered sequence.  By interpreting this simply, we 
find that we may have negative *intervals*.  To be more precise, our 
summation should read $\sum_{i=1}^{N-1} (y_{i+1} + y{i}) * (\Delta x_i) 
/ 2$, where $\Delta x_i = |x_{i+1} - x{i}|$.  This is certainly the 
intent of the usual mathematical construction.  Note that negative 
intervals are not necessary to achieve negative "areas".

This data should integrate to negative "area":
x = [1,  2,  3,  4,  5]
y = [0, -1, -2, -1,  0]
(-4)

This data should integrate to positive "area", even though the x data is 
decreasing (and we have negative intervals):
x = [5, 4, 3, 2, 1]
y = [2, 2, 2, 2, 2]
(currently, get -8, but should get +8).

It might seem reasonable to require a user to always make their x 
sequence monotonically increasing.  However, trapz is silent if this is 
not the case, and at the very least should be modified throw an error 
for non monotonically increasing x sequences.

This issue arises, for example, in the conversion of spectral quantities 
from wavenumber (units of inverse length) to wavelength (units of 
length).  When the spectral quantity is inverted, a monotonically 
increasing sequence becomes a monotonically decreasing sequence, and 
vice versa.  This gives the unexpected result that trapz(<quantity>, 
<spectrum>) works as expected before conversion, but returns values with 
the wrong sign after conversion.

This behavior is even more difficult to recognize if the <quantity> (y-
value) term is not strictly positive, and a user may utilize trapz 
incorrectly with no indication that it is not behaving as expected.