\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}.

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 Scilab 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+function euler_imp(f,a,b,N,yo)+.

On pourra comparer avec des cas connus. Par exemple

\begin{lstlisting}[caption={}]
e = euler_imp(deff('[u] = f(x,y)', 'u = y') ,0,2,20,1);
plot(e(1,:),e(2,:)) , plot(e(1,:),exp(e(1,:)),":r")
\end{lstlisting}

devrait donner:


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


Commentez la fonction:

\begin{lstlisting}[caption={}]
function  M = mystere(f,x,y,xmax,h,liste)
  if x > xmax then M = liste
  else M = mystere(f,x+h,y+h*f(x,y),xmax,h,[liste,y+h*f(x,y)])
  end
endfunction
\end{lstlisting}

Comment l'utiliser pour obtenir le même tracé que précédemment?


Commentez:

\begin{lstlisting}[caption={}]
max(abs((e(2,:) - exp(e(1,:)))))
\end{lstlisting}


\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}[caption={}]
 ans  = 
    0.0036059
\end{lstlisting}

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

\begin{lstlisting}[caption={}]
 ans  =
    0.0000006
\end{lstlisting}

\vfill

\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 Scilab).


\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  Scilab  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
Scilab 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}[caption={}]
t = 0:0.02:1;
uo = 1/5+1D-15;
to = 0;
u = ode(uo,to,t,f);
plot(t,u)
\end{lstlisting}

\newpage

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

\subsection{Systèmes}

Comme  Scilab  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$:

\begin{lstlisting}[caption={}]
deff('du=f(t,u)' , 'du(1) = 3*u(1) - 4*u(2) ; du(2) = u(1) + 3*u(2)');
t0 = 0 ; u0 = [2 ; 0] ; t = 0: 0.1 : 2 ;
z = ode(u0,t0,t,f) ;
plot(t,z(1,:),"ro") // donne u1
\end{lstlisting}

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.





%%% Local Variables: 
%%% TeX-master: "PolyMapleMPstar_12_13"
%%% End: 
