Interpolation and Integration in Python

by Jessica Lu on June 14, 2010

Interpolation

Interpolation: Filter profiles may be reported as transmission vs. wavelength data points. To resample the filter profiles to a different (regularly spaced) wavelength grid, I used scipy's interpolation package.

Interpolation is a used for many astronomical applications. Interpolation is required to combine sub-pixel dithered images or spectroscopy, sample grids of stellar evolution or stellar atmosphere models, calculate extinction from observed extinction curves, and many many more applications. The scipy.interpolate package in python has some nice built-in interpolation functions and I have gathered a few links describing the capabilities (in addition to the documentation). I would recommend using splrep/splev over interp1d for speed.

  • Scipy cookbook on interpolation
  • slideshow outlining the interpolation capabilities in Scipy
  • A thread explaining terms in spline interpolation/smoothing

Integration

Within Scipy, there is an integrate package with several different functions that perform definite or indefinite integrals. I have only played with these briefly; however, I tested both the scipy.integrate.quad() and scipy.integrate.romberg() functions. These all work in roughly the same way by taking a user-defined function, and the upper and lower boundaries of the integral. Here is a brief example:

from scipy import integrate
def myfunc(x, a, b):
return (x**b) + a
 
# These are the arguments that will be passed as a and b to myfunc()
args = (1.0, -2.0)
 
# Integrate myfunc() from 0.5 to 1.5
results = integrate.quad(myfunc, 0.5, 1.5, args)
print 'Integral = ', results[0], ' with error = ', results[1]

The well-established Fortran QUADPACK integration code lies underneath quad() while romberg() appears to be a python only implementation (??). While I would much rather use quad() for its robustness, it only supports user-functions that accept single value floats; while romberg() will pass in vectors of values. For some of my applications, myfunc() is computationally intensive, so there is an order of magnitude speed gain by having a properly vectorized myfunc() used in vector-form by the integration code. Only with romberg() have I matched the speed of IDL’s qpint1d (from Markwardt library, also based on QUADPACK, but vectorized) for identical myfunc().

{ 3 comments… read them below or add one }

1 M June 23, 2010 at 6:13 pm

You might find that the interpolation package (splrep and splint) is faster and more robust than the integration package (quad or romberg) for doing integration, too! If you have a reasonably smooth function (such that the piece-wise polynomial spline is a good model) you can evaluate it at a reasonably small number of points and do the (very fast) spline integration. This can improve speed for a single integration because the number of function evaluations needed to get an accurate spline model may be small compared to the number of evaluations for quad(), but it is particularly useful if you will be integrating the same curve lots of times (i.e., if only the limits of the integration change) or for surface and volume integrals.

Reply

2 Jessica April 28, 2011 at 11:53 am

One of our readers ran into an issue with passing in vector quantities into the integration function when using romberg(). You can find the full discussion of the issue and a suggested solution on the scipy-users mailing list:

http://old.nabble.com/Speedup-with-integrate.romberg%28%29-and-vector-arguments-td31485594.html

Reply

3 bonaventure May 5, 2016 at 10:45 am

How do I integrate an interpolated function in Python. I am interpolate Spectrum data and I want to integrate the approximative function from interpolation. Thank you

Reply

Leave a Comment

Previous post:

Next post: