Partie 1 - Numériques pour équations différentielles ordinaires (EDO)

@c jean-baptiste.apoung@math.u-psud.fr

ATTENTION:

LES EXERCICES 1 et 2 SONT À FAIRE ENTIÈREMENT SUR COPIES

In [229]:
%matplotlib inline
%config InlineBackend.figure_format = 'retina'

import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt
In [230]:
%%javascript
MathJax.Hub.Config({
    TeX: { equationNumbers: { autoNumber: "AMS" } }
});

Exercice 1 : Question de cours

On considère l'équation différentielle ordinaire (edo) suivante où $a > 0$, $ \left \{ \begin{array}{l} x'(t) = - a \, x(t), \, ~ t ~ \in ~ ]0, \infty[,\\ x(0) = x^0 \end{array} \right. $

  • Q-1: Montrer sa solution exacte est $ x(t) = x^0 e^{- a t}$ et quelle vérifie $\displaystyle \lim_{t \rightarrow \infty} x(t) = 0 $.

  • Q-2 : Ecrire un schéma à un pas de votre choix pour la discrétisation de cette edo sur un maillage de pas $\Delta t$.

  • Q-3 : Etudier la convergence du schéma en précisant son ordre exacte de convergence. (On pourra le mettre sous la forme $x_{n+1} = x_n + \Delta t\, \Phi(t_n, x_n, \Delta t)$).

  • Q-4 : Mettre le schéma sous la forme $ x_{n+1} = r(a, \Delta t)\, x_n $. Et déduire son rayon de stabilité asymptotique: c'est-à-dire la plus grande valeur de $\Delta t $ assurant la propriété $\displaystyle \lim_{n\rightarrow \infty} x_n = 0$.

In [ ]:

Exercice 2 :

Soit $f : \mathbb{R} \rightarrow \mathbb{R}$ de classe $C^1$ dont on aimerait approcher une racine en partant d'un point $x^0$.

On pose g(t, x) la fonction définie sur $\mathbb{R}^+ \times \mathbb{R}$ par : $$ g(t, x) = f(x) - e^{-t} f(x^0) $$

Q-1 : Montrer que :

$g(0, x^0) = 0$ et que pour tout $~x\, \in \, \mathbb{R}$, on a $\displaystyle \lim_{t\rightarrow \infty} g(t, x) = f(x)$.

Q-2 : On définie pour tout t, x(t) comme solution de g(t, x(t)) = 0.

Q-2-1 : Montrer qu'on a pour tout $t$, $e^{-t} f(x^0) = f(x(t))$

Q-2-2 : En dérivant g(t,x(t)), montrer que $ x $ est solution de l'équation différentielle ordinaire

$$ (Eq) \left\{ \begin{array}{rll} x'(t) &=& - \frac{f(x(t))}{f'(x(t))} \quad t\in]0,\infty[\\ x(0) &=& x^0. \end{array} \right. $$

Q-2-3 : Montrer que si $\displaystyle \lim_{t\rightarrow \infty} x(t) = x^* $ existe, alors elle est racine de $f$ .

  • En déduire qu'une résolution de l'équation différentielle (Eq) pourra permettre de déterminer une racine de $f$.

Q-3 : Discrétisation de cette équation différentielle.

On se donne un pas de temps $\Delta t > 0$ et on désigne par $x_n$ la valeur approchée de $x(t_n)$ pour $t_n = n \Delta t, n \in \mathbb{N}$.

  • Ecrire le schéma d'Euler explicite pour la résolution numérique de cette équation différentielle.

Q-4 : Comparaison avec la méthode de Newton pour $f(x) = 0$

Q-4-1 : Comparer le schéma d'Euler explicite où $\Delta t = 1$ avec la méthode de Newton appliquée à $f(x) = 0$.

Q-4-2 : Donner un argument en faveur du schéma. (On pourra évoquer le cas des racines multiples )

Q-5 : Un algorithme de résolution

Q-5-1 : Test d'arrêt :

l'équation différentielle étant posée sur un domaine infini, en s'insiprant du test d'arrêt (à savoir $| f(x_k) | < \varepsilon $) dans la méthode de Newton, indiquer comment modifier le schéma d'Euleur explicite ci-dessus afin d'arrêter les calculs.

Q-5-2 : On propose l'algorithme ci-desous , dire s'il remplit les attentes :

On suppose (Eq) réécite sous la forme standard $x'(t) = F(t, x(t))$. Et on désigne par StepEuler la fonction implémentant un pas du schéma d'Euler explicite. A savoir $x_{n+1} = x_n + \Delta t \, F(t_n, x_n)$ est plutôt écrite $x_{n+1} = StepEuler(t_n, x_n, \Delta t, F)$ afin de permettre une généralisation à d'autre schémas.

#######################################################################
Approximation du zero de f par resolution de l' edo 
 x' = F(t, x(t)) où F(t,y) = -f(y)/f'(y)  
#######################################################################
Entree: F, t0, x0, dt, f, tol, iterMax
Sortie: x : valeur approchée de la racine de f
        t : valeur du temps au bout duquel on a pu approcher la solution  
        k : nombre d'itérations nécessaires pour ptoduire 
Initialisation:
--------------
  k = 0
  t = t0
  x = x0
Itérations : 
------------
  Tant que |f(x)| > tol  et  k < iterMax  
      x = StepEuler (t, x, dt, F )
      t = t + dt
      k = k + 1
  retourner x, t, k
#######################################################################
In [ ]:

Exercice 3 : Mise en oeuvre et validations

Q-3-1 : Mise en oeuvre

Mettre en oeuvre l'algorithme Q-5-2 ci-dessus à travers une fonction de prototype

In [285]:
def EDOSolveControle(StepMethod, F, t0, x0, dt, g, tol = 1e-6, maxIter = 500):
    '''
    ENTREES : 
    - StepMethod : est une méthode à un pas par exemple StepEuler, StepHeun, StepPM
                    StepMethod(tn, xn, dt, F)  calcule x_n+1
    - F          : est le second membre de de l'edo
    - t0         : l'instant initial
    - x0         : valeur initiale
    - dt         : le pas de temps
    - g          : la fonction de controle : on approche la racine de g
    - tol        : paramètre du test d'arret i.e  on arrête les calculs dès que  |g(x)| < tol
    - maxIter    : est le nombre maximum d'itérations autorisées
    SORTIES  :
      - x        : valeur approchée de la racine de g
      - t        : l'instant final auquel correspond la solution  x
      - k        : le nombre d'itérations réalisées
    '''
   # COMPLETER

Q-3-2 : Validations

On va essayer d'approcher la racine de l'équation $f(x) = 4 x^3 + 3 x^2 - 1$ sur $[-1, 1]$

Q-3-2-1 :

  • Fournir les $f, f'$ et $F$ respectivement à travers les fonctions Python : f, fp et F(t, x)

  • Et représenter graphiquement $f$ sur $[-1,1]$ et vérifier qu'elle y admet une racine

In [319]:
## METTRE LA REPONSE ICI


## definition des fonctions 



## représentation grapique

Q-3-2-2 : Fournir les fonctions StepEuler, StepHeun, StepPM

implémentant un pas respectivement des schémas d'Euler explicite, de Heun, et du point-milieu.

In [275]:
## METTRE LE CODE ICI

def StepEuler(t,x,dt,f):
   #COMPLETER


def StepPM(t,x,dt,f):
    #COMPLETER

    
def StepHeun(t,x,dt,f):
    #COMPLETER

Q-3-2-3 Evaluation de l'algorithme.

  • Afficher et commenter les résultats fournis par l'algorithme pour les trois méthodes de discrétisation de l'edo dans les cas tests listés ci-dessous. On pourra si nécessaire ajouter d'autres cas pertinents dans la comparaison des méthodes.

Aide : On pourra par exemple afficher les valeurs x,t,k retournées par EDOControl pour la méthode d'Euler explicite et un cas test avec l'instruction print(fmt.format("Euler", dt, x0, x, t, k, f(x))), le formattage fmt étant donné ci-dessous.

In [ ]:
fmt = "{:5s} : dt = {:2.2f}, x0 ={:2.2f} -> x = {:13.11e}  t = {:5.3e}  k = {:3d}  f(x) = {:e}"

Cas : $t_0 = 0, \, x_0 = 1.0, \, \Delta t = 1$

In [317]:
#### METTRE LA REPONSE ICI

t0, x0, dt  = 0,  1.0,  1.
In [ ]:
####  METTRE LES OBSERVATIONS  ICI

Cas : $t_0 = 0, \, x_0 = -1.0, \, \Delta t = 1$

In [315]:
#### METTRE LA REPONSE ICI

t0, x0, dt  = 0,  -1.,  1.
In [ ]:
####  METTRE LES OBSERVATIONS  ICI

Cas : $t_0 = 0, \, x_0 = 0.25, \, \Delta t = 2$

In [316]:
#### METTRE LA REPONSE ICI

t0, x0, dt  = 0,  0.25,  2.
In [ ]:
####  METTRE LES OBSERVATIONS  ICI