%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[a4paper, 10pt]{article}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}

\title{}
\author{}
\date{}

%%------------------------------------------------------------------------------
%% Packages

\usepackage[]{geometry}
\usepackage{amsmath,amsfonts,amsthm}
\usepackage[matha,mathb,mathx]{mathabx} % french mathematical symbols and more
\usepackage[svgnames]{xcolor}
\usepackage{hyperref}
\usepackage{mathtools}
\usepackage{fancyhdr}
\usepackage{lastpage}
\usepackage{enumitem}

\usepackage[french]{babel}

%%------------------------------------------------------------------------------
%% Packages setup

% geometry
\geometry{textwidth = 0.75\paperwidth,
  textheight = 0.75\paperheight}
%defaults : textwidth = 0.7\paperwidth = 1675.8pt, textheight = 0.7\paperheight = 4788pt

% hyperref
\hypersetup{colorlinks=true,
  citecolor=DarkBlue,
  linkcolor=DarkBlue,
  linktoc=page,
  urlcolor=black}

\fancypagestyle{plain}{%
  \fancyhf{}
  \rfoot{\thepage /\pageref{LastPage}}
  \renewcommand{\headrulewidth}{0pt}}

\pagestyle{fancy} % sets page style for the whole document
\fancyhf{} % clear all header and footer fields
\lhead{École Normale Supérieure de Lyon\\Département de mathématiques}
\rhead{2021--2022\\Option B pour l'agrégation}
\rfoot{\thepage /\pageref{LastPage}}
\renewcommand{\footrulewidth}{0.5pt}
\setlength{\headheight}{24.2pt} % required for headers spanning two lines

%%------------------------------------------------------------------------------
%% Custom environments/macros

\theoremstyle{definition}
\newtheorem{defi}{Définition}
\newtheorem{exem}[defi]{Exemple}

\theoremstyle{remark}
\newtheorem{rmrk}[defi]{Remarque}

\theoremstyle{plain}
\newtheorem{thrm}[defi]{Théorème}
\newtheorem{prop}[defi]{Proposition}
\newtheorem{coro}[defi]{Corollaire}
\newtheorem{lemm}[defi]{Lemme}

\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\Sph}{\mathbb{S}}
\newcommand{\mc}{\mathcal}
\newcommand{\eps}{\varepsilon}

\renewcommand{\Im}{\operatorname{Im}}
\renewcommand{\Re}{\operatorname{Re}}
\newcommand{\Id}{\operatorname{Id}}
\DeclareMathOperator{\Tr}{Tr}

\newcommand{\dbra}[1]{\ldbrack #1 \rdbrack}
\DeclarePairedDelimiter{\norm}{\lVert}{\rVert} %% requires mathtools
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\sbra}{[}{]}
\DeclarePairedDelimiter{\pare}{(}{)}
\DeclarePairedDelimiter{\osbral}{]}{]}
\DeclarePairedDelimiter{\osbra}{]}{[}
\DeclarePairedDelimiter{\osbrar}{[}{[}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{center}
  {\Large \bf Résolution numérique d'EDO : méthodes de Runge-Kutta}
\end{center}
\bigskip

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction : méthode du point milieu}

On s'intéresse à l'équation différentielle
\begin{equation}\label{eq:pb_diff}
  y' = f(t,y),\quad t\in I,\,y\in U\subset\R^n
\end{equation}
d'inconnue \(y:I\to \R^n\). On suppose que le problème de Cauchy de condition initiale \(y(t_0)=y_0\) admet une solution maximale unique, comme par exemple dans l'application du théorème de Cauchy-Lipschitz. On note cette solution \(z\), et on fixe \(T>0\) ainsi qu'une discrétisation \(t_0 < t_1 < \dots < t_N = t_0+T\), où \(T\) est choisi de sorte que \(z\) est définie sur \(\sbra{t_0, t_0+T}\). Alors, pour chaque \(k\in\dbra{0,N-1}\), on peut écrire
\[
z(t_{k+1}) - z(t_k) = \int_{t_k}^{t_{k+1}} f(s,z(s))ds,
\]
et l'idée est d'approximer l'intégrale du membre de droite par une méthode de quadrature. Par exemple, en choisissant la méthode des rectanges à gauche, on obitent
\[
z(t_{k+1}) - z(t_k) \simeq (t_{k+1}-t_k) f(t_k,z(t_k)),
\]
c'est la méthode d'Euler explicite. Plutôt que la méthode des rectangles à gauche, on peut par exemple choisir la méthode du point milieu. En notant \(h_k = t_{k+1}-t_k\),
\[
z(t_{k+1}) - z(t_k) \simeq (t_{k+1}-t_k) f\pare*{t_k+\frac{h_k}{2},z\pare*{t_k+\frac{h_k}{2}}},
\]
mais pour cette méthode, on voit qu'il faut d'abord calculer une valeur approchée de \(z(t_k + h_k/2)\), pour laquelle on peut choisir une autre méthode de quadrature, comme par exemple celle des rectangles à gauche (qui correspond à la méthode d'Euler). Dans ce cas, on obtient la formule suivante :
\[
z(t_{k+1}) - z(t_k) \simeq (t_{k+1}-t_k) f\pare*{t_k+\frac{h_k}{2},z(t_k)+\frac{h_k}{2}f(t_k,z(t_k))}.
\]
Cette formule est celle de la méthode dite du point milieu (aussi appelée Runge-Kutta 2). Notons qu'on a utilisé deux méthodes de quadrature, la première d'ordre 1, et la seconde d'ordre 2, il n'est donc pas évident que la méthode du point milieu soit d'ordre supérieur à 1, mais c'est pourtant le cas, comme un calcul explicite le prouve.
\begin{prop}
  La méthode du point milieu est consistante d'ordre 2 si \(f\) est de classe \(\mc C^2\).
\end{prop}
\begin{proof}
  Comme souvent pour ce genre de preuve, il suffit de faire des développements limités et de regrouper les termes. Remarquons d'abord que la solution exacte \(z\) vérifie
  \[
  z''(t) = \frac{d}{dt}f(t,z(t)) = \frac{\partial f}{\partial t}(t,z(t)) + \nabla_y f(t,z(t)) z'(t) = \pare*{\frac{\partial f}{\partial t} + (\nabla_y f) f}(t,z(t)),
  \]
  où \(\nabla_y\) désigne le gradient par rapport à la variable \(y\). Alors
  \begin{align*}
    f\pare*{t_k+\frac{h_k}{2},z(t_k)+\frac{h_k}{2}f(t_k,z(t_k))}
    &= f(t_k,z(t_k)) + \pare*{\frac{h_k}{2}\frac{\partial f}{\partial t} + \frac{h_k}{2}(\nabla_y f) f}(t_k,z(t_k)) + O(h_k^2)\\
    &= z'(t_k) + \frac{h_k}{2}z''(t_k) + O(h_k^2).
  \end{align*}
  Par conséquent, l'erreur de consistance est donnée par
  \begin{align*}
    e_k &= z(t_{k+1}) - z(t_k) - h_kf\pare*{t_k+\frac{h_k}{2},z(t_k)+\frac{h_k}{2}f(t_k,z(t_k))} \\
    &= h_kz'(t_k) + \frac{h_k^2}{2}z''(t_k) - h_k\pare*{ z'(t_k) + \frac{h_k}{2}z''(t_k)} + O(h_k^3)\\
    &= O(h_k^3),
  \end{align*}
  où le dernier \(O\) peut dépendde \(z^{(3)}\) et \(D^2f\).
\end{proof}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Définition des méthodes Runge-Kutta}

Les méthodes de Runge-Kutta constituent une généralisation de la méthode du point milieu : on fixe une subdivision \(0 = c_0 \leq c_1 \leq c_2 \leq \dots \leq c_q = 1\), et on calcule de manière récursive la valeur approchée en \(t_{k,i} = t_k + c_ih_k\) grâce à une formule de quadrature fixée,
\[
\int_0^{c_i} g(u)du = \sum_{j=0}^{i-1} a_{ij} g(c_j),
\]
donc 
\[
z(t_{k,i})=z(t_k)+\int_{t_k}^{t_{k,i}} f(u,z(u))du \simeq h_k\sum_{j=0}^{i-1} a_{ij}f(t_{k,j},z(t_{k,j})).
\]
Notons bien que pour que la méthode reste explicite, il est nécessaire de choisir une méthode de quadrature ne faisant pas intervenir l'indice \(i=j\). La méthode est donc entièrement déterminée par le choix des points \(c_i\) et des coefficients \(a_{ij}\), \(i\in\dbra{1,q}\), \(j\in\dbra{0,i-1}\). On la résume conventionnellement dans un tableau.
\begin{figure}[h]
  \centering
  \begin{minipage}{0.4\textwidth}
    \[
    \begin{array}{c|cccc}
      0=c_0 & & & & \\
      c_1 & a_{1,0} & & & \\
      c_2 & a_{2,0} & a_{2,1} & &  \\
      \vdots & \vdots & \vdots & \ddots & \\ \hline
      1=c_q & a_{q,0} & a_{q,1} & \cdots & a_{q,q-1}
    \end{array}
    \]
  \end{minipage}
  \begin{minipage}{0.45\textwidth}
    \[
    \left\{\begin{aligned}
      p_0 &= f(t_k,y_k) \\
      p_1 &= f(t_k + c_1h_k,\, y_k + h_k a_{1,0}p_0)\\
      \vdots &\\
      p_{q-1} &= f\Big(t_k + c_{q-1}h_k,\, y_k + h_k \sum_{i=0}^{q-2} a_{q-1,i}p_i\Big) \\
      y_{k+1} &= y_k + h_k \sum_{i=0}^{q-1} a_{q,i} p_i
    \end{aligned}\right.
    \]
  \end{minipage}
  \caption{Méthode de Runge-Kutta : tableau de coefficients et algorithme correspondant}
\end{figure}

\begin{exem}
  Le tableau des méthodes d'Euler explicite et du point milieu s'écrivent alors respectivement
  \[
  \begin{array}{c|c}
    0 & \\ \hline
    1 & 1 
  \end{array}
  \qquad\text{et}\qquad
  \begin{array}{c|cc}
    0 & & \\
    1/2 & 1/2 & \\ \hline
    1 & 0 & 1 
  \end{array}
  \]
\end{exem}

\begin{rmrk}
  \begin{itemize}
  \item Pour que la méthode soit d'ordre élevé, il est naturel d'imposer que les méthodes de quadratures intervenant soient toutes d'ordre au moins 0 (donc exactes pour les fonctions constantes), ce qui implique que la somme des coefficients \(a_{i,j}\) de chaque ligne \(i\) est égale à \(c_i\).
  \item Il est possible de choisir, pour certains \(i\), \(c_i=c_{i+1}\) : c'est ce qu'on fait par exemple dans la méthode de Runge-Kutta 4 décrite ci-après. Dans ce cas, \(t_{k,i}=t_{k,i+1}\) et la valeur \(y_{k,i+1}\) est une meilleure approximation de \(z(t_{k,i})\) que ne l'est \(y_{k,i}\).
  \end{itemize}
\end{rmrk}

\begin{exem}\label{ex:rk4}
  La méthode de Runge-Kutta 4 (ou Runge-Kutta classique) est donnée par le tableau
  \[
  \begin{array}{c|cccc}
    0 & & & & \\
    1/2 & 1/2 & & & \\
    1/2 & 0 & 1/2 & & \\
    1 & 0 & 0 & 1 &\\ \hline
    1 & 1/6 & 1/3 & 1/3 & 1/6
  \end{array}
  \]
  Les quadratures utilisées sont successivement les rectangles à gauche, puis à droite, puis le point milieu, et enfin la méthode de Simpson (Newton-Cotes de rang 2).  
\end{exem}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Détermination de l'ordre des méthodes}

On rappelle la proposition suivante, sur une méthode numérique donnée par la fonction \(\Phi\), dont on suppose qu'elle est continue.
\begin{prop}
  La méthode est consistante si, et seulement si,  \(\Phi(t,y,0)=f(t,y)\).
\end{prop}
On peut généraliser ce résultat pour établir la consistance d'ordre \(p\in\N\). L'idée est essentiellement toujours la même, on effectue des développements limités dans la formule de l'erreur de consistance. Pour ce faire, on suppose que \(f\) et \(\Phi\) sont de classe \(\mc C^p\), et on définit la notion de dérivée suivante, qu'on peut voir comme une dérivée particulaire.
\begin{defi}
  Soit \(g:\R\times\R^n\to\R^n\) une fonction de classe \(\mc C^1\). On définit
  \[
  Dg(t,y) = \frac{\partial g}{\partial t}(t,y) + (\nabla_y g(t,y))f(t,y).
  \]
  Alors, en posant \(f^{[j]} = D^jf\), la solution de \eqref{eq:pb_diff} vérifie
  \[
  z^{(j+1)}(t) = f^{[j]}(t,z(t))
  \]
  pour tout \(j\in\dbra{0,p}\).
\end{defi}
\begin{thrm}
  Si \(\Phi\) vérifie
  \begin{equation}\label{eq:cons}
    \frac{\partial^j\Phi}{\partial h}\Phi(t,y,0) = \frac{1}{j+1}f^{[j]}(t,y)
  \end{equation}
  pour tout \(j\in\dbra{0,p-1}\), alors la méthode est d'ordre \(p\).
\end{thrm}
\begin{proof}
  On effectue un développement limité pour \(z\) (par rapport à \(t\)) et pour \(\Phi\) (par rapport à \(h\)) dans la définition de l'erreur de consistance
  \[
  e_k = z(t_{k+1}) - z(t_k) - h_k\Phi(t_k,z(t_k),h_k).
  \]
  D'une part,
  \begin{align*}
    z(t_{k+1}) - z(t_k) &= \sum_{j=1}^{p+1} \frac{h_k^j}{j!}z^{(j)}(t_k) + o(h_k^{p+1})\\
     &= \sum_{j=0}^p \frac{h_k^{j+1}}{(j+1)!}f^{[j]}(t_k,z(t_k)) +  o(h_k^{p+1}),
  \end{align*}
  et d'autre part,
  \begin{align*}
    \Phi(t_k,z(t_k),h_k) = \sum_{j=0}^p\frac{h_k^j}{j!}\frac{\partial^j\Phi}{\partial h^j}(t_k,z(t_k),0) + o(h_k^{p}).
  \end{align*}
  Ainsi, si la relation~\eqref{eq:cons} est vérifiée, on obtient que
  \[
  \abs{e_k} = o(h_k^{p+1}),
  \]
  d'où, en sommant sur \(k\), 
  \[
  \sum_{k=0}^{N-1}\abs{e_k} = O(h_{max}^p).
  \]
\end{proof}

En application directe de ce théorème, on peut montrer que la méthode de Runge-Kutta 4 de l'exemple~\ref{ex:rk4} est d'ordre 4. De manière plus générale, on peut établir des conditions sur les paramètres \(a_{ij}\) et \(c_j\) pour que la méthode de Runge-Kutta correspondante soit d'ordre \(p\), dont certaines sont listées dans le livre de Demailly (p.243). 




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
