%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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éthode d'Euler implicite}
\end{center}
\bigskip

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


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\). Comme précédemment, on suppose que le problème de Cauchy de condition initiale \(y(t_0)=y_0\) admet une solution maximale unique notée \(z\).
Plutôt que d'approximer la différence
\[
z(t_{k+1}) - z(t_k) = \int_{t_k}^{t_{k+1}} f(t,z(t))dt
\]
par la méthode des rectangles à gauche, qui amène à la méthode d'Euler explicite, on peut utiliser la méthode des rectangles à droite, conduisant à la formule suivante
\begin{equation}\label{eq:euler_implicite}
  y_{k+1} = y_k + h_kf(t_{k+1},y_{k+1}).
\end{equation}
C'est la méthode d'Euler dite \emph{implicite}, car elle nécessite la résolution d'une équation non linéaire à chaque étape, par exemple en utilisant l'algorithme de Newton. Elle est en conséquence plus coûteuse que la méthode d'Euler explicite. Les définitions de consistance, de stabilité et de convergence s'étendent naturellement aux méthodes implicites, et on montre, de la même manière que pour Euler explicite, que la méthode implicite est aussi d'ordre 1.
\begin{prop}
  Si \(f\) est de classe \(\mc C^1\), alors la méthode d'Euler implicite est consistante d'ordre 1.
\end{prop}
L'avantage de cette méthode n'est donc pas un ordre élevé, mais, comme on va le voir, c'est sa stabilité pour certaines équations particulières, appelées parfois équation différentielle \emph{raide}. Pour illustrer cette notion avec un exemple, considérons l'équation différentielle très simple
\[
y' = -\lambda y,
\]
avec \(\lambda >0\). Les solutions exactes sont les \(y(t)=y_0e^{-\lambda t}\), qui sont des fonctions qui décroissent rapidement vers 0 lorsque \(t\) croît. Étant donné un pas de temps \(h\), la solution approchée donnée par la méthode d'Euler explicite est
\[
y_{k+1} = y_k -h\lambda y_k,
\]
donc
\[
y_k = y_0(1-h\lambda)^k,
\]
dont l'allure est similaire à une solution exacte que si \(\abs{1-h\lambda}<1\), et donc \(0< h < 2/\lambda\). Si \(\lambda\) est grand, cela implique que le pas \(h\) doit être très petit, ce qui n'est pas toujours possible pour une application pratique. La méthode d'Euler implicite, quant à elle, donne la solution approchée
\[
y_k= \frac{y_0}{(1+h\lambda)^k},
\]
qui converge vers \(0\) inconditionnellement, ce qui est numériquement préférable. Ce genre de comportements très différents entre les schémas implicite et explicite survient typiquement dans des problèmes physiques à plusieurs échelles, comme par exemple plusieurs réactions chimiques simultanées dont certaines sont lentes et d'autres rapides. Dans ce genre de cas, on préférera des méthodes implicites.

Citons pour finir un cas particulier de système où la méthode d'Euler implicite est particulièrement stable.
\begin{defi}
  On dit que le système est dissipatif si pour tout \(t\geq 0\) et \(z_1,z_2\in \R^n\),
  \[
  (f(t,z_1) - f(t,z_2),\, z_1-z_2) \geq 0.
  \]
  De manière équivalente, si \(f\) est de classe \(\mc C^1\),
  \[
  ( \nabla_y f(t,z_1)z_2,\, z_2) \geq 0.
  \]
\end{defi}
Si \(z_1\) et \(z_2\) sont deux solutions de l'équation différentielle~\eqref{eq:pb_diff}, alors
\[
\frac{d}{dt}\abs{z_1-z_2}^2 = 2(f(t,z_1) - f(t,z_2),\, z_1-z_2) \leq 0.
\]
En particulier, tous les points d'équilibre du système sont stables. C'est une condition assez forte sur \(f\), qui n'est vérifiée que dans certains cas particuliers. On peut penser à des systèmes physiques où l'énergie décroît -- elle est dissipée --, comme par exemple le mouvement d'une particule chargée de vitesse \(v\), soumise à un champ électrique \(E\) et un champ magnétique \(B\) constants, ainsi qu'une force de frottement \(-kv\), où \(k \geq 0\). Si l'on normalise sa masse et sa charge, l'équation du mouvement s'écrit
\[
\frac{dv}{dt} = (E + v\times B) - kv,
\]
et alors
\[
 (f(v_1) - f(v_2),\, v_1-v_2) = -k\abs{v_1-v_2}^2 \leq 0.
\]
On dispose dans ce cas particulier du théorème suivant.
\begin{thrm}
  On suppose que \(f\) est de classe \(\mc C^1\), et que le système est dissipatif. Alors 
  \begin{enumerate}
  \item pour tout pas \(h_{max}>0\), la suite \((y_k)_{k\geq 0}\) définie par la méthode d'Euler implicite~\eqref{eq:euler_implicite} est bien définie, et
  \item la méthode est consistante d'ordre 1, et stable avec la constante de stabilité \(S=1\).
  \end{enumerate}
\end{thrm}
\begin{proof}
  1. Commençons par montrer la bonne définition de la suite \((y_k)\). Étant donné \(y_k\), on introduit la fonction \(G(h,z) = z - y_k - hf(t_k+h,z)\). Trouver \(y_{k+1}\) revient alors à démontrer l'existence d'un point \(z\) tel que \(G(h_k,z)=0\). Notons que l'unicité n'est pas obligatoire, mais elle est en tout les cas garantie par le caractère dissipatif du système : on montre en fait que pour tout \(h\geq 0\), \(z\mapsto G(h,z)\) est injective. Si \(z_1,z_2\in\R^n\),
  \begin{align*}
    G(h,z_1) = G(h,z_2) \implies& z_1-z_2 = h(f(t_k+h,z_1)-f(t_k+h,z_2)) \\
    \implies& \abs{z_1-z_2}^2 = h \pare*{f(t_k+h,z_1)-f(t_k+h,z_2),\, z_1-z_2} \leq 0\\
    \implies& z_1=z_2
  \end{align*}
  Soit \(E = \{h^*\geq 0\text{ tel que } \forall h\in\sbra{0,h^*},\exists z\in\R^n,G(h,z)=0\}\) On va montrer que \(E=\R_+\), ce qui prouvera la bonne définition de \((y_k)\). Comme \(G(0,y_k) = 0\), \(E\) est non vide. Par ailleurs, \(G\) est de classe \(\mc C_1\), et
  \[
  \nabla_y G(h,z) = I_n - h\nabla_y f(t_k+h,z),
  \]
  où \(\nabla_y\) désigne le gradient par rapport à la deuxième variable. Ainsi,
  \begin{align*}
    \pare*{G(h,z_1)z_2,\, z_2} &= \abs{z_2}^2 -h\pare*{ \nabla_y f(t_k+h,z_1)z_2,\, z_2} \\
    &\geq \abs{z_2}^2,
  \end{align*}
  et donc \(\nabla_y G(h,z)\) est inversible pour tout \((h,z)\). Par conséquent, le théorème des fonctions implicites assure que si \((h,z)\) vérifie \(G(h,z)=0\), il existe des voisinages ouverts \(V\subset\osbrar{0,+\infty}\) de \(h\) et \(\Omega\) de \((h,z)\), ainsi qu'une fonction \(\phi\) de classe \(\mc C^1\) définie sur \(V\) tels que
  \[
  \pare*{(h,z)\in \Omega \text{ et }G(h,z)=0} \iff  \pare*{h\in V\text{ et }z=\phi(h)}.
  \]
  En particulier, \(E\) est ouvert dans \(\osbrar{0,+\infty}\). Supposons maintenant que \(\sup E = \bar h^* < +\infty\). Par définition, il existe une suite \((h_p)\) croissante et convergeante vers \(h^*\), ainsi qu'une suite \((z_p) \subset\R^n\) telle que \( G(h_p,z_p) = 0 \). Mais alors
  \begin{align*}
    z_p &= y_k + h_pf(t_k+h_p,z_p) \\
    &= y_k + h_p\big(f(t_k+h_p,z_p) - f(t_k+h_p,y_k)\big) + h_pf(t_k+h_p,y_k)
  \end{align*}
  et donc, par dissipativité,
  \begin{align*}
    \abs{z_p-y_k}^2 &\leq h_p( f(t_k+h_p,y_k),\, z_p-y_k) \\
    &\leq h_p\abs{f(t_k+h_p,y_k)} \abs{z_p-y_k},
  \end{align*}
  et on en déduit que la suite \((z_p-y_k)\) est bornée. On conclut en extrayant une sous-suite convergente, et la continuité de \(f\) implique que \(h_*\in E\), et donc que \(E\) est fermé, ce qui montre que \(E\) est nécessairement égal à \(\R_+\) tout entier.

  2. On a déjà montré que la méthode était consistante, il reste donc à montrer qu'elle est stable. Soit \((\eps_k)\) une suite de points dans \(\R^n\), et soient \((y_k)\) et \((\tilde y_k)\) les suites définies par
  \[
  \left\{
  \begin{aligned}
    y_{k+1} &= y_k + h_kf(t_{k+1},y_{k+1}) \\
    \tilde y_{k+1} &= \tilde y_k + h_kf(t_{k+1},\tilde y_{k+1})  + \eps_k.
  \end{aligned}\right.
  \]
  Notons que la suite \((\tilde y_k)\) est bien définie, car les résultats du premier point s'appliquent. On obtient immédiatement
  \begin{align*}
    \abs{y_{k+1} - \tilde y_{k+1}}^2 &\leq (y_k-\tilde y_k-\eps_kn y_{k+1}-\tilde y_{k+1})\\
    &\leq \pare*{\abs{y_k - \tilde y_k}+\abs{\eps_k}}\abs{y_{k+1}-\tilde y_{k+1}},
  \end{align*}
  et donc \( \abs{y_{k+1} - \tilde y_{k+1}} \leq \abs{y_k - \tilde y_k} + \abs{\eps_k}\). Finalement, avec une récurrence immédiate, on en conclut que
  \[
  \max_{0\leq k\leq N-1} \abs{y_k-\tilde y_k} \leq \abs{y_0-\tilde y_0} + \sum_{k=0}^{N-1}\abs{\eps_k},
  \]
  ce qui prouve la stabilité avec constante \(1\).
\end{proof}



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