\chapterimage{Runge}
\chaptertext{
La   résolution  formelle   d'équations  différentielles   s'avère  très
compliquée et limitée: la plupart des problèmes ne peuvent être qu'approchés.

Le très  prolifique mathématicien  suisse Leonard \textsc{Euler}  mit au
point au  XVIII\eme siècle une  méthode de discrétisation  des équations
différentielles que  généralisèrent et  améliorèrent les  Allemands Carl
\textsc{Runge}  (1856-1927) dont  vous  pouvez admirer  le portrait,  et
Martin \textsc{Kutta} (1867-1944). Nous allons donc étudier ces méthodes
de \textsc{Runge-Kutta} très utiles pour l'ingénieur.

Il est à noté que \textsc{Runge} mit également en évidence des fonctions
assez régulières faisant diverger la méthode d'interpolation polynomiale
de \textsc{Lagrange} et donna son nom à ce phénomène.
}

\chapter{Discrétisation d'équations différentielles}


\section{Principe}



On  considère  une  équation  différentielle  ordinaire  (EDO  pour  les
intimes) $u'(t)=f(t,u(t))$. Lorsqu'on ne  connait pas de solution exacte
à  cette EDO,  on  essaye d'en  avoir  une bonne  approximation par  des
méthodes numériques.


Nous  allons par  exemple  \og  observer\fg $u$  à  intervalle de  temps
régulier: notons $h$ cette période d'échantillonnage.

Nous  allons  essayer  d'obtenir  une suite  d'approximations  $u_n$  de
$u(t_n)$ et d'évaluer la  pertinence de cette approximation, le principe
étant résumé sur la figure suivante:


%\begin{figure}
\begin{center}
  \includegraphics{mp-tex-FigDiscretEDO-2}
\end{center}
%\end{figure}


Nous considérerons des fonctions  $f$ $K$-lipschitziennes en chacune des
variables et des fonctions $u$ de  classe $\CR^1$ d'un intervalle I dans
$\bbr$ avec une condition initiale $u(t_0) = y_0$.

Pour une étude  mathématique complète des phénomènes  de convergence, on
pourra étudier l'article \url{http://www.math.ens.fr/~debarre/Exponentielle-Log.pdf
},   lire  le  chapitre   5  de  \textsf{Modélisation  à   l'oral  de
  l'agrégation} de  Laurent \textsc{Dumas}  et les chapitres  4 et  5 de
l'ouvrage \textsf{Analyse  numérique des équations  différentielles} par
Michel \textsc{Crouzeix} et Alain \textsc{Mignot} et surtout \og Solving
Ordinary  Diffrential  Equations\fg{}  de  \textsc{Hairer},  \textsc{Norsett} et \textsc{Wanner} chez Springer.

On notera  $I_n$ l'intervalle  $[t_n,t_{n+1}]$ avec  $t_{n+1}=t_n+h$, la
suite $(t_n)$ constituant une subdivision régulière de l'intervalle I.

On cherche donc à construire une suite finie $(u_n)$ telle que $u_0=y_0$
et  $u_{n+1}=u_n+h\varphi(t_n,u_n)$  (format  explicite  à  un  pas)  ou
$u_{n+1}=u_n+h\varphi(t_{n+1},u_{n+1})$ (format implicite à un pas).


\section{Méthode d'Euler explicite}

Intégrons l'équation  différentielle $u'(t)=f(t,u(t))$  sur l'intervalle
$I_n$:

$$u(t_{n+1})-u(t_n)= \displaystyle\int_{t_n}^{t_{n+1}} f(s,u(s)) \;{\rm d} s $$


L'utilisation de méthodes de quadratures  vues au TP précédent vont nous
permettre de calculer $u_{n+1}$ connaissant $u_n$.

La  méthode des  rectangles à  gauche donne  le schéma  d'\textsc{Euler}
explicite.


  Déterminer une fonction Python qui renvoie la liste des approximations
  des abscisses et des ordonnées de la série de points de la
  solution d'une
  équation différentielle par le schéma d'Euler explicite.

On pourra créer par exemple \verb+euler_imp(f,a,b,N,yo)+.

On pourra comparer avec des cas connus. 

%\begin{figure}
\begin{center}
  \includegraphics[height=5cm]{euler_exp1}
\end{center}
%\end{figure}



\section{Méthode d'Euler modifiée}


Nous avons  vu au  TP précédent  que la  méthode d'intégration  du point
milieu nous faisait gagner un ordre de précision.

Appliquez alors cette méthode à notre problème.

On pourra alors remarquer qu'avec $N=4096$ la méthode d'\textsc{Euler}
a une précision de l'ordre de:

\begin{lstlisting}
 0.0036059
\end{lstlisting}

et avec la méthode d'\textsc{Euler} modifiée:

\begin{lstlisting}
 0.0000006
\end{lstlisting}



\section{Méthode RK4}


On utilise cette fois la méthode de \textsc{Simpson}.
Notons: 

$$S= \frac{h}{6} \left(f(t_n,u(t_n)) + 4 f(t_n+h/2,u(t_n+h/2)) + f(t_{n}+h,u(t_{n}+h))\right) 
 $$

Alors faisons les approximations suivantes:

$$u(t_n)\to u_n = U_0$$

puis:

$$u(t_n+ \frac{h}{2}  ) \to u_n+ \frac{h}{2}  f \left(t_n,u_n\right)=U_1$$
en utilisant le rectangle gauche, ou bien:

$$u(t_n+  \frac{h}{2} )  \to  u_n+ \frac{h}{2}  f \left(t_n+  \frac{h}{2}
  ,u(t_n+ \frac{h}{2} )\right)\to u_n+ \frac{h}{2} f \left(t_n+ \frac{h}{2}
  ,U_1\right)= U_2$$ en utilisant le rectangle droit.

On remplace alors $f(t_n+h/2,u(t_n+h/2))$ par la moyenne de $U_1$ et $U_2$

Enfin:

$$u(t_n+h)\to   u_n+   \frac{h}{2}   \left(  t_n+   \frac{h}{2}   ,u(t_n+
  \frac{h}{2} )\right) \to u_n+ \frac{h}{2} \left( t_n+h/2,U_2\right)=U_3$$ 


Alors:


$$6S \to hf(t_n,u_n)+4h  \frac{f(t_n+h/2,U_1)+f(t_n+h/2,U_2)}{2} +hf(t_n+h,U_3)$$

On note conventionnellement:

$k_1=hf(t_n,u_n)$

$k_2=hf(t_n+h/2,u_n+k_1/2)$

$k_3 = hf(t_n+h/2,u_n+k_2/2)$

$k_4=hf(t_n+h,u_n+k_3)$

Alors $$u_{n+1}=u_n+ \frac{k_1+2k_2+2k_3+k_4}{6} $$


Programmez cette nouvelle méthode. Quelle est alors la précision pour la
fonction habituelle avec $N=4096$?


\section{La boîte de Pandore des approximations}


Encore  une fois,  dès  que  vous tentez  d'approcher  une solution,  de
nombreux  facteurs  peuvent   en  fait  vous  en   éloigner:  une  étude
mathématique sophistiquée doit  alors être envisagée (vous  verrez ça en
master...)

Voici quelques exemples apparemment anodins tirés de la littérature.

\subsection{Problème mal posé mathématiquement}

Que pensez-vous du problème:

$$
\begin{cases}
  u'(t)=2 \sqrt{\left|u(t)\right|}\\
  u(0)=0
\end{cases}
$$

\subsection{Problème mal posé numériquement}

Considérons le problème:

$$
\begin{cases}
  u'(t)=5u(t)-5,\quad t\in [0,50]\\
  u(0)=1
\end{cases}
$$

et   imaginons  une   petite   perturbation  de   la  valeur   initiale:
$u(0)=1+\varepsilon$.      Estimez       $\widetilde{u}(50)-u(50)$      pour
$\varepsilon=10^{-17}$ (de l'ordre de l'epsilon de Python).


\subsection{Problème mal conditionné}

Cette fois, le problème est le suivant:
$$
\begin{cases}
  u'(t) = -150 u(t) - 30\\
  u(0) = 1/5
\end{cases}
$$


Vérifiez qu'il est bien posé mathématiquement et numériquement.

Approchons-le  maintenant  à  l'aide   de  la  méthode  d'\textsc{Euler}
explicite.

Vérifiez  qu'on  obtient  $$u_n= \frac{1}{5}  +  \left(1-  \frac{150}{N}
\right)^n  \left(u_0- \frac{1}{5} \right) 
$$

Qu'obtenez-vous pour $N=50$ et $u_0= \frac{1}{5} +\varepsilon$?


Observez  maintenant  avec nos  fonctions  Python  avec les  différentes
méthodes vues précédemment.


\subsection{Problèmes raides}

Il existe des problèmes qui résistent à toutes les méthodes usuelles. On
trouve par  exemple ces phénomènes en  chimie où deux réactions  ont des
échelles de temps très différentes.

Considérez par exemple:

$$
\begin{cases}
 u'(t) = -1000(u(t) - \exp(-t)) -\exp(-t)\\
 u(0) = 0
\end{cases}
$$

et testez vos fonctions sur $[0,1]$. 

Cependant, vous pouvez utiliser la fonction \verb+ode+ de
Python qui utilise au mieux les méthodes précédentes en adaptant de plus
les pas ce qui permet parfois d'obtenir de bons résultats.

La doc nous dit en effet:

\begin{quote}

{\itshape 
Le solveur lsoda  du package ODEPACK est utilisé par  défaut. Il choisit
automatiquement  entre un  schéma  prédicteur-correcteur  d'Adams et  un
schéma   adapté   au  systèmes   raides   (stiff)   de  type   \og{}Backward
Differentiation Formula\fg{} (BDF).

 Initialement  le  schéma  adapté  aux
systèmes  non raides  est  choisi  puis la  méthode  adaptée est  ensuite
choisie dynamiquement. 
}

\end{quote}

Pour le problème mal conditionné du paragraphe précédent:


\begin{lstlisting}
uo = 1/5+1E-15;
to = 0;
u = ode(uo,to,t,f);
\end{lstlisting}



\section{Systèmes différentiels et ordre 2}

\subsection{Systèmes}

Comme  Python  travaille  sur  des  vecteurs,  on  peut  considérer  des
fonctions dans $\bbr^2$ de manière naturelle.


Supposons que nous voulions résoudre le système:

$$
\begin{cases}
u_1'(t) = 3u_1(t) - 4u_2(t)\\
u'_2(t) = u_1(t) + 3u_2(t)\\
u_1(0) = 2\\
u_2(0) = 0
\end{cases}
$$

On utilise  la même  syntaxe que  précédemment...mais avec  une fonction
$u : \bbr\to \bbr^2$.


Comparez avec la solution attendue.


\subsection{Ordre 2}


On considère à présent l'équation bien connue du pendule pesant amorti:
$y''+y'+\sin(y)=0$ avec $y(0)=0$ et $y'(0)=1$.


Posez $u = \tr (y,y')$ et essayez de revenir à un système
différentiel que vous résoudrez sur $[0,10]$.

Vous tracerez le portrait de phase.


\section{CCP MP : Sujet 0 de SI 2015}

{\itshape 
Après  25 questions,  voici l'éternelle,  l'inévitable question  sur la  méthode
d'Euler... Apparemment, pour un mécanicien,  le rôle de l'informatique se résume
à donner des solutions approchées d'une équation différentielle...

Nous sommes donc à la page 14 du sujet :-)
On introduit l'inévitable équa diff.
}

[On obtient] l'équation différentielle suivante, liant la vitesse de rotation du
moteur $ ω_m(t)$ à la vitesse de consigne $ ω_{\rm const}(t)$:

$$\frac{1}{ω_0^2}·\frac{{\rm d}^2  ω_m(t)}{{\rm d}t^2}+\frac{2m}{ω_0}·\frac{{\rm
    d}ω_m(t)}{{\rm d}t}+ω_m(t)=K·ω_{\rm const}(t) $$


Dans la suite,  on s'intéressera uniquement à la réponse  indicielle unitaire du
système  dans les  conditions de  Heaviside, ce  qui revient  à prendre  $ ω_{\rm
  const}(t)=u(t)$ avec $ ω_m(0)=0$ et $ ω'_m(0)=0$.

En posant le vecteur $ \mathbf{Y}$ tel que $ \mathbf{Y}(t)=
\begin{pmatrix}
  ω_m(t)\\
  \frac{{\rm d}ω_m(t)}{{\rm d}t}
\end{pmatrix}
$, l'équation différentielle à résoudre peut se mettre sous la forme:

$$\frac{{\rm     d}\mathbf{Y}(t)}{{\rm     d}t}=\mathbf{F}(t,\mathbf{Y}(t))\quad
\text{avec}\quad \mathbf{F}(t,\mathbf{Y}(t))=
\begin{pmatrix}
  \frac{{\rm d}ω_m(t)}{{\rm d}t}\\
  Kω_0^2-ω_0^2ω_m(t)-2mω_0\frac{{\rm d}ω_m(t)}{{\rm d}t}
\end{pmatrix}
$$

La réponse  $ ω_m(t)$ recherchée  sur l'intervalle $ [0,T_{\max}]$  sera obtenur
par...la méthode d'\textsc{Euler} explicite ;-)

Le pas  de temps,  noté \lstinline{pas}, sera  choisi constant.  L'intervalle de
temps    discrétisé     est    alors    représenté    par     le    tableau    $
T=[t_0=0,t_1,...,t_{N-1}=T_{\max}]$.

Pour chaque pas de temps, une valeur  approchée $ \mathbf{Y}_i$ de la solution $
\mathbf{Y}(t_i)$ de  l'équation différentielle est recherchée.  L'ensemble des $
\mathbf{Y}_i$ représente  $ N$ vecteurs  de dimension  2, qui seront  stockés en
mémoire sous la forme du tableau:

$$\mathbf{SY}=
\begin{pmatrix}
  ω_m(0)&ω'_m(0)\\
  ω_m(t_1)&ω'_m(t_1)\\
  \vdots & \vdots\\
  ω_m(T_{\max{}})&ω'_m(T_{\max{}})\\
\end{pmatrix}
$$

{Et voici enfin la première question...}

\begin{enumerate}
\item Écrire un  algorithme ou une fonction permettant de  calculer le tableau $
  T$.
\item Écrire une fonction \lstinline{f1(ti,yi)} qui prend en arguments la valeur
  du  temps discrétisé  $ i$  et la  valeur du  vecteur $  \mathbf{Y}$ au  temps
  discrétisé    \lstinline{ti}    et   qui    retourne    la    valeur   de    $
  \mathbf{F}(t,\mathbf{Y}(t))$.

\item  Donner  la  relation  de  récurrence qui  lie  $  \mathbf{Y}_{i+1}$  à  $
  \mathbf{Y}_i$ et à $ \mathbf{F}(t,\mathbf{Y}(t))$  en fonction du pas de temps
  \lstinline{pas}.

\item Écrire une fonction \lstinline{EulerExplicite(Yini,h,Tmax,F)} qui prend en
  arguments \lstinline{Yini}, un  tableau de dimension 2  contenant la condition
  initiale de $ \mathbf{Y}(t)$, \lstinline{h}  le pas de temps, \lstinline{Tmax}
  l'instant  final  du calcul,  et  \lstinline{F}  la  fonction du  problème  de
  \textsc{Cauchy}.

Cette fonction renverra le tableau \lstinline{SY}. L'appel à cette fonction dans
le    programme     se    fera     avec    la    commande     \lstinline{SY = EulerExplicite(Y0,pas,tmax,f1)}.

\item Si le pas de temps es divisé par un facteur 10, comment évolue l'erreur de
  calcul?

\item Donner  la complexité  de cette  méthode pour $  T_{\max{}}$ fixé  et indiquer
  comment évolue  le temps de  calcul quand  le pas de  temps est divisé  par un
  facteur 10.

On suppose  que la  quantité de  mémoire nécessaire pour  réaliser le  calcul se
limite  au stockage  de la  matrice \textbf{SY}  et du  vecteur \textbf{T}.  Ces
éléments sont  représentés en  mémoire sous  forme de  tableaux de  flottants en
double précision.

\item Déterminer  le nombre d'octets  nécessaire en mémoire pour  réaliser cette
  simulation numérique avec nombre de pas de temps $ N=\num{10000}$.

\item \textbf{Convergence  de l'algorithme}.  Avant d'appliquer  l'algorithme de
  recherche  du temps  de réponse  à  $ 5\%$,  il  est nécessaire  de savoir  si
  l'algorithme  de simulation  converge. Après  calcul, la  solution $  ω_m$ est
  stockée dans la variable \textbf{W}.

On suppose que  la solution converge si toutes  les valeurs de $ ω_m$  pour $ t$
appartenant  à $  [0.9  T_{\max{}},T_{\max{}}]$ ne  s'écartent pas  de  la valeur  finale
calculée, $ ω_m(T_{\max{}})$, de plus de $ 0,1\%$.

Écrire une  fonction \lstinline{TestConvergence(t,w)}  qui prend en  argument le
tableau  des temps  \lstinline{t} et  la solution  \lstinline{w} et  qui renvoie
\lstinline{True} si le critère est vérifié et \lstinline{False} sinon.

\item \textbf{Détermination du temps de réponse à $ \mathbf{5\%}$ du système}

Écrire une  fonction \lstinline{CalculT5} qui renvoie  le temps de réponse  à $
5\%$ si la convergence est vérifiée et $ -1$ s'il n'est pas vérifié. Les entrées
et sorties de cette fonction seront clairement définies.

La programmation et la simulation de l'algorithme précédent, dans les conditions
de l'étude, ont permis de déterminer une temps de réponse de $ \num{0,045}s$.















\end{enumerate}


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "PolyIPT2014_15"
%%% End: 


