Etant donné un ensemble de n+1 points (obtenu, par exemple, à la suite d'une expérience), on cherche un polynôme de degré n qui passe par tous ces points.
Soit donc une série de n+1 points (xi,yi) avec i compris entre 0 et n, on souhaite trouver les coefficients d'un polynôme d'interpolation de la forme :
P(x) = anxn + an-1xn-1 + ... + a1x + a0
Ce qui revient en fait à résoudre un système de n+1 équations à n+1 inconnues.
L'objectif de ce billet est de décrire, à l'aide de tableaux, les différentes étapes permettant de simplifier progressivement notre système d'équations initial, pour finalement aboutir à un système triangulaire facile à résoudre.
On montrera en particulier l'intérêt de la différence finie et de la différence divisée dans ce type de problème, et comment les utiliser pour trouver les valeurs des coefficients du polynôme, solutions du système de départ.
On proposera également à la fin une implémentation de l'algorithme en Python avec une représentation graphique d'un polynôme d'interpolation passant par 3 points.
Note importante : on obtiendra donc un polynôme exprimé sous une forme développée et réduite, la forme usuelle, contrairement par exemple aux polynômes d'interpolation de Newton ou de Lagrange.
Les notations utilisées pour représenter les différences finies sont adaptées aux variables indicées xi et yi.
II. Interpolation polynomiale pour des points uniformément répartis
Soit une série de n+1 points (x0,y0), (x1,y1), ..., (xn,yn), et on recherche un polynôme de la forme P(x) = anxn + ... + a1x + a0 qui passe par ces points.
Ceci revient à résoudre le système à n+1 équations :
P(x0) = y0
P(x1) = y1
...
P(xn) = yn
Qui peut également s'écrire avec f(x) = P(x) :
f(x0) = y0
f(x1) = y1
...
f(xn) = yn
Ou encore :
anx0n + ... + a1x0 + a0 = y0
anx1n + ... + a1x1 + a0 = y1
...
anxnn + ... + a1xn + a0 = yn
Les valeurs de xi et de yi étant connues, il s'agit donc de trouver les valeurs des coefficients ai dans ce système de n+1 équations à n+1 inconnues.
On notera la différence entre deux yi consécutifs :
Δyi = yi+1 - yi
Et la différence entre deux f(xi) successifs :
Δ[f](xi) = f(xi+1) - f(xi)
A l'ordre 2 :
Δ2yi = Δyi+1 - Δyi
Δ2[f](xi) = Δ[f](xi+1) - Δ[f](xi)
Plus généralement à l'ordre k :
Δkyi = Δk-1yi+1 - Δk-1yi
Δk[f](xi) = Δk-1[f](xi+1) - Δk-1[f](xi)
Note importante : On peut faire le lien avec la différence finie. Si h est l'intervalle entre deux xi successifs, alors la différence avant en xi s'écrit :
Δh[f](xi) = f(xi+h) - f(xi) = f(xi+1) - f(xi)
On n'a pas besoin de mentionner l'intervalle h dans la notation, étant donné que dans notre cas xi et yi sont des variables indicées.
D'autre part, la différence avant à l'ordre k et au point x de la fonction f définie par f(x) = anxn + an-1xn-1 + ... + a0 , ne conserve que les coefficients de degré supérieur ou égal à k (an, an-1, ..., ak).
Les coefficients de degré inférieur à k sont éliminés.
Cette propriété va nous permettre d'éliminer progressivement les coefficients en passant successivement de f(xi) à Δ[f](xi), puis à Δ2[f](xi), etc.
L'idée va être en fait de prendre la différence entre deux yi et deux f(xi) consécutifs, pour transformer les séries de yi et f(xi) , en les différences Δyi et Δ[f](xi), et pour passer d'un système de n+1 équations, à un système de n équations.
On va répéter cette procédure à l'ordre n, jusqu'à obtenir une équation à une seule inconnue, et au final un système triangulaire qui peut être résolu facilement.
II-A. Interpolation linéaire entre 2 points
L’interpolation linéaire est la méthode la plus simple pour estimer la valeur prise par une fonction continue entre deux points déterminés (interpolation). Elle consiste à utiliser pour cela la fonction affine, de la forme f(x) = ax + b, passant par les deux points déterminés.
Prenons un exemple simple :
Soient 2 points M0(1,3) et M1(2,5), on souhaite déterminer les coefficients a et b du polynôme P(x) = ax + b qui passe par ces points.
Ordre 0 : système à 2 équations
On utilise au départ la formule f(x)=ax + b, pour exprimer f(x0) et f(x1) en fonction de a et b.
Avec :
f(xi) = axi + b
On obtient ainsi :
f(x0) = f(1) = a + b
f(x1) = f(2) = 2a + b
Puis, en utilisant l'égalité f(xi) = yi, on a logiquement le système de 2 équations :
a + b = 3
2a + b = 5
Ordre 1 : équation finale
On calcule la différence entre les deux yi consécutifs :
Δy0 = y1 - y0 = 2
Et pour les deux f(xi) successifs, on a :
Δ[f](x0) = f(x1) - f(x0) = a
Comme f(xi)= yi, on en déduit que Δ[f](xi) = Δyi, et que :
a = 2
Finalement, en prenant la 1re équation de chacun des systèmes à l'ordre 0 et 1, on obtient le système triangulaire :
a + b = 3
a = 2
Représentation graphique du polynôme d'interpolation obtenu :
II-B. Interpolation parabolique pour 3 points uniformément répartis
L’interpolation parabolique consiste à rechercher un polynôme de la forme P(x) = ax2 + bx + c qui passe par 3 points donnés.
On posera comme précédemment :
f(x) = P(x)
Ordre 0 : système de 3 équations
On utilise au départ la formule f(x)=ax2 + bx + c, pour exprimer f(x0), f(x1) et f(x2) en fonction de a, b et c.
Avec :
f(xi) = axi2 + bxi + c
On obtient ainsi :
f(x0) = f(1) = a + b + c
f(x1) = f(2) = 4a + 2b + c
f(x2) = f(3) = 9a + 3b + c
Puis, en utilisant l'égalité f(xi) = yi, on a logiquement le système de 3 équations :
a + b + c = 4
4a + 2b + c = 11
9a + 3b + c = 22
Ordre 1 : système à 2 équations
On transforme les yi et f(xi), en Δyi et Δ[f](xi). Ceci va nous permettre d'éliminer c et d'obtenir ensuite un système de 2 équations à 2 inconnues.
On calcule pour cela les différences entre deux yi consécutifs :
Δy0 = y1 - y0 = 11 - 4 = 7
Δy1 = y2 - y1 = 22 - 11 = 11
Et pour f(xi) :
Δ[f](x0) = f(2) - f(1) = (4a + 2b + c) - (a + b + c ) = 3a + b
Δ[f](x1) = f(3) - f(2) = (9a + 3b + c) - (4a + 2b + c ) = 5a + b
Comme f(xi) = yi, on en déduit que Δ[f](xi) = Δyi, et que :
3a + b = 7
5a + b = 11
Ordre 2 : équation finale
On transforme les différences d'ordre 1, Δyi et Δ[f](xi), en différences d'ordre 2, Δ2yi et Δ2[f](xi).
Ceci nous permet d'éliminer b et d'obtenir une équation à 1 inconnue :
2a = 4
Finalement, en prenant la 1re équation de chacun des systèmes à l'ordre 0, 1 et 2, on obtient le système d'équations triangulaire :
a + b + c = 4
3a + b = 7
2a = 4
Représentation graphique du polynôme d'interpolation obtenu :
III. Interpolation polynomiale généralisée
Pour réaliser une interpolation polynomiale de points pas nécessairement uniformément répartis, on va utiliser la différence divisée entre 2 points.
Par exemple, les différences divisées d'ordre 0, 1 et 2 sont définies comme suit :
On obtient ainsi le schéma récursif suivant :
Les différences divisées interviennent dans les polynômes de Newton.
Contrairement aux différences finies vues précédemment, les différences divisées tiennent compte de la distance de 2 points par rapport à l'axe des abscisses.
Cette propriété va ainsi nous permettre de déterminer un polynôme de degré n passant par une série de points pas nécessairement uniformément répartis.
On notera les différences divisées entre deux yi et deux f(xi) consécutifs :
[yi, yi+1] = (yi+1 - yi) / (xi+1 - xi)
f[xi, xi+1] = (f(xi+1) - f(xi)) / (xi+1 - xi)
Ou plus généralement à l'ordre k :
[yi, ..., yi+k] = ([yi+1, ..., yi+k] - [yi, ..., yi+k-1]) / (xi+k - xi)
f[xi, ..., xi+k] = (f[xi+1, ..., xi+k] - f[xi, ..., xi+k-1]) / (xi+k - xi)
Note importante : la différence divisée à l'ordre k notée f[xi, ..., xi+k], avec f définie par f(x) = anxn + an-1xn-1 + ... + a0, ne conserve que les coefficients de degré supérieur ou égal à k (an, an-1, ..., ak).
Les coefficients de degré inférieur à k sont éliminés.
Cette propriété va nous permettre d'éliminer progressivement les coefficients en passant successivement de f(xi) à f[xi, xi+1], etc.
L'idée va être de prendre la différence divisée entre deux yi et deux f(xi) consécutifs, pour transformer les yi et f(xi), en les différences [yi, yi+1] et f[xi, xi+1], et pour passer d'un système de n+1 équations, à un système de n équations.
On va répéter cette procédure à l'ordre n, jusqu'à obtenir une équation à une seule inconnue, et au final un système triangulaire qui peut être résolu facilement.
III-A. Interpolation parabolique généralisée
On recherche un polynôme de la forme P(x) = ax2 + bx + c qui passe par 3 points donnés pas nécessairement répartis uniformément.
On posera à nouveau :
f(x) = P(x)
Ordre 0 : système à 3 équations
De la même façon que dans l'exemple précédent, à partir de la série de 3 points (1,4), (2,11) et (4,37), on obtient le système d'équations :
a + b + c = 4
4a + 2b + c = 11
16a + 4b + c = 37
Ordre 1 : système à 2 équations
On transforme les yi et f(xi), en les différences divisées [yi, yi+1] et f[xi, xi+1]. Ceci va nous permettre d'éliminer c et d'obtenir un système de 2 équations à 2 inconnues.
On évalue pour cela les différences divisées pour deux yi consécutifs :
[y0, y1] = (y1 - y0) / (x1 - x0) = (11 - 4) / (2 - 1) = 7
[y1, y2] = (y2 - y1) / (x2 - x1) = (37 - 11) / (4 - 2) = 13
Même chose pour les f(xi) :
f[x0, x1] = (f(x1) - f(x0)) / (x1 - x0) = ((4a + 2b + c) - (a + b + c)) / 1 = 3a + b
f[x1, x2] = (f(x2) - f(x1)) / (x2 - x1) = ((16a + 4b + c ) - (4a + 2b + c )) / 2 = 6a + b
Comme f(xi) = yi, on en déduit que f[xi, xi+1] = [yi, yi+1], et que :
3a + b = 7
6a + b = 13
Ordre 2 : équation finale
On transforme de la même manière les différences divisées d'ordre 1, [yi, yi+1] et f[xi, xi+1], en les différences d'ordre 2, [yi, yi+1, yi+2] et f[xi, xi+1, xi+2].
On utilise pour cela les formules des différences divisées à l'ordre 2 :
[yi, yi+1, yi+2] = ([yi+1, yi+2] - [yi, yi+1]) / (xi+2 - xi)
f[xi, xi+1, xi+2] = (f[xi+1, xi+2] - f[xi, xi+1]) / (xi+2 - xi)
Ceci nous permet d'éliminer b et d'obtenir une équation à 1 inconnue :
a = 2
Le système d'équations final s'écrit donc :
a + b + c = 4
3a + b = 7
a = 2
Note importante : on remarque que le coefficient de a dans la dernière équation est égal à 1, même chose pour le coefficient de b dans la 2e équation, et de c dans la 1re.
Représentation graphique du polynôme d'interpolation :
IV. Implémentation en Python
La fonction interpolation_polynomiale prend en arguments 2 listes x et y, représentant en fait la série de points (xi,yi), et renvoie la liste des coefficients du polynôme recherché.
Déroulé de la fonction :
- Transformation des 2 listes x et y, en 2 listes df et dy représentant les membres de gauche et de droite des équations du système initial ;
- Parcours des ordres de k=1 à k=n ;
- --- Parcours simultané des éléments des listes dy et df ;
- ------ Pour chaque paire d'éléments, calcul des différences divisées entre 2 valeurs successives de dy et de df ;
- ------ Copie des différences dans les 2 listes dy et df qui représentent le système d'équations à l'ordre k ;
- --- Ajout de la 1re équation du système d'ordre k au système final ;
- Résolution du système triangulaire final « en remontant » les valeurs des variables dans les équations.
| Code Python : | Sélectionner tout |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 | def interpolation_polynomiale(x,y): n=len(x)-1 # nombre de valeurs de la série - 1, il correspond au degré du polynôme, exemple pour n=2 : a2x^2 + a1x + a0, où a2, a1 et a0 représentent les coefficients dy = y # on copie la liste y dans dy : différences divisées à l'ordre 0 df = [] # on initialise la liste permettant d'enregistrer les sommes de termes correspondant aux membres de gauche des équations du système initial for i in range(n+1): # on parcourt les indices des x fi = somme_termes(x[i],n) # on remplace x par sa valeur dans le polynôme p(x) de degré n pour obtenir une somme de termes sous la forme [4, 2, 1] -> 4a2 + 2a1 + a0 df = df + [fi] # on ajoute la somme de termes fi à la liste des sommes df : différences divisées à l'ordre 0 systeme_triangulaire = [(df[0],dy[0])] # on initialise la liste permettant d'enregistrer le système triangulaire final : ajout de la 1re équation for k in range(1,n+1): # on parcourt les ordres suivants des différences divisées : ordre 1, 2, ..., n for i in range(n-k+1): # on parcourt les indices de x, dy et df jusqu'à n-k h = (x[i+k] - x[i]) # on détermine la différence entre xi+k et xi pour le calcul de la différence divisée à l'ordre k dyi = (dy[i+1] - dy[i])/h # on évalue la différence divisée à l'ordre k : [yi, ..., yi+k]) dfi = difference_divisee(df[i],df[i+1],h) # on détermine la différence divisée à l'ordre k : f[xi, ..., xi+k]). dy[i] = dyi # on remplace l'élément d'indice i de la liste dy par la valeur de la différence divisée dyi df[i] = dfi # on remplace l'élément d'indice i de la liste df par la valeur de la différence divisée dfi systeme_triangulaire = systeme_triangulaire + [(df[0],dy[0])] # on ajoute à la liste la 1re équation du système d'ordre k. a = [] # initialisation de la liste des coefficients remontés # on résout notre système triangulaire « en remontant » les valeurs des coefficients dans les équations for eq in reversed(systeme_triangulaire): # on parcourt la liste des équations dans l'ordre inverse ai = eval_coefficient(eq,a) # on évalue le coefficient ai a = a + [ai] # on ajoute la valeur du coefficient à la liste return a # on renvoie la liste des coefficients |
Exemple de test avec une série de 3 points :
| Code Python : | Sélectionner tout |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 | import matplotlib.pyplot as plt # module permettant de générer le graphique def correct_value(v): # corrige l'erreur éventuelle dans les résultats des opérations en virgule flottante return round(v,5) def somme_termes(x,n): s=[] # initialisation de la liste représentant la somme de termes for i in reversed(range(n+1)): # parcours des n termes du polynôme coef = (x**i) # évaluation du coefficient de ai (degré i) s = s + [coef] # ajout du coefficient de ai à la liste return s # on renvoie la liste représentant la somme de termes - membre de gauche de l'équation : 2a2 + a1 + a0 -> [2, 1, 1] def difference_divisee(f1,f2,h=1): # on initialise la liste permettant d'enregistrer les coefficients des ai, après calcul de la différence divisée entre les sommes de termes : # exemple pour h=1 : ((4a2 + 2a1 + a0) - (a2 + a1 + a0)) / h = 3a2 + a1-> [3, 1] dfi=[] for c1,c2 in zip(f1[:-1],f2[:-1]): # on parcourt simultanément les 2 listes f1 et f2 jusqu'à l'avant-dernier élément : le dernier est éliminé coef=(c2-c1)/h # calcul de la différence divisée pour chaque paire de termes: on soustrait 2 à 2 les coefficients des variables, puis on divise le résultat par l'intervalle h entre les xi: dfi = dfi + [coef] # on ajoute le résultat à la liste return dfi # on retourne la liste représentant la somme de termes : 3a2 + a1 -> [3, 1] def eval_coefficient(eq,a): fi = eq[0] # membre de gauche de l'équation correspondant à la somme de termes : [4, 2, 1] -> (4a2 + 2a1 + a0) yi = eq[1] # membre de droite de l'équation : nombre réel s=0 # initialisation de la somme de termes for j in range(len(fi)-1): # on parcourt les termes de la somme situés avant le terme de ai c = a[j] # copie de la valeur du coefficient d'indice j et de degré n-j s = s + fi[j]*c # on ajoute à la somme le produit du coefficient c et de fi[j] ai = (yi-s) # calcul de la valeur de ai return correct_value(ai) # on renvoie la valeur du coefficient ai def interpolation_polynomiale(x,y): n=len(x)-1 # nombre de valeurs de la série - 1, il correspond au degré du polynôme, exemple pour n=2 : a2x^2 + a1x + a0, où a2, a1 et a0 représentent les coefficients dy = y # on copie la liste y dans dy : différences divisées à l'ordre 0 df = [] # on initialise la liste permettant d'enregistrer les sommes de termes correspondant aux membres de gauche des équations du système initial for i in range(n+1): # on parcourt les indices des x fi = somme_termes(x[i],n) # on remplace x par sa valeur dans le polynôme p(x) de degré n pour obtenir une somme de termes sous la forme [4, 2, 1] -> 4a2 + 2a1 + a0 df = df + [fi] # on ajoute la somme de termes fi à la liste des sommes df : différences divisées à l'ordre 0 systeme_triangulaire = [(df[0],dy[0])] # on initialise la liste permettant d'enregistrer le système triangulaire final : ajout de la 1re équation for k in range(1,n+1): # on parcourt les ordres suivants des différences divisées : ordre 1, 2, ..., n for i in range(n-k+1): # on parcourt les indices de x, dy et df jusqu'à n-k h = (x[i+k] - x[i]) # on détermine la différence entre xi+k et xi pour le calcul de la différence divisée à l'ordre k dyi = (dy[i+1] - dy[i])/h # on évalue la différence divisée à l'ordre k : [yi, ..., yi+k]) dfi = difference_divisee(df[i],df[i+1],h) # on détermine la différence divisée à l'ordre k : f[xi, ..., xi+k]). dy[i] = dyi # on remplace l'élément d'indice i de la liste dy par la valeur de la différence divisée dyi df[i] = dfi # on remplace l'élément d'indice i de la liste df par la valeur de la différence divisée dfi systeme_triangulaire = systeme_triangulaire + [(df[0],dy[0])] # on ajoute à la liste la 1re équation du système d'ordre k. a = [] # initialisation de la liste des coefficients remontés # on résout notre système triangulaire « en remontant » les valeurs des coefficients dans les équations for eq in reversed(systeme_triangulaire): # on parcourt la liste des équations dans l'ordre inverse ai = eval_coefficient(eq,a) # on évalue le coefficient ai a = a + [ai] # on ajoute la valeur du coefficient à la liste return a # on renvoie la liste des coefficients def eval_polynome(a,x): s=0; n=len(a)-1 for i in range(n+1): s = s + a[i]*(x**(n-i)) return s # on teste la fonction pour une série de 3 points (-3, 17) (-1, 3), (3, 23). On cherche donc un polynôme de degré 2. x = [-3, -1, 3] # liste des valeurs d'abscisse y = [17, 3, 23] # liste des valeurs d'ordonnée # trace les points dont les abscisses x et les ordonnées y sont fournies en arguments. plt.scatter(x, y, color='blue') a = interpolation_polynomiale(x,y) # renvoie la liste des coefficients [an, ..., a1, a0] n= len(x) - 1 for i in range(n+1): # parcourt de la liste des coefficients obtenus print("a" + str(n-i) + " = " + str(a[i])) # affiche la valeur du coefficient d'indice i et de degré n-i pas = (x[-1]-x[0])/(n*10) # intervalle entre 2 valeurs d'abscisse x0 = x[0] # 1re valeur d'abscisse x = [] # initialisation de la liste des abscisses f = [] # initialisation de la liste destinée à enregistrer les valeurs successives du polynôme p(x) = f(x) i=0 # on génère la nouvelle liste des abscisses x et des ordonnées f pour tracer la courbe du polynôme d'interpolation (avec un intervalle entre les x plus petit). for i in range(n*10+1): # parcourt les n*10 indices (0, 1, 2, ...) xi = x0 + pas*i # copie de la valeur d'abscisse correspondant à l'indice i x = x + [xi] # ajout de la valeur d'abscisse à la liste fi = eval_polynome(a,xi) # évaluation du polynôme pour xi f = f + [fi] # ajout de la valeur de fi à la liste f # trace la courbe qui relie les points dont les abscisses x et les ordonnées f sont fournies en arguments. Polynôme utilisé P(x) = f(x) = 2*x^2 + x + 2 plt.plot(x, f, color='red') # copie du titre du graphique plt.title("Polynôme d'interpolation P(x) = 2x^2 + x + 2") plt.show() # affiche la figure à l'écran |
Résultat obtenu :
a2 = 2.0
a1 = 1.0
a0 = 2.0
Le code affiche ensuite le graphique représentant le polynôme d'interpolation passant par les 3 points :
Note importante : il faut préciser que l'erreur d'interpolation peut s'amplifier démesurément, dans le cas particulier où les points sont uniformément répartis, et lorsqu'on augmente le nombre de points pour un intervalle [x0 , xn] donné. On parle dans ce cas de phénomène de Runge.
V. Conclusion
Après avoir posé notre système d'équations initial, nous avons compris comment le réduire progressivement à l'aide des différences finies et des différences divisées, pour finalement obtenir les valeurs des coefficients du polynôme recherché.
Ceci nous a permis ensuite d'implémenter une fonction en Python qui évalue les coefficients d'un polynôme à partir de 2 listes de valeurs représentant la série de points (xi,yi).
Sources :
https://fr.wikipedia.org/wiki/Interp..._lin%C3%A9aire
https://fr.wikipedia.org/wiki/Interpolation_polynomiale
https://fr.wikipedia.org/wiki/Diff%C3%A9rence_finie
https://fr.wikipedia.org/wiki/Diff%C..._divis%C3%A9es
https://fr.wikipedia.org/wiki/Interpolation_newtonienne