Résolution numérique d’équations différentielles ordinaires
Cours ENPC — Pratique du calcul scientifique
Introduction
Objectif
$$
% % %
%
$$
\[ \definecolor{gris_C}{RGB}{96,96,96} \definecolor{blanc_C}{RGB}{255,255,255} \definecolor{pistache_C}{RGB}{176,204,78} \definecolor{pistache2_C}{RGB}{150,171,91} \definecolor{jaune_C}{RGB}{253,235,125} \definecolor{jaune2_C}{RGB}{247,208,92} \definecolor{orange_C}{RGB}{244,157,84} \definecolor{orange2_C}{RGB}{239,119,87} \definecolor{bleu_C}{RGB}{126,151,206} \definecolor{bleu2_C}{RGB}{90,113,180} \definecolor{vert_C}{RGB}{96,180,103} \definecolor{vert2_C}{RGB}{68,141,96} \definecolor{pistache_light_C}{RGB}{235,241,212} \definecolor{jaune_light_C}{RGB}{254,250,224} \definecolor{orange_light_C}{RGB}{251,230,214} \definecolor{bleu_light_C}{RGB}{222,229,241} \definecolor{vert_light_C}{RGB}{216,236,218} \]
But de ce chapitre
Étudier les méthodes numériques pour résoudre les problèmes de Cauchy : \[ \left \{ \begin{aligned} & \mathbf{\boldsymbol{x}}'(t) = \mathbf{\boldsymbol{f}}\bigl(t, \mathbf{\boldsymbol{x}}(t)\bigr), \\ & \mathbf{\boldsymbol{x}}(t_0) = \mathbf{\boldsymbol{x}}_0. \end{aligned} \right. \tag{1}\] où \(\mathbf{\boldsymbol{f}}\colon \mathbf R\times \mathbf R^n \to \mathbf R^n\).
Applications
Mécanique lagrangienne ;
Mécanique céleste ;
Circuits électriques ;
Équations aux dérivées partielles dépendant du temps : \[ \left \{ \begin{aligned} \partial_t u &= \partial_x^2 u \\ u(0, x) &= u_0(x). \end{aligned} \right. \xrightarrow{\text{discrétisation}} \left \{ \begin{aligned} \mathbf{\boldsymbol{x}}' &= \mathsf A \mathbf{\boldsymbol{x}}. \\ \mathbf{\boldsymbol{x}}(0) &= \mathbf{\boldsymbol{x}}_0. \end{aligned} \right. \] \({\color{green} \rightarrow}\) omniprésentes en sciences.
Remarque : équations différentielles d’ordre supérieur
Les problèmes de Cauchy d’ordre supérieur peuvent se ramener à la forme Équation 1.
Exemple : particule chargée (avec masse et charge 1) dans un champ électrique : \[ \mathbf{\boldsymbol{x}}''(t) = \mathbf{\boldsymbol{E}}\bigl(\mathbf{\boldsymbol{x}}(t)\bigr). \] En introduisant \(\mathbf{\boldsymbol{v}} = \mathbf{\boldsymbol{x}}'\), on réécrit cette équation sous la forme \[ \left \{ \begin{aligned} & \mathbf{\boldsymbol{v}}'(t) = \mathbf{\boldsymbol{E}}\bigl(\mathbf{\boldsymbol{x}}(t)\bigr), \\ & \mathbf{\boldsymbol{x}}'(t) = \mathbf{\boldsymbol{v}}(t). \end{aligned} \right. \]
Remarque : un cas particulier utile
Lorsque \(\mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{f}}(t)\), la solution est une intégrale : \[ \mathbf{\boldsymbol{x}}(t) = \mathbf{\boldsymbol{x}}_0 + \int_{t_0}^{t} \mathbf{\boldsymbol{f}}(t) \, \mathrm dt. \] \({\color{green} \rightarrow}\) utile pour deviner l’ordre de convergence des méthodes numériques
Analyse du problème continu
Définition : solution locale, maximale, globale
Une fonction \(\mathbf{\boldsymbol{x}} \colon I \to \mathbf R^n\) est une solution locale si \(t_0 \in I\) et \[ \left \{ \begin{aligned} & \mathbf{\boldsymbol{x}}'(t) = \mathbf{\boldsymbol{f}}\bigl(t, \mathbf{\boldsymbol{x}}(t)\bigr) \qquad \forall t \in I \\ & \mathbf{\boldsymbol{x}}(t_0) = \mathbf{\boldsymbol{x}}_0. \end{aligned} \right. \] Une solution locale \(\mathbf{\boldsymbol{x}} \colon I \to \mathbf R^n\) est maximale si elle ne se prolonge pas à un intervalle plus grand, et globale si \(I = \mathbf R\).
Exemple
Considérons l’équation \[ \left \{ \begin{aligned} & x'(t) = x(t)^2, \\ & x(0) = 1. \end{aligned} \right. \] La solution \(x_* \colon (-\infty, 1) \to \mathbf R\) donnée par \[ \mathbf{\boldsymbol{x}}_*(t) = \frac{1}{1-t}. \] est maximale mais pas globale.
Théorème : existence et unicité
Supposons que
La fonction \(\mathbf{\boldsymbol{f}}\) est continue ;
Pour tout ensemble borné \(D \subset \mathbf R\times \mathbf R^n\) contenant \((t_0, \mathbf{\boldsymbol{x}}_0)\), il existe \(L(D)\) tel que \[ \forall \bigl((t, \mathbf{\boldsymbol{x}}_1), (t, \mathbf{\boldsymbol{x}}_2)\bigr) \in D \times D, \qquad {\lVert { \mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}_1) - \mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}_2) } \rVert} \leqslant L {\lVert {\mathbf{\boldsymbol{x}}_1 - \mathbf{\boldsymbol{x}}_2} \rVert}. \]
Alors
Il existe une unique solution maximale \(\mathbf{\boldsymbol{x}}_* \colon (T_-, T_+) \to \mathbf R^n\) ;
Explosion : Si \(T_+\) est fini, alors \(\lim_{t \to T_+} {\lVert {\mathbf{\boldsymbol{x}}(t)} \rVert} = \infty\), et de même pour \(T_-\).
\(~\)
Dorénavant, \(t_0 = 0\) pour simplifier.
Problème jouet
Problème jouet
Pour l’illustration numérique, on considère l’équation \[ \left\{ \begin{aligned} x'(t) &= \alpha \bigl(x(t) - \sin(t) \bigr) + \cos(t) \\ x(0) &= x_0. \end{aligned} \right. \]
La solution exacte est donnée par
\[x(t) = \sin(t) + e^{\alpha t} x_0.\]
Les solutions convergent si \(\alpha < 0\), et divergent si \(\alpha > 0\).
Problème jouet
Problème jouet
Pour l’illustration numérique, on considère l’équation \[ \left\{ \begin{aligned} x'(t) &= \alpha \bigl(x(t) - \sin(t) \bigr) + \cos(t) \\ x(0) &= x_0. \end{aligned} \right. \]
La solution exacte est donnée par
\[x(t) = \sin(t) + e^{\alpha t} x_0.\]
Les solutions convergent si \(\alpha < 0\), et divergent si \(\alpha > 0\).
Problème jouet
Problème jouet
Pour l’illustration numérique, on considère l’équation \[ \left\{ \begin{aligned} x'(t) &= \alpha \bigl(x(t) - \sin(t) \bigr) + \cos(t) \\ x(0) &= x_0. \end{aligned} \right. \]
La solution exacte est donnée par
\[x(t) = \sin(t) + e^{\alpha t} x_0.\]
Les solutions convergent si \(\alpha < 0\), et divergent si \(\alpha > 0\).
Les méthodes d’Euler explicite et implicite
La méthode d’Euler explicite
La méthode
La méthode d’Euler explicite est l’itération \[ \mathbf{\boldsymbol{x}}_{n+1} = \mathbf{\boldsymbol{x}}_{n} + \Delta \mathbf{\boldsymbol{f}}(t_n, \mathbf{\boldsymbol{x}}_n), \] avec \(\Delta\) le pas de temps.
Remarque : lien avec la quadrature numérique
Lorsque \(\mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{f}}(t)\), la méthode effectue l’approximation \[ \int_0^{n \Delta} \mathbf{\boldsymbol{f}}(t) \, \mathrm dt \approx \sum_{i=0}^{n-1} \Delta \mathbf{\boldsymbol{f}}(i \Delta). \] \({\color{green} \rightarrow}\) C’est la règle du point gauche pour l’intégration numérique.
\({\color{green} \rightarrow}\) L’erreur est en \(\mathcal O(\Delta)\) dans ce cas.
Question : est-ce vrai en général ?
La méthode d’Euler explicite : interprétation géométrique
Problème exemple
On considère le problème jouet du début : \[ \left\{ \begin{aligned} x'(t) &= \alpha \bigl(x(t) - \sin(t) \bigr) + \cos(t) \\ x(0) &= 0, \end{aligned} \right. \] avec \(\alpha = 0.15\).
La méthode d’Euler explicite : illustration de la convergence
Estimation numérique de la vitesse de convergence
Estimation de l’ordre de convergence
Question : comment estimer l’ordre \(p\) ? \[ \max_n |x(t_n) - x_n| = \mathcal O(\Delta^{\color{blue} p}) \]
Réponse : supposons que \[ u(Δ) = C \Delta^{\color{blue} p}. \] Alors \[ \log\bigl(u(\Delta)\bigr) = \log C + {\color{blue} p} \log(\Delta) \] Autrement dit, \(\log(u)\) est une fonction affine de \(\log(\Delta)\) !
\({\color{green} \rightarrow}\) Permet de calculer \(p\) par régression linéaire.
La méthode d’Euler explicite : analyse d’erreur
Théorème de convergence
Supposons que
il existe une solution unique \(\mathbf{\boldsymbol{x}}(t) \colon [0, T] \to \mathbf R^n\), de classe \(C^2\) avec \[ M := \sup_{t \in [0,T]} {\lVert {\mathbf{\boldsymbol{x}}''(t)} \rVert}. \]
la fonction \(\mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}})\) est globalement lipschitzienne en \(\mathbf{\boldsymbol{x}}\) de constante \(L\).
Alors l’estimation d’erreur suivante est vérifiée : \[\begin{equation} \forall n, \quad {\lVert {\mathbf{\boldsymbol{x}}(t_n) - \mathbf{\boldsymbol{x}}_n} \rVert} \leqslant \frac{{\color{blue} \Delta} M}{2} \left( \frac{\mathop{\mathrm{e}}^{L t_n} - 1}{L} \right), \end{equation}\]
Preuve (dimension 1 pour simplifier)
Soit \(e_n = x(t_n) - x_n\). Par Taylor, il existe \(\tau_n \in [t_n, t_{n+1}]\) tel que \[ e_n = e_{n-1} + \Delta \Bigl( f\bigl(t_{n-1}, x(t_{n-1})\bigr) - f\bigl(t_{n-1}, x_{n-1}\bigr) \Bigr) + \frac{\Delta^2}{2} x''(\tau_n). \]
Par la continuité lipschitzienne, \[ {\lvert {e_{n}} \rvert} \leqslant(1 + \Delta L) {\lvert {e_{n-1}} \rvert} + \frac{M \Delta^2}{2}. \]
En itérant cette inégalité, on déduit \[ \begin{aligned} |e_{n}| &\leqslant(1 + \Delta L) \Bigl( (1 + \Delta L) {\lvert {e_{n-2}} \rvert} + \frac{M \Delta^2}{2} \Bigr) + \frac{M \Delta^2}{2}\\ &\leqslant\dotsc \leqslant(1 + \Delta L)^{n} {\lvert {e_{0}} \rvert} + \sum_{i=1}^{n} (1 + \Delta L)^{n-i} \frac{M \Delta^2}{2}. \end{aligned} \]
En utilisant la formule des séries géométriques, on a \[ {\lvert {e_{n}} \rvert} \leqslant(1 + \Delta L)^{n} {\lvert {e_{0}} \rvert} + \frac{(1+\Delta L)^{n} - 1}{\Delta L} \left( \frac{\Delta^2 M}{2} \right). \]
La méthode d’Euler explicite : instabilité numérique
Définition : instabilité numérique
Une méthode est numériquement instable lorsque la solution numérique diverge mais pas la solution exacte.
Problème exemple
Considérons le problème jouet du début \[ \left\{ \begin{aligned} x'(t) &= f(t, x) := \alpha \bigl(x(t) - \sin(t) \bigr) + \cos(t) \\ x(0) &= x_0. \end{aligned} \right. \]
Pour \(\alpha < 0\), les solutions exactes convergent avec l’échelle de temps \(|\alpha|^{-1}\).
Pour \(\alpha \ll -1\), présence de deux échelles de temps très séparées.
\({\color{green} \rightarrow}\) échelle courte \(|\alpha|^{-1}\), et échelle longue \(2\pi\) (période de \(\sin\))
\({\color{green} \rightarrow}\) l’équation est dite raide.
On observe une instabilité numérique si \(|1 + \Delta \alpha| > 1\). Pourquoi ?
\({\color{green} \rightarrow}\) en général, la stabilité numérique dépend de l’échelle la plus courte.
\(\alpha = -10\), \(\Delta = .18\)
\(~\)
La méthode d’Euler explicite : instabilité numérique
Définition : instabilité numérique
Une méthode est numériquement instable lorsque la solution numérique diverge mais pas la solution exacte.
Problème exemple
Considérons le problème jouet du début \[ \left\{ \begin{aligned} x'(t) &= f(t, x) := \alpha \bigl(x(t) - \sin(t) \bigr) + \cos(t) \\ x(0) &= x_0. \end{aligned} \right. \]
Pour \(\alpha < 0\), les solutions exactes convergent avec l’échelle de temps \(|\alpha|^{-1}\).
Pour \(\alpha \ll -1\), présence de deux échelles de temps très séparées.
\({\color{green} \rightarrow}\) échelle courte \(|\alpha|^{-1}\), et échelle longue \(2\pi\) (période de \(\sin\))
\({\color{green} \rightarrow}\) l’équation est dite raide.
On observe une instabilité numérique si \(|1 + \Delta \alpha| > 1\). Pourquoi ?
\({\color{green} \rightarrow}\) en général, la stabilité numérique dépend de l’échelle la plus courte.
\(\alpha = -10\), \(\Delta = .22\)
\(~\)
Stabilité absolue
Définition : stabilité absolue
Considérons l’équation modèle \[ \left\{ \begin{aligned} x'(t) &= \lambda x(t), \\ x(0) &= 1. \end{aligned} \right. \] Un schéma numérique est dit absolument stable si \[ \lvert x_n \rvert \xrightarrow[n \to \infty]{} 0. \] La région de stabilité absolue est par définition \[ \mathcal A := \{ z \in \mathbf C: \text{${\lvert {x_n} \rvert} \to 0$ lorsque $\Delta \lambda = z$} \} \subset \mathbf C. \]
Stabilité absolue de la méthode d’Euler explicite
L’itération d’Euler explicite pour l’équation modèle s’écrit \[ x_{n+1} = x_n + \Delta \lambda x_n = (1 + \Delta \lambda) x_n. \] Le schéma est absolument stable ssi \(|1 + \Delta \lambda| < 1\).
xlims, ylims = (-5, 1), (-2, 2)
x = LinRange(xlims..., 200)
y = LinRange(ylims..., 200)
stable = @. clamp(abs(1 + x' + im*y), 0, 2)
Plots.contourf(x, y, stable, c=cgrad(:Pastel1_3, rev=true),
aspect_ratio=:equal, levels=[0., 1., 2.],
xlabel=L"\Re(\Delta \lambda)",
ylabel=L"\Im(\Delta \lambda)",
xlims=xlims, ylims=ylims, bg=:transparent)
Plots.contour!(x, y, stable, color=:red,
levels=[0., 1., 2.], colorbar=:none)
Plots.vline!([0], color=:gray)
Plots.hline!([0], color=:gray)
Stabilité absolue : illustration pour l’équation modèle
Problème exemple
Considérons l’équation \[ \left\{ \begin{aligned} x'(t) &= - x(t), \\ x(0) &= 1. \end{aligned} \right. \] Stabilité absolue d’Euler explicite si \(\Delta < 2\).
\(\Delta = \frac{1}{2}\)
\(~\)
Stabilité absolue : illustration pour l’équation modèle
Problème exemple
Considérons l’équation \[ \left\{ \begin{aligned} x'(t) &= - x(t), \\ x(0) &= 1. \end{aligned} \right. \] Stabilité absolue d’Euler explicite si \(\Delta < 2\).
\(\Delta = 2\)
\(~\)
Stabilité absolue : illustration pour l’équation modèle
Problème exemple
Considérons l’équation \[ \left\{ \begin{aligned} x'(t) &= - x(t), \\ x(0) &= 1. \end{aligned} \right. \] Stabilité absolue d’Euler explicite si \(\Delta < 2\).
\(\Delta = 2.5\)
\(~\)
Application principale de la stabilité absolue : EDP
Exemple : équation de la chaleur
Discrétisation spatiale de l’équation de la chaleur \[ \left \{ \begin{aligned} &\partial_t u = \partial_x^2 u, \\ &u(t, 0) = u(t, 1) = 0, \end{aligned} \right. \] par la méthode des différences finies conduit à l’équation différentielle \[ \mathbf{\boldsymbol{u}}' := \mathsf A \mathbf{\boldsymbol{u}}, \qquad \mathsf A := \frac{1}{h^2} \begin{pmatrix} -2 & 1 \\ 1 & -2 & 1 \\ & 1 & -2 & 1 \\ & & \ddots & \ddots & \ddots & \\ & & & 1 & -2 & 1 \\ & & & & 1 & -2 \\ \end{pmatrix}, \] où \(h\) est le pas spatial. La méthode d’Euler explicite pour ce système s’écrit \[ \mathbf{\boldsymbol{u}}_{n+1} = \mathbf{\boldsymbol{u}}_{n} + \Delta \mathsf A \mathbf{\boldsymbol{u}}_{n} = (\mathsf I + \Delta \mathsf A) \mathbf{\boldsymbol{u}}_{n}. \] \({\color{green} \rightarrow}\) stabilité numérique si \(\Delta \lambda \in \mathcal A\), c.-à-d. \(|1 + \Delta \lambda|< 1\), pour toute valeur propre \(\lambda\) de \(\mathsf A\).
\({\color{green} \rightarrow}\) puisque \(\max_i |\lambda_i| \propto \frac{1}{h^2}\), la stabilité requiert \(\Delta = \mathcal O(h^2)\)…
La méthode d’Euler implicite
La méthode
La méthode d’Euler implicite est l’itération \[ \mathbf{\boldsymbol{x}}_{n+1} = \mathbf{\boldsymbol{x}}_{n} + \Delta \mathbf{\boldsymbol{f}}({\color{red}t_{n+1}}, {\color{red}\mathbf{\boldsymbol{x}}_{n+1}}), \] avec \(\Delta\) le pas de temps.
\({\color{green} \rightarrow}\) Équation non linéaire à résoudre à chaque étape !
\({\color{green} \rightarrow}\) Une itération à point fixe est l’approche la plus simple pour résoudre cette équation : \[ {\color{magenta}\mathbf{\boldsymbol{x}}_{n+1}^{(k+1)}} = \mathbf{\boldsymbol{x}}_{n} + \Delta \mathbf{\boldsymbol{f}}(t_n, {\color{magenta}\mathbf{\boldsymbol{x}}_{n+1}^{(k)}}). \]
Remarque : lien avec la quadrature numérique
Lorsque \(\mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{f}}(t)\), la méthode effectue l’approximation \[ \int_0^{n \Delta} \mathbf{\boldsymbol{f}}(t) \, \mathrm dt \approx \sum_{i=1}^{n} \Delta \mathbf{\boldsymbol{f}}(i \Delta). \] \({\color{green} \rightarrow}\) C’est la règle du point droit pour l’intégration numérique.
\({\color{green} \rightarrow}\) L’erreur est en \(\mathcal O(\Delta)\) dans ce cas.
Proposition : stabilité absolue
La région de stabilité de la méthode d’Euler implicite est \[ \mathcal A = \left\{ z \in \mathbf C: \frac{1}{|1 - \Delta \lambda|} < 1 \right\} \]
\({\color{green} \rightarrow}\) Stable pour tout \(\Re(\Delta \lambda) < 0\), «stabilité inconditionnelle»
Stabilité absolue : illustration pour l’équation modèle
Problème exemple
Considérons l’équation \[ \left\{ \begin{aligned} x'(t) &= - x(t), \\ x(0) &= 1. \end{aligned} \right. \] Stabilité absolue pour tout \(\Delta > 0\) !
\(\Delta = 2.5\)
\(~\)
Méthodes à un pas d’ordre élevé
Méthodes de Taylor explicites
Idée motivante
En dérivant l’équation \[ \mathbf{\boldsymbol{x}}'(t) = \mathbf{\boldsymbol{f}}\bigl(t, \mathbf{\boldsymbol{x}}(t)\bigr), \] on a accès aux dérivées d’ordre supérieur de \(\mathbf{\boldsymbol{x}}\), par exemple \[ \begin{aligned} \mathbf{\boldsymbol{x}}'' &= \mathbf{\boldsymbol{f}}_2\bigl(t, \mathbf{\boldsymbol{x}}(t)\bigr), \qquad \mathbf{\boldsymbol{f}}_2(t, \mathbf{\boldsymbol{x}}) := (\partial_t + \mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}) \cdot \nabla_{\mathbf{\boldsymbol{x}}}) \mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}) \\ \mathbf{\boldsymbol{x}}''' &= \mathbf{\boldsymbol{f}}_3\bigl(t, \mathbf{\boldsymbol{x}}(t)\bigr), \qquad \mathbf{\boldsymbol{f}}_3(t, \mathbf{\boldsymbol{x}}) := (\partial_t + \mathbf{\boldsymbol{f}}(t, \mathbf{\boldsymbol{x}}) \cdot \nabla_{\mathbf{\boldsymbol{x}}}) \mathbf{\boldsymbol{f}}_2(t, \mathbf{\boldsymbol{x}}), \\ &\dots \end{aligned} \] Ceci motive les soi-disant méthodes de Taylor : \[ \mathbf{\boldsymbol{x}}_{n+1} = \mathbf{\boldsymbol{x}}_n + \Delta \mathbf{\boldsymbol{f}}(t_n, \mathbf{\boldsymbol{x}}_n) + \frac{\Delta^2}{2} \mathbf{\boldsymbol{f}}_2(t_n, \mathbf{\boldsymbol{x}}_n) + \dotsc + \frac{\Delta^p}{p!} \mathbf{\boldsymbol{f}}_p(t_n, \mathbf{\boldsymbol{x}}_n). \]
Remarques
Notez que pour \(p = 1\), c’est la méthode d’Euler explicite ;
L’erreur est en \(\mathcal O(\Delta^p)\) sous des hypothèses appropriées.
Par construction, exacte lorsque \(\mathbf{\boldsymbol{x}}(t)\) est un polynôme de degré \(p\).
Proposition : stabilité absolue
La région de stabilité absolue de la méthode de Taylor est \[ \mathcal A = \left\{ z \in \mathbf C: \left|1 + z + \dotsc + \frac{z^p}{p!} \right| < 1 \right\} \]
Méthodes de Runge-Kutta
Définition : méthode de Runge-Kutta explicite
Une méthode de Runge-Kutta explicite à \(p\) étages est de la forme \[ \begin{aligned} \mathbf{\boldsymbol{x}}_{n+1} &= \mathbf{\boldsymbol{x}}_n + \Delta \sum_{i=1}^p b_i \mathbf{\boldsymbol{k}}_i \\ \mathbf{\boldsymbol{k}}_1 &= \mathbf{\boldsymbol{f}}(t_n, \mathbf{\boldsymbol{x}}_n), \\ \mathbf{\boldsymbol{k}}_2 &= \mathbf{\boldsymbol{f}}\bigl(t_n + c_2 \Delta, \mathbf{\boldsymbol{x}}_n+\Delta(a_{21} \mathbf{\boldsymbol{k}}_1)\bigr), \\ \mathbf{\boldsymbol{k}}_3 &= \mathbf{\boldsymbol{f}}\bigl(t_n + c_3 \Delta, \mathbf{\boldsymbol{x}}_n+\Delta(a_{31} \mathbf{\boldsymbol{k}}_1 + a_{32} \mathbf{\boldsymbol{k}}_2)\bigr), \\ &\;\;\vdots \\ \mathbf{\boldsymbol{k}}_p &= \mathbf{\boldsymbol{f}}\left(t_n + c_p \Delta, \mathbf{\boldsymbol{x}}_n + \Delta \sum_{j = 1}^{p-1} a_{pj} \mathbf{\boldsymbol{k}}_j\right), \end{aligned} \]
\({\color{green} \rightarrow}\) Les dérivées de \(\mathbf{\boldsymbol{f}}\) ne sont pas nécessaires !
Remarques
Construction fastidieuse et non discutée ici ;
Les plus utilisées en pratique ;
L’erreur est en \(\mathcal O(\Delta^p)\) pour des coefficients appropriés ;
Exemple : méthode de Heun
La méthode de Heun est une méthode de Runge-Kutta d’ordre 2 et s’écrit \[ \mathbf{\boldsymbol{x}}_{n+1} = \mathbf{\boldsymbol{x}}_n + \frac{\Delta}{2}\mathbf{\boldsymbol{f}}(t_n, \mathbf{\boldsymbol{x}}_n) + \frac{\Delta}{2} \mathbf{\boldsymbol{f}} \bigl( t_n + \Delta, \mathbf{\boldsymbol{x}}_n + \Delta \mathbf{\boldsymbol{f}}(t_n, \mathbf{\boldsymbol{x}}_n) \bigr). \]
Proposition : même région de stabilité que la méthode de Taylor
La région de stabilité absolue d’une méthode de Runge-Kutta d’ordre \(p\) est \[ \mathcal A = \left\{ z \in \mathbf C: \left|1 + z + \dotsc + \frac{z^p}{p!} \right| < 1 \right\} \]
Pour aller plus loin
Parmi les sujets non abordés ici :
Méthodes multi-pas ;
Méthodes de Runge-Kutta implicites ;
Contrôle adaptatif du pas de temps ;
Schémas symplectiques pour les systèmes hamiltoniens ;
…
