#!/usr/bin/env python
# coding: utf-8

# In[1]:


import numpy as np
import scipy as sp

import matplotlib.pyplot as plt
from scipy.integrate import odeint
from numpy.linalg import inv


# $\newcommand{\R}{\mathbb{R}}$
# $\newcommand{\N}{\mathbb{N}}$
# 
# #  TP 3 : Schémas explicites pour les EDO
# 
# L'objectif de ce TP est de mettre en oeuvre différents schémas numériques pour approcher la solution d'une équation différentielle ordinaire. Nous allons considérer différents schémas numériques, dont certains  nous n'avons pas encore vus en cours. L'objectif est de comprendre ce que c'est une solution approchée donnée par une méthode numérique, mais aussi de comprendre des notions comme l' erreur entre la solution exacte et la solution approchée, la  convergence d'une méthode numérique et l' ordre de précision d'une méthode numérique
# 
# 
# 
# Dans ce TP on considère des équations différentielles ordinaires, ou des systèmes d'équations différentielles ordinaires, de la forme
# \begin{equation}
# \label{EDO}
# \left\{\begin{aligned}
#   y'(t) & =  f(t,y(t)), \\
#   y(t_0) & =  y_0,
# \end{aligned} \right.
# \tag{EDO}
# \end{equation}
# 
# où $f$ : $I \times \R^n \longrightarrow \R^n$ est une
# fonction dans les conditions du théorème de Cauchy-Lipschitz, avec $I\subseteq\R$ un intervalle ouvert de $\R$ et où $(t_0,y_0)\in I\times\R^n$ est donné. Le problème de Cauchy (EDO) admet alors une unique solution maximale définie dans un intervalle ouvert $J\subseteq I$.
# 
# 
# On souhaite calculer une solution approchée de ce problème dans un intervalle de la forme $[t_0,t_f]=[t_0,t_0+T]\subseteq J$, avec $T>0$. Pour ce faire, on se donne $N\in\N$ et on considère une discrétisation de $[t_0,t_0+T]$ de pas $h=\frac TN,$ donnée par les points 
# $$
# t_0< t_1< \cdots <t_N=t_0+T,\ \ \ t_n=nh,\ n=0,\dots,N,
# $$
# et l'on utilise un schéma numérique afin de calculer une approximation $y^n$ de $y(t_n)$ pour $n=0,\dots,N$. 
# 
# Les schémas explicites qui suivent sont fréquemment utilisés pour discrétiser le problème (EDO). On pose $y^0=y_0,\,$ et, pour $n=0,\dots,N-1,$ 
# 
# 
# - **Euler explicite** :   $y^{n+1}=y^{n}+h f(t_n,y^{n})$ où   $f^{n}=f(t_n,y^{n})$ ;
# 
# -  **Heun** :  $y^{n+1}=y^{n}+\frac{h}{2}\big(f(t_n,y^{n})+f\big(t_n+h,y^{n}+h f(t_n,y^{n})\big)\big)$ où $p_1=f(t_n,y^{n})$ et $p_2=f(t_n+h,y^{n}+h p_1)$ ;
# 
# 
# 

# ### Exercice 1
# **Q1)** Écrire une fonction python de la forme `euler_exp(fct, t0, tf, y0, h)`
# 
# prenant en argument  `fct` la fonction $f(t,y)$ de (EDO), les extrémités `t0` et `tf` de l'intervalle de temps $[t_0,t_f]=[t_0,t_0+T],$ la donnée initiale `y0` et le pas de temps `h`. Cette fonction devra retourner deux tableaux :  
# 
# - `[t_0,\, t_1,\, ...,\, t_N]`, tableau  unidimensionnel de taille $(N+1)\times 1$ représentant la subdivision de l'intervalle $[t_0,t_f]$ de pas $h$ considérée, 
# 
# - `[y^0,\, y^1,\, ...,\, y^N]` tableau de taille $(N+1)\times n$ représentant la solution approchée aux instants $t_n,\ n=0,\dots,N.$

# In[2]:


# <completer>




# **Q2)** Testez d'abord votre fonction `euler_exp` sur un exemple simple d'une EDO scalaire. Vous pouvez par exemple considérer le problème
# \begin{equation}\label{edoexp}
# \begin{cases}
# y'(t)=y(t),\\
# y(0)=1,
# \end{cases}
# \end{equation}
# dont la solution est $y(t)=e^t.$ 
# 
# **Q2.a)** Tracer sur une même fenêtre graphique :
# 
# -  La solution exacte sur l'intervalle $[0,1]$ discrétisé avec un pas de temps de $10^{-4}$ (autrement dit, en remplaçant l'intervalle $[0,1]$ par l'ensemble discret correspondant aux points d'une sub-division de $[0,1]$ de pas $10^{-4}$) ;
# 
# - Les solutions approchées sur le même intervalle obtenues avec la méthode d'Euler avec des pas de temps $h$ égaux à $1/4,\ 1/10,\ 1/100.$ Observer que lorsque le pas diminue, la solution approchée est  *plus proche* de la solution exacte.

# In[3]:


# <completer>



# **Q2.b)** Donner les valeurs approchées de $y(1)=e$ obtenues avec $h=1/10,\ 1/100,\ 1/1000$, et donner les erreurs entre la valeur exacte de $e$ et les valeurs approchées obtenues avec ces pas de temps. Observer la décroissance de l'erreur lorsque le pas de temps diminue.
# 
# **Q2.c)** Tracer, sur une autre figure, la différence en valeur absolue entre la solution exacte et la solution approchée, en fonction du temps **(i.e. $|y(t_n)-y^n|$ en fonction de $t_n$)**, en utilisant des discrétisations de $[0,1]$ avec des pas de temps $h$ égaux à $1/10,\ 1/100,\ 1/1000.$ **Vous devez alors tracer trois courbes et observer que lorsque le pas diminue, cette différence s'approche de plus en plus de 0.**

# In[47]:


# <completer>


# **Q3)** Testez votre fonction `euler_exp` sur le modèle logistique
# $$
# (P_1)\ \ \ \ 
# \begin{cases}
# y'(t)=c y (1 - \frac{y}{b}),\\
# y(0)=a,
# \end{cases}
# $$
# dont la solution exacte est
# $$
# y(t) = \frac{b}{1 + \frac{b-a}{a} e^{-ct}},
# $$
# avec les données $c=1$, $b=2$, $a=0.1$, dans l'intervalle $[0,15]$, avec un pas $h=0.2$. Tracer sur la même fenêtre la solution exacte et la solution approchée, obtenue avec le pas $h$ et avec un pas égal à $\frac h2$.

# In[48]:


# <completer>


# **Q4)** Testez ensuite votre fonction dans le cas vectoriel $n>1$ ($n=2$) sur le problème
# $$
# (P_2)\ \ \ \
# \begin{cases}
# y''(t)  =    -y(t) + \cos(t)  \\
# y(0)  =   5,\ \ y'(0)  =   1,  
# \end{cases}
# $$
# dont la solution exacte est 
# 
# $$
# y(t) = \frac{1}{2} \sin(t) t  + 5 \cos(t) + \sin(t),
# $$
# 
# dans l'intervalle $[0,15]$, avec un pas  $h=0.2$.\\
# Pour ce faire, il faut écrire l'équation d'ordre 2 de $(P_2)$ sous la forme d'un système de deux équations d'ordre 1 dans les nouvelles variables $u(t)=y(t)$ et $v(t)=y'(t)$. On se ramènera alors à la résolution d'une équation de la forme
# 
# $$
# X'=F(t,X),\ \ \ \textrm{avec}\ X=(u,v)=(y,y')^T.
# $$
# 
# Représenter à nouveau la solution exacte et la solution approchée dans une même fenêtre graphique.
# 
# **Remarque :** La solution $y$ de $(P_2)$ correspond à la première composante du vecteur $X$ ci-dessus. Votre fonction `euler_ex`! retournera dans ce cas, si vous avez respecté la structure conseillée, un tableau de taille $(N+1)\times2$, $N$ étant le nombre de points de la discrétisation. Ce tableau donne les valeurs approchées de $X$ au points de la discrétisation considérée, la première colonne correspondant à la première composante de $X$, la seconde à la seconde composante de $X$. La solution approchée de $(P_2)$ que l'on cherche correspond alors à la première colonne de ce tableau.

# In[49]:


# <completer>


# ### Exercice 2 : Étude de l'erreur et convergence du schéma d'Euler explicite
# 
#  Dans cet exercice on essaie de comprendre ce que c'est l'erreur associé à un schéma numérique et ce que c'est un schéma numérique convergent.
#  
# On se donne un problème de Cauchy de la forme (EDO) et une méthode numérique pour approcher la solution de (EDO) dans un intervalle de la forme $[t_0,t_0+T]$. Pour un certain pas de temps $h$ fixé (ou pour un certain nombre de points de la discrétisation $N$ fixé), l'erreur globale entre la solution approchée associée à une discrétisation de pas $h$ de l'intervalle $[t_0,t_0+T]$ et la solution exacte est donnée par :
# 
# \begin{equation*}
# E_h =  \max_{k=0,\cdots,N}( |y(t_k)-y^{k}|). 
# \end{equation*}
# 
# Ci dessus, $y(t_k)$ est la solution exacte à l'instant $t_k$ et $y^{k}$ la valeur approchée de $y(t^k),$ donnée par le schéma numérique.
# 
# **Remarque 1 :**  l'erreur globale $E$ dépend du pas $h$, ou, de manière équivalente, du nombre de points $N$ de la discrétisation. Pour des valeurs de $N$ différentes, ou de manière équivalente pour des valeurs de $h$ différentes, les  discrétisations de l'intervalle $[t_0,t_0+T]$ sont différentes (elles ont un nombre de points différent), et les solutions approchées sont différentes. On s'attend à que, lorsque $N$ augmente, ou de manière équivalente lorsque $h$ diminue, l'erreur $E_h$ diminue, puisque l'on considère dans ce cas de plus en plus de points dans la discrétisation de l'intervalle $[t_0,t_0+T]$ et les approximations que l'on a faites pour construire le schéma numérique deviennent alors de plus en plus précises.
# 
# **Remarque 2 :** Parfois dans la littérature on ne spécifie pas la dépendance de $E$ par rapport à $h$, mais on doit toujours garder à l'esprit cette dépendance et que l'erreur est donc une fonction de $h$

# Considérons le problème
# $$
# (P_3)\ \ \ \
# \begin{cases}
# y'(t)=\frac{\cos(t)-y(t)}{1+t},\\
# y(0)=-\frac14, 
# \end{cases}
# $$
# dont la solution exacte est 
# $$
# y(t) = \frac{\sin(t)-1/4}{1 + t}.
# $$
# 
# **Q1)** Calculez les solutions approchées de $(P_3)$ obtenues avec le schéma d'Euler explicite, avec $h=1/2^s$ pour $s = 1,2,...,8$ ; représentez dans la même figure la différence en valeur absolue entre la solution exacte et la solution approchée, en fonction du temps, pour chaque valeur de $h$.

# In[50]:


# <completer>


# **Q2)** Calculer, pour chaque valeur de $h=1/2^s,\ s = 1,2,...,8,$ l'erreur globale $E_h$ correspondante. Vérifiez que $E_h$ diminue et s'approche de $0$ lorsque $h$ diminue.
# 
# **On dira qu'une méthode numérique converge si $\displaystyle{\lim_{h\to0}E_h}=0$.**
# 
# Représentez ensuite, en échelle logarithmique, l'erreur en fonction du pas de temps $h$, autrement dit, représentez $\log(E_h)$ en fonction de $\log(h)$. Vous devez obtenir des points qui sont à peu près alignés sur une droite de pente 1. Vérifiez graphiquement que c'est le cas, en estimant la pente de la droite passant au plus prêt des points (ou en représentant une droite de pente $1$ qui passe par un des points et en vérifiant que tous les points sont à peu près sur cette droite).
# 
# *Ceci signifie que $\log(E_h) \sim C+\log(h)$ et donc que $E_h\sim \widetilde{C}h$, pour certaines constantes $C$ et $\widetilde{C}$. On dit que la méthode d'Euler explicite est d' ordre 1 : c'est l'ordre de la puissance de $h$ dans cette relation. On a donc que l'erreur globale $E_h$ tend vers 0 comme $h$. L'ordre d'une méthode donne une indication sur sa vitesse de convergence. Une méthode d'ordre $p$ est une méthode dont l'erreur globale tend vers $0$ comme $h^p$. Donc plus l'ordre est élevé, plus la méthode converge plus vite*
# 
# **Remarque :** pour étudier numériquement l'ordre d'une méthode, on utilise souvent l'échelle logarithmique pour tracer l'erreur en fonction du pas de discrétisation $h$. La pente de la droite obtenue donne l'ordre $p$ de la méthode : si $E_h \sim Ch^p$ alors $\log(E_h)\sim \log(C) + p\log(h)$.

# In[51]:


# <completer>


# ### Exercice 3 : Un schéma d'ordre 2
# 
# **Q1)** À l'image de ce que vous avez fait pour le schéma d'Euler explicite, écrire une fonction python de la forme `Heun(fct, t0, tf, y0, h)` prenant les mêmes arguments que la fonction `euler_exp` et retournant les mêmes tableaux. 

# In[52]:


# <completer>


# **Q2)** Pour le problème $(P_1),$ comparez les solutions approchées obtenues avec la méthode d'Euler explicite et avec la méthode de Heun avec le même pas de temps. Pour ce faire, représenter dans la même fenêtre graphique la solution exacte et les solutions approchées obtenues avec chacune des méthodes. Quelle méthode semble être plus précise ? Pour le confirmer, représenter dans une autre figure la différence, en valeur absolue, entre la solution exacte et la solution approchée, pour les deux méthodes.

# In[53]:


# <completer>


# **Q3)** Reprenez l'exercice **2** pour la méthode de Heun. Vérifiez graphiquement que cette méthode est d'ordre $2,$ en représentant l'erreur globale en fonction du pas de discrétisation en échelle logarithmique et en estimant la pente de la droite passant au plus prêt des points correspondant à cette représentation ou en représentant une droite de pente $2$ qui passe par un des points.

# In[54]:


# <completer>


# In[ ]:





# In[ ]:





# In[ ]: