\chapterimage{ingenieur}
\chaptertext{}

\chapter{Ingénieur, quel beau métier!}

\section{Pour en finir avec les résolutions numériques d'équas diff}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\modeexercice


\subsection{Position du problème}

Soit $U$ un ouvert de $ℝ×ℝ^n$ et $f\colon U\longrightarrow ℝ^n$ une
application continue sur l'ouvert $U$.

Alors pour $(t_0,x_0)∈U$, on appelle \emph{solution du problème de
\textsc{Cauchy}} :
\[
\begin{cases}
x'(t) = f(t, x(t))\\
x(t_0) = x_0 \qquad {\mathcal{C}}
\end{cases}
\]
tout couple $(I,x)$ où $I$ est un intervalle contenant $t_0$ et $x\colon
I\longrightarrow ℝ^n$ une solution sur $I$ de l'équation différentielle
$x'(t)= f(t, x(t))$ telle que $ x(t_0) = x_0$.

L'équation différentielle qui apparaît dans $(\mathcal{C})$ sera appelée équation
différentielle scalaire si $n=1$, et système différentiel sinon.


Si l'application $f\colon U\longrightarrow ℝ^n$ est de \emph{classe}
$\mathcal{C}^1$ sur l'ouvert $U$, alors, d'après le théorème de
\textsc{Cauchy-Lipschitz},
pour $(t_0,x_0)∈U$, il existe une unique solution maximale au problème
de \textsc{Cauchy} $(\mathcal{C})$.



Dans la plupart des cas, on ne sait pas résoudre explicitement le problème de
Cauchy $(\mathcal{C})$ ; d'où la nécessité de mettre au point des méthodes
numériques de résolution approchée d'un tel problème.

%---------------------------
\subsection{Méthode d'\textsc{Euler}}
%---------------------------

La méthode d'\textsc{Euler} est la méthode la plus simple pour résoudre numériquement une
équation différentielle. Elle présente un réel intérêt théorique, puisque elle
peut être utilisée pour démontrer le théorème de \textsc{Cauchy-Lipschitz}. Toutefois,
son intérêt pratique est limité par sa faible précision.

L'idée d'\textsc{Euler} consiste à utiliser l'approximation :
\[
x'(t) ≈ \frac{x(t+h) - x(t)}{h}
\quad \text{ pour } h \text{ petit }
\]
Autrement dit,
\[
x(t+h) ≈ x(t) + h ⋅ x'(t) = x(t) + h ⋅ f(t, x(t))
\]
Partant du point $(t_0, x_0)$, on suit alors la droite de pente $f(t_0, x_0)$
sur l'intervalle de temps $[t_0, t_0+h]$.


On pose alors :
\[
\begin{cases}{t_1 = t_0 + h}\\
{x_1 = x_0 + h ⋅ f(t_0, x_0)}
\end{cases}
\]

De nouveau, partant du point $(t_1, x_1)$, on suit alors la droite de pente $f(t_1, x_1)$
sur l'intervalle de temps $[t_1, t_1+h]$, et ainsi de suite.

On construit ainsi une suite de points de la manière suivante :
\[
∀k∈⟦0, N⟧, \quad
\begin{cases}{t_{k+1} = t_k + h}\\
{x_{k+1} = x_k + h ⋅ f(t_k, x_k)}
\end{cases}
\]
La ligne brisée joignant les points $\bigl\{(t_k, x_k) | k∈⟦0,N⟧\bigr\}$ donnera une
solution approchée de notre équation différentielle.

Il faut bien comprendre qu'à chaque étape, on repart dans la direction d'une
solution exacte de l'équation différentielle, mais qui n'est pas celle qui est
solution du problème de \textsc{Cauchy} initial. Sur le graphique ci-contre, on trace une
solution approchée pour la condition initiale $x(0)=1$ et les courbes intégrales
de notre équation qui passent par les points $(t_k, x_k)$.


%\begin{figure}
\begin{center}
   \includegraphics[scale=0.65]{ode_exponentielle2}
\end{center}
%\end{figure}


Pour juger de la qualité d'une méthode (ou schéma) numérique de résolution d'équations
différentielles, il faut prendre en compte plusieurs critères.

\begin{itemize}
    \item
        \index{erreur de consistance}
        L'\emph{erreur de consistance} donne l'ordre de grandeur de l'erreur
        effectué à chaque pas. Par exemple, dans le cas de la méthode
        d'\textsc{Euler},
        l'erreur de consistance mesure l'erreur qu'entraîne le fait d'approcher
        le nombre dérivé par un taux d'accroissement.
        La sommation des erreurs de consistance à chaque pas donnera l'ordre de
        grandeur de l'erreur globale (sous des hypothèses de régularité pour la
        fonction $f$).
        À l'aide de la formule de \textsc{Taylor-Lagrange}, on peut montrer que dans le
        cas de la méthode d'\textsc{Euler}, l'erreur de consistance est dominée par $h^2$
        ; on dit alors que cette méthode est d'ordre 1.
        Ainsi, d'un point de vue théorique,  plus le pas est petit, meilleure
        sera l'approximation.

        \index{schéma numérique!consistant}
        Un schéma numérique est dit \emph{consistant} si la somme des erreurs de
        consistance tend vers 0 quand le pas $h$ tend vers 0 ; ce qui est le cas
        de la méthode d'\textsc{Euler} (en gros, car $Nh^2∼h\to0$ où $N$ désigne le
        nombre d'itérations).

    \item
        \index{schéma numérique!stable}
        La \emph{stabilité} contrôle la différence entre deux solutions
        approchées correspondant à deux conditions initiales voisines : un
        schéma est dit \emph{stable} si un petit écart entre les conditions
        initiales $x(t_0)=x_0$ et $\tilde{x}(t_0)=x_0 + ε$ et de petites erreurs
        d'arrondi dans le calcul récurrent des $\tilde{x}_k$ provoquent une erreur
        finale $\bigl|x_k-\tilde{x}_k\bigr|$ contrôlable.
        La méthode d'\textsc{Euler} est une méthode stable.

    \item
        \index{schéma numérique!convergent}
        Un schéma est dit \emph{convergent} lorsque l'erreur globale (qui
        est le maximum des écarts entre la solution exacte et la solution
        approchée) tend vers 0 lorsque le pas tend vers 0. En fait, pour que le
        schéma soit convergent, il suffit que la méthode soit consistante
        \emph{et} stable.

    \item Enfin pour la mise en application du schéma, il faut aussi prendre en
        compte l'influence des erreurs d'arrondi. En effet, afin de minimiser
        l'erreur globale théorique, on pourrait être tenté d'appliquer la
        méthode d'\textsc{Euler} avec un pas très petit, par exemple de l'ordre de
        $10^{-16}$, mais ce faisant, outre que le temps de calcul deviendrait
        irréaliste, les erreurs d'arrondi feraient diverger la solution
        approchée très rapidement, puisque pour des flottants Python de l'ordre de
        $10^{-16}$, les calculs ne sont plus du tout exacts !

        En pratique, il faut prendre $h$ assez petit (pour que la méthode
        converge assez rapidement du point de vue théorique), mais pas trop
        petit non plus (pour que les erreurs d'arrondis ne donnent pas lieu à
        des résultats incohérents, et que les calculs puissent être effectués en
        un temps fini). La question de trouver de manière théorique
        le pas optimal peut s'avérer un problème épineux. Dans un premier temps,
        on peut se contenter de faire des tests pratiques comme dans l'exemple
        suivant.
\end{itemize}



Dans le script suivant, on applique la méthode d'\textsc{Euler} sur l'intervalle
$[0, 1]$ avec la condition initiale $x(0)=1$.

\begin{pythoncode}
from math import exp


print("{:^10} | {:^12} | {:^12} | {:^13}".format('h','x(t)','exp(t)','erreur'))
print('-'*57)
for i in range(1, 8):
    h, t, x  = 10**(-i), 0, 1
    while t < 1:
        t, x = t + h, x * (h + 1)
    print("{0:>10g} | {1:>.10f} | {2:>.10f} | {3:> .10f}"
            .format(h, x, exp(t),  exp(t) - x))
\end{pythoncode}


\begin{pythoncode}
    h      |     x(t)     |    exp(t)    |    erreur
---------------------------------------------------------
       0.1 | 2.8531167061 | 3.0041660239 |  0.1510493178
      0.01 | 2.7048138294 | 2.7182818285 |  0.0134679990
     0.001 | 2.7169239322 | 2.7182818285 |  0.0013578962
    0.0001 | 2.7184177414 | 2.7185536702 |  0.0001359288
     1e-05 | 2.7182954199 | 2.7183090114 |  0.0000135915
     1e-06 | 2.7182804691 | 2.7182818285 |  0.0000013594
     1e-07 | 2.7182819660 | 2.7182820996 |  0.0000001336
     1e-08 | 2.7182817983 | 2.7182818347 |  0.0000000363
     1e-09 | 2.7182820738 | 2.7182818299 | -0.0000002438
\end{pythoncode}

Pour un pas $h\lesssim 10^{-9}$, le temps de calcul devient trop important, de
plus, on voit poindre l'effet des erreurs d'arrondis, puisque la solution
approchée est devenu supérieure à l'exponentielle, ce qui est mathématiquement
faux, en vertu de l'inégalité :
\[
∀n∈ℕ, \quad \left(1+\frac{1}{n}\right)^n ≤ e
\]

En général, il ne suffit pas qu'un schéma numérique soit convergent pour qu'il
donne de bons résultats sur n'importe quelle équation différentielle.
Encore faut-il que le problème soit mathématiquement bien posé (en
particulier, les hypothèses du théorème de \textsc{Cauchy-Lipschitz} doivent être
vérifiées), qu'il soit numériquement bien posé (continuité suffisamment bonne
par rapport aux conditions initiales) et qu'il soit bien conditionné (temps de
calcul raisonnable).



Considérons par exemple le problème de \textsc{Cauchy} :
\[
\begin{cases}
{x(1)=1}\\
{x'(t) = 3\frac{x(t)}{t} - \frac{5}{t^3}}
\end{cases}
\]
dont la solution est la fonction ???%$φ\colon t\mapsto 1/t^2$.



\begin{center}
    \includegraphics[scale=0.9]{mal_pose}
\end{center}

On constate qu'ici la solution approchée s'écarte assez rapidement de la
solution exacte $φ$. En effet, la forme générale des solutions de l'équation
différentielle précédente est donnée par $t\mapsto λx^3 + 1/x^2$.  La condition
$x(1)=1$ impose $λ=0$, mais dès qu'on s'écarte de la solution $φ$ on suit
\og tangentiellement\fg{} des courbes intégrales qui comportent un terme en
$λx^3$ et donc qui diffère notablement de la solution.  Donc le problème est
\emph{mal posé}. Ici le schéma numérique n'est pas en cause.

Soit à résoudre à présent le problème de \textsc{Cauchy} :
\[
\begin{cases}
{x(1)=1}\\
{x'(t) =100\,(\sin(t)-x)}
\end{cases}
\]
dont la solution est la fonction $φ\colon t\mapsto \frac1{10001} (-100\cos t +
10000\sin t + 100\exp(-100⋅t)) ≈\sin t$.
Que se passe-t-il si on résout ce problème avec un pas de l'ordre de $0.02$ ?


\begin{center}
    \includegraphics[scale=1]{mal_conditionne}
\end{center}
Avec un tel pas, la solution approchée oscille en s'éloignant de plus en plus de la
solution exacte.
En effet, le schéma numérique est donné par :
\[
x_{k+1} = x_k + 100\,h\, \left(\sin t_k - x_k\right)
= (1 - 100h) x_k + 100h \sin t_k
\]
Donc une erreur de $ε_k$ sur $x_k$ aura les répercussions suivantes sur le calcul
des termes ultérieurs :
\[
ε_{k+1} = (1-100h)ε_k \implies ε_{k+n} = (1-100h)^n ε_k
\]
Ainsi donc, tant que $1-100h≥-1$, c'est-à-dire tant que $h\lesssim 0.002$, une petite
erreur sur l'un des termes aura des répercussions allant en s'amplifiant.

Bien que le problème soit numériquement bien posé, il est en fait  \emph{mal
conditionné}.




%---------------------------
\subsection{Méthode de \textsc{Runge-Kutta} d'ordre 4}
%---------------------------



Si, comme nous l'avons déjà dit, la méthode d'\textsc{Euler} présente un intérêt
théorique, on préfère en pratique des méthodes d'ordre plus élevé.
Parmi la multitude des schémas numériques (méthodes à un pas comme celle de
\textsc{Taylor}, méthodes à pas
multiples comme celles d'\textsc{Adams-Bashforth}, d'\textsc{Adams-Moulton}, de
prédiction-correction) celle qui présente le meilleur rapport
précision/complexité est certainement celle de \textsc{Runge-Kutta} d'ordre 4.

En voici le schéma pour $∀k∈⟦0,N⟧$,
\[
\left\lbrace{}
\begin{array}{rcl}
t_{k+1} &=& t_k + h \\
x_{k+1} &=& x_k + \frac16 \left(k_1+2k_2+2k_3+k_4\right)
\end{array}\right.
\quad \text{ où }\quad
\left\lbrace{}
\begin{array}{rcl}
    k_1 &=& h⋅f(t_n, x_n) \\
    k_2 &=& h⋅f\left(t_n + \frac{h}2, x_n + \frac{k_1}2\right) \\
    k_3 &=& h⋅f\left(t_n + \frac{h}2, x_n + \frac{k_2}2\right) \\
    k_4 &=& h⋅f\left(t_n + h, x_n + k_3\right)
\end{array}\right.
\]
Les coefficients qui apparaissent dans ces formules mystérieuses sont
judicieusement ajustés pour obtenir une méthode d'ordre 4, sans pour autant
avoir à calculer les dérivées successives de $f$ (comme c'est le cas dans les
méthodes de \textsc{Taylor}), ou à recourir à une formule de récurrence d'ordre au moins
2 pour définir $x_k$ (comme c'est le cas dans les méthodes à pas multiples).




Définissez  une  fonction  adaptée  à  ce schéma.  Vous  devriez  observer  pour
l'exponentielle


\[
\begin{cases}{x'(t) = x(t)}\\
{x(0) = 1}
\end{cases}
\qquad{$\mathcal{C}$}
\]

\begin{pythoncode}
    h      |   erreur euler    |    erreur rk4
----------------------------------------------------
       0.1 |  1.5104931784e-01 |  2.5338874514e-06
      0.01 |  1.3467999038e-02 |  2.2464341498e-10
     0.001 |  1.3578962232e-03 |  2.2204460493e-14
    0.0001 |  1.3592881559e-04 | -2.6645352591e-13
     1e-05 |  1.3591551171e-05 | -5.2180482157e-12
     1e-06 |  1.3591611072e-06 |  2.1464163780e-11
\end{pythoncode}



On constate que pour $h=0.001$, la méthode d'\textsc{Euler} ne donne que trois décimales
significatives du nombre $e$, alors que la méthode de \textsc{Runge-Kutta} en donne
quatorze ! On remarque également qu'il est inutile d'appliquer la méthode de
\textsc{Runge-Kutta} avec un pas $h\gtrsim 0.001$, puisque dans ce cas les erreurs
d'arrondi prennent le dessus par rapport au gain théorique de précision.


\subsection{Système}

Soit à représenter les courbes intégrales du système autonome suivant, appelé
oscillateur de van der Pol :
\[
\begin{cases}
\frac{\dd x}{\dd t} = y\\
\frac{\dd y}{\dd t} = \frac y2 -x -y^3
\end{cases}
\]

Les variables $x$,
$k1$,... $k4$ représentent à présent des listes. Or une instruction comme
$x+k1/2$ renverra un message d'erreur. En effet, la division d'une liste par un
entier n'est pas définie ; de plus, l'addition de deux listes effectue la
concaténation des listes, et non l'addition des listes composante par
composante.

Il faut alors travailler non plus avec des listes, mais avec des objets
supportant la vectorisation, comme c'est le cas des tableaux de la librairie de
tierce partie \verb?NumPy? ; en outre, cette option a l'avantage de réduire de
manière significative les temps de calcul.


%\begin{figure}
\begin{center}
 \includegraphics[scale=.65]{syst_autonome}
\end{center}
%\end{figure}


%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\subsubsection{Un exemple d'équation différentielle linéaire d'ordre 2}
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

On rappelle qu'une équation différentielle linéaire d'ordre $n$ peut toujours
s'écrire comme un système différentiel d'ordre 1. Illustrons ceci dans le cas
$n=2$ :
en notant :
\[
A(t)  = \begin{pmatrix} 0 & 1 \\ -b(t) & -a(t) \end{pmatrix},\quad
X(t)=\begin{pmatrix} x(t)\\ x'(t) \end{pmatrix},\quad
C(t)  = \begin{pmatrix} 0  \\ c(t)  \end{pmatrix}
\]
on a les équivalences

\begin{align*}
    ∀t∈I\quad x''+a(t)x'+b(t)x = c(t)
& \iff ∀t∈I\quad \begin{cases}
    x' =        x'    \\
    x''= -bx - ax' + c 
  \end{cases}
  \\[1ex]
    &\iff ∀t∈I\quad
    \frac{\dd}{\dd t}\begin{pmatrix} x\\ x' \end{pmatrix}
    = \begin{pmatrix} 0 & 1 \\ -b & -a \end{pmatrix}
    ⋅\begin{pmatrix} x\\ x' \end{pmatrix}
    + \begin{pmatrix} 0  \\ c  \end{pmatrix} \\[1ex]
    &\iff ∀t∈I\quad X' = A(t)⋅X +C(t)
\end{align*}

\index{équation!de \textsc{Bessel} }
Considérons, par exemple,  l'équation différentielle du 2\ieme{} ordre dite de
\textsc{Bessel}  :
\[
t^2 x''(t) + t x'(t) + \left(t^2-α^2\right) x(t) = 0
\iff
X'(t) = A(t)⋅X(t)
\enskip \text{ avec }\enskip
A(t) = \begin{pmatrix}
    0 & 1 \\
    \tfrac{α^2-t^2}{t^2} & -\tfrac{1}{t}
\end{pmatrix}
\]


\begin{center}
    \includegraphics[scale=1]{bessel}
\end{center}


\subsection{Lien avec les méthodes de quadrature}



Version objet:


\begin{pythoncode}
class Integration(object):
    def __init__(self, a, b, n):
        self.a, self.b, self.n = a, b, n
        self.poids_pivots = self.quadrature()

    def quadrature(self):
        raise NotImplementedError

    def integrale(self, f):
        return sum(w * f(x) for w, x in self.poids_pivots)

class Rectangles(Integration):
    def quadrature(self):
        a, b, n = self.a, self.b, self.n
        h = (b-a) / n
        return [(h, a + i*h) for i in range(n)]

class Trapezes(Integration):
    def quadrature(self):
        a, b, n = self.a, self.b, self.n
        h = (b-a) / n
        return [(h, a + (i+0.5)*h) for i in range(n)]

class Simpson(Integration):
    def coeff(self, i, N):
        if i == 0 or i == N:
            return 1
        elif i % 2 == 1:
            return 4
        else:
            return 2

    def quadrature(self):
        a, b, n = self.a, self.b, self.n
        if n % 2 != 1:
            n += 1
        h = (b-a) * 0.5 / n
        return [(h/3*self.coeff(i, 2*n), a + i*h) for i in range(2*n+1)]

class Romberg(object):
    def __init__(self, a, b, m):
        self.a, self.b, self.m = a, b, m

    def integrale(self, f):
        a, b, m = self.a, self.b, self.m
        A = [[(b-a) * (f(a)+f(b)) * 0.5]]
        for k in range(1, m+1):
            h = (b-a) / 2**k
            Ap = h * sum([f(a + j*h) for j in range(1, 2**k, 2)])
            A.append([0.5 * A[k-1][0] + Ap])
        for t in range(1, m+1):
            for k in range(t, m+1):
                A[k].append((4**t * A[k][t-1] - A[k-1][t-1]) / (4**t-1))
        return A[m][m]
\end{pythoncode}

   Représenter la spirale de \textsc{Cornu} définie par
\[
x(t) = ∫_{0}^{t} \cos\left(u^2\right) \dd u
\quad \text{et} \quad
y(t) = ∫_{0}^{t} \sin\left(u^2\right) \dd u
\]


\section{Recherches}





\begin{exercice}[Lynx et lapins]


  \begin{multicols}{2}
    
\begin{figure}
\begin{center}
  \includegraphics[height=4cm]{lokta}
\end{center}
\caption{Alfred \textsc{Lokta}\\US (1880-1949)}
\end{figure}

\begin{figure}
\begin{center}
  \includegraphics[height=4cm]{volterra}
\end{center}
\caption{Vito Volterra\\It. (1860-1940)}
\end{figure}

  \end{multicols}

Vous connaissez ce problème par c\oe ur...



$$    \left\{\begin{array}{ccc}\dfrac{\mathrm{d}x(t)}{\mathrm{d}t} &=& x(t)\ \Big(\alpha - \beta y(t)\Big) \\ \dfrac{\mathrm{d}y(t)}{\mathrm{d}t} &=& -y(t) \ \Big(\gamma - \delta x(t)\Big) \end{array}\right. 
$$

\begin{itemize}
\item     $t$ est le temps ;
   \item  $x(t)$ est l'effectif des proies en fonction du temps ;
  \item   $y(t)$ est l'effectif des prédateurs en fonction du temps ;
    \item        les       dérivées        $\mathrm{d}x(t)/\mathrm{d}t$       et
      $\mathrm{d}y(t)/\mathrm{d}t$ représentent la  variation des populations au
      cours du temps.

    \end{itemize}
    

Les paramètres suivants caractérisent les interactions entre les deux espèces :
\begin{itemize}
\item     $\alpha$, taux de reproduction des proies (constant, indépendant du nombre de prédateurs) ;
  \item     $\beta$, taux de mortalité des proies dû aux prédateurs rencontrés ;
  \item   $\gamma$, taux de mortalité des prédateurs (constant, indépendant du nombre de proies) ;
   \item  $\delta$, taux de reproduction des prédateurs en fonction des proies rencontrées et mangées ;
\end{itemize}



% \begin{pythoncode}
% def voltera(t,V):
%     x,y = V[0], V[1]
%     return np.array([x*(3 - 2*y), y*(-4 + x)]) 
% \end{pythoncode}


%Prenez par exemple $α=3, β=-2, γ=-1, δ=4$.

Vous connaissez la méthode d'Euler par c\oe ur...

% \begin{pythoncode}
% def EulerV(F,V0,h,n):
%     t = 0
%     V = V0
%     X , Y , T = [V[0]] , [V[1]], [0]
%     for i in range(n):
%         V = V + h*F(t, V)
%         t += h
%         T.append(t)
%         X.append(V[0])
%         Y.append(V[1])
%     plt.plot(T,X)
%     plt.plot(T,Y)
%     # plt.plot(X,Y)
% \end{pythoncode}

Tracez $x$ et $y$ en fonction de $t$ puis $y$ en fonction de $x$.

\end{exercice}



\begin{exercice}[Deux corps dans l'espace...]
  

Bla bla Principe Fondamental de la Dynamique bla bla bla:

$$
r''-rθ'^2=-\frac{GMt}{r^2}+\frac{Ar}{M_0-At}\qquad rθ''+2r'θ'=\frac{Arθ'}{M_0-At}
$$


Mais  comment faire!?  Mais  comment  appliquer la  méthode  d'Euler?  Y a  deux
inconnues et en plus y a des dérivées secondes...


On pose $u_1=r$, $u_2=r'$, $u_3=θ$, $u_4=θ'$ et tout s'arrange.


\end{exercice}



\begin{exercice}[Centrale 2015 avec une solution exacte du problème des 3 corps]
 \textit{Sur une idée de \href{http://rouxph.blogspot.fr/}{Philippe \textsc{Roux}}}

% Voici le PFD appliqué à trois corps célestes:

% $${d^2\over dt^2}X_i =\sum_{j\neq i}-{{\cal G}m_j(X_i-X_j)\over |X_i-X_j|^3}=\sum_{j\neq i} F_g(X_i-X_j,m_j) ~~~~i=1,2,3$$
% avec 6 conditions initiales $X_i(0), {d\over dt}X_i(0)$. 

% $$
% U=\begin{pmatrix}X_1\\ {d\over dt}X_1\\ X_2\\ {d\over dt}X_2\\ X_3\\{d\over dt}X_3 \end{pmatrix} 
% ~~~~{d\over dt}U =F(U)=\begin{pmatrix} 
% U_2\\ 
% F_g(U_1-U_3,m_2)+F_g(U_1-U_5,m_3)\\ 
% U_4\\ 
% F_g(U_3-U_1,m_1)+F_g(U_3-U_5,m_3)\\ 
% U_6\\ 
% F_g(U_5-U_3,m_2)+F_g(U_5-U_1,m_1)\\ 
%  \end{pmatrix}
% $$

% Alors:


% $${d\over dt}U_{2i-1}= {d\over dt}X_i=U_{2i}\qquad {d\over dt}U_{2i}= {d^2\over dt^2}X_i=F_{G,i}(X_1,\dots,X_N)$$

% Dans   cet  \href{http://arxiv.org/abs/math/0011268}{article},   des  conditions
% initiales permettent  d'obtenir une solution  périodique pour des corps  de même
% masse.

% Si on se place  dans le cas d'un mouvement plan les termes  $U_i$ sont des vecteurs
% de $\bbr^2$ contenant les positions $x_i,y_i$ et les vitesses $vx_i,vy_i$ des 3 corps :

% $$
% U_1=\begin{pmatrix}x_1\\y_1\end{pmatrix},~~U_3=\begin{pmatrix}x_2\\y_2\end{pmatrix},~~U_5=\begin{pmatrix}x_3\\y_3\end{pmatrix}~~~~  
% U_2=\begin{pmatrix}vx_1\\vy_1\end{pmatrix},~~U_4=\begin{pmatrix}vx_2\\vy_2\end{pmatrix},~~U_6=\begin{pmatrix}vx_3\\vy_3\end{pmatrix}
% $$



Nous allons terminer la résolution du problème des N corps de Centrale 2015 avec
$N=3$ et les conditions initiales suivantes:

$$
x_1=−x_2=0.97000436~~~y_1=-y_2=−0.24308753,~~~x3=y_3=0
$$

$$
vx_3=−2vx_1=−2vx_2=−0.93240737~~~vy_3=−2vy_1=−2vy_2=−0.86473146
$$


avec des unités telles que les masses et la constante $\cal G$ valent 1.

%A-t-on un schéma périodique? Quelle est la période?

On rappelle que la PFD donne $$\ve{F}_{k/j}=G \frac{m_jm_k}{r_{jk}^3}\ve{P_jP_k}$$

avec $r_{ij}$ la distance entre les corps $j$ et $k$.

On  utilisera  des  tableaux  de  \verb+numpy+  pour  effectuer  des  opérations
vectorielles.

\begin{pythoncode}
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

G = 1
m = 1
position = np.array([[0.97000436,-0.24308753],[-0.97000436,0.24308753],[0,0]])
#position2 = np.array([[0.99000436,-0.24308753],[-0.97000436,0.24308753],[0,0]])
vitesse = np.array([[ 0.46620368,  0.43236573],[ 0.46620368,  0.43236573],[-0.93240737,-0.86473146]])

def norme(v) :
    """ norme d'un vecteur"""
    return 

def force2(m1,p1, m2, p2) :
    """ force exercée par le corps 2 sur le corps 1"""
    p1p2 =  
    return 

def force(j,m,pos) :
    """ somme des forces exercées sur le corps j
        on pourra utiliser np.sum(array,axis=0)"""
    return

def gravit(pos) :
    """renvoie le vecteur dont les composantes sont les accélérations de chaque corps pour une position donnée """
    return 

def verlet(f, pos0, vit0, h, n) :
    """ on utilise le schéma de Verlet vu en cours vectorialisé"""

def phase(f, pos0, vit0, h, n, no) :
    """ trace les mouvements des trois corps"""
\end{pythoncode}


On obtient une courbe fermée:


%\begin{figure}
\begin{center}
  \includegraphics[height=5cm]{../Python/TroisCorps0}
\end{center}
%\end{figure}

qui ne l'est plus si une valeur initiale est légèrement perturbée (1\%):


%\begin{figure}
\begin{center}
  \includegraphics[height=5cm]{../Python/TroisCorps1}
\end{center}
%\end{figure}

et si une valeur initiale est un peu plus perturbée (10\%):


%\begin{figure}
\begin{center}
  \includegraphics[height=5cm]{../Python/TroisCorps2}
\end{center}
%\end{figure}


\end{exercice}















\begin{exercice}[CCP MP : sujet 0 de SΙ 2015]
{\itshape 
Après  25 questions,  voici l'éternelle,  l'inévitable question  sur la  méthode
d'Euler...

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{exercice}


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