Comparaison de 3 méthodes

Interpolation de Lagrange en utilisant

Rappel

Soit \(f\) une fonction et soient \(n+1\) réels distincts \((x_i)_{0\leq i\leq n}\) (dans l'ensemble de définition de \(f\)). On cherche à évaluer le polynôme, noté \(p\), d'interpolation de Lagrange de la fonction aux points \(x_i\).

Les différences divisées et la formule barycentrique sont grosso modo équivalentes en nombre d'opérations. Avec Python et

le but est de faire une fonction interp(x,f,xx) où xx est un vecteur (tableau dimension 1) qui renvoie les valeurs du polynôme en xx sous forme d'un vecteur.

Différences divisées

\[\begin{split}p(x)= & f[x_0] + f[x_0,x_1] (x-x_0) + \ldots \\ & + f[x_0,\ldots,x_n](x-x_0)\cdots (x-x_{n-1})\end{split}\]

Avec \(f[x_i,\ldots,x_{i+p}]\) les différences divisées

\[f[x_i]=f(x_i)\]\[f[x_i,\ldots,x_{i+p}]=\frac{f[x_{i+1},\ldots,x_{i+p}] - f[x_i,\ldots,x_{i+p-1}]}{x_{i+p}-x_i}\]

Évaluation via la méthode de Horner

Difficulté : aucune, il s'agit juste d'utiliser la vectorisation et un tableau de dimension 1 pour les différences divisées

Formule barycentrique

Vient des polynômes de Lagrange avec une réécriture pour obtenir un calcul en \(O(n^2)\).

\[\begin{split}w_i=\prod_{\substack{0\leq j\leq n\\j\neq i}}\frac{1}{x_i-x_j} \quad\text{pour tout $i\in\{0,\ldots,n\}$}\end{split}\]\[\forall x\notin \lbrace x_0,\ldots, x_n\rbrace \quad p(x)=\dfrac{\displaystyle \sum_{i=0}^n \dfrac{w_i}{x-x_i} y_i}{\displaystyle \sum_{i=0}^n \dfrac{w_i}{x-x_i}}.\]

Difficulté : gestion des valeurs interdites, vectorisation

Le code 1

scipy.interpolate.lagrange(x,y)

retourne le polynôme cherché, c'est un type poly1D

In [1]: import numpy as np

In [2]: import scipy.interpolate as interp

In [3]: x =np.linspace(-5,5,5)

In [4]: type(interp.lagrange(x,np.cos(x)))
 Out[4]: numpy.lib.polynomial.poly1d

In [5]: print(interp.lagrange(x,np.cos(x)))
         4             3          2
0.01384 x + 1.518e-18 x - 0.3747 x + 6.418e-17 x + 1

Le code 2

import numpy as np
import matplotlib.pyplot as plt
import sys
def diffdiv(x,y):
	n=np.size(x)
	if n != np.size(y):
	    sys.exit("erreur de taille")
	df=y.copy()
	for i in range(n):
	    for j in range(n-1,i,-1):
		    df[j]=(df[j]-df[j-1])/(x[j]-x[j-i-1])
	return df

def horner(x,df,xx):
    n=np.size(x)
    if n != np.size(df):
        sys.exit("erreur de taille")
    s=df[n-1]
    for i in range(n-1,0,-1):
        s=df[i-1]+s*(xx-x[i-1])
    return s

Le code 3

def cw(x):
    n , w =np.size(x), np.ones(np.size(x))
    for i in range(n):
        for j in range(n):
            if i!=j:
                w[i]=(x[i]-x[j])*w[i]
    return 1./w
def bary(x,y,xx):
    n, m =np.size(x) , np.size(xx)
    nu , de = np.zeros(m) , np.zeros(m)
    pb=(de!=0)   #valeurs à problèmes
    w=cw(x)	
    for i in range(n):
        dx=xx-x[i]
        pb=pb+(dx==0)
        dx[dx==0]=1
        nu=nu+y[i]*w[i]/dx
        de=de+w[i]/dx
    s=nu/de   # division vectorisée terme à terme
    for j in np.transpose(np.nonzero(pb)):  #traitement des valeurs 
        s[j]=y[np.transpose(np.nonzero((x-xx[j])==0))]  # à problèmes
    return s

Expérience code 1

Si \(f\) a toute ses dérivées successives uniformément bornées sur l'intervalle de définition alors quand le pas de la subdivision équidistante tend vers 0 le polynôme d'interplation de Lagrange converge uniformément vers la fonction.

\[n=20\]

(Source code)

_images/deux_inter-1.png

Expérience code 1 n=40

import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate as itg
x, xx =np.linspace(-5,5,40), np.linspace(-5,5,1000)
plt.plot(xx,itg.lagrange(x,np.cos(x))(xx),'b')
plt.show()

(Source code)

_images/deux_inter-2.png

Expérience code 2 n=40

In [1]: x, xx =np.linspace(-5,5,40), np.linspace(-5,5,1000)

In [2]: plt.plot(xx,horner(x,diffdiv(x,np.cos(x)),xx),'b') ;
_images/diffdiv1.png

Expérience code 2 n=65

In [1]: x, xx =np.linspace(-5,5,65), np.linspace(-5,5,1000)

In [2]: plt.plot(xx,horner(x,diffdiv(x,np.cos(x)),xx),'b') ;
_images/diffdiv2.png

Expérience code 2 n=100

In [1]: x, xx =np.linspace(-5,5,100), np.linspace(-5,5,1000)

In [2]: plt.plot(xx,horner(x,diffdiv(x,np.cos(x)),xx),'b') ;
_images/diffdiv3.png

Expérience code 3 n=65

In [1]: x, xx = np.linspace(-5,5,65), np.linspace(-5,5,1000)

In [2]: plt.plot(xx,bary(x,np.cos(x),xx),'g') ;
_images/bary1.png

Expérience code 3 n=100

In [1]: x, xx = np.linspace(-5,5,100), np.linspace(-5,5,1000)

In [2]: plt.plot(xx,bary(x,np.cos(x),xx),'g') ;
_images/bary2.png

Conclusions