Scipy Simpson integration of a positive quantity comes out negative

Question

The following numerical integration of y values that are positive evaluates to a negative integration value using the builtin provided Simpson integration method of Scipy. This is much to my surprise. I'm a beginner on this. Any help would be appreciated.

>>> import numpy as np 
>>> from scipy.integrate import simpson 
>>> LOS = np.array([1.00000000e-06, 8.40496052e-06, 7.06433614e-05, 5.93754663e-04, 4.99048451e-03, 4.19448253e-02, 3.52544600e-01, 2.96312345e+00, 2.49049356e+01, 2.09325001e+02]) 
>>> ne  = np.array([1.89989649e-04, 1.89989583e-04, 1.89989021e-04, 1.89984303e-04, 1.89944615e-04, 1.89609139e-04, 1.86662296e-04, 1.56895875e-04, 4.13547470e-05, 3.17491868e-06]) 
>>> simpson(ne, LOS) 
-0.006258317426859777
>>> np.trapz(ne, LOS) 
0.006795909558151828

I expected a positive integral value with Simpson but what I get is negative.

Answer 1

Your data are a little strange for integration, because overwhelmingly the last few data points will have most of the influence over the integral and everything else will have no effect, since you have log scale.

This means that the Simpson integral is very sensitive to the choice of even; quoting the doc:

If there are an even number of samples, N, then there are an odd number of intervals (N-1), but Simpson’s rule requires an even number of intervals. The parameter ‘even’ controls how this is handled. [...]

‘first’: Use Simpson’s rule for the first N-2 intervals with a trapezoidal rule on the last interval.

import numpy as np
from scipy.integrate import simpson

LOS = np.array([
    1.00000000e-06, 8.40496052e-06, 7.06433614e-05, 5.93754663e-04, 4.99048451e-03, 4.19448253e-02, 3.52544600e-01,
    2.96312345e+00, 2.49049356e+01, 2.09325001e+02])
ne = np.array([
   1.89989649e-04, 1.89989583e-04, 1.89989021e-04, 1.89984303e-04, 1.89944615e-04, 1.89609139e-04, 1.86662296e-04,
   1.56895875e-04, 4.13547470e-05, 3.17491868e-06])

order = np.argsort(LOS)
print(simpson(y=ne, x=LOS, even='first'))
print(np.trapz(y=ne, x=LOS))

0.006355138100371379
0.006795909558151828