CODEWITHSUNDEEP
pythonnumpyscipy
Shu Liao
Shu Liao

Reputation: 145

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)]

Upvotes: 7

Views: 7640

Answers (2)

Warren Weckesser
Warren Weckesser

Reputation: 114911

Instead of interp1d, you could use scipy.interpolate.InterpolatedUnivariateSpline. That interpolator has the method integral(a, b) that computes the definite integral.

Here's an example:

import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline
import matplotlib.pyplot as plt


# Create some test data.
x = np.linspace(0, np.pi, 21)
np.random.seed(12345)
y = np.sin(1.5*x) + np.random.laplace(scale=0.35, size=len(x))**3

# Create the interpolator.  Use k=1 for linear interpolation.
finterp = InterpolatedUnivariateSpline(x, y, k=1)

# Create a finer mesh of points on which to compute the integral.
xx = np.linspace(x[0], x[-1], 5*len(x))

# Use the interpolator to compute the integral from 0 to t for each
# t in xx.
qq = [finterp.integral(0, t) for t in xx]

# Plot stuff
p = plt.plot(x, y, '.', label='data')
plt.plot(x, y, '-', color=p[0].get_color(), label='linear interpolation')
plt.plot(xx, qq, label='integral of linear interpolation')
plt.grid()
plt.legend(framealpha=1, shadow=True)
plt.show()

The plot:

plot

Upvotes: 5

Stelios
Stelios

Reputation: 5521

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.

Upvotes: 5

Related Questions