16.6 Analyse numérique des équations différentielles

16.6.1 Méthode d’Euler

Soit f : J × E → E de classe {C}^{1} , soit {t}_{o} dans J , {y}_{o} dans E et φ : {J}_{0} → E la solution maximale vérifiant la condition φ({t}_{o}) = {y}_{o}. {J}_{0} est un intervalle contenu dans J et contenant {t}_{o}, et dire que la solution est maximale, c’est dire que φ ne peut se prolonger en une solution de l’équation différentielle sur un intervalle strictement plus grand que {J}_{0}. Notre but est de trouver une approximation de la fonction φ .

Pour cela soit h un nombre réel suffisamment petit et t dans {J}_{0} tel que t + h appartienne encore à {J}_{0}. Alors φ(t + h) est peu différent de φ(t) + hφ'(t). Mais φ'(t) = f(t,φ(t)). Donc φ(t + h) est peu différent de φ(t) + hf(t,φ(t)) = y + f(t,y) si y = φ(t). Ceci nous amène à définir pour un nombre réel h donné

Nous espérons bien entendu que la fonction {φ}_{h} ainsi définie sera une approximation de φ pour h petit. Nous allons montrer que c’est effectivement le cas, tout au moins sur un segment [a,b] contenu dans {J}_{0} et dans le cas simple où E = ℝ (bien que le résultat reste valable si E est un espace vectoriel normé).

Remarquons tout d’abord que puisque f est de classe {C}^{1} et que φ'(t) = f(t,φ(t)) , φ' est dérivable et

\begin{eqnarray*} φ''(t)& =&{ ∂f \over ∂t} (t,φ(t)) + φ'(t){ ∂f \over ∂y} (t,φ(t)) %& \\ & =&{ ∂f \over ∂t} (t,φ(t)) + f(t,φ(t)){ ∂f \over ∂y} (t,φ(t))%& \\ \end{eqnarray*}

Soit α > 0 et soit K = \{(t,y) ∈ [a,b] × E\mathrel{∣}φ(t) − α ≤ y ≤ φ(t) + α\}.

K est un compact qui contient le graphe de φ. La fonction (t,y)\mathrel{↦}{ ∂f \over ∂t} (t,y) + f(t,y){ ∂f \over ∂y} (t,y) est continue sur ce compact, donc bornée. Soit

M ={\mathop{ sup}}_{(t,y)∈K}|{ ∂f \over ∂t} (t,y) + f(t,y){ ∂f \over ∂y} (t,y)|

Alors on a pour tout t dans [a,b] |φ''(t)|≤ M. La formule de Taylor-Lagrange nous donne alors φ(t + h) = φ(t) + hφ'(t) +{ {h}^{2} \over 2} φ''(ξ), d’où

\left |{ φ(t + h) − φ(t) \over h} − f(t,φ(t))\right |≤ M{ |h| \over 2}

D’autre part la fonction (t,y)\mathrel{↦}{ ∂f \over ∂y} (t,y) est continue sur K donc bornée. Soit A ={\mathop{ sup}}_{(t,y)∈K}|{ ∂f \over ∂y} (t,y)|. Alors on a si (t,y) ∈ K et (t,y') ∈ K, f(t,y) − f(t,y') = (y − y'){ ∂f \over ∂y} (t,z) pour un certain z appartenant à [y,y'], donc |f(t,y) − f(t,y')|≤ A|y − y'| .

Définissons alors une suite {({t}_{i})}_{i∈ℕ} par {t}_{i} = {t}_{o} + ih, une suite {({y}_{i})}_{i∈ℕ} par {y}_{i+1} = {y}_{i} + hf({t}_{i},{y}_{i}). Nous allons mesurer l’erreur {e}_{i} = |{y}_{i} − φ({t}_{i})|. Comme nous ne sommes malheureusement pas sûrs que les couples ({t}_{i},{y}_{i}) appartiennent à K nous allons définir une fonction g : [a,b] × ℝ → ℝ par

g(x,y) = \left \{ \cases{ f(t,y) &si φ(t) − α ≤ y ≤ φ(t) + α \cr f(t,φ(t) − α)&si y < φ(t) − α \cr f(t,φ(t) + α)&si y > φ(t) + α)) } \right .

et une suite {({z}_{i})}_{i∈ℕ} par {z}_{i+1} = {z}_{i} + hg({t}_{i},{z}_{i}). Il est clair que {z}_{i} = {y}_{i} tant que ({t}_{i},{y}_{i}) appartient à K. Posons {ε}_{i} = |{z}_{i} − φ({t}_{i})|. La fonction g est continue sur [a,b] × ℝ et vérifie |g(t,y) − g(t,y')|≤ A|y − y'| pour tout t ∈ [a,b] et tous y,y' ∈ ℝ d’après (2). On a alors

\begin{eqnarray*}{ ε}_{i+1}& =& |{z}_{i+1} − φ({t}_{i+1})| = |{z}_{i} + hg({t}_{i},{z}_{i}) − φ({t}_{i} + h)|%& \\ & ≤& |{z}_{i} − φ({t}_{i})| + |h||g({t}_{i},{z}_{i}) − g({t}_{i},φ({t}_{i}))| %& \\ & & +|φ({t}_{i}) + hg({t}_{i},φ({t}_{i})) − φ({t}_{i} + h)| %& \\ & ≤& |{z}_{i} − φ({t}_{i})| + A|h||{z}_{i} − φ({t}_{i})| %& \\ & & +|h|\left |{ φ({t}_{i} + h) − φ({t}_{i}) \over h} − f({t}_{i},φ({t}_{i}))\right | %& \\ \end{eqnarray*}

car g(t,φ(t)) = f(t,φ(t)). On obtient donc en utilisant (1) {ε}_{i+1} ≤ (1 + A|h|){ε}_{i} + M{ |h{|}^{2} \over 2} . Comme {ε}_{o} = {e}_{o} = 0, on a donc par récurrence

\begin{eqnarray*}{ ε}_{i}& ≤& M{ |h{|}^{2} \over 2} (1 + (1 + A|h|) + ... + {(1 + A|h|)}^{i−1})%& \\ & =& M{ |h{|}^{2} \over 2} { {(1 + A|h|)}^{i} − 1 \over |h|A} %& \\ \end{eqnarray*}

soit, puisque 1 + x ≤ {e}^{x},

\begin{eqnarray*}{ ε}_{i}& ≤ M&{ |h| \over 2} { {e}^{Ai|h|}− 1 \over A} ≤|h|{ M \over 2A} ({e}^{A|t−ti|}− 1)%& \\ & ≤ & |h|{ M \over 2A} ({e}^{A(b−a)} − 1) %& \\ \end{eqnarray*}

On voit donc que pour h assez petit, on a pour tout i tel que {t}_{i} appartienne à [a,b], {ε}_{i} ≤ α, donc {z}_{i} = {y}_{i}, et donc {ε}_{i} = {e}_{i}, avec une erreur {e}_{i} ≤|h|{ M \over 2A} ({e}^{A(b−a)} − 1).

Soit maintenant x ∈]{t}_{i},{t}_{i+1}[. Considérons la fonction affine g qui vérifie g({t}_{i}) = φ({t}_{i}) et g({t}_{i+1}) = φ({t}_{i+1}). Soit h(t) = φ(t) − g(t) − μ{ (t−{t}_{i})(t−{t}_{i+1}) \over 2} μ est choisi de telle sorte que h(x) = 0. Deux applications du théorème de Rolle à la fonction h qui s’annule en {t}_{i}, x et {t}_{i+1} montrent qu’il existe ξ tel que h''(ξ) = 0. Or h''(ξ) = φ''(ξ) − μ puisque g'' = 0. On a donc en écrivant que h(x) = 0,

φ(x) − g(x) = φ''(ξ){ (x − {t}_{i})(x − {t}_{i+1}) \over 2}

soit

|φ(x) − g(x)|≤ M{ {({t}_{i+1} − {t}_{i})}^{2} \over 8} = M{ {h}^{2} \over 8}

puisque l’on a vu que pour tout t dans [a,b], |φ''(t)|≤ M. Or g et {φ}_{h} sont affines sur [{t}_{i},{t}_{i+1}] et donc

\begin{eqnarray*} |g(x) − {φ}_{h}(x)|& ≤& \mathop{max}(|g({t}_{i}) − {φ}_{h}({t}_{i})|,|g({t}_{i+1}) − {φ}_{h}({t}_{i+1})|)%& \\ & =& \mathop{max}({e}_{i},{e}_{i+1}) ≤|h|{ M \over 2A} ({e}^{A(b−a)} − 1) %& \\ \end{eqnarray*}

On en déduit donc que

\mathop{∀}x ∈ [a,b], |φ(x) − {φ}_{h}(x)|≤|h|{ M \over 2A} ({e}^{A(b−a)} − 1) + M{ {h}^{2} \over 8}

(en fait on n’a vu cette majoration que sur [{t}_{o},b], mais il suffit de changer h en − h pour avoir le même résultat sur [a,{t}_{o}], c’est pour cela que volontairement nous avons laissé les valeurs absolues partout). Les fonctions {φ}_{h} convergent donc uniformément vers φ sur [a,b] quand h tend vers 0, avec une majoration du type

\mathop{∀}x ∈ [a,b],|φ(x) − {φ}_{h}(x)|≤ B|h|

Dans la pratique il faut faire attention aux accumulations d’erreurs d’arrondis. L’erreur sur le calcul de {y}_{n} est de l’ordre de , où ε est la précision de calcul de l’ordinateur. On en déduit que l’erreur sur le calcul de {φ}_{h}(x) est de l’ordre de grandeur de { b−a \over h} ε. Soit une erreur totale du type Bh +{ b−a \over h} ε. Une étude simple de cette fonction montre que l’erreur totale est minimale pour des fonctions usuelles quand h est de l’ordre de \sqrt{ε}. On préférera prendre une valeur de h un peu trop grande, plutôt que trop petite. Il faut se méfier également du temps de calcul qui croit rapidement si l’on prend des valeurs de h trop petites.

16.6.2 Méthode de Runge et Kutta

La méthode est inspirée de la même idée que celle de la méthode d’Euler, mais on améliore l’approximation faite. Dans la méthode d’Euler nous prenions pour approximation de φ(t + h) l’expression φ(t) + hφ'(t) = φ(t) + hf(t,φ(t)). Ici on prendra comme approximation de φ(t + h) l’expression φ(t) + hk(t,h)k(t,h) est défini en posant

\begin{eqnarray*}{ k}_{1}(t,h)& =& f(t,φ(t)), %& \\ {k}_{2}(t,h)& =& f(t +{ h \over 2} ,φ(t) +{ h \over 2} {k}_{1}(t,h)),%& \\ {k}_{3}(t,h)& =& f(t +{ h \over 2} ,φ(t) +{ h \over 2} {k}_{2}(t,h)),%& \\ {k}_{4}(t,h)& =& f(t + h,φ(t) + h{k}_{3}(t,h)) %& \\ \end{eqnarray*}

et enfin

k(t,h) ={ 1 \over 6} ({k}_{1}(t,h) + 2{k}_{2}(t,h) + 2{k}_{3}(t,h) + {k}_{4}(t,h)).

On définit donc notre suite {y}_{i} par la relation de récurrence {y}_{i+1} = {y}_{i} + hk(i) où l’on a posé

\begin{eqnarray*}{ k}_{1}(i)& =& f({t}_{i},{y}_{i}), %& \\ {k}_{2}(i)& =& f({t}_{i} +{ h \over 2} ,{y}_{i} +{ h \over 2} {k}_{1}(i)), %& \\ {k}_{3}(i)& =& f({t}_{i} +{ h \over 2} ,{y}_{i} +{ h \over 2} {k}_{2}(i)), %& \\ {k}_{4}(i)& =& f({t}_{i} + h,{y}_{i} + h{k}_{3}(i)), %& \\ k(i)& =&{ 1 \over 6} ({k}_{1}(i) + 2{k}_{2}(i) + 2{k}_{3}(i) + {k}_{4}(i))%& \\ \end{eqnarray*}

pour i ≥ 0 , pour i < 0 on change h en − h.

On montre alors par un calcul pénible (à base de formule de Taylor) que

\left |{ φ(t+h)−φ(t) \over h} − k(t,h)\right |≤ M{ |h{|}^{4} \over 2} et la même démonstration que dans la méthode d’Euler montre qu’il existe une constante C telle que

\mathop{∀}x ∈ [a,b], |φ(x) − {φ}_{h}(x)|≤ C|h{|}^{4} + D|h{|}^{2}

Le terme en {h}^{2} provient en fait de l’interpolation linéaire. Aux points {t}_{i} l’erreur est en fait en C|h{|}^{4}. On obtient ainsi une convergence beaucoup plus rapide que dans la méthode d’Euler. L’étude de l’accumulation des erreurs montre que le meilleur h possible (pour les points {t}_{i}) est de l’ordre de \root{5}\of{ε}ε est la précision de l’ordinateur. On pourra par exemple prendre un h de l’ordre de 1{0}^{−2} ou 1{0}^{−3}.

16.6.3 Equations différentielles d’ordre supérieur

Il suffit de rappeler qu’une équation différentielle d’ordre p du type {y}^{(p)} = f(t,y,y',...,{y}^{(p−1)}) se ramène à un système différentiel

\left \{\array{ {y}_{1}' & = {y}_{2} \cr &\mathop{\mathop{…}}\cr {y}_{ p−1}'& = {y}_{p} \cr {y}_{p}' & = f(t,{y}_{1},{y}_{2},...,{y}_{p}) } \right .

en posant {y}_{1} = y,{y}_{2} = y',...,{y}_{p} = {y}^{(p−1)}. On appliquera donc l’une des deux méthodes précédentes à ce système.