TD3 - Révisions

Auteur·rice

Gloria FACCANONI

Date de publication

25 mars 2024

1 Exercice : Étude d’un problème mathématiquement mal posé

Considérons le problème de Cauchy \( \begin{cases} y'(t)=y^2(t)+t^2, & t\in[0;1-\varepsilon],\\ y(0)=1. \end{cases} \)

Soit \(\varepsilon=0.2\). Utiliser la méthode d’Euler explicite puis la méthode RK4 pour calculer la solution approchée avec un pas \(h=0.1\) et comparer les deux résultats en \(t=1\). Même exercice avec \(h=0.05\) et \(h=0.01\).

Répéter avec \(\varepsilon=0.1\) puis \(\varepsilon=0\).

Commenter les résultats (hint: tracer le champ de vecteurs).

1.1 Correction

Code

%reset -f
%matplotlib inline

from matplotlib.pylab import *
# rcdefaults()
rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})
from scipy.optimize import fsolve

# ODE
phi = lambda t,y : y**2+t**2
    
# SCHEMAS
def EE(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        uu.append( uu[i] + h*phi(tt[i],uu[i]) )
    return uu

def EI(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        uu.append( fsolve(lambda x : -x + uu[i] + h*phi(tt[i+1],x) , uu[i] )[0] )
    return uu

def RK4(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        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
y0     = 1

# DISCRETISATION
for epsilon in [0.2, 0.1, 0]:
    print(f"Interval [0,{</span><span class="dv">1</span><span class="op">-</span>epsilon<span class="sc">}]")
    for h in [0.1,0.05,0.01]:
        N = int((1-epsilon)/h)
        tt = linspace(0,1-epsilon,N+1) # ne pas utiliser arange, pb avec les arrondis
        uu_EE = EE(phi,tt,y0)
        uu_EI = EI(phi,tt,y0)
        uu_RK4 = RK4(phi,tt,y0)
        h = tt[1]-tt[0]
        print(f"\tAvec h={</span>h<span class="sc">} on a \t EE({</span>tt[<span class="op">-</span><span class="dv">1</span>]<span class="sc">})={</span>uu_EE[<span class="op">-</span><span class="dv">1</span>]<span class="sc">:g},\t EI({</span>tt[<span class="op">-</span><span class="dv">1</span>]<span class="sc">})={</span>uu_EI[<span class="op">-</span><span class="dv">1</span>]<span class="sc">:g},\t RK4({</span>tt[<span class="op">-</span><span class="dv">1</span>]<span class="sc">})={</span>uu_RK4[<span class="op">-</span><span class="dv">1</span>]<span class="sc">:g}")

Interval [0,0.8]
    Avec h=0.1 on a      EE(0.8)=3.50783,    EI(0.8)=4.99976,    RK4(0.8)=5.84201
    Avec h=0.05 on a     EE(0.8)=4.20131,    EI(0.8)=9.99988,    RK4(0.8)=5.84811
    Avec h=0.01 on a     EE(0.8)=5.34276,    EI(0.8)=6.56639,    RK4(0.8)=5.84862
Interval [0,0.9]
    Avec h=0.1 on a      EE(0.9)=4.80232,    EI(0.9)=4.99976,    RK4(0.9)=14.0218
    Avec h=0.05 on a     EE(0.9)=6.4606,     EI(0.9)=9.99988,    RK4(0.9)=14.2712
    Avec h=0.01 on a     EE(0.9)=10.8041,    EI(0.9)=28.1192,    RK4(0.9)=14.3048
Interval [0,1]
    Avec h=0.1 on a      EE(1.0)=7.18955,    EI(1.0)=4.99976,    RK4(1.0)=735.099
    Avec h=0.05 on a     EE(1.0)=12.3209,    EI(1.0)=9.99988,    RK4(1.0)=175863
    Avec h=0.01 on a     EE(1.0)=90.7555,    EI(1.0)=50.0038,    RK4(1.0)=inf

/tmp/ipykernel_56231/1239910182.py:10: RuntimeWarning: overflow encountered in double_scalars
  phi = lambda t,y : y**2+t**2

Le problème présuppose qu’une solution existe sur l’intervalle \([0;1]\) alors qu’en fait aucune solution n’existe. La théorie générale dit qu’une solution unique existe pour un certain intervalle contenant \(t=0\), mais elle ne nous dit pas jusqu’où s’étend cet intervalle. Nous pouvons montrer qu’une solution existe pour \(t>\pi/4\). Cependant, la solution devient non bornée quelque part entre \(\pi/4\) et \(1\). En effet, la solution cherchée \(y\) est bornée par les deux fonctions suivantes: \( z(t)\le y(t) \le w(t), \qquad t\in[0;1] \) où \( \begin{cases} z'(t)=z^2(t), & t\in[0;1],\\ z(0)=1; \end{cases} \qquad \begin{cases} w'(t)=w^2(t)+1, & t\in[0;1],\\ w(0)=1. \end{cases} \) On trouve \(z(t)=\frac{1}{1-t}\) qui est définie sur \([0;1[\) et \(w(t)=\tan\left(t+\frac{\pi}{4}\right)\) qui est définie sur \([0;\pi/4[\).

Code

figure(figsize=(10,7))

# quiverplot
# define a grid and compute direction at each point
g1  = linspace(0,1,21)
g2  = linspace(0,3,21)
T,Y = meshgrid(g1,g2)         # create a grid
DT, DY = 1, phi(T,Y)   # compute growth rate on the grid
M = sqrt(DT**2+DY**2)         # norm growth rate 
M[ M==0 ] = 1                 # avoid zero division errors 
quiver(T,Y, DT/M, DY/M, M, pivot='mid')
# streamplot(T,Y, DT/M, DY/M, color=M)

grid()
xlabel('t')
ylabel('y')
title('Champ des pentes');

Cet exemple montre l’importance de l’interaction entre la théorie et le calcul numérique. Si vous vous contentez de calculer numériquement des solutions sans savoir qu’une solution existe, vous pouvez penser que vous avez réussi alors que vous êtes en train de calculer quelque chose qui n’existe pas. Prouvez l’existence et l’unicité avant de calculer. La théorie vient en premier. Par exemple, ce serait bien de connaître la taille de l’intervalle d’existence de la solution avant de calculer quoi que ce soit, mais ce n’est pas toujours possible. De même, il n’est pas évident, en regardant l’équation, de voir qu’il y a un problème pour \(t \to 1\). Les difficultés que nous avons rencontrées avec le calcul numérique suggèrent qu’il pourrait y avoir un problème théorique. Je ne dirais pas que la théorie ou le calcul viennent nécessairement en premier. Je dirais qu’il faut itérer entre les deux, en commençant par l’approche la plus facile à mettre en œuvre. Les résultats théoriques sont plus satisfaisants lorsqu’ils sont disponibles, mais la théorie ne nous en dit souvent pas autant que ce que nous aimerions savoir. En outre, les gens font des erreurs dans les calculs théoriques, tout comme dans les calculs numériques. Il est préférable que les travaux théoriques et numériques se valident mutuellement.

Ce problème montre l’importance de se préoccuper de l’existence et de l’unicité, mais les méthodes théoriques ne sont pas les seules méthodes pour explorer l’existence. Une bonne pratique numérique, c’est-à-dire l’essai de plus d’une taille de pas ou de plus d’une méthode numérique, est également précieuse. Quoi qu’il en soit, le problème montre que sans une certaine diligence - théorique ou numérique - vous pourriez naïvement calculer une “solution” approximative là où il n’existe aucune solution.

Source : https://www.johndcook.com/blog/2009/08/11/approximating-a-solution-that-doesnt-exist/
Source : (page 438 et 447) Boyce and DiPrima Elementary Differential Equations and Boundary Value Problems, seventh edition

2 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} \) 1. Calculer la solution exacte. 2. 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\). 3. Même exercice avec la méthode d’Euler explicite. 3. Commenter les résultats.

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

Code

%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())

EDO

Solution generale

En prenant en compte la CI on trouve

Solution particuliere

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = 3 y{\left(t \right)} - 4 e^{- t}\)

\(\displaystyle y{\left(t \right)} = C_{1} e^{3 t} + e^{- t}\)

\(\displaystyle \left\{ C_{1} : \epsilon\right\}\)

\(\displaystyle y{\left(t \right)} = \epsilon e^{3 t} + e^{- t}\)

Code

%matplotlib inline
from matplotlib.pylab import *
# rcdefaults()
rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})

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));

2.2 Correction : solution approchée (avec RK41 et EE)

Code

%reset -f
%matplotlib inline

from matplotlib.pylab import *
# rcdefaults()
rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})

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))

3 Exercice : A-stabilité du schéma d’Euler expicite dans le cas d’un système

Soit le problème de Cauchy: \(\begin{cases} y''(t)+1001y'(t)+1000y(t)=0, & \forall t \in \mathbb{R},\\ y(0)=1,\\ y'(0)=0. \end{cases}\) 1. Écrire le problème comme un système de 2 EDO d’ordre 1 et en calculer la solution. 1. Écrire le schéma le schéma d’Euler explicite pour ce problème. 4. Pour quelles valeurs de \(h\) le schéma est A-stable?

3.1 Correction

Soit \(\mathbf{z}(t)=(y(t),y'(t))^T\) alors \(y''(t)+1001y'(t)+1000y(t)=0\) se réécrit \( \mathbf{z}'(t)= -\begin{pmatrix} 0&-1\\ 1000&1001 \end{pmatrix} \mathbf{z}(t) \) Cette matrice a pour valeurs propres \(\lambda_1=1\) et \(\lambda_2=1000\) et pour vecteurs propres associés \(\mathbf{v}_1=(-1,1)^T\) et \(\mathbf{v}_2=(0,1000)^T\) (ci dessous le calcul formel des valeurs et vecteurs propes). Ainsi \( \mathbf{z}(t)=c_1\mathbf{v}_1e^{-\lambda_1t}+c_2\mathbf{v}_2e^{-\lambda_2t} \) et \( y(t)=-c_1e^{-t}-c_2e^{-1001t}. \) Il s’agit d’un problème stiff!

Code

%reset -f
import sympy as symb
symb.init_printing()

symb.var('t')
M = symb.Matrix(((0,-1), (1000,1001)))
display(M)

p=M.charpoly(t)
display(symb.factor(p))

P, D = M.diagonalize()
display(P,D)

display(symb.simplify(P*D*P**-1))

\(\displaystyle \left[\begin{matrix}0 & -1\\1000 & 1001\end{matrix}\right]\)

\(\displaystyle \operatorname{PurePoly}{\left( t^{2} - 1001 t + 1000, t, domain=\mathbb{Z} \right)}\)

\(\displaystyle \left[\begin{matrix}-1 & -1\\1 & 1000\end{matrix}\right]\)

\(\displaystyle \left[\begin{matrix}1 & 0\\0 & 1000\end{matrix}\right]\)

\(\displaystyle \left[\begin{matrix}0 & -1\\1000 & 1001\end{matrix}\right]\)

3.2 Correction

Soit \(h>0\) fixé et \(t_n=nh\) pour tout \(n\in\mathbb{Z}\). Le schéma d’Euler explicite pour le système donné construit la suite \( \begin{aligned} \mathbf{u}_{n+1} &= \mathbf{u}_n -h \begin{pmatrix} 0&-1\\ 1000&1001 \end{pmatrix} \mathbf{u}_n \\ &= \begin{pmatrix} 1 & h\\ -1000h & 1-1001h \end{pmatrix} \mathbf{u}_n \\ &= \begin{pmatrix} 1 & h\\ -1000h & 1-1001h \end{pmatrix}^{n+1} \mathbf{u}_0,\\ &= \begin{pmatrix} -1 & -1\\ 1000 & 1 \end{pmatrix} \begin{pmatrix} 1-1000h & 0\\ 0 & 1-h \end{pmatrix}^{n+1} \begin{pmatrix} -1 & -1\\ 1000 & 1 \end{pmatrix}^{-1} \mathbf{u}_0\\ &= \begin{pmatrix} -1 & -1\\ 1000 & 1 \end{pmatrix} \begin{pmatrix} (1-1000h)^{n+1} & 0\\ 0 & (1-h)^{n+1} \end{pmatrix} \begin{pmatrix} -1 & -1\\ 1000 & 1 \end{pmatrix}^{-1} \mathbf{u}_0, \qquad \forall n \in \mathbb{N}. \end{aligned} \)

Code

symb.var('h')
M = symb.Matrix(((1,h), (-1000*h,1-1001*h)))
display(M)

P, D = M.diagonalize()
display(P , D)

\(\displaystyle \left[\begin{matrix}1 & h\\- 1000 h & 1 - 1001 h\end{matrix}\right]\)

\(\displaystyle \left[\begin{matrix}-1 & -1\\1000 & 1\end{matrix}\right]\)

\(\displaystyle \left[\begin{matrix}1 - 1000 h & 0\\0 & 1 - h\end{matrix}\right]\)

3.3 Correction

Le schéma est A-stable si \(\lim_{n\to+\infty}\mathbf{u}_{n+1}=\mathbf{0}\). Il faut donc que \(\begin{cases} |1-1000h|<1\\ |1-h|<1 \end{cases}\) d’où la condition \(h<\frac{2}{1000}\).

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

4.1 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} \)

Code

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]))

u(4.50)=0.0003302503231810

5 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)\).
1. 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).
2. 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.
3. 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.

5.1 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)\):

Code

from matplotlib.pylab import *
# rcdefaults()
rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})

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

5.2 Correction étude numérique

Code

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)

5.2.1 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}\)

Code

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');

Euler explicite : |energie totale (t=Tfinal) - Energie totale (t=0)|= 3.1764943624320936

5.2.2 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:

Code

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');

Euler symplectique : |energie totale (t=Tfinal)-Energie totale (t=0)|= 0.030048843024607974

5.2.3 Méthodes d’Euler modifiée

Code

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');

Euler modifie : |energie totale (t=Tfinal)-Energie totale (t=0)|= 0.02457609735191002

5.2.4 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} \)

Code

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');

Heun : |energie totale (t=Tfinal)-Energie totale (t=0)|= 0.025432751256537878

5.2.5 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} \)

Code

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');

AB2 : |energie totale (t=Tfinal)-Energie totale (t=0)|= 0.07960924681985027

5.2.6 Méthode AB3

\( \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_2=x_1+\frac{h}{2}\left( 3\varphi_1(t_1,x_1,y_1) - \varphi_1(t_{0},x_{0},y_{0})\right),\\ y_2=y_1+\frac{h}{2}\left( 3\varphi_2(t_1,x_1,y_1) - \varphi_2(t_{0},x_{0},y_{0})\right),\\(0.5em] x_{n+1}=x_n+\frac{h}{12}\left( 23\varphi_1(t_n,x_n,y_n) - 16\varphi_1(t_{n-1},x_{n-1},y_{n-1}) + 5\varphi_1(t_{n-2},x_{n-2},y_{n-2})\right),\\ y_{n+1}=y_n+\frac{h}{12}\left( 23\varphi_2(t_n,x_n,y_n) - 16\varphi_2(t_{n-1},x_{n-1},y_{n-1}) + 5\varphi_2(t_{n-2},x_{n-2},y_{n-2})\right). \end{cases} \)

Code

def AB3(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]))
    q.append(q[1]+h*(3*phi1(tt[1],q[1],p[1])-phi1(tt[0],q[0],p[0]))/2. )
    p.append(p[1]+h*(3*phi2(tt[1],q[1],p[1])-phi2(tt[0],q[0],p[0]))/2.)
    for i in range(2,len(tt)-1):
        qtilde1 = phi1( tt[i],q[i],p[i] )
        qtilde2 = phi1( tt[i-1],q[i-1],p[i-1] )
        qtilde3 = phi1( tt[i-2],q[i-2],p[i-2] )
        ptilde1 = phi2( tt[i],q[i],p[i] )
        ptilde2 = phi2( tt[i-1],q[i-1],p[i-1] )
        ptilde3 = phi2( tt[i-2],q[i-2],p[i-2] )
        q.append( q[i] + (23*qtilde1-16*qtilde2+5*qtilde3)*h/12 )
        p.append( p[i] + (23*ptilde1-16*ptilde2+5*ptilde3)*h/12 )
    return [q,p]
[q_AB3, p_AB3] = AB3(phi1,phi2,tt,y0)
Ec_AB3 = [p**2/2. for p in p_AB3]
Ep_AB3 = [-math.cos(q) for q in q_AB3]
Et_AB3 = [Ec_AB3[i]+Ep_AB3[i] for i in range(len(Ec_AB3))]
print( 'AB3 : |energie totale (t=Tfinal)-Energie totale (t=0)|=', abs(Et_AB3[-1]-Et_AB3[0]))

figure(figsize=(18,5))
subplot(131)
plot(tt,q_AB3,tt,p_AB3)
xlabel('t')
legend(['x(t)','y(t)'])
title('AB3 - x(t) et y(t)')
subplot(132)
plot(q_AB3,p_AB3)
xlabel('x(t)')
ylabel('y(t)')
title('AB3 - y(x)')
subplot(133)
plot(tt,Ec_AB3,tt,Ep_AB3,tt,Et_AB3)
xlabel('t')
legend(['Cinetique','Potentielle','Totale'])
title('AB3 - Energies');

AB3 : |energie totale (t=Tfinal)-Energie totale (t=0)|= 0.051119920140796915

5.2.7 Méthode RK4

\( \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] \check x_n=x_n+\frac{h}{2}\varphi_1(t_n+h/2,\tilde x_n,\tilde y_n),\\ \check y_n=y_n+\frac{h}{2}\varphi_2(t_n+h/2,\tilde x_n,\tilde y_n),\\(0.5em] \hat x_n=x_n+h\varphi_1(t_{n+1},\check x_n,\check y_n),\\ \hat y_n=y_n+h\varphi_2(t_{n+1},\check x_n,\check y_n),\\(0.5em] x_{n+1}=x_n+\frac{h}{6}\left( \varphi_1(t_n,x_n,y_n) + 2\varphi_1(t_{n}+h/2,\tilde x_n,\tilde y_n) + 2\varphi_1(t_{n}+h/2,\check x_n,\check y_n) + \varphi_1(t_{n+1},\hat x_n,\hat y_n)\right),\\ y_{n+1}=y_n+\frac{h}{6}\left( \varphi_2(t_n,x_n,y_n) + 2\varphi_2(t_{n}+h/2,\tilde x_n,\tilde y_n) + 2\varphi_2(t_{n}+h/2,\check x_n,\check y_n) + \varphi_2(t_{n+1},\hat x_n,\hat y_n)\right). \end{cases} \)

Code

def RK4(phi1,phi2,tt,y0):
    h = tt[1]-tt[0]
    q = [y0[0]]
    p = [y0[1]]
    for i in range(len(tt)-1):
        qtilde = q[i] + 0.5*h * phi1( tt[i], q[i], p[i] )
        ptilde = p[i] + 0.5*h * phi2( tt[i], q[i], p[i] )
        qv = q[i] + 0.5*h * phi1( tt[i]+h*0.5, qtilde, ptilde )
        pv = p[i] + 0.5*h * phi2( tt[i]+h*0.5, qtilde, ptilde )
        qhat = q[i] + h * phi1( tt[i+1], qv, pv )
        phat = p[i] + h * phi2( tt[i+1], qv, pv )
        q.append( q[i]+h/6*( phi1(tt[i],q[i],p[i]) + 2*phi1(tt[i]+0.5*h,qtilde,ptilde)+ 2*phi1(tt[i]+0.5*h,qv,pv)+ phi1(tt[i+1],qhat,phat) ) )
        p.append( p[i]+h/6*( phi2(tt[i],q[i],p[i]) + 2*phi2(tt[i]+0.5*h,qtilde,ptilde)+ 2*phi2(tt[i]+0.5*h,qv,pv)+ phi2(tt[i+1],qhat,phat) ) )
    return [q,p]
[q_RK4, p_RK4] = RK4(phi1,phi2,tt,y0)
Ec_RK4 = [p**2/2. for p in p_RK4]
Ep_RK4 = [-math.cos(q) for q in q_RK4]
Et_RK4 = [Ec_RK4[i]+Ep_RK4[i] for i in range(len(Ec_RK4))]
print ('RK4 : |energie totale (t=Tfinal)-Energie totale (t=0)|=', abs(Et_RK4[-1]-Et_RK4[0]))

figure(figsize=(18,5))
subplot(131)
plot(tt,q_RK4,tt,p_RK4)
xlabel('t')
legend(['x(t)','y(t)'])
title('RK4 - x(t) et y(t)')
subplot(132)
plot(q_RK4,p_RK4)
xlabel('x(t)')
ylabel('y(t)')
title('RK4 - y(x)')
subplot(133)
plot(tt,Ec_RK4,tt,Ep_RK4,tt,Et_RK4)
xlabel('t')
legend(['Cinetique','Potentielle','Totale'])
title('RK4 - Energies');

RK4 : |energie totale (t=Tfinal)-Energie totale (t=0)|= 8.32124615898211e-05

Retour au sommet

1 Exercice : Étude d’un problème mathématiquement mal posé

1.1 Correction

2 Exercice : Étude d’un problème numériquement mal posé

2.1 Correction : solution exacte (avec SymPy)

2.2 Correction : solution approchée (avec RK41 et EE)

3 Exercice : A-stabilité du schéma d’Euler expicite dans le cas d’un système

3.1 Correction

3.2 Correction

3.3 Correction

4 Exercice : EDO d’ordre 2

4.1 Correction

5 Exercice : Étude théorique et numérique du mouvement d’un pendule

5.1 Correction étude théorique

5.2 Correction étude numérique

5.2.1 Méthode d’Euler explicite

5.2.2 Méthode d’Euler explicite/implicite

5.2.3 Méthodes d’Euler modifiée

5.2.4 Méthodes de Heun

5.2.5 Méthode AB2

5.2.6 Méthode AB3

5.2.7 Méthode RK4

2.1 Correction : solution exacte (avec `SymPy`)