Méthodes numériques

Y. Morel
Télécharger le document au format pdf   pdf icon     

Résolution d'une équation de la forme $f(x)=0$



Résolution par dichotomie

La récherche d'une solution approchée à l'équation $ f(x)=0$ sur l'intervalle initial $ [a;b]$ se fait en découpant à chaque itération l'intervalle en deux sous intervalles de même longueur, puis en sélectionnant celui contenant la solution.

Ce nouvel intervalle est alors redécoupé en deux...et ainsi de suite jusqu'à ce que la longueur de l'intervalle dans lequel on recherche la solution soit inférieure à la précision $ \varepsilon$ souhaitée.

$\displaystyle \fbox{\begin{minipage}[t]{6cm} Initialisation de $a$, $b$\ et $\... ... \hspace*{\ProgIndent} Fin Si \\ Fin Tant que \\ Afficher $x$ \end{minipage}}$

Méthode de Newton

C'est une méthode de descente, c'est-à-dire que l'approximation trouvée à chaque itération est meilleure que celle trouvée à l'étape précédente.

En partant d'une valeur $x$ , on "descend" donc pour la méthode de Newton le long de la tangente (dont le coefficient directeur est la dérivée de $ f$ en $ x$ : $ f'(x)$ ) à la courbe représentative de $ f$ jusqu'à l'axe des abscisses. On utilise alors l'abscisse de ce point comme nouvelle valeur pour $ x$ , jusqu'à atteindre la précision $ \varepsilon$ souhaitée.

Initialisation de $x$\ \\ $x_{new}... ...} $x_{new}=x-f(x)/f'(x)$\\ Fin Tant que \\ Afficher $x_{new}$

Méthode de la sécante

Comme on peut ne pas connaître explicitement la fonction $ f$ (pas de formule algébrique donnant $ f(x)$ ), où qu'il peut être plus ou moins fastidieux de déterminer sa dérivée, on peut penser à utiliser une approximation plus facile à calculer pour $ f'(x)$ .
Dans la méthode de la sécante, on approxime la tangente à la courbe de $ f$ par une sécante; cela revient à approximer la dérivée $ f'(x_i)$ dans la méthode de Newton par le taux d'accroissement: $ \displaystyle f'(x_i)\simeq \dfrac{f(x_i)-f(x_{i-1}}{x_i-x_{i-1}}$ .
L' algorithme devient:
$Initialisation de $x_1$\ et $x_2... ...ght)/\left(x_2-x_1\right)$\ \\ Fin Tant que \\ Afficher $x_3$



Calcul numérique d'intégrales


On cherche à calculer l'intégrale $ \displaystyle \int_a^b f(t)\,dt$ .
On discrétise l'intervalle $ [a;b]$ en $ n+1$ points: $ t_0=a$ ; $ t_1$ ; $ t_2$ ; ..., $ t_n=b$ , régulièrement, c'est-à-dire que le pas $ h_i=t_{i+1}-t_i=h$ est constant.
On a alors,
$\int_a^b f = \sum_{i=0}^{n-1} \int_{t_i}^{t_{i+1}} f \,.$

A l'aide d'un changement de variable, on ramène les $ n$ intégrales de cette somme à l'intervalle  $ [0;1]$ :
$\int_{t_i}^{t_{i+1}} f =h\int_0^1 g(t)\,dt$     avec $\quad g(t)=f(x_i+th) \,.$

Il reste à calculer cette intégrale.


Méthode des rectangles


On approxime "grossièrement" l'aire sous la courbe de la fonction par l'aire d'un rectangle.
Rectangle à gauche: $\int_0^1 g \simeq g(0)$
Rectangle à gauche 		Rectangle à droite: $\int_0^1 g \simeq g(1)$
Rectangle à gauche


Méthode du point milieu


$\int_0^1 g \simeq g(1/2)$ Rectangle à gauche

Méthode des trapèzes


$\int_0^1 g \simeq \dfrac{1}{2}\Bigl( g(0)+g(1) \Bigr)$ Rectangle à gauche

Méthode de Simpson


On interpole la fonction $ g$ par un polynôme $ P$ du second degré aux points $ (0;g(0))$ , $ (1/2;g(1/2))$ , et $ (1;g(1))$ (polynôme d'interpolation).
On approxime alors l'intégrale de $ g$ par celle du polynôme d'interpolation, c'est-à-dire l'aire sous la courbe par l'aire sous la parabole.
\begin{array}{ll} \displaystyle \int_0^1 g &\displaystyle \si... ...space{0.2cm}\\ &=\dfrac{1}{6}\left[2a+3b+6c\right] \end{array} \begin{pspicture}(0.5,-0.8)(1.5,3.2) %\psset{xunit=1.5cm,yunit=1.4cm}%,algebrai... ...0.5 2.2 mul 2 exp mul add 2 add) \rput(0.5,-0.4){$\frac{1}{2}$} \end{pspicture}

Il reste à déterminer les coefficients $ a$ , $ b$ et $ c$ du trinôme du second degré. On écrit pour cela les relations d'interpolations aux trois points:
\begin{array}{ll} \left\{\begin{array}{ll} P(0)=c=g(0) \\ P(... ...0)-b=2g(1)+2g(0)-4g(\frac{1}{2}) \end{array}\right. \end{array}

En remplaçant ces valeurs dans l'expression de l'intégrale, on obtient:
$\int_0^1 g \simeq \dfrac{1}{6}\left[2a+3b+c\right] =\dfrac{1}{6}\left[g(1)+4g\left(\frac{1}{2}\right)+g(0)\right] $

Plus généralement:

Toutes les méthodes précédentes consistent à interpoler la fonction par un polynôme, puis à remplacer le calcul de l'intégrale de celle-ci par l'intégrale du polynôme d'interpolation.
Pour les méthodes des rectangles, on interpole la fonction en un seul point avec un polynôme constant (à gauche, milieu ou droite de l'intervalle); la méthode des trapèzes utilise l'interpolation de la fontion par un polynôme du premier degré ( $ P(x)=ax+b$ ) aux deux extrémités de l'intervalle; enfin, la méthode de Simpson utilise un polynôme d'interpolation de degré 2.
On peut construire ainsi des méthodes d'intégration numérique d'ordre quelconque, en augmentant le degré du polynôme d'interpolation. Il s'agit plus généralement des méthodes de Newton-Cotes.


Méthode de Gauss

Les méthodes numériques précédentes utilisent toutes une discrétisation régulière de l'intervalle d'intégration. On peut penser, au lieu de découper cet intervalle régulièrement, à choisir des points de manière plus adaptée.
C'est ce à quoi aboutissent les méthodes de Gauss pour l'intégration numérique ...


Résolution numérique d'EDO


On considère l'équation différentielle ordinaire (EDO) du premier ordre
$y'=f(t,y) \quad,\quad t\in[a;b]\,.$


EDO et intégration numérique

Après discrétisation de l'intervalle $ [a;b]$ , et en intégrant sur l'intervalle $ [t_i;t_{i+1}]$ , on obtient:
$\int_{t_i}^{t_{i+1}}y'(t)\,dt=y(t_{i+1})-y(t_i)=\int_{t_i}^{t_{i+1}} f(t,y)\,dt $

soit,
$\displaystyle y_{i+1}=y_i+\int_{t_i}^{t_{i+1}} f(t,y)\,dt\,. $

Définir un schéma d'intégration numérique pour une EDO revient alors à utiliser une méthode d'intégration numérique (rectangle, trapèzes, Simpson, ...).


Méthodes d'Euler ( $ \simeq$ 1770)


La méthode des rectangles à gauche fournit la méthode d'Euler:
$y_{i+1}=y_i+\int_{t_i}^{t_{i+1}} f(t,y)\,dt \simeq y_i+h f(t_i,y_i)$

On peut utiliser aussi la méthode des rectangles à droite. On aboutit alors au schéma d'Euler implicite, ou rétrograde:
$y_{i+1}=y_i+\int_{t_i}^{t_{i+1}} f(t,y)\,dt \simeq y_i+h f(t_{i+1},y_{i+1})$

C'est un schéma implicite (expression de l'inconnue $ y_{i+1}$ en fonction d'elle-même, $ y_{i+1}$ ).
On résout cette équation, dont la seule inconnue est $ y_{i+1}$ par une méthode numérique de résolution d'équation (dichotomie, tangente, Newton, ...).
Ce schéma est donc un peu plus compliqué, et un peu plus lourd à mettre en \oeuvre; néanmoins, il a l'avantage d'être inconditionnellement stable (au contraire de la méthode d'Euler directe).
Avec la même idée, on peut penser à utiliser une intégration numpérique par la méthode des trapèzes. La méthode d'Euler obtenue est souvent alors appelée méthode d'Euler modifiée:
$y_{i+1}=y_i+\int_{t_i}^{t_{i+1}} f(t,y)\,dt \simeq y_i+\dfrac{h}{2}\Bigl(f(t_i,y_i)+f(t_{i+1},y_{i+1})\Bigr) $

Formule de Taylor et différences finies

Les méthodes de différences finies reposent sur la formule de taylor: au voisinage de $ t_i$ , on a à l'ordre $ n+1$ ,
$y(t)=y(t_i)+(t-t_i)y'(t_i)+\dfrac{(t-t_i)^2}{2}y''(t_i) +\dots+ \dfrac{(t-t_i)^n}{n!}y^{(n)}(t_i) +O\left((t-t_i)^{n+1}\right) $

d'où, en $ t_{i+1}$ , et pour $ h=t_{i+1}-t_i$ petit,
$y_{i+1}=y_i+hy'(t_i)+\dfrac{h^2}{2}y''(t_i) +\dfrac{h^3}{6}y^{(3)}(t_i) +\dots +\dfrac{h^n}{n!}y^{(n)}(t_i) +O\left(h^{n+1}\right) $


$ \blacksquare$ Au premier ordre, on a ainsi, $ y_{i+1}=y_i+hy'(t_i)+O\left(h^2\right) \iff y'(t_i)=\dfrac{y_{i+1}-y_i}{h}+O\left(h\right)$ , d'où l'approximation $ y'(t_i)\simeq \dfrac{y_{i+1}-y_i}{h}$ qui est utilisé par la méthode d'Euler. C'est une approximation de la dérivée d'ordre $ 1$ , c'est-à-dire pour laquelle l'erreur commise tend vers 0 comme le pas $ h$ .


$ \blacksquare$ Au second ordre:
\begin{array}{lll} \mbox{en } t=t_{i_+1} , \ &y_{i+1}&=y_i+h... ...i-hy'(t_i)+\dfrac{h^2}{2}y''(t_i)+O\left(h^3\right) \end{array}

soit, en soustrayant ces deux relations:
$y_{i+1}-y_{i-1}=2hy'(t_i)+O\left(h^3\right)$

d'où,
$y'(t_i)=\dfrac{y_{i+1}-y_{i-1}}{2h}+O\left(h^2\right)$

C'est un schéma centré, fournissant une approximation de la dérivée d'ordre 2.


$ \blacksquare$ Ordre supérieur:
on peut aussi par exemple, et pour plus de précision, prendre en compte les points $ y_{i-2}$ et $ y_{i+2}$ dans le calcul de $ y'(t_i)$ .
\begin{array}{ll} y_{i+2} &=y_i+(2h)y'(t_i)+\dfrac{(2h)^2}{2}... ...)+\dfrac{(-2h)^4}{24}y^{(4)}(t_i)+O\left(h^5\right) \end{array}

soit,
$\begin{array}{rll} y_{i+2} &=y_i+2hy'(t_i)+2h^2y'... ...(t_i)+0\tm h^3y^{(3)}(t_i)+0\tm h^4 y^{(4)}(t_i) +o\left(h^4\right) \end{array}$

d'où,
$y'(t_i)=\dfrac{-y_{i+2}+8y_{i+1}-8y_{i-1}+y_{i-2}}{12h} +O\left(h^4\right)$

C'est un schéma centré, fournissant une approximation de la dérivée $ y'(t_i)$ d'ordre 4.
On peut aussi, au lieu d'utiliser les valeurs $ y_{i+2}$ et $ y_{i-2}$ , utiliser les points "milieux" $ y_{i+\frac{1}{2}}$ et $ y_{i-\frac{1}{2}}$ , ou encore que les valeurs "passées" $ y_{i-1}$ , $ y_{i-2}$ , $ y_{i-3}$ , ...fournissant ainsi un schéma explicite, donc plus aisé à mettre en \oeuvre.
La méthode des différences finies permet de discrétiser la plupart des équations différentielles ordinaires et des équations aux dérivées partielles de la physique assez simplement.


Différences finies et méthode d'Euler (1768)

Approximation de la dérivée $ y'(t_i)$ par la différence finie:
$y'(t_i)\simeq \dfrac{y(t_{i+1})-y(t_i)}{t_{i+1}-t_i} =\dfrac{y_{i+1}-y_i}{h}$

d'où la relation de récurrence:
$f(t_i,y_i)=y'(t_i)\simeq \dfrac{y_{i+1}-y_i}{h} \quad\iff\quad y_{i+1}=y_i+hf(t_i,y_i)$

Il s'agit donc d'une méthode d'ordre $ 1$ .


Méthodes de Runge Kutta (RK) ( $ \simeq 1900$)


Les méthodes de Runge et Kutta repose fondamentalement sur la formule de Taylor. L'idée de Runge et Kutta est d'écrire une solution approchée $ \tilde{y}$ de $ y$ où interviennent des évaluations de la fonction $ f$ (et non pas de ses dérivées, qui pourraient être fastidieuses à calculer), de manière à ce que cette solution coïncide avec le développement en série de Taylor à un ordre souhaité.
Runge-Kutta 1 (RK1): La formule de Taylor à l'ordre 1 s'écrit
$y(t+h)=y(t)+hy'(t)+O(h^2)$

Soit, en utilsant l'équation différentielle: $ y'=f(t;y)$ ,
$y(t+h)\simeq y(t)+hf(t;y)$

Elle coïncide donc avec la méthode d'Euler.
Runge-Kutta 2 (RK2): La formule de Taylor à l'ordre 2 s'écrit:
$y(t+h)=y(t)+hy'(t)+\dfrac{h^2}{2}y''(t)+O(h^3)$

dans laquelle on peut détailler à l'aide de l'équation différentielle:
$\left\{\begin{array}{ll} y'(t)=f(t;y) \vspace{0.2cm}\\ y''(t)=... ...(t)) =\partial_1 f(t;y(t))+f(t;y(t))\,\partial_2 f(t;y(t)) \end{array}\right.$

On substitue alors ces deux expressions dans le développement de Taylor:
$y(t+h)=y(t)+hf(t;y)+ \dfrac{h^2}{2}\Bigl[ \partial_1 f(t;y)+f(t;y)\,\partial_2 f(t;y) \Bigr]+O(h^3)$

L'idée de Runge et Kutta est d'introduire une combinaison linéaire dans le calcul de $ y(t+h)$ entre $ k_0=f(t;y)=y'(t)$ , la dérivée en $ t$ , et $ k_1=f(t+\alpha h,y+\beta hk_0)$ , la dérivée un "peu plus loin":
$y(t+h)=y(t)+h\Bigl[ Ak_0+Bk_1\Bigr]$

Ensuite, en appliquant le développement de Taylor de $ k_1$ , puis comparant à celui de $ y(t+h)$ précédent, on espère obtenir des valeurs pour les coefficients $ A$ , $ B$ , $ \alpha$ et $ \beta$ . Ce développement s'écrit:
\begin{array}{ll} k_1&=f(t+\alpha h,y+\beta hk_0) \vspace{0.... ...artial_1 f(t;y) +\beta hf(t;y)\,\partial_2 f(t;y) \end{array}

En reportant ces expressions de $ k_0$ et $ k_1$ dans la combinaison linéaire souhaitée de $ y(t+h)$ , on obtient:
$y(t+h)=y(t)+(A+B)h f(t;y) +B\alpha h^2 \,\partial_1 f(t;y) +B\beta h^2 f(t;y)\,\partial_2 f(t;y)$

puis, en identifiant avec la série de Taylor à l'ordre 2, on trouve que les coefficients doivent vérifier
$\left\{\begin{array}{ll} A+B&=1 \vspace{0.2cm}\\ B\alpha&=\dfrac{1}{2} \vspace{0.2cm}\\ B\beta&=\dfrac{1}{2} \end{array}\right.$

Ainsi, lorsque les paramètres $A$ , $ B$ , $ \alpha$ et $ \beta$ vérifient ces relations, la combinaison linéaire proposée aura la même ordre d'approximation que la série de Taylor, c'est-à-dire un ordre 2.
Ce système est surdéterminé: seulement 3 équations pour 4 inconnues, et admet donc une infinité de solutions.
On peut fixer par exemple arbitrairement $ A=0$ , et alors, $ B=1$ , $ \alpha=\beta=\frac{1}{2}$ , et le schéma correspondant s'écrit:
$y_{i+1}=y_i+h k_1$    avec, $\quad\left\{\begin{array}{ll} k_0&=f(t_i;y_i) \vspace{0.2cm}\\ k_1&=f\left(x_i+\frac{1}{2}h;y_i+\frac{h}{2}k_0\right) \end{array}\right.$

Runge-Kutta 4 (RK4): L'idée de la méthode reste la même que précédemment. On écrit $ y_{i+1}$ comme une combinaison linéaire:
$y_{i+1}=y_i+h\Bigl( A k_0 +B k_1 +C k_2 + D k_2\Bigr)$

avec
$\left\{\begin{array}{ll} k_0=f(t_i;y_i) \vspace{0.2cm}\\ k_1=f... ...\right)\vspace{0.2cm}\\ k_3=f\left(t_i+h; y_i+hk_3\right) \end{array}\right.$

En utilisant le développement de Taylor à l'ordre 4:
$y(t+h)=y(t)+hy'(t)+\dfrac{h^2}{2}y''(t) +\dfrac{h^3}{6}y'''(t) +\dfrac{h^4}{24}y^{(4)}(t) +O(h^5)$

et en identifiant les coefficients, après écrit les développements de Taylor au même ordre de $ k_1$ , $ k_2$ et $ k_3$ , on trouve $ A=D=\frac{1}{6}$ , et $ B=C=\frac{2}{6}$ , d'où le schéma RK4:
$y_{i+1}=y_i+\dfrac{h}{6}\Bigl( k_0 + 2k_1 + 2k_2 + k_2\Bigr)$

avec les valeurs de $ k_0$ , $ k_1$ , $ k_2$ et $ k_3$ précédentes.

Voir aussi


LongPage: h2: 5 - h3: 12