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.

Your `LOS` is very much a logarithmic scale. Are all of your points of equal importance? — Reinderien, Dec 28 '22 at 19:37
Somewhat related: https://stackoverflow.com/questions/36803745/python-simpsons-rule-negative-answer-for-positive-area-under-the-curve — Reinderien, Dec 28 '22 at 20:26
Thanks for the inputs. Looks like LOS needs to be either sampled linearly or in log space with a higher resolution for Simpson to properly function. — curious, Dec 30 '22 at 03:38

Reinderien · Answer 1 · 2022-12-28T20:31:56.167

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

Scipy Simpson integration of a positive quantity comes out negative

1 Answers1