Fonction intégrale d'un polynôme interpolant en python

Question

J'ai un ensemble de points, et sur ces points j'ai utilisé scipy pour calculer le polynôme interpolant. Je souhaite avoir la primitive de cette fonction

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

Je peux facilement utiliser scipy pour calculer la valeur de l'intégrale sur un intervalle, en utilisant

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

Ce que je veux au lieu est d'avoir la primitive de self.p(). J'ai trouvé sympy, mais je n'ai pas de version analytique de mon polynôme d'interpolation.

Que feriez-vous en ces occasions?

Answer 1

En supposant que vous utilisez un interpolateur piecewise (par opposition à une interpolation polynomiale global), il y a plusieurs façons dont vous pouvez le faire avec scipy:

Méthode 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 

Deux manies de UnivariateSpline: premièrement, utiliser s=0 pour l'interpolation (par opposition à ajustement des moindres carrés). Deuxièmement, faites attention à l'extrapolation pour les hors-limites. Par défaut, UnivariateSpline extrapole les valeurs hors limites (ceci peut être contrôlé dans le constructeur), mais .integral suppose que la spline est à zéro hors limites.

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

Méthode 2: splev, splrep et attelle. Ceci est équivalent à l'utilisation de UnivariateSpline, seule l'interface est un peu différente. Voir les documents pour plus de détails.

Méthode 3: interp1d.

Sous le capot, interp1d utilise également b-splines (sauf si vous demandez kind = 'linear' ou 'nearest'), mais les routines d'évaluation sont différentes. interp1d construit un appelable, qui peut ensuite être transmis à un intégrateur généraliste.

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) 

Encore une fois, méfiez-vous des valeurs hors limites: Le comportement du résultat construit par interp1d n'est pas très utile (même si elle est contrôlable dans une certaine mesure).

Answer 2

Ce interp1d method semble simplement vous donner un objet fonction pour la fonction d'interpolation qui vous permet d'évaluer des valeurs de fonction arbitraires.

En regardant la documentation je ne vois aucune interface à la représentation interne. Je suppose que c'est une somme de polynômes cubiques définis par morceaux. Dans ce cas votre primitive serait une somme de polynômes quadratiques définis par morceaux. Cela vous serait-il vraiment utile?

calcul Outre les cannelures directement (comme suggéré dans les commentaires par askewchan), vous pouvez essayer d'utiliser les valeurs de la fonction avec approximate_taylor_polynomial (doc) rééchantillonner et obtenir poly1d (doc) objets sur chaque sous-intervalle, puis utilisez poly1d.integ (doc) sur chaque sous-intervalle pour obtenir les coefficients d'une primitive.