[FIXED] scipy.cumtrapz and np.trapz give different answers

Issue

For some reason using numpy's trapz and scipy's cumtrapz yields different solutions.

x = np.linspace(-2, 4, num=20)
y = (x)
y_int = integrate.cumtrapz(y, x, initial=0)
y_tr = np.trapz(y, x, axis = 0) 
display(y_tr)
display(y_int)

The final value given by trapz rule is twice that given by cumtrapz.

-5.551115123125783e-17
array([ 0.00000000e+00, -3.98891967e-01, -7.53462604e-01, -1.06371191e+00,
       -1.32963989e+00, -1.55124654e+00, -1.72853186e+00, -1.86149584e+00,
       -1.95013850e+00, -1.99445983e+00, -1.99445983e+00, -1.95013850e+00,
       -1.86149584e+00, -1.72853186e+00, -1.55124654e+00, -1.32963989e+00,
       -1.06371191e+00, -7.53462604e-01, -3.98891967e-01, -2.77555756e-16])

Is there a reason for this?

Solution

It seems the reason is the way memory is allocated and floating operations are carried out.

Both -2.77555756e-16 and -5.551115123125783e-17 are essentially 0. Running the same code by changing the limit from -2 to 4 gives the correct answer:

x = np.linspace(-2, 4, num=20)
y = (x)
y_int = integrate.cumtrapz(y, x, initial=0)
y_tr = np.trapz(y, x, axis = 0) 
display(y_tr)
display(y_int)

Giving:

6.000000000000001

array([ 0.        , -0.58171745, -1.06371191, -1.44598338, -1.72853186,
       -1.91135734, -1.99445983, -1.97783934, -1.86149584, -1.64542936,
       -1.32963989, -0.91412742, -0.39889197,  0.21606648,  0.93074792,
        1.74515235,  2.65927978,  3.67313019,  4.7867036 ,  6.        ])

The result is 6 in each case.

Answered By - R Walser