from IPython.display import display, Latex
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
import sys #only needed to determine Python version number
print('Python version ' + sys.version)
M62_TD3 : Révisions¶
Exercice : construction d'un schéma¶
Considérons le problème de Cauchy:
trouver $y \colon [t_0,T]\subset \mathbb{R} \to \mathbb{R}$ tel que $$\begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in [t_0,T],\\ y(t_0) = y_0. \end{cases}$$ Supposons que l'on ait montré l'existence d'une unique solution $y$.
On subdivise l'intervalle $[t_0,T]$ en $N$ intervalles de longueur $h=(T-t_0)/N=t_{n+1}-t_n$. Pour chaque nœud $t_n=t_0 + nh$ ($1 \le n \le N$) on cherche la valeur inconnue $u_n$ qui approche $y(t_n)$. L'ensemble de $N+1$ valeurs $\{u_0 = y_0, u_1,\dots, u_{N} \}$ représente la solution numérique.
Dans cette exercice on va construire des nouveaux schémas numériques basés sur l'intégration de l'EDO $y'(t)=\varphi(t,y(t))$ entre $t_n$ et $t_{n+2}$: $$ y(t_{n+2})=y(t_n)+\int_{t_n}^{t_{n+2}} \varphi(t,y(t))\mathrm{d}t. $$
- En utilisant la formule de quadrature du point milieu pour approcher le membre de droite écrire un schéma numérique explicite permettant de calculer $u_{n+2}$ à partir de $u_{n+1}$ et $u_{n}$. Notons que ce schéma a besoin de deux valeurs initiales; on posera alors $u_0=y_0$ et $u_{1}$ sera approché par une prédiction d'Euler progressive.
- En utilisant la formule de quadrature de Cavalieri-Simpson pour approcher le membre de droite écrire un schéma numérique implicite permettant de calculer $u_{n+2}$ à partir de $u_{n+1}$ et $u_{n}$. Notons que ce schéma a besoin de deux valeurs initiales; on posera alors $u_0=y_0$ et $u_{1}$ sera approché par une prédiction d'Euler progressive.
- Proposer une modification du schéma au point précédent pour qu'il devient explicite.
Rappels: une formule de quadrature interpolatoire approche l'intégrake d'une fonction $f$ par l'intègrale dun polynôme qui interpole $f$ en des points donnés: $\int_a^b f(t)\mathrm{d}t\approx\int_a^b p(t)\mathrm{d}t$
- la formule de quadrature du point milieu considère comme points d'interpolation le point $\frac{a+b}{2}$ ainsi $$ \int_a^b f(t)\mathrm{d}t\approx\int_a^b f\left(\frac{a+b}{2}\right)\mathrm{d}t = (b-a)f\left(\frac{a+b}{2}\right) $$
- la formule de Simpson considère comme points d'interpolation les points $\left\{a,\frac{a+b}{2},b\right\}$ ainsi $$ \int_a^b f(t)\mathrm{d}t\approx \frac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right] $$
Correction¶
Si on utilise la formule de quadrature du point milieu sur l'intervalle $[t_n;t_{n+2}]$, ie $$ \int_{t_n}^{t_{n+2}} \varphi(t,y(t))\mathrm{d}t\approx 2h \varphi\left(t_{n+1},y(t_{n+1})\right) $$ on obtient $$\begin{cases} u_0=y(t_0)=y_0,\\ u_1\approx y(t_1)\\ u_{n+2}=u_{n}+2h \varphi(t_{n+1},u_{n+1})& n=0,1,\dots N-2 \end{cases}$$ où $u_{1}$ est une approximation de $y(t_1)$. Nous pouvons utiliser une prédiction d'Euler progressive pour approcher $u_1$. Nous avons construit ainsi un nouveau schéma explicite à 2 pas: $$\begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_{0},u_{0}),\\ u_{n+2}=u_{n}+2h \varphi(t_{n+1},u_{n+1})& n=0,1,\dots N-2 \end{cases}$$
Correction¶
Si on utilise la formule de quadrature de Cavalieri-Simpson sur l'intervalle $[t_n;t_{n+2}]$, ie $$ \int_{t_n}^{t_{n+2}} \varphi(t,y(t))\mathrm{d}t\approx \frac{h}{3} \left( \varphi(t_{n},y(t_{n}))+4\varphi(t_{n+1},y(t_{n+1}))+\varphi(t_{n+2},y(t_{n+2})) \right), $$ et avec une prédiction d'Euler explicite pour approcher $u_1$, nous obtenons un nouveau schéma implicite à 3 pas: $$\begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_{0},u_{0}),\\ u_{n+2}=u_{n}+\frac{h}{3} \left( \varphi(t_{n},u_{n})+4\varphi(t_{n+1},u_{n+1})+\varphi(t_{n+2},u_{n+2}) \right)& n=0,1,\dots N-2 \end{cases}$$
Correction¶
Pour éviter le calcul implicite de $u_{n+2}$, nous pouvons utiliser une technique de type predictor-corrector et remplacer le $u_{n+2}$ dans le terme $\varphi(t_{n+2},u_{n+2})$ par exemple par $\tilde u_{n+2}=u_{n+1}+h \varphi(t_{n+1},u_{n+1})$. Nous avons construit ainsi un nouveau schéma explicite à 3 pas: $$\begin{cases} u_0=y(t_0)=y_0,\\ u_1=u_0+h\varphi(t_{0},u_{0}),\\ u_{n+2}=u_{n}+\frac{h}{3} \left( \varphi(t_{n},u_{n})+4\varphi(t_{n+1},u_{n+1})+\varphi(t_{n+2},u_{n+1}+h \varphi(t_{n+1},u_{n+1}) \right)& n=0,1,\dots N-2 \end{cases}$$
Exercice : EDO d'ordre 2¶
Considérons le problème de Cauchy
$$
\begin{cases}
y''(t)= ty(t), &t\in[0;4.5]\\
y(0)=0.355028053887817,\\
y'(0)=-0.258819403792807.
\end{cases}
$$
Transformer l'EDO d'ordre 2 en un système de deux EDO d'ordre 1.
Calculer la valeur approchée de $y$ en $t=4.5$ obtenue avec le schéma RK4 avec un pas de discrétisation $h=0.001$.
Correction¶
On appelle $x(t)=y(t)$ et $z(t)=y'(t)$, alors on a $x'(t)=y'(t)=z(t)$ et $z'(t)=y''(t)=ty(t)=tx(t)$ donc
$$ \begin{cases} x'(t)= z(t)\\ z'(t)=tx(t)\\ x(0)=0.355028053887817,\\ z(0)=-0.258819403792807. \end{cases} $$def RK4(phi1,phi2,tt,y0,yp0):
uu1 = [y0]
uu2 = [yp0]
for i in range(len(tt)-1):
k1_1 = phi1( tt[i] , uu1[i] , uu2[i] )
k1_2 = phi2( tt[i] , uu1[i] , uu2[i] )
k2_1 = phi1( tt[i]+h/2 , uu1[i]+h*k1_1/2, uu2[i]+h*k1_2/2 )
k2_2 = phi2( tt[i]+h/2 , uu1[i]+h*k1_1/2, uu2[i]+h*k1_2/2 )
k3_1 = phi1( tt[i]+h/2 , uu1[i]+h*k2_1/2, uu2[i]+h*k2_2/2 )
k3_2 = phi2( tt[i]+h/2 , uu1[i]+h*k2_1/2, uu2[i]+h*k2_2/2 )
k4_1 = phi1( tt[i+1] , uu1[i]+h*k3_1 , uu2[i]+h*k3_2 )
k4_2 = phi2( tt[i+1] , uu1[i]+h*k3_1 , uu2[i]+h*k3_2 )
uu1.append( uu1[i] + (k1_1+2.0*k2_1+2.0*k3_1+k4_1)*h/6 )
uu2.append( uu2[i] + (k1_2+2.0*k2_2+2.0*k3_2+k4_2)*h/6 )
return [uu1,uu2]
phi1 = lambda t,y1,y2 : y2
phi2 = lambda t,y1,y2 : t*y1
t0=0.0
y0=0.355028053887817
yp0=-0.258819403792807
h=0.001
tt=[t0+i*h for i in range(4501)]
[uu1,uu2]=RK4(phi1,phi2,tt,y0,yp0)
print("u(%1.2f)=%1.16f"%(tt[-1],uu1[-1]))
Exercice : Étude d'un problème numériquement mal posé¶
Considérons le problème de Cauchy $$ \begin{cases} y'(t)=3y(t)-4e^{-t}, & t\in[0;1],\\ y(0)=1+\varepsilon. \end{cases} $$
- Calculer la solution exacte.
- Utiliser la méthode RK4 pour calculer la solution approchée avec un pas $h=\frac{1}{10}$ et comparer à la solution exacte pour trois valeurs des $\varepsilon$, à savoir $\varepsilon=1$, $\varepsilon=\frac{1}{10}$ et $\varepsilon=0$.
- Même exercice avec la méthode d'Euler explicite.
- Commenter les résultats.
Correction : solution exacte (avec SymPy)¶
La solution exacte est $y(t)=\varepsilon e^{3t}+e^{-t}$.
Si $\varepsilon=0$, la solution exacte devient donc $y(t)=e^{-t}$ mais le problème est mal conditionné car toute (petite) erreur de calcul a le même effet qu'une perturbation de la condition initiale: on "réveil" le terme dormant $e^{3t}$.
%reset -f
import sympy as sym
sym.init_printing(pretty_print=True)
sym.var('t,epsilon')
y=sym.Function('y')
edo=sym.Eq( sym.diff(y(t),t), 3*y(t)-4*sym.exp(-t) )
print("EDO")
display(edo)
print("\nSolution generale")
soln=sym.dsolve(edo,y(t))
display(soln.expand())
print("\nEn prenant en compte la CI on trouve")
t0 = 0
y0 = 1+epsilon
constants = sym.solve([soln.rhs.subs(t,t0) - y0])
display(constants)
print("\nSolution particuliere")
sym.var('C1')
solp = soln.subs(constants)
display(solp.expand())
%matplotlib inline
from matplotlib.pylab import *
func = sym.lambdify((t,epsilon),solp.rhs,'numpy')
tt=linspace(0,1,101)
figure(figsize=(18,5))
EPSILON=[1,1.e-1,0]
for i,epsi in enumerate(EPSILON):
subplot(1,3,i+1)
plot(tt,func(tt,epsi))
title("$\epsilon=$"+str(epsi));
Correction : solution approchée (avec RK41 et EE)¶
%reset -f
%matplotlib inline
from matplotlib.pylab import *
phi = lambda t,y : 3*y-4*exp(-t)
sol_exacte = lambda t,y0 : (y0-1)*exp(3*t)+exp(-t)
def EE(phi,tt,y0):
h = tt[1]-tt[0]
uu = [y0]
for i in range(N):
uu.append( uu[i] + h*phi(tt[i],uu[i]) )
return uu
def RK4(phi,tt,y0):
h = tt[1]-tt[0]
uu = [y0]
for i in range(N):
k1 = phi( tt[i], uu[i] )
k2 = phi( tt[i]+h/2 , uu[i]+h*k1/2 )
k3 = phi( tt[i]+h/2 , uu[i]+h*k2/2 )
k4 = phi( tt[i+1] , uu[i]+h*k3 )
uu.append( uu[i] + (k1+2*k2+2*k3+k4)*h/6 )
return uu
# INITIALISATION
t0 = 0
tfinal = 1
# DISCRETISATION
N = 10
h = (tfinal-t0)/N
tt = [ t0+i*h for i in range(N+1) ]
EPSILON=[1,1.e-1,0]
figure(figsize=(18,5))
for i,epsilon in enumerate(EPSILON):
y0 = 1+epsilon
yy = [sol_exacte(t,y0) for t in tt]
uu_RK4 = RK4(phi,tt,y0)
uu_EE = EE(phi,tt,y0)
subplot(1,3,i+1)
plot(tt,yy,'b-',tt,uu_RK4,'y:D',tt,uu_EE,'m:o')
legend(['Exacte','RK4','EE']);
title("$\epsilon=$"+str(epsilon))
Exercice : Étude théorique et numérique du mouvement d'un pendule¶
Considérons le problème du mouvement d'un pendule de masse $m$ suspendu au point $O$ par un fil non pesant de longueur $\ell$.
On se propose d'étudier l'évolution au cours du temps $t\ge0$ de l'angle $\vartheta(t)$ autour du point $0$.
L'angle $\vartheta(t)$ est mesuré par rapport à une verticale passante par $0$.
Considérons pour simplifier seulement la force de gravité $g$.
La fonction $\vartheta$ est alors solution de l'EDO d'ordre 2 $$ \vartheta''(t)=-\frac{g}{\ell}\sin(\vartheta(t)). $$ Si on introduit le paramètre $\omega$ tel que $\omega^2=\dfrac{g}{\ell}$, l'équation devient $$ \vartheta''(t)+\omega^2\sin(\vartheta(t))=0. $$ Pour simplicité on pose $\omega=1$ et on note $q=\vartheta$ et $p=\vartheta'$. On peut alors transformer l'EDO d'ordre 2 en un système de deux EDO d'ordre 1: $$ \begin{cases} q'(t)=p(t),\\ p'(t)=-\sin(q(t)). \end{cases} $$ Considérons l'énergie totale du système: $$ E(q,p)=\underbrace{-\cos(q)}_{\text{Énergie potentielle}}+\underbrace{\frac{1}{2}p^2}_{\text{Énergie cinétique}}. $$
Étude théorique
- Montre analytiquement que l'énergie totale reste constante au cours du temps, i.e. que $\frac{\mathrm{d}}{\mathrm{d}t}E(q(t),p(t))=0$.
- Tracer avec Python les courbes de niveau de la fonction $(q,p)\mapsto E(q,p)$, i.e. les ensembles $$ \mathcal{N}_\kappa =\{(q,p)\ |\ E(q,p)=\kappa\} =\{(q,p)\ |\ p^2=2(\kappa+\cos(q))\}. $$ pour $(q,p)\in[-\pi;\pi]\times\mathbb{R}$.
- Lorsque $\vartheta$ est petit on peut considérer l'approximation $\sin(\vartheta)\simeq\vartheta$. Calculer la solution exacte de cette équation approchée.
Étude numérique
Dans une simulation numérique on aimerait que la propriété de conservation de l'énergie soit préservée aussi bien que possible. Nous allons regarder ce qui se passe avec certaines méthodes vues en cours.On notera $x_n\approx q(t_n)=\vartheta(t_n)$ et $y_n\approx p(t_n)=\vartheta'(t_n)$.
À chaque instant $t_n$, on calculera les valeurs approchées de l'énergie cinétique $E_c(t_n)=\frac{1}{2}y_n^2$, de l'énergie potentielle $E_p(t_n)=-\cos(x_n)$ et de l'énergie totale $E(t_n)=E_c(t_n)+E_p(t_n)$.
Choisissons $h=t_{n+1}-t_n=0.1$ et $(x_0,y_0)=(\pi/4,0)$.- Dans un premier temps on se propose d'appliquer la méthode d'Euler explicite à la résolution du système.
- Tracer $(t_n,x_n)$ et $(t_n,y_n)$ sur un même graphique;
- sur un autre graphique représenter les trois énergies simultanément en fonction du temps;
- tracer enfin $(x_n,y_n)$ sur un autre graphique et vérifier que la solution numérique tourne vers l’extérieur (i.e. l'énergie augmente au cours du temps).
- Considérons le schéma obtenu en écrivant le schéma d'Euler explicite pour la première équation et le schéma d'Euler implicite pour la seconde.
- Montrer que ce schéma peut s'écrire sous forme parfaitement explicite;
- tracer $(t_n,x_n)$ et $(t_n,y_n)$ sur un même graphique;
- sur un autre graphique représenter les trois énergies simultanément en fonction du temps;
- tracer enfin $(x_n,y_n)$ sur un autre graphique. Que peut-on constater? Commenter les résultats.
- On se propose d'appliquer les méthodes d'Euler modifiée, de Heun, AB2, AB3 et RK4 à la résolution du système. Pour chaque méthode:
- Tracer $(t_n,x_n)$ et $(t_n,y_n)$ sur un même graphique;
- sur un autre graphique représenter les trois énergies simultanément en fonction du temps;
- tracer enfin $(x_n,y_n)$ sur un autre graphique. Que peut-on constater? Commenter les résultats.
- Dans un premier temps on se propose d'appliquer la méthode d'Euler explicite à la résolution du système.
Correction étude théorique¶
Il s'agit d'un système de deux EDO scalaires d'ordre 1 du type $$ \begin{cases} q'(t)=\varphi_1(t,q(t),p(t)),\\ p'(t)=\varphi_2(t,q(t),p(t)), \end{cases} \qquad\text{avec}\qquad \begin{cases} \varphi_1(t,q(t),p(t))=p(t),\\ \varphi_2(t,q(t),p(t))=-\sin(q(t)). \end{cases} $$
L'énergie totale reste constante au cours du temps: \begin{align*} \frac{\mathrm{d}}{\mathrm{d}t}E(q(t),p(t)) &= \sin(q(t))q'(t)+p(t)p'(t) {}= \sin(q(t))p(t)+p(t)p'(t) \\ &= p(t)\left( \sin(q(t))+p'(t) \right) {}= p(t)\left( \sin(q(t))-\sin(q(t)) \right)=0. \end{align*}
Courbes de niveau de la fonction $(q,p)\mapsto E(q,p)$:
from matplotlib.pylab import *
Etot = lambda q,p : -cos(q)+p**2/2
q,p = meshgrid(linspace(-8*pi,8*pi,201),linspace(-3,3,201))
E = Etot(q,p)
figure(1,figsize=(20,5))
graphe=contour(q,p,E,20)
clabel(graphe,inline=1);
- $\vartheta''(t)+\omega^2\vartheta(t)=0$ est une EDO d'ordre 2 homogène à coefficients constants dont la solution est $\vartheta(t)=\vartheta(0)\cos(\omega t)+\frac{\vartheta'(0)}{\omega}\sin(\omega t)$.
Correction étude numérique¶
t0 = 0
tfinal = 8*pi
y0 = [pi/4,0]
N = 100
tt = linspace(t0,tfinal,N+1)
#h = tt[1]-tt[0]
phi1 = lambda t,q,p : p
phi2 = lambda t,q,p : -sin(q)
Méthode d'Euler explicite¶
$$\begin{cases} x_0=q(0),\\ y_0=p(0),\\[0.5em] x_{n+1}=x_n+h\varphi_1(t_n,x_n,y_n),\\ y_{n+1}=y_n+h\varphi_2(t_n,x_n,y_n). \end{cases}$$def euler_progressif(phi1,phi2,tt,y0):
h = tt[1]-tt[0]
q = [y0[0]]
p = [y0[1]]
for i in range(len(tt)-1):
q.append(q[i]+h*phi1(tt[i],q[i],p[i]))
p.append(p[i]+h*phi2(tt[i],q[i],p[i]))
return [q,p]
[q_ep, p_ep] = euler_progressif(phi1,phi2,tt,y0)
Ec_ep = [p**2/2 for p in p_ep]
Ep_ep = [-math.cos(q) for q in q_ep]
Et_ep = [Ec_ep[i]+Ep_ep[i] for i in range(len(Ec_ep))]
print ('Euler explicite : |energie totale (t=Tfinal) - Energie totale (t=0)|=', abs(Et_ep[-1]-Et_ep[0]))
figure(figsize=(18,5))
subplot(131)
plot(tt,q_ep,tt,p_ep)
xlabel('t')
legend(['x(t)','y(t)'])
title('Euler progressif - x(t) et y(t)')
subplot(132)
plot(q_ep,p_ep)
xlabel('x(t)')
ylabel('y(t)')
title('Euler progressif')
subplot(133)
plot(tt,Ec_ep,tt,Ep_ep,tt,Et_ep)
xlabel('t')
legend(['Cinetique','Potentielle','Totale'])
title('Euler progressif - Energies');
Méthode d'Euler explicite/implicite¶
Schéma obtenu en écrivant le schéma d'Euler explicite pour la première équation et le schéma d'Euler implicite pour la seconde:
$$ \begin{cases} x_0=q(0),\\ y_0=p(0),\\[0.5em] x_{n+1}=x_n+h\varphi_1(t_n,x_n,y_n),\\ y_{n+1}=y_n+h\varphi_2(t_{n+1},x_{n+1},y_{n+1}). \end{cases} $$Dans notre cas $\varphi_2(t_{n+1},x_{n+1},y_{n+1})=-\sin(x_{n+1})$ donc le schéma est parfaitement explicite:
def euler_symplectique(phi1,phi2,tt,y0):
h = tt[1]-tt[0]
q = [y0[0]]
p = [y0[1]]
for i in range(len(tt)-1):
q.append(q[i]+h*phi1(tt[i],q[i],p[i]))
p.append(p[i]+h*phi2(tt[i+1],q[i+1],p[i])) # il faudrait p[i+1] (schema implicite) mais comme dans notre cas phi2 ne depend pas de "p", le schema est explicite
return [q,p]
[q_es, p_es] = euler_symplectique(phi1,phi2,tt,y0)
Ec_es = [p**2/2 for p in p_es]
Ep_es = [-math.cos(q) for q in q_es]
Et_es = [Ec_es[i]+Ep_es[i] for i in range(len(Ec_es))]
print ('Euler symplectique : |energie totale (t=Tfinal)-Energie totale (t=0)|=', abs(Et_es[-1]-Et_es[0]))
figure(figsize=(18,5))
subplot(131)
plot(tt,q_es,tt,p_es)
xlabel('t')
legend(['x(t)','y(t)'])
title('Euler symplectique - x(t) et y(t)')
subplot(132)
plot(q_es,p_es)
xlabel('x(t)')
ylabel('y(t)')
title('Euler symplectique - y(x)')
subplot(133)
plot(tt,Ec_es,tt,Ep_es,tt,Et_es)
xlabel('t')
legend(['Cinetique','Potentielle','Totale'])
title('Euler symplectique - Energies');
Méthodes d'Euler modifiée¶
$$ \begin{cases} x_0=q(0),\\ y_0=p(0),\\[0.5em] \tilde x_n=x_n+\frac{h}{2}\varphi_1(t_n,x_n,y_n),\\ \tilde y_n=y_n+\frac{h}{2}\varphi_2(t_n,x_n,y_n),\\[0.5em] x_{n+1}=x_n+h\varphi_1(t_n,\tilde x_n,\tilde y_n),\\ y_{n+1}=y_n+h\varphi_2(t_{n},\tilde x_{n},\tilde y_{n}). \end{cases} $$def euler_modifie(phi1,phi2,tt,y0):
h = tt[1]-tt[0]
q = [y0[0]]
p = [y0[1]]
for i in range(len(tt)-1):
qtilde = phi1( tt[i], q[i], p[i] )
ptilde = phi2( tt[i], q[i], p[i] )
q.append( q[i]+h*phi1(tt[i]+h/2.,q[i]+qtilde*h/2.,p[i]+ptilde*h/2.) )
p.append( p[i]+h*phi2(tt[i]+h/2.,q[i]+qtilde*h/2.,p[i]+ptilde*h/2.) )
return [q,p]
[q_em, p_em] = euler_modifie(phi1,phi2,tt,y0)
Ec_em = [p**2/2 for p in p_em]
Ep_em = [-math.cos(q) for q in q_em]
Et_em = [Ec_em[i]+Ep_em[i] for i in range(len(Ec_em))]
print ('Euler modifie : |energie totale (t=Tfinal)-Energie totale (t=0)|=', abs(Et_em[-1]-Et_em[0]))
figure(figsize=(18,5))
subplot(131)
plot(tt,q_em,tt,p_em)
xlabel('t')
legend(['x(t)','y(t)'])
title('Euler modifie - x(t) et y(t)')
subplot(132)
plot(q_em,p_em)
xlabel('x(t)')
ylabel('y(t)')
title('Euler modifie - y(x)')
subplot(133)
plot(tt,Ec_em,tt,Ep_em,tt,Et_em)
xlabel('t')
legend(['Cinetique','Potentielle','Totale'])
title('Euler modifie - Energies');
Méthodes de Heun¶
$$ \begin{cases} x_0=q(0),\\ y_0=p(0),\\[0.5em] \tilde x_n=x_n+h\varphi_1(t_n,x_n,y_n),\\ \tilde y_n=y_n+h\varphi_2(t_n,x_n,y_n),\\[0.5em] x_{n+1}=x_n+\frac{h}{2}\left( \varphi_1(t_n,x_n,y_n) + \varphi_1(t_{n+1},\tilde x_n,\tilde y_n)\right),\\ y_{n+1}=y_n+\frac{h}{2}\left( \varphi_2(t_n,x_n,y_n) + \varphi_2(t_{n+1},\tilde x_n,\tilde y_n)\right). \end{cases} $$def heun(phi1,phi2,tt,y0):
h = tt[1]-tt[0]
q = [y0[0]]
p = [y0[1]]
for i in range(len(tt)-1):
qtilde = phi1( tt[i], q[i], p[i] )
ptilde = phi2( tt[i], q[i], p[i] )
q.append( q[i]+h*0.5*( phi1(tt[i],q[i],p[i]) + phi1(tt[i+1],q[i]+h*qtilde,p[i]+h*ptilde) ) )
p.append( p[i]+h*0.5*( phi2(tt[i],q[i],p[i]) + phi2(tt[i+1],q[i]+h*qtilde,p[i]+h*ptilde) ) )
return [q,p]
[q_heun, p_heun] = heun(phi1,phi2,tt,y0)
Ec_heun = [p**2/2 for p in p_heun]
Ep_heun = [-math.cos(q) for q in q_heun]
Et_heun = [Ec_heun[i]+Ep_heun[i] for i in range(len(Ec_heun))]
print ('Heun : |energie totale (t=Tfinal)-Energie totale (t=0)|=', abs(Et_heun[-1]-Et_heun[0]))
figure(figsize=(18,5))
subplot(131)
plot(tt,q_heun,tt,p_heun)
xlabel('t')
legend(['x(t)','y(t)'])
title('Heun - x(t) et y(t)')
subplot(132)
plot(q_heun,p_heun)
xlabel('x(t)')
ylabel('y(t)')
title('Heun - y(x)')
subplot(133)
plot(tt,Ec_heun,tt,Ep_heun,tt,Et_heun)
xlabel('t')
legend(['Cinetique','Potentielle','Totale'])
title('Heun - Energies');
Méthode AB2¶
$$ \begin{cases} x_0=q(0),\\ y_0=p(0),\\[0.5em] x_1=x_0+h\varphi_1(t_0,x_0,y_0),\\ y_1=y_0+h\varphi_2(t_0,x_0,y_0),\\[0.5em] x_{n+1}=x_n+\frac{h}{2}\left( 3\varphi_1(t_n,x_n,y_n) - \varphi_1(t_{n-1},x_{n-1},y_{n-1})\right),\\ y_{n+1}=y_n+\frac{h}{2}\left( 3\varphi_2(t_n,x_n,y_n) - \varphi_2(t_{n-1},x_{n-1},y_{n-1})\right). \end{cases} $$def AB2(phi1,phi2,tt,y0):
h = tt[1]-tt[0]
q = [y0[0]]
p = [y0[1]]
q.append(q[0]+h*phi1(tt[0],q[0],p[0]))
p.append(p[0]+h*phi2(tt[0],q[0],p[0]))
for i in range(1,len(tt)-1):
qtilde1 = phi1( tt[i],q[i],p[i] )
qtilde2 = phi1( tt[i-1],q[i-1],p[i-1] )
ptilde1 = phi2( tt[i],q[i],p[i] )
ptilde2 = phi2( tt[i-1],q[i-1],p[i-1] )
q.append( q[i] + (3*qtilde1-qtilde2)*h/2 )
p.append( p[i] + (3*ptilde1-ptilde2)*h/2 )
return [q,p]
[q_AB2, p_AB2] = AB2(phi1,phi2,tt,y0)
Ec_AB2 = [p**2/2. for p in p_AB2]
Ep_AB2 = [-math.cos(q) for q in q_AB2]
Et_AB2 = [Ec_AB2[i]+Ep_AB2[i] for i in range(len(Ec_AB2))]
print ('AB2 : |energie totale (t=Tfinal)-Energie totale (t=0)|=', abs(Et_AB2[-1]-Et_AB2[0]))
figure(figsize=(18,5))
subplot(131)
plot(tt,q_AB2,tt,p_AB2)
xlabel('t')
legend(['x(t)','y(t)'])
title('AB2 - x(t) et y(t)')
subplot(132)
plot(q_AB2,p_AB2)
xlabel('x(t)')
ylabel('y(t)')
title('AB2 - y(x)')
subplot(133)
plot(tt,Ec_AB2,tt,Ep_AB2,tt,Et_AB2)
xlabel('t')
legend(['Cinetique','Potentielle','Totale'])
title('AB2 - Energies');