CODEWITHSUNDEEP
pythonnumpyscipyintegral
asdf
asdf

Reputation: 698

Integral function of a interpolating polynomial in python

I have a set of points, and on those points I have used scipy in order to calculate the interpolating polynomial. I wish to have the primitive of that function 

self.p=interpolate.interp1d(self.__discreteDistribution['x'], self.__discreteDistribution['y'], kind='cubic')

I can easily use scipy to calculate the value of the integral over an interval, using 

integrate.quad(self.p, 0, max)

What I want instead is to have the primitive of self.p(). I have found sympy, but I do not have an analytical version of my interpolating polynomial.

What would you do in these occasion?

Upvotes: 0

Views: 2830

Answers (2)

ev-br
ev-br

Reputation: 26070

Assuming that you're using a piecewise interpolator (as opposed to a global polynomial interpolation), there are several ways you can do it with scipy:

Method 1: UnivariateSpline.

In [1]: import numpy as np

In [2]: x = np.arange(8)

In [3]: y = x

In [4]: from scipy.interpolate import interp1d

In [5]: from scipy.interpolate import interp1d, UnivariateSpline

In [6]: spl = UnivariateSpline(x, y, s=0)

In [7]: spl.<TAB>
spl.antiderivative        spl.get_coeffs            spl.roots
spl.derivative            spl.get_knots             spl.set_smoothing_factor
spl.derivatives           spl.get_residual          
spl.ext                   spl.integral              

In [8]: spl.integral(0, 1)
Out[8]: 0.5000000000000001

Two quirks of UnivariateSpline: first, use s=0 for interpolation (as opposed to least-square fitting). Second, watch out for extrapolation for out-of-bounds. By default, UnivariateSpline extrapolates out-of-bounds values (this can be controlled in constructor), but .integral assumes the spline is zero out of bounds.

In [9]: spl.integral(-1, 1)
Out[9]: 0.5000000000000001

Method 2: splev, splrep and splint.

In [13]: from scipy.interpolate import splev, splint, splrep

In [14]: tck = splrep(x, y, s=0)

In [15]: splint(0, 1, tck)
Out[15]: 0.5000000000000001

This is equivalent to using UnivariateSpline, only the interface is a bit different. See the docs for details.

Method 3: interp1d.

Under the hood, interp1d also uses b-splines (unless you request kind='linear' or 'nearest'), but the evaluation routines are different. interp1d constructs a callable, which can then be fed to a general-purpose integrator.

In [18]: from scipy.interpolate import interp1d

In [19]: interp = interp1d(x, y, kind='cubic')

In [20]: from scipy.integrate import quad

In [21]: quad(interp, 0, 1)
Out[21]: (0.5000000000000024, 5.5511151231258095e-15)

Again, watch out for out-of-bounds values: The behavior of the result constructed by interp1d is not very useful (even though it's controllable to an extent).

Upvotes: 5

mvw
mvw

Reputation: 5105

That interp1d method seems to just give you a function object for the interpolation function which allows you to evaluate arbitrary function values.

Looking at the documentation I see no interface to the internal representation. I guess it is a sum of piecewise defined cubic polynomials. In that case your primitive would be a sum of piecewise definied quadratic polynomials. Would that really be useful to you?

Aside calculating the splines directly (as suggested in the comments by askewchan) you could try to use the function values with approximate_taylor_polynomial (doc) to resample and get poly1d (doc) objects on each subinterval and then use poly1d.integ (doc) on each subinterval to get the coefficients of a primitive.

Upvotes: 1

Related Questions