Interpolation, dérivation et integration numérique
Interpolation
Supposons que l'on dispose d'un certain nombre de points, notés xi, auxquels on associe des valeurs yi, qui peuvent être par exemple des résultats de mesures. On souhaite déterminer une fonction qui permette de calculer la réponse y à un entrée x non comprise dans les xi de départ.
1) Lorsque les points xi sont "peu nombreux", on peut chercher à construire une fonction passant exactement par les yi. Dans cette optique, l'idée la plus simple est de chercher un polynôme p de degré n qui vérifie, pour des couples de réels (xi, pi) donnés,
p(xi) = pi.
La solution la plus évidente est d'écrire
p(x) = a0 + ... + aixi,
puis d'identifier les coefficients ai en écrivant les conditions p(xi) = pi pour chaque xi. On est amené à former la matrice, dite de Vandermonde, associée aux points xi. Il faut ensuite l'inverser pour déterminer p(x) ...
Une solution plus intéressante consiste à former la base des polynômes de Lagrange, définis par
qk(x) = (x-x0)...(x-xk-1)(x-xk+1)...(x-xn) / [(xk-x0)....(xk-xk-1)(xk-xk+1)...(xk-xn)]
On a qk(xj) = 1 si k = j et 0 sinon, et le polynôme p(x) = p0q0(x) + ... + pnqn(x) satisfait alors les conditions voulues.
Les problèmes d'interpolation de Lagrange font intervenir les valeurs de polynômes en certains points, mais ne tiennent pas compte des dérivées. Si les valeurs des dérivées sont données, on parle d'interpolation d'Hermite. Si l'on donne par exemple les valeurs d'une fonction et de sa dérivée en deux points x0 et x1, on peut construire un polynôme de degré 3 à partir d'une base de quatre polynômes vérifiant des conditions dont l'écriture est inspirée de celles de Lagrange.
L'interpolation d'une fonction pose des problèmes liés à l'amplitude des
oscillations engendrées par un polynôme de degré élevé. On peut limiter
cet inconvénient en réalisant des interpolations par intervalles. Avec
Lagrange, on aura une discontinuité de la dérivée première à chaque borne
d'intervalles, et de la dérivée seconde avec Hermite.
2) Lorsque les xi sont "nombreux", le degré des polynômes d'interpolation augmente et avec lui le risque de fortes oscillations ne correspondant pas forcément à la fonction que l'on souhaite approximer. On peut alors choisir de construire des approximations par morceaux de faible degré, les différents "morceaux" étant indépendants les uns des autres ou bien liés par des conditions de continuité de pente voire de courbure (c'est le principe des splines), ou encore des approximations globales ne passant pas exactement par les yi mais minimisant l'erreur entre les yi et les f(xi) (c'est le principe des méthodes de moindre carré).
- Les splines cubiques sont des polynômes de degré 3 permettant d'obtenir une continuité en pente et en courbure et qui constituent un bon compromis entre complexité et précision. Sur chaque intervalle, la spline est déterminée par les quatre coefficients du polynôme, ce qui porte à 4n le nombre d'inconnues pour les n-1 intervalles. L'égalité des valeurs de splines "voisines" sur les n-1 points intérieurs conduit à écrire 2(n-1) conditions. L'égalité des pentes et des courbures en ces n-1 points conduit à écrire 2(n-1) conditions supplémentaires. Les valeurs imposées en x1 et xn amènent 2 autres équations. Il manque donc deux équations pour pouvoir déterminer tous les coefficients, équations qui sont obtenues en spécifiant les valeurs des dérivées secondes en x1 et xn.
- Approximation linéaire par les moindres carrés : si une droite y = ax +
b doit minimiser la somme quadratique S des erreurs (yi-a-bxi),
on peut déterminer a et b en écrivant que dS/da et dS/db sont nuls, ce qui
conduit à résoudre un système de la forme
a.n + b.Σi=1,n xi = Σi=1,n yi
a.Σi=1,n xi + b. Σi=1,n xi2
= Σi=1,n xiyi
Dérivation numérique
On peut partir de l'opérateur de différence première progressive pour calculer la dérivée d'une fonction f en un point x0
f '(x0) = limh->0 [(f(x0+h)-f(x0)) / h]
mais également de ses homologues rétrograde et centrée :
f '(x0) = limh->0 [(f(x0)-f(x0-h))
/ h]
et
f '(x0) = limh->0 [(f(x0+h)-f(x0-h))
/ 2h]
On utilise numériquement une valeur de h "suffisamment" petite pour que
l'approximation soit assimilable à un passage à la limite.
Il faut ajouter une précision, cependant, sur la ... précision d'un calcul
numérique de dérivée. Si n est la précision relative du calculateur,
l'erreur absolue |x-xapp| obtenue sur l'évaluation xapp
du nombre x est de l'ordre de nx. Donc l'erreur absolue commise sur
l'évaluation de la dérivée en x0 est de l'ordre de 2n |f(x0)|
/ h. On voit donc que, si la précision relative du calculateur n'est pas
suffisante, un h trop petit aura des conséquences sur la précision de
l'approximation effectuée. C'est l'erreur d'arrondi.
Si les differences finies progressive et rétrograde sont la source d'une
erreur de troncature d'ordre 1 en h, la différence finie centrée est la
source d'une erreur d'ordre 2 en h.
On peut obtenir des approximations plus précises de la dérivée en un point, par exemple à l'aide de l'extrapolation de Richardson. Il est facile de montrer, par un développement limité à l'ordre 5, que l'on a une approximation de la différence finie centrée avec
f'(x0) + f'''(x0).h2 / 24 + O(h4)
En réécrivant ce développement avec h/2 au lieu de h, on arrive par combinaison linéaire des deux développements à :
f'(x0) = [8f(x0+h/4)-8f(x0-h/4)+f(x0-h/2)-f(x0+h/2)] / 3h
Le principe de l'extrapolation de Richardson est d'ailleurs valable dans bon nombre de cas. Supposons que l'on dispose d'une approximation numérique qapp(h) d'une certaine quantité qexa inconnue. L'approximation est fonction du paramètre numérique h. Supposons que l'approximation soit d'ordre n. On a donc
qexa = qapp(h) + cnhn + cn+1hn+1 + ...
Les constantes cn dépendent de la méthode numérique utilisée. Remplaçons h par h/2 dans l'équation précédente. On a alors
qexa = qapp(h/2) + cn(h/2)n + cn+1(h/2)n+1 + ...
En multipliant cette équation par 2n et en lui soustrayant la première, on peut faire disparaître le terme d'ordre n :
(2n-1)qexa = 2nqapp(h/2) - qapp(h) - cn+1hn+1/2 - ...
soit
qexa = [2nqapp(h/2) - qapp(h)] / (2n-1) + O(hn+1)
L'extrapolation de Richardson permet donc de gagner au moins un ordre de
convergence, voire davantage si cn+1 = 0 dès le départ.
Intégration numérique
Les méthodes d'intégration numérique s'appuient sur une approximation de
la fonction à intégrer par des fonctions plus simples, le plus souvent des
polynômes. C'est en faisant varier le degré de ces polynômes que l'on
obtient la plupart des méthodes présentées ci-dessous.
Méthodes de Newton-Cotes : trapèzes, rectangles, Simpson ...
Les méthodes de Newton-Cotes rassemblent les techniques d'intégration qui consistent à approximer une fonction par un polynôme, puis à intégrer ledit polynôme sur l'intervalle voulu. La plus simple est la méthode des trapèzes, où l'on choisit un polynôme du premier degré prenant les mêmes valeurs que la fonction aux bornes de l'intervalle. Cela revient à faire une somme pondérée de valeurs de la fonction en des points à déterminer, puisque l'on approche ∫a->b f(x) par
[(f(a) + f(b))(b-a) / 2]
Il est évidemment intéressant de découper l'intervalle principal, et de répéter la méthode des trapèzes dans chaque sous-intervalle. Si l'on note h = (b-a)/n la longueur de chacun, on peut montrer que l'erreur globale est donnée par (b-a)/12 f''(u) h2, pour u dans [a,b]. La méthode des trapèzes est donc d'ordre 2.
Le degré de précision d'une méthode d'intégration étant donné par la valeur maximale de n pour laquelle ladite méthode intègre exactement un polynôme de degré inférieur ou égal à n, il est évident que la méthode des trapèzes est de degré 1.
La méthode du rectangle est dans l'idée aussi simple que celle des trapèzes, puisqu'elle approche ∫a->b f(x) par
(b-a).f((a+b)/2)
On utilise ici la valeur de la fonction en le point milieu de chaque (sous-)intervalle, et l'on suppose que la fonction prend une valeur constante sur cet intervalle. On peut la ranger dans les méthodes de Newton-Cotes en considérant que cela consiste à approximer la fonction par un polynôme de degré 0 ...
Si l'on approxime la fonction par un polynôme de degré 2, on tombe sur les méthodes de Simpson. Il faut évidemment trois points dans chaque intervalle : les deux extrémités et le point milieu.. On a alors
∫a->b f(x) ~ (b-a) (f(a) + 4f((a+b)/2) + f(b)) / 3
On l'appelle formule de Simpson 1/3, du fait du coefficient qui apparaît dans la pondération. C'est en fait une moyenne des méthodes des trapèzes et des rectangles, qui associe aux trois points les poids respectifs 1/3, 4/3 et 1/3.
On peut montrer que la méthode de Simpson 1/3 composée (c'est-à-dire appliquée après décomposition en sous-intervalles) est d'ordre 4 en h, et de degré 3.
En utilisant un polynôme de degré 3, on tombe sur la méthode de Simpson 3/8 :
∫a->b f(x) ~ 3/8 (f(a) + 3f((2a+b)/3) + 3f((a+2b)/3) + f(b)) (b-a)
Cette variante n'est pourtant pas plus performante (en termes d'erreur et de degré) que sa cousine 1/3, qu'on lui préfère donc généralement puisque moins coûteuse.
La formule de Boole part d'un polynôme du quatirème degré, et on a l'approximation suivante :
∫a->b f(x) ~ 2/45 (7f(a) + 32f((3a+b)/4) + 12f((a+b)/2) + 32f((3a+b)/4) + 7f(b)) (b-a)
Cette approximation est d'ordre 6 et de degré 5.
Méthode de Romberg
La méthode de Romberg est basée sur une utilisation de la méthode des trapèzes combinée à l'extrapolation de Richardson. En effet, le terme d'erreur de la méthode des trapèzes ne contient que des puissances paires de h. On gagnera donc deux ordres de convergence au lieu d'un seul à chaque extrapolation. Notons T1,i le résultat obtenu avec 2i-1 intervalles. Pour passer de T1,i à T1,i+1, on doit doubler le nombre d'intervalles, donc diviser h par 2. On définit alors, grâce à l'extrapolation de Richardson :
T2,i = (22T1,i+1 - T1,i) / (22-1)
puis, de la même manière
T3,i = (24T2,i+1 - T2,i) / (24-1)
T4,i = (26T3,i+1 - T3,i) / (26-1)
etc ...
On peut alors définir un triangle de la forme
T1,1 T1,2 T1,3 T1,4 ...
T2,1 T2,2 T2,3 ...
T3,1 T3,2 ...
T4,1 ...
Chaque ligne est plus précise de deux ordres de convergence que la précédente. Et sur chaque ligne la précision augmente lorsque l'on se dirige vers la droite puisque la taille des sous-intervalles est divisée par 2 à chaque pas.
Quadratures de Gauss
L'idée de l'intégration de Gauss-Legendre est de placer les points xi de manière optimale (et non régulière comme dans les méthodes de Newton-Cotes) et de déterminer les poids wi associés de manière à ce que la quadrature soit exacte pour des polynômes de degré aussi grand que possible.
Ainsi, en raisonnant sur l'intervalle [-1,1] (que l'on peut évidemment ramener à n'importe quel intervalle [a.b] par changementde variable), on va déterminer les termes inconnus xi et wi de l'approximation ∫-1=>1 f(x) ~ w1f(x1) + w2f(x2) en appliquant cette quadrature aux fonctions f(x) = 1, f(x) = x, f(x) = x2 et f(x) = x3. On obtient un système non-linéaire mais résoluble analytiquement, qui donne les deux points +/- sqrt(1/3) et les poids associés 1 et 1. La formule de Gauss à deux points est de degré 3. Plus généralement, on peut montrer que la formule de Gauss-Legendre à n points est exacte pour des polynômes de degré 2n-1, et que les points d'intégration sont les racines des polynômes de Legendre définis par L0(x) = 1, L1(x) = x et la formule de récurrence :
(n+1)Ln+1(x) = (2n+1)xLn(x) - nLn-1(x)