TP 2 - Implémentation des schémas “classiques” et étude de la convergence – M64 – Data Science & Scientific Computing 6 : Schémas numériques pour les problèmes de Cauchy

TP 2 - Implémentation des schémas “classiques” et étude de la convergence

Auteur·rice

Gloria Faccanoni

Date de publication

11 février 2026

NotePlan de la séance de TP

Cette séance de TP comporte 4 exercices :

  • Exercice 1 : Implémentation et test des schémas d’Euler explicite et implicite sur un exemple avec solution exacte. Votre travail consiste à ajouter les autres schémas classiques.

  • Exercice 2 : Étude de la convergence empirique des schémas d’Euler explicite et implicite sur le même exemple. Votre travail consiste à ajouter l’étude de la convergence des autres schémas implémentés précédemment.

  • Exercice 3 : Application de la méthode d’Euler explicite pour approcher une intégrale elliptique (aucune solution exacte n’est disponible dans ce cas).

  • Exercice 4 : Implémentation des schémas classiques appliqués à un système de deux EDO d’ordre 1. Bien que la solution exacte soit disponible pour ce problème, l’accent sera mis sur la conservation d’un invariant, souvent lié à un principe de conservation physique tel que l’énergie totale du système.

Calculons la solution exacte en utilisant le module sympy:

Code 
import sympy as sym
sym.init_printing()

# Variables
t = sym.Symbol('t')
y = sym.Function('y')

# Equation différentielle
edo = sym.Eq( sym.diff(y(t),t) , -4*y(t)+t**2 )
display(edo)

# Conditions initiales
t0 = 0
y0 = 1
display(sym.Eq(y(t0),y0))

# Méthode directe
solpar = sym.dsolve(edo, y(t), ics={y(t0):y0}) 
display(solpar)

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

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

\(\displaystyle y{\left(t \right)} = \frac{t^{2}}{4} - \frac{t}{8} + \frac{1}{32} + \frac{31 e^{- 4 t}}{32}\)

On définit la solution exacte à utiliser pour estimer les erreurs:

Code 
sol_exacte = sym.lambdify(t,solpar.rhs,'numpy')

On commence par importer les modules numpy et matplotlib et la fonction fsolve du module scipy.optimize pour résoudre l’équation implicite présente dans le schéma implicite :

Code 
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams.update({'font.size': 12}) # plt.rcParams.update(rcParamsDefault)
from scipy.optimize import fsolve # pour les schémas implicites

On initialise le problème de Cauchy

Code 
# variables globales
t0     = 0
tfinal = 1
y0     = 1

On définit l’équation différentielle : phi est une lambda function qui contient la fonction mathématique \(\varphi(t, y)=-4y+t^2\) dépendant des variables \(t\) et \(y\).

Code 
phi = lambda t,y : -4*y+t**2

On écrit les schémas numériques : les valeurs \([u_0,u_1,\dots,u_{N}]\) pour chaque méthode sont contenues dans le vecteur uu.

Schéma d’Euler explicite :

\( \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi(t_n,u_n)& n=0,1,2,\dots N-1 \end{cases} \)

Code 
def EE(phi, tt, y0):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt) - 1):
        uu[n+1] = uu[n] + h * phi(tt[n], uu[n])
    return uu
Code 
# !!!!!! Ajouter ici le schéma d'Euler modifié : en markdown ET en python

Schéma d’Euler modifié:

\( \begin{cases} u_0=y_0,\\ \tilde u = u_n+\frac{h}{2}\varphi(t_n,u_n),\\ u_{n+1}=u_n+h\varphi\left(t_n+\frac{h}{2},\tilde u\right)& n=0,1,2,\dots N-1 \end{cases} \)

qu’on peut réécrire sous la forme

\( \begin{cases} u_0=y_0,\\ k_1 = \varphi(t_n,u_n),\\ u_{n+1}=u_n+h\varphi\left(t_n+\frac{h}{2}, u_n+h\frac{k_1}{2}\right)& n=0,1,2,\dots N-1 \end{cases} \)

Code 
def EM(phi,tt, y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        uu[n+1] = uu[n]+h*phi( tt[n]+h/2 , uu[n]+h*k1/2 ) 
    return uu
Code 
# !!!!!! Ajouter ici le schéma de Heun : en markdown ET en python

Schéma de Heun:

\( \begin{cases} u_0=y_0,\\ \tilde u = u_n+h\varphi(t_n,u_n)\\ u_{n+1}=u_n+\frac{h}{2}\Bigl(\varphi(t_n,u_n)+\varphi(t_{n+1},\tilde u)\Bigr)& n=0,1,2,\dots N-1 \end{cases} \)

qu’on peut réécrire sous la forme

\( \begin{cases} u_0=y_0,\\ k_1 = \varphi(t_n,u_n),\\ k_2 = \varphi(t_{n+1},u_n+hk_1),\\ u_{n+1}=u_n+h\frac{k_1+k_2}{2}& n=0,1,2,\dots N-1 \end{cases} \)

Code 
def heun(phi,tt, y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        k2 = phi( tt[n+1], uu[n] + h*k1  )
        uu[n+1] = uu[n] + h*(k1+k2)/2 
    return uu

Schéma d’Euler implicite :

\( \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi(t_{n+1},u_{n+1})& n=0,1,2,\dots N-1 \end{cases} \)

Attention : \(u_{n+1}\) est solution de l’équation \(x=u_n+h\varphi(t_{n+1},x)\), c’est-à-dire un zéro de la fonction (en générale non linéaire)

\( x\mapsto -x+u_n+h\varphi(t_{n+1},x) \)

Code 
def EI(phi, tt, y0):
    h = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt) - 1):
        temp = fsolve(lambda x: -x + uu[n] + h * phi(tt[n + 1], x), uu[n])
        uu[n+1] = temp[0]
    return uu
Code 
# !!!!!! Ajouter ici le schéma de  Crank-Nicolson : en markdown ET en python

Schéma de Crank-Nicolson :

\( \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+\frac{h}{2}\Bigl(\varphi(t_n,u_n)+\varphi(t_{n+1},u_{n+1})\Bigr)& n=0,1,2,\dots N-1 \end{cases} \)

avec \(u_{n+1}\) zéro de la fonction \(x\mapsto -x+u_n+\frac{h}{2}(\varphi(t_n,u_n)+\varphi(t_{n+1},x))\).

Code 
def CN(phi,tt, y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt)-1):
        temp = fsolve(lambda x: -x+uu[n]+0.5*h*( phi(tt[n+1],x)+phi(tt[n],uu[n]) ), uu[n])
        uu[n+1] = temp[0]
    return uu
Code 
# !!!!!! Ajouter ici le schéma RK1_M : en markdown ET en python

Schéma RK1_M:

\( \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi\left(\dfrac{t_n+t_{n+1}}{2},\dfrac{u_n+u_{n+1}}{2}\right)& n=0,1,2,\dots N-1 \end{cases} \)

avec \(u_{n+1}\) zéro de la fonction \(x\mapsto -x+u_n+h\varphi\left(\dfrac{t_n+t_{n+1}}{2},\dfrac{u_n+x}{2}\right)\).

Code 
def RK1_M(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt)-1):
        temp = fsolve(lambda x: -x+uu[n]+h*phi( (tt[n]+tt[n+1])/2,(uu[n]+x)/2 ), uu[n])
        uu[n+1] = temp[0]
    return uu

On introduit la discrétisation : les nœuds d’intégration \([t_0,t_1,\dots,t_{N}]\) sont contenus dans le vecteur tt qui contient \(N+1\) points espacés de \(h=\frac{t_N-t_0}{N}\).

On construit la solution exacte discrète et les solutions approchées en ces points :

Code 
N = 8 
tt = np.linspace(t0,tfinal,N+1)

# yy      = np.array([sol_exacte(t) for t in tt], dtype=object)
yy      = sol_exacte(tt)  # si vectorisée
uu_ee   = EE(phi, tt, y0)
uu_ei   = EI(phi, tt, y0)
uu_em   = EM(phi, tt, y0)
uu_rk1m = RK1_M(phi, tt, y0)
uu_cn   = CN(phi, tt, y0)
uu_heun = heun(phi, tt, y0)

On compare les graphes de la solution exacte discrète et des solutions approchées et on affiche le maximum de l’erreur :

\( e = \max_{n\in[0,N]} \left| y(t_n)-u_n \right| = ||\mathbf{y}-\mathbf{u}||_{\infty} \)

Code 

plt.figure(figsize=(24, 14))

rows, cols, pos = 2, 3, 0

pos +=1 
plt.subplot(rows, cols, pos)
plt.plot(tt, yy, 'b-')
plt.plot(tt, uu_ee, 'r-D')
# err_ee = max( [abs(uu_ee[i] - yy[i]) for n in range(N)] )
err_ee = np.linalg.norm( uu_ee-yy, np.inf)
plt.title(f'Euler explicite\nmax(|erreur|)= {err_ee:1.10f}')
plt.grid()

pos +=1 
plt.subplot(rows, cols, pos)
plt.plot(tt, yy, 'b-', tt, uu_ei, 'r-D')
# err_ei = max( [abs(uu_ei[i] - yy[i]) for n in range(N)] )
err_ei = np.linalg.norm(uu_ei-yy, np.inf)
plt.title(f'Euler implicite\nmax(|erreur|)={err_ei:1.2f}')
plt.grid();


pos +=1 
plt.subplot(rows, cols, pos)
plt.plot(tt, yy, 'b-', tt, uu_em, 'r-D')
# err_em = max( [abs(uu_em[i] - yy[i]) for n in range(N)] )
err_em = np.linalg.norm(uu_em-yy, np.inf)
plt.title(f'Euler modifié\nmax(|erreur|)={err_em:1.5f}')
plt.grid()

pos +=1 
plt.subplot(rows, cols, pos)
plt.plot(tt, yy, 'b-', tt, uu_rk1m, 'r-D')
# err_rk1m = max( [abs(uu_rk1m[i] - yy[i]) for n in range(N)] )
err_rk1m = np.linalg.norm(uu_rk1m-yy, np.inf)
plt.title(f'RK1_m\nmax(|erreur|)={err_rk1m:1.5f}')
plt.grid()

pos +=1 
plt.subplot(rows, cols, pos)
plt.plot(tt, yy, 'b-', tt, uu_cn, 'r-D')
# err_cn = max( [abs(uu_cn[i] - yy[i]) for n in range(N)] )
err_cn = np.linalg.norm(uu_cn-yy,np.inf)
plt.title(f'Crank Nicolson\nmax(|erreur|)={err_cn:1.5f}')
plt.grid()

pos +=1 
plt.subplot(rows, cols, pos)
plt.plot(tt, yy, 'b-', tt, uu_heun, 'r-D')
# err_heun = max( [abs(uu_heun[i] - yy[i]) for n in range(N)] )
err_heun = np.linalg.norm(uu_heun-yy, np.inf)
plt.title(f'Heun\nmax(|erreur|)={err_heun:1.5f}')
plt.grid();

Pour chaque schéma, on calcule la solution approchée avec différentes valeurs de \(h_k=1/N_k\) et on sauvegarde les valeurs de \(h_k\) dans le vecteur H.

Pour chaque valeur de \(h_k\), on calcule le maximum de la valeur absolue de l’erreur et on sauvegarde toutes ces erreurs dans le vecteur err_schema de sort que err_schema[k] contient \(e_k=\max_{i=0,\dots,N_k}|y(t_i)-u_{i}|\).

Code 
H = []
err_ee = []
err_ei = []
# !!!!!! Ajouter ici les initialisations des nouveaux schémas
err_em = []
err_rk1m = []
err_cn = []
err_heun = []

for k in range(7):

    N = 2**(k + 3)
    
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)  # si vectorisée
    
    uu_ee   = EE(phi, tt, y0)
    uu_ei   = EI(phi, tt, y0)
    # !!!!!! Ajouter ici l'appel aux nouveaux schémas
    uu_em   = EM(phi, tt, y0)
    uu_rk1m = RK1_M(phi, tt, y0)
    uu_cn   = CN(phi, tt, y0)
    uu_heun = heun(phi, tt, y0)
    
    H.append(h)
    
    err_ee.append(np.linalg.norm(uu_ee-yy,np.inf))
    err_ei.append(np.linalg.norm(uu_ei-yy,np.inf))
    # !!!!!! Ajouter ici le calcul des erreurs
    err_em.append(np.linalg.norm(uu_em-yy,np.inf))
    err_rk1m.append(np.linalg.norm(uu_rk1m-yy,np.inf))
    err_cn.append(np.linalg.norm(uu_cn-yy,np.inf))
    err_heun.append(np.linalg.norm(uu_heun-yy,np.inf))

Pour afficher l’ordre de convergence on utilise une échelle logarithmique, i.e. on représente \(\ln(h)\) sur l’axe des abscisses et \(\ln(\text{err})\) sur l’axe des ordonnées.
En effet, si \(\text{err}=Ch^p\) alors \(\ln(\text{err})=\ln(C)+p\ln(h)\).
En échelle logarithmique, \(p\) représente donc la pente de la ligne droite \(\ln(\text{err})\).

Code 
plt.figure(figsize=(10,7))
# plt.plot(np.log(H), np.log(err_ee), 'r-o', label='Euler Explicite')
# plt.plot(np.log(H), np.log(err_ei), 'g-+', label='Euler Implicite')
plt.loglog(H, err_ee, 'r-o', label='Euler Explicite')
plt.loglog(H, err_ei, 'g-+', label='Euler Implicite')
# !!!!!! Ajouter ici les graphes des nouveaux schémas
plt.loglog(H, err_em,   'm-x', label='Euler Modife')
plt.loglog(H, err_rk1m, 'b-^', label='RK1_m')
plt.loglog(H, err_cn,   'y-d', label='Crank Nicolson')
plt.loglog(H, err_heun, 'c-v', label='Heun')
plt.xlabel(r'$h$')
plt.ylabel(r'$e$')
plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left')
plt.grid(True)

Pour estimer l’ordre de convergence on doit estimer la pente de la droite qui relie l’erreur au pas \(k\) à l’erreur au pas \(k+1\) en echelle logarithmique. Pour estimer la pente globale de cette droite (par des moindres carrés) on peut utiliser la fonction polyfit (du module numpy).

Code 
# Calcul des ordres de convergence = pente des droites de regression
# ==================================================================
order_ee = np.polyfit(np.log(H), np.log(err_ee), 1)[0]
order_ei = np.polyfit(np.log(H), np.log(err_ei), 1)[0]
# !!!!!! Ajouter ici les pentes des nouveaux schémas
order_em   = np.polyfit(np.log(H), np.log(err_em), 1)[0]
order_rk1m = np.polyfit(np.log(H), np.log(err_rk1m), 1)[0]
order_cn   = np.polyfit(np.log(H), np.log(err_cn), 1)[0]
order_heun = np.polyfit(np.log(H), np.log(err_heun), 1)[0]


# Print des ordres de convergence
# ===============================

# Version avec print
# print (f'Euler progressif : ordre = {order_ee:1.2f}')
# print (f'Euler regressif  : ordre = {order_ei:1.2f}')
# print (f'Euler modifié  : ordre = {order_em:1.2f}')
# print (f'RK1_M  : ordre = {order_rk1m:1.2f}')
# print (f'Cranck Nicolson  : ordre = {order_cn:1.2f}')
# print (f'Heun  : ordre = {order_heun:1.2f}')

# Version avec tabulate pour une affichage ameliore
from tabulate import tabulate
T = []
T.append( ["Schéma", "Ordre"] )
T.append( ["Euler progressif", order_ee] )
T.append( ["Euler regressif", order_ei] )
T.append( ["Euler modifié", order_em] )
T.append( ["RK1_M", order_rk1m] )
T.append( ["Cranck Nicolson", order_cn] )
T.append( ["Heun", order_heun] )
display(tabulate(T,headers="firstrow",tablefmt="html",floatfmt=".2f"))

# Version avec pandas pour une affichage ameliore
# import pandas as pd
# schemes = [
#     ("Euler progressif", order_ee),
#     ("Euler regressif", order_ei),
#     ("Euler modifié", order_em),
#     ("RK1_M", order_rk1m),
#     ("Cranck Nicolson", order_cn),
#     ("Heun", order_heun)
# ]
# T = pd.DataFrame(schemes, columns=["Schéma", "Ordre"])
# display(T)
Schéma Ordre
Euler progressif 1.05
Euler regressif 0.96
Euler modifié 2.08
RK1_M 2.00
Cranck Nicolson 2.01
Heun 2.08

Rappel sur la dérivée d’une fonction définie par une intégrale :

\( F(x)=\displaystyle\int_{u(x)}^{v(x)}g(x,t)\mathrm{d}t \qquad\leadsto\qquad F'(x)=\displaystyle g\left(x,v(x)\right)v'(x)-g\left(x,u(x)\right)u'(x)+ \int_{u(x)}^{v(x)} \frac{\partial}{\partial x} g(x,t)\mathrm{d}t. \)

La fonction \(f\) est solution du problème de Cauchy

\( \begin{cases} y'(t)= \sqrt{1-\frac{1}{4}\sin^2(t)},\\ y(0)=0. \end{cases} \)

En effet, \(f(0)=0\) et si on intègre l’EDO entre \(0\) et \(x\) on a

\( f(x)=\int_0^x \sqrt{1-\frac{1}{4}\sin^2(t)}\mathrm{d}t=\int_0^x y'(t)\mathrm{d}t=y(x)-y(0)=y(x). \)

Pour estimer la fonction \(f\) en \(\frac{\pi}{6}\) et \(\frac{\pi}{2}\), nous devons nous assurer que ces valeurs fassent partie des nœuds de discrétisation. Nous avons posé

  • \(t_0=0\) le point initial,
  • \(t_N=\frac{\pi}{2}\) le point final,
  • \(t_i\) les nœuds de discrétisation définis par \(t_i=t_0+ih\) avec le pas de discrétisation \(h=\frac{t_N-t_0}{N}=\frac{\pi}{2N}\).

Nous voulons que \(\frac{\pi}{6}\) soit l’un des nœuds \(t_i\), ce qui revient à chercher un indice \(i\) tel que \(t_i = \frac{\pi}{6}\).

En remplaçant \(t_i\) par sa définition, nous avons \(ih = \frac{\pi}{6}\), soit encore \(i \frac{\pi}{2N} = \frac{\pi}{6}\). On obtient \(i = \frac{2N}{6} = \frac{N}{3}\).

Cela montre que \(\frac{\pi}{6}\) sera un nœud de discrétisation si \(i = \frac{N}{3}\) est un entier. Autrement dit, \(N\) doit être un multiple de \(3\) pour garantir que \(i\) soit entier. Par ailleurs, \(t_N = \frac{\pi}{2}\) est automatiquement inclus comme nœud terminal. Ainsi, pour inclure \(\frac{\pi}{6}\) et \(\frac{\pi}{2}\) parmi les nœuds de discrétisation, il suffit de choisir \(N\) comme un multiple de \(3\).

On calcule enfin la solution approchée et on vérifie que \(f(\pi/6)\simeq0.517881934859938\) et \(f(\pi/2)\simeq1.46746220933943\) et on affiche le graphe de la fonction approchée.

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


def EE(phi, tt, y0):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt) - 1):
        uu[n+1] = uu[n] + h * phi(tt[n], uu[n])
    return uu

# EDO
phi = lambda t,y : np.sqrt(1-0.25*(np.sin(t))**2)

# CI
t0, y0 = 0, 0

# DISCRETISATION
tfinal = np.pi/2
k = 100
N = 3*k # ainsi t_N=pi/2=3pi/6 et t_k=pi/6
tt = np.linspace(t0,tfinal,N+1)

# CALCUL DE LA SOLUTION APPROCHÉE
uu = EE(phi,tt, y0)

print(f"f(pi/{np.pi/tt[k]}) = {uu[k]:g}") 
print(f"f(pi/{np.pi/tt[-1]}) = {uu[-1]:g}") 

# AFFICHAGE, NON DEMANDÉ
plt.plot(tt,uu)
plt.plot([tt[0],tt[k],tt[k]],[uu[k],uu[k],0],'m--o');
plt.plot([tt[0],tt[-1],tt[-1]],[uu[-1],uu[-1],0],'c--o');
f(pi/6.0) = 0.517965
f(pi/2.0) = 1.46781

Et sympy ? Il dit bien qu’il ne peut pas calculer cette intégrale et qu’il s’agit d’une intégrale elliptique de seconde espèce.

Code 
import sympy as sym
sym.init_printing()

t   = sym.Symbol('t')
y   = sym.Function('y')
edo = sym.Eq( sym.diff(y(t),t) , sym.sqrt(1-sym.Rational(1,4)*(sym.sin(t))**2) )
display(edo)
# solgen = sym.dsolve(edo)
# display(solgen)

t  = sym.Symbol('t')
f = sym.Integral(sym.sqrt(1-sym.Rational(1,4)*(sym.sin(t))**2))
display(f)
f.doit()

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = \sqrt{1 - \frac{\sin^{2}{\left(t \right)}}{4}}\)

\(\displaystyle \int \sqrt{1 - \frac{\sin^{2}{\left(t \right)}}{4}}\, dt\)

\(\displaystyle E\left(t\middle| \frac{1}{4}\right)\)

Il est cependant possible d’afficher la valeur “exacte” de l’intégrale elliptique de seconde espèce en utilisant la fonction elliptic_e du module sympy.

Code 
display( sym.elliptic_e(np.pi/6,1/4) ) # si on utilise sym.pi, il ne donne pas le résultat float
display( sym.elliptic_e(np.pi/2,1/4) ) 

# help(sym.elliptic_e)

\(\displaystyle 0.517881934859938\)

\(\displaystyle 1.46746220933943\)

Partie 1 : Invariant

Pour tout \(t\) on a

\( I'(t)=\left(x^2(t) + y^2(t)\right)' = 2x'(t)x(t)+2y'(t)y(t)=-2x(t)y(t)+2x(t)y(t)=0 \)

donc \(I(t)=I(0)=0^2+1^2=1\).

Partie 2 : Approximations

Code 
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams.update({'font.size': 12}) 
# from scipy.optimize import fsolve # non nécessaire car on résout l’équation implicite à la main

import sympy as sym
sym.init_printing()
Code 
t0 = 0
tfinal = 4*np.pi


N = 192 #*100
tt = np.linspace(0, tfinal, N+1)  # np.arange(t0, tfinal+h/2, h)


# Notation vectorielle
# ====================

# CI (2 conditions)
zz0 = np.array([1,0])

# EDO (2 equations) 
# phi : RxR^2 -> R^2
phi = lambda t,zz : np.array( [ -zz[1], zz[0] ] )

Pour simplifier les affichages, on définit une fonction affichage qui affiche les graphes de \(u_n\) et \(w_n\) en fonction de \(t_n\), le graphe de \(J_n\) et enfin le portrait de phase \((u_n,w_n)\). Il suffira alors d’appeler cette fonction pour chaque schéma.

Code 
def affichage(zz, JJ, tt, zz0, phi, schema):
    
    # zz[:, 0] = x(t) [vecteur uu]
    # zz[:, 1] = y(t) [vecteur ww]
    
    plt.figure(figsize=(24, 7))

    # Graphe de x(t) et y(t)
    plt.subplot(1, 3, 1)
    plt.plot(tt, zz[:, 0], label='x(t)')
    plt.plot(tt, zz[:, 1], label='y(t)')
    plt.grid()
    plt.xlabel('t')
    plt.legend()
    plt.title(f'{schema} - x(t) et y(t)')

    # Graphe de l'invariant JJ
    plt.subplot(1, 3, 2)
    plt.plot(tt, JJ)
    plt.grid()
    plt.xlabel('t')
    plt.title(f'{schema} - Invariant')

    # Diagramme de phase y(x)
    plt.subplot(1, 3, 3)
    u_min, u_max = min(zz[:, 0]), max(zz[:, 0])
    w_min, w_max = min(zz[:, 1]), max(zz[:, 1])
    Y1, Y2 = np.meshgrid(np.linspace(u_min, u_max, 20), np.linspace(w_min, w_max, 20))
    
    # Calcul du champ de vecteurs
    V = np.array([phi(0, [y1, y2]) for y1, y2 in zip(Y1.flatten(), Y2.flatten())])
    V1 = V[:, 0].reshape(Y1.shape)
    V2 = V[:, 1].reshape(Y2.shape)
    r = np.sqrt(V1**2 + V2**2)
    
    # Affichage du champ et des trajectoires
    plt.quiver(Y1, Y2, V1 / r, V2 / r)  
    plt.plot(zz[:, 0], zz[:, 1], 'r')
    plt.scatter(zz[0, 0], zz[0, 1], c='g', s=100)  # point initial
    plt.xlabel('x(t)')
    plt.ylabel('y(t)')
    plt.title(f'{schema} - y(x)')
    plt.axis('equal')
    plt.show()

Euler explicite

\( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+h\varphi_1(t_n,u_n,w_n),\\ w_{n+1}=w_n+h\varphi_2(t_n,u_n,w_n) \end{cases} \)

Code 
def EE(phi,tt,zz0):
    nb_eq = len(zz0) # nombre d'équations
    N = len(tt)-1 # N+1 points, donc N intervalles
    zz = np.zeros((N+1,nb_eq)) # colonnes = variables, lignes = temps
    zz[0] = zz0
    h = tt[1]-tt[0]
    for n in range(N):
        zz[n+1] = zz[n]+h*phi(tt[n],zz[n])
    return zz


# Version qui retourne un tableau 1D si zz0 est scalaire
# def EE(phi, tt, zz0):
#     zz0 = np.atleast_1d(zz0)  # Convertir zz0 en tableau (même si scalaire)
#     nb_eq = len(zz0)  # Nombre d'équations
#     N = len(tt) - 1  # N+1 points, donc N intervalles
#     zz = np.zeros((N + 1, nb_eq))  # Colonnes = variables, lignes = temps
#     zz[0] = zz0
#     h = tt[1] - tt[0]
#     for n in range(N):
#         zz[n + 1] = zz[n] + h * phi(tt[n], zz[n])
#     if nb_eq == 1:
#         return zz.flatten()  # Retourner un tableau 1D si scalaire
#     return zz

On remarque que

\( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n-hw_n,\\ w_{n+1}=w_n+hu_n \end{cases} \)

On a

\( J_{n+1} =u_{n+1}^2+w_{n+1}^2 =(u_n-hw_n)^2+(w_n+hu_n)^2 =(1+h^2)(u_n^2+w_n^2) =(1+h^2)J_n \)

soit encore

\( J_n=(1+h^2)^nJ_0=(1+h^2)^n. \)

On voit que

  • l’invariant n’est conservé que si \(h=0\), ce qui est impossible,
  • pour \(h\) fixé, \(\lim\limits_{n\to+\infty}J_n=+\infty\) (comportement pour \(t\to\infty\)),
  • pour \(n\) fixé, l’erreur \(|J_n-I_n|\) diminue comme \(h^2\) (convergence sur un intervalle fixé).
Code 
zz_ep = EE(phi,tt,zz0) # N+1 lignes, 2 colonnes
uu_ep, ww_ep = zz_ep.T
JJ_ep = uu_ep**2+ww_ep**2
affichage(zz_ep, JJ_ep, tt, zz0, phi, 'Euler Explicite')

Euler implicite

\( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+h\varphi_1(t_{n+1},u_{n+1},w_{n+1})=u_n-hw_{n+1},\\ w_{n+1}=w_n+h\varphi_2(t_{n+1},u_{n+1},w_{n+1})=w_n+hu_{n+1} \end{cases} \)

En résolvant le système linéaire on obtient une écriture explicite

\( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=\dfrac{u_n-hw_{n}}{1+h^2},\\ w_{n+1}=\dfrac{w_n+hu_{n}}{1+h^2}. \end{cases} \)

Code 
def EI(phi,tt,zz0): 
    nb_eq = len(zz0) # nombre d'équations
    N = len(tt)-1 # N+1 points, donc N intervalles
    zz = np.zeros((N+1,nb_eq))
    zz[0] = zz0
    h = tt[1]-tt[0]
    for n in range(N):
        zz[n+1,0] = (zz[n,0]-h*zz[n,1])/(1+h**2)
        zz[n+1,1] = (zz[n,1]+h*zz[n,0])/(1+h**2)
    return zz
Code 
#  VERSION AVEC fsolve (non demandé)
from scipy.optimize import fsolve

def EI_bis(phi,tt,zz0):
    nb_eq = len(zz0) # nombre d'équations
    N = len(tt)-1 # N+1 points, donc N intervalles
    zz = np.zeros((N+1,nb_eq))
    zz[0] = zz0
    h = tt[1]-tt[0]
    for n in range(N):
        system = lambda z : z - zz[n] - h*phi(tt[n+1],z)
        zz[n+1] = fsolve( system , zz[n] ) 
    return zz


# Version qui retourne un tableau 1D si zz0 est scalaire
# def EI_bis(phi, tt, zz0):
#     zz0 = np.atleast_1d(zz0)  # Convertir zz0 en tableau (même si scalaire)
#     nb_eq = len(zz0)  # Nombre d'équations
#     N = len(tt) - 1  # N+1 points, donc N intervalles
#     zz = np.zeros((N + 1, nb_eq))
#     zz[0] = zz0
#     h = tt[1] - tt[0]
#     for n in range(N):
#         system = lambda z: z - zz[n] - h * phi(tt[n + 1], z)
#         zz[n + 1] = fsolve(system, zz[n])
#     if nb_eq == 1:
#         return zz.flatten()  # Retourner un tableau 1D si scalaire
#     return zz

On a

\( J_{n+1} =u_{n+1}^2+w_{n+1}^2 =\frac{(u_n-hw_n)^2+(w_n+hu_n)^2}{(1+h^2)^2} =(u_n^2+w_n^2)\frac{1+h^2}{(1+h^2)^2} =\frac{1+h^2}{(1+h^2)^2}J_n =\frac{1}{1+h^2}J_n \)

soit encore

\( J_n=\left(\frac{1}{1+h^2}\right)^nJ_0=\left(\frac{1}{1+h^2}\right)^n. \)

On voit que

  • l’invariant n’est conservé que si \(h=0\), ce qui est impossible,
  • pour \(h\) fixé, \(\lim\limits_{n\to+\infty}J_n=0^+\),
  • pour \(n\) fixé, l’erreur \(|J_n-I_n|\) diminue comme \(h^2\).
Code 
# zz_er = EI(phi,tt,zz0)
zz_er = EI_bis(phi,tt,zz0)
uu_er, ww_er = zz_er.T
JJ_er = uu_er**2+ww_er**2
affichage(zz_er, JJ_er, tt, zz0, phi, 'Euler Implicite')

Crank Nicolson

\( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+\frac{h}{2}\Big(\varphi_1(t_{n},u_{n},w_{n})+\varphi_1(t_{n+1},u_{n+1},w_{n+1})\big)=u_n-\frac{h}{2}w_{n}-\frac{h}{2}w_{n+1},\\ w_{n+1}=w_n+\frac{h}{2}\Big(\varphi_2(t_{n},u_{n},w_{n})+\varphi_2(t_{n+1},u_{n+1},w_{n+1})\big)=w_n+\frac{h}{2}u_{n}+\frac{h}{2}u_{n+1}. \end{cases} \)

En résolvant le système linéaire on obtient une écriture explicite :

\( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=\dfrac{(4-h^2)u_n-4hw_{n}}{4+h^2},\\ w_{n+1}=\dfrac{(4-h^2)w_n+4hu_{n}}{4+h^2}. \end{cases} \)

De plus, on a

\( \begin{aligned} J_{n+1} &=u_{n+1}^2+w_{n+1}^2 \\ &=\frac{((4-h^2)u_n-4hw_n)^2+((4-h^2)w_n+4hu_n)^2}{(4+h^2)^2} \\ &=\frac{(4-h^2)^2u_n^2+16h^2w_n^2-8(4-h^2)hu_nw_n+(4-h^2)^2w_n^2+16h^2u_n^2+8(4-h^2)hu_nw_n}{(4+h^2)^2} \\ &=\frac{\big((4-h^2)^2+16h^2\big)(u_n^2+w_n^2)}{(4+h^2)^2} \\ &=\frac{(4+h^2)^2(u_n^2+w_n^2)}{(4+h^2)^2} \\ &=(u_n^2+w_n^2)=J_n \end{aligned} \)

On voit que l’invariant est conservé pour tout \(h\)

Vérifions ces calculs avec sympy:

Code 
dt = sym.Symbol('h') # pour ne pas redefinir h comme un symbole, on écrit dt mais on affiche h
u_n = sym.Symbol('u_n')
u_np1 = sym.Symbol('u_{n+1}')
w_n = sym.Symbol('w_n')
w_np1 = sym.Symbol('w_{n+1}')

eq1 = sym.Eq( u_np1 , u_n-dt/2*w_n-dt/2*w_np1 )
eq2 = sym.Eq( w_np1 , w_n+dt/2*u_n+dt/2*u_np1 )
sol = sym.solve([eq1,eq2],[u_np1,w_np1],dict=True)[0]
display(sol)

J_np1 = u_np1**2+w_np1**2
J_np1.subs(sol).simplify()

\(\displaystyle \left\{ u_{n+1} : \frac{- h^{2} u_{n} - 4 h w_{n} + 4 u_{n}}{h^{2} + 4}, \ w_{n+1} : \frac{- h^{2} w_{n} + 4 h u_{n} + 4 w_{n}}{h^{2} + 4}\right\}\)

\(\displaystyle u_{n}^{2} + w_{n}^{2}\)

Code 
def cn(phi,tt,zz0):
    nb_eq = len(zz0) # nombre d'équations
    N = len(tt)-1 # N+1 points, donc N intervalles
    zz = np.zeros((N+1,nb_eq))
    zz[0] = zz0
    h = tt[1]-tt[0]
    for n in range(N):
        zz[n+1,0] = ((4-h**2)*zz[n,0]-4*h*zz[n,1])/(4+h**2)
        zz[n+1,1] = ((4-h**2)*zz[n,1]+4*h*zz[n,0])/(4+h**2)
    return zz
Code 
#  VERSION AVEC fsolve (non demandé)
from scipy.optimize import fsolve

def cn_bis(phi,tt,zz0):
    nb_eq = len(zz0) # nombre d'équations
    N = len(tt)-1 # N+1 points, donc N intervalles
    zz = np.zeros((N+1,nb_eq))
    zz[0] = zz0
    h = tt[1]-tt[0]
    for n in range(N):
        system = lambda z : z - zz[n] - h/2*( phi(tt[n],zz[n]) + phi(tt[n+1],z) )
        zz[n+1] = fsolve( system , zz[n] ) 
    return zz

# Version qui retourne un tableau 1D si zz0 est scalaire
# def cn_bis(phi, tt, zz0):
#     zz0 = np.atleast_1d(zz0)  # Convertir zz0 en tableau (même si scalaire)
#     nb_eq = len(zz0)  # Nombre d'équations
#     N = len(tt) - 1  # N+1 points, donc N intervalles
#     zz = np.zeros((N + 1, nb_eq))
#     zz[0] = zz0
#     h = tt[1] - tt[0]
#     for n in range(N):
#         system = lambda z: z - zz[n] - h / 2 * (phi(tt[n], zz[n]) + phi(tt[n + 1], z))
#         zz[n + 1] = fsolve(system, zz[n])
#     if nb_eq == 1:
#         return zz.flatten()  # Retourner un tableau 1D si scalaire
#     return zz
Code 
# zz_cn = cn(phi,tt,zz0)
zz_cn = cn_bis(phi,tt,zz0)
uu_cn, ww_cn = zz_cn.T
JJ_cn = uu_cn**2+ww_cn**2
affichage(zz_cn, JJ_cn, tt, zz0, phi, 'Crank Nicolson')

Avec odeint du module scipy.integrate:

Code 
from scipy.integrate import odeint

zz_odeint = odeint(phi,zz0,tt,tfirst=True)
uu_odeint, ww_odeint = zz_odeint.T

JJ_odeint = uu_odeint**2+ww_odeint**2
affichage(zz_odeint, JJ_odeint, tt, zz0, phi, 'odeint')

Bonus : avec solve_ivp du module scipy.integrate:

Code 
from scipy.integrate import solve_ivp
Code 
sol_RK45 = solve_ivp(phi,[t0,tfinal],zz0,t_eval=tt, method='RK45')
uu_RK45, ww_RK45 = sol_RK45.y

JJ_RK45 = uu_RK45**2+ww_RK45**2

affichage(sol_RK45.y.T, JJ_RK45, tt, zz0, phi, 'RK45')

Code 
sol_Radau = solve_ivp(phi,[t0,tfinal],zz0,t_eval=tt, method='Radau')
uu_Radau, ww_Radau = sol_Radau.y

JJ_Radau = uu_Radau**2+ww_Radau**2

affichage(sol_Radau.y.T, JJ_Radau, tt, zz0, phi, 'Radau')

Code 
sol_BDF = solve_ivp(phi,[t0,tfinal],zz0,t_eval=tt, method='BDF')
uu_BDF, ww_BDF = sol_BDF.y

JJ_BDF = uu_BDF**2+ww_BDF**2

affichage(sol_BDF.y.T, JJ_BDF, tt, zz0, phi, 'BDF')

Code 
sol_LSODA = solve_ivp(phi,[t0,tfinal],zz0,t_eval=tt, method='LSODA')
uu_LSODA, ww_LSODA = sol_LSODA.y

JJ_LSODA = uu_LSODA**2+ww_LSODA**2

affichage(sol_LSODA.y.T, JJ_LSODA, tt, zz0, phi, 'LSODA')

Partie 3 : Solution exacte

La solution exacte est

\(\begin{align*} x(t)&=\cos(t),\\ y(t)&=\sin(t). \end{align*}\)

Cette solution est périodique de période \(2\pi\) et on a bien

\( I(t)=x^2(t) + y^2(t) = 1 \quad \forall t. \)

Vérifions-le avec le module sympy:

Code 
t = sym.Symbol('t')
x = sym.Function('x')
y = sym.Function('y')

t_0 = 0
x_0 = zz0[0]
y_0 = zz0[1]


edo1 = sym.Eq( sym.diff(x(t),t) , -y(t) )
edo2 = sym.Eq( sym.diff(y(t),t) , x(t) )
display(edo1)
display(edo2)


solpar_1, solpar_2 = sym.dsolve([edo1,edo2],[x(t),y(t)], ics={x(t_0): x_0, y(t_0): y_0})


# AFFICHAGE
func_1 = sym.lambdify(t,solpar_1.rhs,'numpy')
func_2 = sym.lambdify(t,solpar_2.rhs,'numpy')

tt = np.linspace(0,4*np.pi,193)
xx = func_1(tt)
yy = func_2(tt)
JJ = xx**2+yy**2

affichage(np.array([xx,yy]).T, JJ, tt, zz0, phi, 'Solution exacte')

\(\displaystyle \frac{d}{d t} x{\left(t \right)} = - y{\left(t \right)}\)

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = x{\left(t \right)}\)