cond_initiales.tex

\section{Problèmes aux conditions initiales}

\begin{frame}[t]
  \frametitle{Motivation}
  \framesubtitle{Le modèle Lorenz 63}
  
  \begin{columns}
    \column{6cm}
     
  \pdfpcmovie{\includegraphics[width=6cm]{Images/lorenz63.jpg}}{Animations/lorenz_attractor.mp4}

    \column{6cm}

\begin{block}
  {Le modèle de Saltzman-Lorenz}

Contexte : modèle de la convection atmosphérique, fortement tronqué
\begin{eqnarray*}
\frac{\mathop{d}x}{\mathop{d}t} &=& \sigma ( y - x ) \\
\frac{\mathop{d}y}{\mathop{d}t} &=& x ( \rho - z ) -  y \\
\frac{\mathop{d}z}{\mathop{d}t} &=& x y - \beta z
\end{eqnarray*}

\end{block}
  \end{columns}

\end{frame}

\begin{frame}
\frametitle{Exemples}
\begin{exampleblock}{Décroissance radioactive}
$\ddt{x} = -\lambda x$, \quad avec $x(t_0)=x_0$. \newline
Posons $\hat{x}=x/x_0$ et $\hat{t}=(t-t_0)\lambda$, on a 
\begin{equation}
\dd{\hat{x}}{\hat{t}} = - \hat {x} \quad\text{avec}\quad \hat{x}(0)=1.
\label{eq:rad}
\end{equation}
\end{exampleblock}
\begin{exampleblock}{Equation logistique}
\begin{equation}
\ddt{y}= \lambda y ( 1-y)
\end{equation}
\end{exampleblock}
\begin{exampleblock}{Oscillateur harmonique}
\begin{equation}
\ddt{p}=-\pd{H}{q} \, ; \, \ddt{q}=\pd{H}{p} \quad\text{avec\ } H=\frac{1}{2}(p^2+q^2)
\end{equation}
et $p(0)=p_0$ et $q(0)=q_0$
\end{exampleblock}
\end{frame}
\begin{frame} \frametitle{4 étapes pour un schéma en différencies finies}
De façon générale, le système s'écrit:
\begin{equation}
\ddt{\vec{x}} = \vec{F}(\vec{x}, t)
\end{equation}
\begin{enumerate}
\item Définir les points de discrétisation dans le temps
\item Exprimer le système sous la forme d'une différence finie
\item Laisser tomber l'erreur de troncature
\item Étudier les propriétés : stabilité, convergence, et autres (voir plus loin)
\end{enumerate}
\end{frame}
\begin{frame}\frametitle{1. Discrétisation dans le temps}
\begin{tikzpicture}
 \draw[->](-1, 0) -- ( 7, 0);
 \foreach \i in {0,1,2,3,4,5,6} 
{ \node  at (\i, 0) {$\bullet$};  } 
\foreach \i in {0,1, 2} { \draw (\i, -0.5) node {$t_\i$}; } 
\foreach \i in {6}{ \draw (\i, -0.3) node {$t_n$}; }
\end{tikzpicture}

Dans la suite de ce cours, nous allons poser : 
\begin{itemize}
\item $t_i$ :  le temps discrétisé \\
\item $x(t_i)$ :  la solution exacte au temps $t_i$ \\
\item $x_i$ :  la solution du système en diff. finies au temps $t_i$ \\
\end{itemize}
\begin{example}
\begin{align*}
t_0 &= t(0)  & 
t_i &= t_0 + i\dt
\end{align*}
\end{example}
\end{frame}

\begin{frame}
\frametitle {2. Exprimer le système sous forme de différence finie}
Posons : $F_i := F(x_i, t_i)$. Le problème est d'établir une relation itérative entre $x_i$ et $x_{i+1}$ en utilisant uniquement de l'information disponible aux $t_i$. On 
recourt au développement de Taylor. 
\medskip
On peut faire:
\begin{align}
x(t_{i+1}) &= x(t_i) + F_i k + k\cdot \tau , & \quad \tau &= \dddt {x} (x(t_i),t_i) {k \over 2} + \mathcal{O}(k^2) \label{eq:d1} 
\end{align} ou bien:
\begin{align}
x(t_{i})   &= x(t_{i+1}) - F_{i+1} k + k\cdot\tau, &\quad \tau &= \dddt {x} (x(t_{i+1}),t_{i+1}) {k\over 2} + \mathcal{O}(k^2) \label{eq:d2}
\end{align}

\end{frame}

\begin{frame}
\frametitle{\ldots et laisser tomber les termes de troncature}

On crée le système aux différences finies en laissant tomber les termes de troncature:

\begin{definition}
Le schéma est dit \textit{consistent} si l'itération converge vers l'expression
du système originel à la limite $k\rightarrow 0$. 
\end{definition}

\begin{center}
\colorbox{yellow!20}
{ \parbox{\textwidth}
{
\begin{tabular}{llll} 
\hline
Euler avant     & (\ref{eq:d1}) & $x_{i+1} = x_i + F_i \dt$ &  $\tau = \mathcal{O}(k)$ \\ 
Euler arrière   & (\ref{eq:d2}) & $x_{i+1} = x_i + F_{i+1} \dt$ & $\tau = \mathcal{O}(k)$ \\ 
Euler centré   & (\eqref{eq:d1} + \eqref{eq:d2}) / 2 &  $x_{i+1} = x_i + F_{i+i/2} k$ & $\tau = \mathcal{O}(k^2)$  \\
\hline
\end{tabular}
où $F_{i+i/2} := {1 \over 2} (F_{i}  + F_{i+1} ) $.
}}
\end{center}

\end{frame}

\begin{frame} \frametitle{Example : décroissance radioactive}

\begin{block}{Euler avant} 
$x_{i+1} = x_{i} - \dt x_{i} = (1-k) x_{i} \Longrightarrow \color{red}{x_i = (1-\dt)^i}x_0$
\end{block}

\begin{block}{Euler arrière} 
$x_{i+1} = x_{i} - \dt x_{i+1}$ : schéma implicite !

$x_{i+1} = (1+\dt)^{-1} x_{i} \Longrightarrow \color{red}{ x_i = (1+\dt)^{-i}}x_0$
\end{block}

\begin{block}{Euler centré} 
$x_{i+1} = x_{i} - {\dt \over 2}(x_{i}+x_{i+1})$ : schéma (semi-implicite)

$x_{i+1} = \frac{1-k/2}{1+k/2} x_{i} \quad \Longrightarrow \quad \color{red}{x_i = \left(\frac{1-k/2}{1+k/2}\right)^i}x_0$
\end{block}
\end{frame}

\begin{frame}
\includegraphics{radio}
\end{frame}

\begin{frame}
\includegraphics{error}
\end{frame}


\begin{frame} \frametitle{Etude des propriétés. I. Convergence}
\begin{block}{Definition} {$\forall t>t_0\, : \, \lim_{k\rightarrow 0} x_i = x(t_i)$} \end{block}
\begin{exampleblock}{Euler avant, radioactivité : (rappel: $x(t_i) = e^{-t}$)}
\begin{align*}
  \log x_{i} &= i \log(1-k) = i ( - k - k^2/2 + \mathcal{O}(k^3) ) = - t - tk/2 + \mathcal{O}(k^3),  \\
 \Rightarrow x_{i} &= e^{-t (1+k/2)} + \mathcal{O}(k^2).
\end{align*}
\end{exampleblock}
\begin{exercice}
Dans cet exemple, l'Euler avant sous-estime la quantité de produit radioactif. Montrer que l'Euler arrière le surestime, et que l'Euler centré est correct au premier ordre.
\end{exercice}
\end{frame}
\begin{frame}
\frametitle{Ordre de grandeur d'approximation du schéma numérique}
\begin{block}{Règle générale}
  Si, étant connu $x_i$, l'approximation sur $x_{i+1}$ est de l'ordre de k$\tau$, alors l'erreur sur la solution à un temps $t$ sera (si $t$ n'est pas trop grand) de l'ordre de $T\cdot\tau$. 
\end{block}
En effet: à chaque pas de temps on accumule une erreur de l'ordre de $k \tau $, et $T/k$ pas de temps sont nécessaires pour intégrer sur une période $T$. 
\end{frame}

\begin{frame} \frametitle{Étude des propriétés. II. stable}
  \begin{block}{Définition} Le système est dit \emph{A-stable} si les $x_i$ tendent vers 0 (quand $t\rightarrow \infty$) dans une problème où cette solution tend exponentiellement vers 0. \end{block}
  \begin{block}{Définition} Le système est dit \emph{L-stable} si les $x_i$ tendent vers 0 (quand $k\rightarrow \infty$) dans une problème où cette solution tend exponentiellement vers 0. \end{block}
\begin{exampleblock}{Euler avant, radioactivité : }
stable si $|(1-\dt)|\leq1, \quad \Rightarrow \dt\leq2$. 
Le schéma est dit  \emph{conditionnellement A-stable}.
\end{exampleblock}
\begin{exercice}
Déterminer les conditions de stabilité de  l'Euler arrière et l'Euler centré.
\end{exercice}
\end{frame}

\begin{frame} \frametitle{Etude des propriétés. III. monotone}
\begin{block}{Definition} si $x(t_i) < x(t_j)$, alors $x_i < x_j$. \end{block}
\begin{exampleblock}{Euler avant, radioactivité : }
Monotone si $\dt<1$.
Le schéma est dit  \emph{conditionnellement monotone}.
\end{exampleblock}
\begin{exercice}
Déterminer la condition de monotonicité pour l'Euler-centré, et montrer que l'Euler arrière est inconditionnellement monotone. 
\end{exercice}
\end{frame}


\begin{frame} \frametitle{Etude des propriétés. IV. Conserve l'énergie}
S'illustre pour l'oscillateur harmonique. Soit $\vec{x}=(p, q)\T$, 
le problème se formule comme $\ddt {\vec{x}} = M\vec{x}$ où 
$M:=\left( \begin{array}{cc} 0 &-1\\1 &0 \end{array}  \right)$. 

Le problème discrétisé s'écrit
\begin{equation*}
\vec{x}_{i+1} = M_d \vec{x}_i\quad \text{avec}
\end{equation*}
Euler avant : $M_d=(1+M\dt)$; Euler arrière : $M_d=(1-M\dt)^{-1}$ ; 
Euler centré : $M_d=(1-Mk/2)^{-1}(1+M\dt/2)$.

Sachant que le Hamiltonien est ici $\frac{1}{2}x\T x$, l'énergie
augmente de $\vec{x}_i$ à $\vec{x}_{i+1}$ d'un facteur ${M_d}\T M_d$.
Cette matrice doit donc être égale à la matrice identité dans un schéma conservatif. 

\begin{exercice}
Vérifier que le schéma Euler avant crée de l'énergie; que l'Euler arrière en détruit;
et que seul l'Euler centré est conservatif. Illustrer ces propriétés en intégrant le système.
\end{exercice}
\end{frame}

\begin{frame}
\begin{numerictip}
L'opération $M^{-1} x$  (inverse d'une matrice $\times$ vecteur, ou inverse d'une matrice $\times$ matrice) s'effectue efficacement au moyen du package \texttt{LAPACK}, implémenté dans la plupart des langages interprêtés destinés au calcul scientifique:
\medskip
\begin{tabular}{ll}
\hline
matlab / octave : & \texttt{M  $\setminus$ x} \\
R  : & \texttt{ solve(M, x) } \\
numpy / scipy : & \texttt{ np.solve(M, x) }  \\
\hline
\end{tabular}
\end{numerictip}
\end{frame}


\begin{frame}
\includegraphics{oscillateur_euler}
\end{frame}

\begin{frame}
\frametitle{Schéma leapfrog}
\begin{block}{Idée}
Mimer le schéma centré, mais selon un schéma explicite à deux pas de temps
\begin{equation}
x_{i+1} = x_{i-1} + 2\dt F(x_i)
\end{equation}
\end{block}

nécessite d'initialiser $x_0$ et $x_1$.

\begin{exercice}
Démontrer que le schéma leapfrog est instable dans le problème de décroissance radioactive. 
Indice : tenter une solution $x_i = \alpha^i$ et chercher les solutions pour $\alpha$. 
Le schéma est instable si $|\alpha|>1$. 
\end{exercice}
\end{frame}

\begin{frame}
\includegraphics{radio_leapfrog}
\end{frame}

\begin{frame}
\frametitle{Schéma de Heun (RK2)}
\begin{block}{Idée}
Approximer $F_{i+1}$ par $\tilde{F}_{i+1} := F(x_i + k F(x_i), t_{i+1})$. Peut se réecrire comme:
\begin{align*}
\tilde{x}_{i+1} &= x_i + \dt F_i \\
x_{i+1} &= x_i + {k \over 2}\left(F\left(x_i, t_i \right) + F\left(\tilde{x}_{i+1}, t_{i+1}\right)\right),
\end{align*}
qui fait apparaître la notion \emph{prédicteur - correcteur}.
\end{block}
%
\begin{exercice}
Développer $\tilde{F}_{i+1}$ en série autour de $F(x_i)$ et montrer que la méthode de Heun est
convergente à l'ordre 2 en k.
\end{exercice}
\end{frame}
\begin{frame}
\frametitle{Schéma de Runge-Kutta d'ordre 4 (RK4)}
\begin{align*}
k_1 &= k\cdot F_i \\
k_2 &= k\cdot F(x_i + k_1/2, t_i + k/2) \\
k_3 &= k\cdot F(x_i + k_2/2, t_i + k/2) \\
k_4 &= k\cdot F(x_i + k_3, t_i + k) \\
x_{i+1} &= x_i + {1\over 6}\left( k_1 + 2 k_2 + 2 k_3 + k_4 \right)
\end{align*}
\begin{block}{propriétés}
\begin{itemize}
\item Explicite
\item Conditionnellement stable
\item Troncation d'ordre 5; solution convergente à l'ordre 4.
\end{itemize}
\end{block}
\end{frame}
%

\begin{frame}
\includegraphics{error_withRK4}
\end{frame}

\begin{frame}
\begin{exercice}
L'équation logistique :
\begin{equation}
\ddt {x} = \lambda x ( 1-x)\, \quad x{0}=\alpha 
\end{equation}
a pour solution exacte:
\begin{equation}
x(t) = \frac{\alpha}{\alpha +   (1 - \alpha ) e^{-\lambda t }}.
\end{equation}
Comparer les erreurs des schémas explicites Euler avant, Heun et RK4 pour différents pas de temps comme sur la figure précédente. 

Le schéma Euler arrière peut s'établir analytiquement pour l'équation logistique, mais cela implique de résoudre, à chaque pas de temps, une équation du second degré,\ldots qui a donc deux solutions. Laquelle des deux solutions retenir?
\end{exercice}
\end{frame}

\begin{frame}
\frametitle{Erreurs d'arrondi (erreur machine)}
\begin{numerictip}  
Dans l'opération $a+b$, $|log_2(a)-log_2(b)|$  doit être plus grand que $m$, où
$m$ est le nombre de bits réservés pour la mantisse:
\begin{center}
\begin{tabular}{l|l|l}
type & m & $ log_{10} (2^m)$\\
\hline
float16 & 10 &  $3.0 $\\
float32 & 23 & $6.9 $\\
float64 & 52 & $15.6$\\
float128& 64 & $19.2$\\
float256 & 112 & $33.7$ \\
\hline
\end{tabular}
\end{center}
\end{numerictip}
\end{frame}

\begin{frame}
\includegraphics{rounding_error}
\end{frame}


\begin{frame} \frametitle{Schémas symplectiques}
Reprenons l'oscillateur harmonique, et considérons :
\begin{exampleblock}{Le schéma Verlet}
\begin{align*}
\tilde{p}_{i+1} &= p_i - \dt  q_i \\
q_{i+1} &= q_i + {\dt\over 2} (p_i + \tilde{p}_{i+1} ) \\
p_{i+1} &= p_i - {\dt\over 2} (q_i + q_{i+1} ) 
\end{align*}
\end{exampleblock}
\end{frame}

\begin{frame}
\includegraphics{oscillateur_energy1}
\end{frame}

\begin{frame}
\includegraphics{oscillateur_energy2}
\end{frame}

\begin{frame}
\frametitle{Méthode symplectique = respecte la géométrie du problème}
\begin{columns}
\column{0.45\textwidth}
\begin{tikzpicture}
 \draw [-open triangle 60] (-0.5, 0) -- (3, 0) node [below] {$q$}; 
 \draw [-open triangle 60] ( 0, -0.5) -- (0, 3) node [left] {$p$}; 

 \draw [-open triangle 60] ( 0.5, 0.7) node [below] {$x^a_{i}$} -- (1.3, 2.7) node [left] {$x^b_{i}$}; 
 \draw [-open triangle 60] ( 0.5, 0.7) -- (2.5, 0.9) node [right] {$x^c_{i}$}; 
 \draw [dashed] ( 2.5, 0.9) -- (3.3, 2.9) -- (1.3, 2.7) ; 
\end{tikzpicture}

\column{0.45\textwidth}
\begin{tikzpicture}
 \draw [-open triangle 60] (-0.5, 0) -- (3, 0) node [below] {$q$}; 
 \draw [-open triangle 60] ( 0, -0.5) -- (0, 3) node [left] {$p$}; 

 \draw [-open triangle 60] ( 0.5, 0.7) node [below] {$x^{a}_{i+1}$} -- (1.3, 2.2) node [left] {$x^{b}_{i+1}$}; 
 \draw [-open triangle 60] ( 0.5, 0.7) -- (3.0, 0.9) node [right] {$x^{c}_{i+1}$}; 
 \draw [dashed] ( 3.0, 0.9) -- (3.8, 2.4) -- (1.3, 2.2) ; 
\end{tikzpicture}
\end{columns}
Le schéma symplectique conserve l'aire des parallélogrammes (dans une métrique donnée)
\begin{theorem}
Soit $p_{i+1}=f(p_i, q_i)$ et $q_{i+1}=g(p_i, q_i)$. Alors le schéma est symplectique ssi
$\left| \begin{array}{cc}f_p & f_q \\ g_p & g_q \end{array}\right|=1, \ \forall{p, q}.$
\end{theorem}
\end{frame}
\begin{frame}
\begin{exercice}
Montrer que la méthode de Heun n'est pas symplectique pour l'oscillateur Harmonique; mais que la méthode de Verlet l'est bien.
\end{exercice}
\end{frame}