Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

scipy -- how to integrate a linearly interpolated function?

I have a function which is an interpolation of a relative large set of data. I use linear interpolation interp1d so there are a lot of non-smooth sharp point like this. The quad function from scipy will give warning because of the sharp points. I wonder how to do the integration without the warning?

Thank you!


Thanks for all the answers. Here I summarize the solutions in case some others run into the same problem:

  1. Just like what @Stelios did, use points to avoid warnings and to get a more accurate result.
  2. In practice the number of points are usually larger than the default limit(limit=50) of quad, so I choose quad(f_interp, a, b, limit=2*p.shape[0], points=p) to avoid all those warnings.
  3. If a and b are not the same start or the end point of the data set x, the points p can be chosen by p = x[where(x>=a and x<=b)]
like image 831
Shu Liao Avatar asked Jun 28 '17 20:06

Shu Liao


1 Answers

quad accepts an optional argument, called points. According to the documentation:

points : (sequence of floats,ints), optional

A sequence of break points in the bounded integration interval where local difficulties of the integrand may occur (e.g., singularities, discontinuities). The sequence does not have to be sorted.

In your case, the "difficult" points are exactly the x-coordinates of the data points. Here is an example:

import numpy as np 
from scipy.integrate import quad
np.random.seed(123)

# generate random data set 
x = np.arange(0,10)  
y = np.random.rand(10)

# construct a linear interpolation function of the data set 
f_interp = lambda xx: np.interp(xx, x, y)

Here is a plot of the data points and f_interp: enter image description here

Now calling quad as

quad(f_interp,0,9)

return a series of warnings along with

(4.89770017785734, 1.3762838395159349e-05)

If you provide the points argument, i.e.,

quad(f_interp,0,9, points = x)

it issues no warnings and the result is

(4.8977001778573435, 5.437539505167948e-14)

which also implies a much greater accuracy of the result compared to the previous call.

like image 142
Stelios Avatar answered Sep 24 '22 20:09

Stelios