TD4 - Révisions – M64 – Data Science & Scientific Computing 6 : Schémas numériques pour les problèmes de Cauchy

TD4 - Révisions

Auteur·rice

Gloria Faccanoni

Date de publication

14 février 2026

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 les méthodes d’Euler implicite, de Crank-Nicolson et 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).

  • 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

Correction

Code 
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

# Configuration des paramètres d'affichage
plt.rcParams.update({'font.size': 14})

# Définition de l'équation différentielle
phi = lambda t,y : y**2 + t**2

# Schéma d'Euler Explicite (EE)
def euler_explicite(phi, tt, y0):
    h = tt[1] - tt[0]
    uu = np.empty_like(tt)  
    uu[0] = y0  
    for i in range(len(tt) - 1):
        uu[i + 1] = uu[i] + h * phi(tt[i], uu[i])  
    return uu

# Schéma d'Euler Implicite (EI)
def euler_implicite(phi, tt, y0):
    h = tt[1] - tt[0]
    uu = np.empty_like(tt)  
    uu[0] = y0  
    for i in range(len(tt) - 1):
        uu[i + 1] = fsolve(lambda x: -x + uu[i] + h * phi(tt[i + 1], x), uu[i])[0]
    return uu

# Schéma de Crank-Nicolson (CN)
def crank_nicolson(phi, tt, y0):
    h = tt[1] - tt[0]
    uu = np.empty_like(tt)  
    uu[0] = y0  
    for i in range(len(tt) - 1):
        f = lambda x: -x + uu[i] + h * (phi(tt[i], uu[i]) + phi(tt[i + 1], x)) / 2
        uu[i + 1] = fsolve(f, uu[i])[0]
    return uu

# Schéma de Runge-Kutta d'ordre 4 (RK4)
def runge_kutta_4(phi, tt, y0):
    h = tt[1] - tt[0]
    uu = np.empty_like(tt)  
    uu[0] = 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[i + 1] = uu[i] + (k1 + 2 * k2 + 2 * k3 + k4) * h / 6
    return uu

# Initialisation des paramètres
t0, tfinal, y0 = 0, 1, 1

# Création d'un tableau pour stocker tous les résultats
results = []


# Création de la figure avec 3 lignes et 3 colonnes
fig, axes = plt.subplots(3, 3, figsize=(18, 18))  


for row, epsilon in enumerate([0.2, 0.1, 0]):
    for col, h in enumerate([0.1, 0.05, 0.01]):
        N = int((1 - epsilon) / h)
        tt = np.linspace(0, 1 - epsilon, N + 1)  # Utilisation de linspace pour éviter les erreurs d'arrondi
        uu_EE = euler_explicite(phi, tt, y0)
        uu_EI = euler_implicite(phi, tt, y0)
        uu_CN = crank_nicolson(phi, tt, y0)
        uu_RK4 = runge_kutta_4(phi, tt, y0)
        
        # Enregistrer les résultats dans le tableau
        results.append([epsilon, h, uu_EE[-1], uu_EI[-1], uu_CN[-1], uu_RK4[-1]])

        # Tracer les solutions dans le sous-graphe correspondant
        if row == 0:
            axes[row, col].set_title(f"$h={h}$", fontsize=14)
        axes[row, col].plot(tt, uu_EE, '*-', label=f'EE', color='r')
        axes[row, col].plot(tt, uu_EI, 'o-', label=f'EI', color='g')
        axes[row, col].plot(tt, uu_CN, 's-', label=f'CN', color='m')
        axes[row, col].plot(tt, uu_RK4, 'D-', label=f'RK4', color='b')
        axes[row, col].set_xlabel('$t$', fontsize=12)
        axes[row, col].set_ylabel('$y(t)$', fontsize=12)
        axes[row, col].set_xlim([0, 1])
        axes[row, col].grid(True)
        axes[row, col].legend(loc='best')

    axes[row, 0].set_ylabel(f'tfinal=$={1 - epsilon}$', fontsize=16)  # Ajouter l'étiquette de epsilon

plt.tight_layout()  # Pour un espacement optimal entre les sous-figures
plt.show()

# Afficher tous les résultats dans un tableau
from tabulate import tabulate

headers = ["tfinal", "h", "EE(tfinal)", "EI(tfinal)", "CN(tfinal)", "RK4(tfinal)"]
table = []

# Remplir le tableau avec les résultats, en ajoutant un séparateur après chaque changement de tfinal
for result in results:
    table.append([1 - result[0], result[1], result[2], result[3], result[4], result[5]])

table_str = tabulate(table, headers=headers, tablefmt="html", floatfmt=".2f")
display(table_str)

/usr/lib/python3/dist-packages/scipy/optimize/_minpack_py.py:177: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last five Jacobian evaluations.
  warnings.warn(msg, RuntimeWarning)
/usr/lib/python3/dist-packages/scipy/optimize/_minpack_py.py:177: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last ten iterations.
  warnings.warn(msg, RuntimeWarning)
/tmp/ipykernel_254831/4035068944.py:9: RuntimeWarning: overflow encountered in scalar power
  phi = lambda t,y : y**2 + t**2

tfinal h EE(tfinal) EI(tfinal) CN(tfinal) RK4(tfinal)
0.80 0.10 3.51 5.00 7.45 5.84
0.80 0.05 4.20 10.00 6.08 5.85
0.80 0.01 5.34 6.57 5.86 5.85
0.90 0.10 4.80 5.00 10.00 14.02
0.90 0.05 6.46 10.00 20.01 14.27
0.90 0.01 10.80 28.12 14.44 14.30
1.00 0.10 7.19 5.00 10.00 735.10
1.00 0.05 12.32 10.00 20.00 175862.67
1.00 0.01 90.76 50.00 100.00 inf

Le problème est que nous avons présupposé que le problème est mathématiquement bien posé, c’est-à-dire qu’une unique solution existe sur l’intervalle d’étude \([0;1]\). Il se trouve qu’aucune solution n’existe sur l’intervalle \([0;1]\).

La théorie générale permet de montrer 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 au moins \(t>\pi/4\). Cependant, la solution devient non bornée quelque part entre \(\pi/4\) et \(1\). Nous pouvons le prouver en observant que la solution cherchée \(y\) est bornée par les deux fonctions \(z\) et \(w\) suivantes

\( z(t)\le y(t) \le w(t), \qquad t\in[0;1] \)

vérifiant les problèmes de Cauchy

\( \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[\),
  • \(w(t)=\tan\left(t+\frac{\pi}{4}\right)\) qui est définie sur \([0;\pi/4[\).

Puisque \(z(t)\to+\infty\) quand \(t\to 1^-\), la solution cherchée \(y\) devient non bornée quelque part entre \(\pi/4\) et \(1\).

Même sans connaître la solution exacte, nous pouvons étudier le comportement sur l’intervalle \([0;1]\) en affichant le champ de vecteurs :

Code 
plt.figure(figsize=(10,7))

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

plt.grid()
plt.xlabel('t')
plt.ylabel('y')
plt.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 une solution approchée avec un schéma et un seul pas \(h\) sans savoir qu’une solution unique existe, vous pouvez penser que vous avez réussi alors que vous êtes en train de calculer quelque chose qui n’existe pas.
  • Lorsque cela est possible, il faut essayer d’abord de prouver que le problème est mathématiquement bien posé avant d’en calculer une solution approchée. 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. Bien-sûr, ce n’est pas toujours possible. Dans notre exemple il n’est pas évident, en regardant l’équation, de voir qu’il y a un problème pour \(t \to 1\). L’affichage du champ de vecteurs peut aider à voir ce genre de problème.
  • Une bonne pratique numérique, c’est-à-dire l’essai de plus d’une taille de pas et/ou de plus d’une méthode numérique, est également précieuse. Les difficultés que nous avons rencontrées avec les comparaisons suggèrent qu’il pourrait y avoir un problème théorique sous-jacent qu’il faut clarifier.

On ne peut pas dire que la théorie ou le calcul viennent nécessairement en premier, 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. 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.

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 cinq valeurs des \(\varepsilon\), à savoir \(\varepsilon=\pm 1\), \(\varepsilon=\pm \frac{1}{10}\) et \(\varepsilon=0\).
  3. Même exercice avec la méthode d’Euler explicite.
  4. Commenter les résultats.

Correction : solution exacte

La solution exacte est \(y(t)=\varepsilon e^{3t}+e^{-t}\) (ci-dessous, le calcul de la solution exacte avec sympy).

Code 
import sympy as sym
sym.init_printing(pretty_print=True)

t = sym.Symbol('t', real=True, positive=True)
epsilon = sym.Symbol(r'\varepsilon', real=True, positive=True)
y = sym.Function('y')

edo = sym.Eq( sym.diff(y(t),t), 3*y(t)-4*sym.exp(-t) )
display(edo)

t0 = 0
y0 = 1+epsilon
display( sym.Eq(y(t).subs(t,t0), y0) )

sol = sym.dsolve(edo,y(t), ics={y(t0):y0})
display(sol.expand())

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

\(\displaystyle y{\left(0 \right)} = \varepsilon + 1\)

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

Si \(\varepsilon=0\), la solution exacte devient donc \(y(t)=e^{-t}\) mais le problème est mal conditionné : toute (petite) erreur de calcul a le même effet qu’une perturbation de la condition initiale. En d’autres termes, on “réveil” le terme dormant \(e^{3t}\).

Affichons l’allure de la solution exacte dans les cinq cas \(\varepsilon=\pm1\), \(\varepsilon=\pm\frac{1}{10}\) et \(\varepsilon=0\). On remarque que la solution exacte est très sensible à la condition initiale.

Code 
import numpy as np
import matplotlib.pyplot as plt

# rcdefaults()
plt.rcParams.update({'font.size': 14})

func = sym.lambdify((t,epsilon),sol.rhs,'numpy')

tt = np.linspace(0,1,101)
plt.figure(figsize=(9,5))
EPSILON=[1,1.e-1,0,-1.e-1,-1]
for i,epsi in enumerate(EPSILON):
    plt.plot( tt, func(tt,epsi),
         #label=r"$\epsilon=$"+str(epsi), 
         label=r'$y_0=$'+str(1+epsi),
         linestyle=['-','--',':','-.'][i//2])

plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') 
plt.grid(); plt.xlabel('$t$'); plt.ylabel('$y(t)$');

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

Code 
import numpy as np
import matplotlib.pyplot as plt
# plt.rcdefaults()
plt.rcParams.update({'font.size': 16})

phi = lambda t,y : 3*y-4*np.exp(-t)
sol_exacte = lambda t,y0 : (y0-1)*np.exp(3*t)+np.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,-1.e-1,-1]

plt.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)
    plt.subplot(1,5,i+1)
    plt.plot(tt,yy,'b-',tt,uu_RK4,'y:D',tt,uu_EE,'m:o')
    plt.legend(['Exacte','RK4','EE']);
    plt.title(r"$\epsilon=$"+str(epsilon))

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

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

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.
  2. Écrire le schéma le schéma d’Euler explicite pour ce problème.
  3. Pour quelles valeurs de \(h\) le schéma est A-stable?

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

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

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

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

  1. É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)\in[-\pi;\pi]\times\mathbb{R}\ |\ E(q,p)=\kappa\} =\{(q,p)\in[-\pi;\pi]\times\mathbb{R}\ |\ p^2=2(\kappa+\cos(q))\}. \)
    • Lorsque \(\vartheta\) est petit on peut considérer l’approximation \(\sin(\vartheta)\simeq\vartheta\). Calculer la solution exacte de cette équation approchée.

  2. É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.

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{aligned} \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{aligned}\)

  • Courbes de niveau de la fonction \((q,p)\mapsto E(q,p)\):

Code 
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams.update({'font.size': 16})

Etot = lambda q,p : -np.cos(q)+p**2/2

q,p = np.meshgrid(np.linspace(-8*np.pi,8*np.pi,201),np.linspace(-3,3,201))
E = Etot(q,p)

plt.figure(1,figsize=(20,5))
graphe = plt.contour(q,p,E,20)
plt.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

Code 
t0 = 0
tfinal = 8*np.pi
y0 = [np.pi/4,0]

N = 100 
tt = np.linspace(t0,tfinal,N+1)
#h = tt[1]-tt[0]
    
phi1 = lambda t,q,p : p
phi2 = lambda t,q,p : -np.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}\)

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 = [-np.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]))

plt.figure(figsize=(18,5))
plt.subplot(131)
plt.plot(tt,q_ep,tt,p_ep)
plt.xlabel('t')
plt.legend(['x(t)','y(t)'])
plt.title('Euler progressif - x(t) et y(t)')
plt.subplot(132)
plt.plot(q_ep,p_ep)
plt.xlabel('x(t)')
plt.ylabel('y(t)')
plt.title('Euler progressif')
plt.subplot(133)
plt.plot(tt,Ec_ep,tt,Ep_ep,tt,Et_ep)
plt.xlabel('t')
plt.legend(['Cinetique','Potentielle','Totale'])
plt.title('Euler progressif - Energies');
Euler explicite : |energie totale (t=Tfinal) - Energie totale (t=0)|= 3.1764943624320936

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 = [-np.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]))

plt.figure(figsize=(18,5))
plt.subplot(131)
plt.plot(tt,q_es,tt,p_es)
plt.xlabel('t')
plt.legend(['x(t)','y(t)'])
plt.title('Euler symplectique - x(t) et y(t)')
plt.subplot(132)
plt.plot(q_es,p_es)
plt.xlabel('x(t)')
plt.ylabel('y(t)')
plt.title('Euler symplectique - y(x)')
plt.subplot(133)
plt.plot(tt,Ec_es,tt,Ep_es,tt,Et_es)
plt.xlabel('t')
plt.legend(['Cinetique','Potentielle','Totale'])
plt.title('Euler symplectique - Energies');
Euler symplectique : |energie totale (t=Tfinal)-Energie totale (t=0)|= 0.030048843024607974

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

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 = [-np.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]))

plt.figure(figsize=(18,5))
plt.subplot(131)
plt.plot(tt,q_em,tt,p_em)
plt.xlabel('t')
plt.legend(['x(t)','y(t)'])
plt.title('Euler modifie - x(t) et y(t)')
plt.subplot(132)
plt.plot(q_em,p_em)
plt.xlabel('x(t)')
plt.ylabel('y(t)')
plt.title('Euler modifie - y(x)')
plt.subplot(133)
plt.plot(tt,Ec_em,tt,Ep_em,tt,Et_em)
plt.xlabel('t')
plt.legend(['Cinetique','Potentielle','Totale'])
plt.title('Euler modifie - Energies');
Euler modifie : |energie totale (t=Tfinal)-Energie totale (t=0)|= 0.02457609735191002

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 = [-np.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]))

plt.figure(figsize=(18,5))
plt.subplot(131)
plt.plot(tt,q_heun,tt,p_heun)
plt.xlabel('t')
plt.legend(['x(t)','y(t)'])
plt.title('Heun - x(t) et y(t)')
plt.subplot(132)
plt.plot(q_heun,p_heun)
plt.xlabel('x(t)')
plt.ylabel('y(t)')
plt.title('Heun - y(x)')
plt.subplot(133)
plt.plot(tt,Ec_heun,tt,Ep_heun,tt,Et_heun)
plt.xlabel('t')
plt.legend(['Cinetique','Potentielle','Totale'])
plt.title('Heun - Energies');
Heun : |energie totale (t=Tfinal)-Energie totale (t=0)|= 0.025432751256537878

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 = [-np.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]))

plt.figure(figsize=(18,5))
plt.subplot(131)
plt.plot(tt,q_AB2,tt,p_AB2)
plt.xlabel('t')
plt.legend(['x(t)','y(t)'])
plt.title('AB2 - x(t) et y(t)')
plt.subplot(132)
plt.plot(q_AB2,p_AB2)
plt.xlabel('x(t)')
plt.ylabel('y(t)')
plt.title('AB2 - y(x)')
plt.subplot(133)
plt.plot(tt,Ec_AB2,tt,Ep_AB2,tt,Et_AB2)
plt.xlabel('t')
plt.legend(['Cinetique','Potentielle','Totale'])
plt.title('AB2 - Energies');
AB2 : |energie totale (t=Tfinal)-Energie totale (t=0)|= 0.07960924681985027

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 = [-np.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]))

plt.figure(figsize=(18,5))
plt.subplot(131)
plt.plot(tt,q_AB3,tt,p_AB3)
plt.xlabel('t')
plt.legend(['x(t)','y(t)'])
plt.title('AB3 - x(t) et y(t)')
plt.subplot(132)
plt.plot(q_AB3,p_AB3)
plt.xlabel('x(t)')
plt.ylabel('y(t)')
plt.title('AB3 - y(x)')
plt.subplot(133)
plt.plot(tt,Ec_AB3,tt,Ep_AB3,tt,Et_AB3)
plt.xlabel('t')
plt.legend(['Cinetique','Potentielle','Totale'])
plt.title('AB3 - Energies');
AB3 : |energie totale (t=Tfinal)-Energie totale (t=0)|= 0.051119920140796915

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 = [-np.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]))

plt.figure(figsize=(18,5))
plt.subplot(131)
plt.plot(tt,q_RK4,tt,p_RK4)
plt.xlabel('t')
plt.legend(['x(t)','y(t)'])
plt.title('RK4 - x(t) et y(t)')
plt.subplot(132)
plt.plot(q_RK4,p_RK4)
plt.xlabel('x(t)')
plt.ylabel('y(t)')
plt.title('RK4 - y(x)')
plt.subplot(133)
plt.plot(tt,Ec_RK4,tt,Ep_RK4,tt,Et_RK4)
plt.xlabel('t')
plt.legend(['Cinetique','Potentielle','Totale'])
plt.title('RK4 - Energies');
RK4 : |energie totale (t=Tfinal)-Energie totale (t=0)|= 8.32124615898211e-05