CM 1 - Introduction à l’approximation numérique d’EDO

Auteur·rice

Gloria FACCANONI

Date de publication

2 mars 2025

1 Position du problème

Problème de Cauchy:

trouver une fonction \(y \colon I\subset \mathbb{R} \to \mathbb{R}\) définie sur un intervalle \(I\) telle que

\( \begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in I=]t_0,T[,\\ y(t_0) = y_0, \end{cases} \)

avec \(y_0\) une valeur donnée et supposons que l’on ait montré l’existence et l’unicité d’une solution \(y\) pour \(t\in I\).

Nous disposons donc de :

une valeur initiale \(y_0\) de \(y(t)\) à un certain temps \(t_0\),
l’EDO elle-même qui relie le changement de \(y\) par rapport à \(t\) à une fonction \(\varphi(t, y)\), c’est-à-dire \(y'(t) = \varphi(t, y(t))\).

Exemple : rechercher la fonction \(t\mapsto y(t)\) pour \(t\) entre \(t_0 = 0\) et \(T = 2\) telle que \(y'(t) = -2y(t)\) avec \(y_0 = y(t_0) = 17\). La solution exacte existe, est unique, et est \(y(t) = 17 \exp(-2t)\).

On ne peut expliciter des solutions analytiques que pour des équations différentielles ordinaires très particulières.

Dans certains cas, on ne peut exprimer la solution que sous forme implicite. C’est le cas par exemple de l’EDO \(y'(t)=\dfrac{y(t)-t}{y(t)+t}\) dont les solutions vérifient la relation implicite

\(\frac{1}{2}\ln(t^2+y^2(t))+\arctan\left( \frac{y(t)}{t} \right)=C, \)

où \(C\) est une constante arbitraire.

Dans d’autres cas, on ne parvient même pas à représenter la solution sous forme implicite. C’est le cas par exemple de l’EDO \(y'(t)=e^{-t^2}\) dont les solutions ne peuvent pas s’écrire comme composition de fonctions élémentaires.

Lorsqu’on ne connait pas de solution exacte à un problème de Cauchy qui admet une et une seule solution, on essaye d’en avoir une bonne approximation par des méthodes numériques.

2 Solution exacte discrète VS solution approchée

L’idée la plus simple pour résoudre de manière approchée un problème de Cauchy est de discrétiser l’intervalle de temps avec un pas \(h\) et d’approcher la dérivée temporelle sur chaque sous-intervalle de longueur \(h\).

Pour \(h>0\) soit \(t_n\equiv t_0+nh\) avec \(n=0,1,2,\dots,N\) une suite de \(N+1\) nœuds de \(I\) induisant une discrétisation de \(I\) en \(N\) sous-intervalles \(I_n=[t_n;t_{n+1}]\) chacun de longueur \(h=\frac{T-t_0}{N}>0\) (appelé le pas de discrétisation).

Pour chaque nœud \(t_n\), on cherche la valeur inconnue \(u_n\) qui approche la valeur exacte \(y_n\equiv y(t_n)\).

L’ensemble de \(N+1\) valeurs \(\{t_0, t_1=t_0+h,\dots , t_{N}=T \}\) représente les points de la discrétisation.
L’ensemble de \(N+1\) valeurs \(\{y_0, y_1,\dots , y_{N} \}\) représente la solution exacte discrète.
L’ensemble de \(N+1\) valeurs \(\{u_0 = y_0, u_1,\dots , u_{N} \}\) représente la solution numérique.

Code

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">14</span>})

def euler_explicit(f, y0, t_values):
    h = t_values[1] - t_values[0]
    y_values = np.zeros_like(t_values)
    y_values[0] = y0
    for i in range(len(t_values) - 1):
        y_values[i+1] = y_values[i] + h * f(t_values[i], y_values[i])
    return y_values

def plot_solutions(n, show_exact_line=True, show_exact_points=True, show_approx_points=True):
    t0 = 0
    tn = 2
    y0 = 17
    phi = lambda t, y : -2 * y
    exact_solution = lambda t : 17 * np.exp(-2 * t)

    # Calcul des solutions
    t_exact = np.linspace(t0, tn, 100)  # pour la solution exacte continue
    y_exact = exact_solution(t_exact)
    t_exact_discrete = np.linspace(t0, tn, n+1)  # pour les solutions exacte et approchée discrètes
    y_exact_discrete = exact_solution(t_exact_discrete)
    y_approx = euler_explicit(phi, y0, t_exact_discrete)  

    plt.figure(figsize=(14, 6))
    plt.xlabel('t')
    plt.ylabel('y')
    plt.title('Comparaison des solutions')

    if show_exact_line:
        plt.plot(t_exact, y_exact, label='Solution exacte')

    if show_exact_points:
        plt.scatter(t_exact_discrete, y_exact_discrete, color='red', label='Solution exacte discrète')
        for i in range(len(y_exact_discrete)):
            plt.annotate('$y_{'</span><span class="op">+</span><span class="bu">str</span>(i)<span class="op">+</span><span class="st">'}$', (t_exact_discrete[i], y_exact_discrete[i] + 0.5), color='red')

    if show_approx_points:
        plt.plot(t_exact_discrete, y_approx, 'g--')
        plt.scatter(t_exact_discrete, y_approx, color='green', label='Solution approchée')
        for i in range(len(y_approx)):
            plt.annotate('$u_{'</span><span class="op">+</span><span class="bu">str</span>(i)<span class="op">+</span><span class="st">'}$', (t_exact_discrete[i]-0.04, y_approx[i] - 0.75), color='green')

    plt.legend()
    plt.grid(True)
    plt.show()


# si interactif
# from ipywidgets import interact, widgets
# interact(plot_solutions, n=widgets.IntSlider(min=5, max=20, step=1, value=10), 
#          show_exact_line=True, show_exact_points=False, show_approx_points=False);

# si non interactif
plot_solutions(10, show_exact_line=True, show_exact_points=True, show_approx_points=True);

3 Construction élémentaire des méthodes d’Euler explicite et implicite

Une méthode classique, la méthode d’Euler explicite (ou progressive, de l’anglais forward), est obtenue en considérant l’équation différentielle en chaque nœud \(t_n\) et en remplaçant la dérivée exacte \(y'(t_n)\) par le taux d’accroissement

\( \varphi(t_n,y(t_n))=y'(t_n)=\frac{y(t_{n}+h)-y(t_n)}{h}+\mathscr{O}(h). \)

Cela permet de construire une solution numérique par une suite récurrente:

Schéma EE:

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

On peut donner une intérprétation géométrique de la méthode d’Euler explicite : en un point donné \(t_n\), on avance d’un pas \(h\) le long de la droite tangente à la solution (inconnue).

De même, en utilisant le taux d’accroissement

\( \varphi(t_{n+1},y(t_{n+1}))=y'(t_{n+1})\simeq\frac{y(t_{n+1})-y(t_n)}{h} \)

pour approcher \(y'(t_{n+1})\), on obtient la méthode d’Euler implicite (ou rétrograde, de l’anglais backward)

Schéma EI:

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

Ces deux méthodes sont dites à un pas: pour calculer la solution numérique \(u_{n+1}\) au nœud \(t_{n+1}\), on a seulement besoin des informations disponibles au nœud précédent \(t_n\). Plus précisément,

pour la méthode d’Euler progressive, \(u_{n+1}\) ne dépend que de la valeur \(u_n\) calculée précédemment,
pour la méthode d’Euler rétrograde, \(u_{n+1}\) dépend aussi “de lui-même” à travers la valeur de \(\varphi(t_{n+1},u_{n+1})\).

C’est pour cette raison que la méthode d’Euler progressive est dite explicite tandis que la méthode d’Euler rétrograde est dite implicite.
Les méthodes implicites sont plus coûteuses que les méthodes explicites car, si la fonction \(\varphi\) est non linéaire, un problème non linéaire doit être résolu à chaque temps \(t_{n+1}\) pour calculer \(u_{n+1}\). Néanmoins, nous verrons que les méthodes implicites jouissent de meilleures propriétés de stabilité que les méthodes explicites.

4 Implémentation des schémas d’Euler explicite et implicite

Voyons un exemple complet.

Exemple: considérons le problème de Cauchy

\( \begin{cases} y'(t) = 2ty(t), &\forall t \in I=[0,1],\\ y(0) = 1. \end{cases} \)

Sachant que la solution est \(y(t)=e^{t^2}\), estimer la qualité du schéma en affichant la plus grande erreur commise (en valeur absolue) et varier le pas de discrétisation \(h\) pour voir son impact sur l’erreur.

On commence par importer

la fonction fsolve du module scipy.optimize pour résoudre les équations implicites présentes dans le schéma implicite.
le module matplotlib et le module numpy:
- soit en important séparemment matplotlib.pyplot et numpy (avec les alias classiques);
- soit en important juste matplotlib.pylab, ce qui importe à la fois matplotlib et numpy (avec par exemple l’alias commun pyl ou sans alias si on utilise *).

Rappels: numpy rédéfinit toutes les fonctions mathématiques du module math et ces fonctions sont vectorisées (e.g. on pourra écrire directement yy=sin(xx) avec xx une liste/tuple/array au lieu d’écrire yy=[sin(x) for x in xx]).

Code

%reset -f
%matplotlib inline

from scipy.optimize import fsolve # pour les schémas implicites

#from matplotlib.pylab import *

import numpy as np
import matplotlib.pyplot as plt
# plt.rcdefaults()
plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})

On initialise le problème de Cauchy

Code

t0     = 0
tfinal = 1
y0     = 1

On définit l’équation différentielle : phi est une fonction python qui contient la fonction mathématique \(\varphi\colon [t_0,T]\times\mathbb{R}\to\mathbb{R}\) définie par \(\varphi(t, y)=2ty\). Bien noter que \(\varphi\) dépend de deux variables: \(t\) et \(y\). Dans l’EDO on a \(\varphi(t,y(t))\): il s’agit d’une fonction d’une seule variable obtenue par composition de la fonction \(\varphi\) avec la fonction \(t\mapsto y(t)\).

Code

phi = lambda t, y : 2*y*t

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

Code

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

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 progressif : \(\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 euler_progressif(phi,tt,y0):
    h  = tt[1]-tt[0]
    N  = len(tt)-1 # len(tt) = N+1 car tt contient N+1 points
    uu = np.zeros_like(tt) # array de même taille que tt et qui ne contient que des zéros
    uu[0] = y0
    for n in range(N): # range(N) = 0,1,...,N-1                  
        uu[n+1] = uu[n]+h*phi(tt[n],uu[n])
    return uu

Schéma d’Euler régressif : \(\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éral non linéaire) \(x\mapsto -x+u_n+h\varphi(t_{n+1},x)\)
la fonction fsolve du module scipy.optimize requiert deux paramètres : une fonction et un point de départ. Elle renvoie une liste avec la valeur du zéro approché (c’est une liste, car fsolve résout un système, ici avec une seule équation). Un choix naturel pour le point de départ est la valeur précédente \(u_n\) ou la valeur prédite par la méthode explicite \(u_n+h\varphi(t_n,u_n)\).

Code

def euler_regressif(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.zeros_like(tt) # array de même taille que tt et qui ne contient que des zéros
    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] # fsolve renvoie un array qui ne contient qu'un seul élément
    return uu

On calcule les solutions approchées:

Code

uu_ep = euler_progressif(phi,tt,y0)
uu_er = euler_regressif(phi,tt,y0)

Comme on la connait, on définit la solution exacte pour calculer les erreurs:

Code

sol_exacte = lambda t : y0*np.exp(t**2)
yy = sol_exacte(tt) # = [sol_exacte(t) for t in tt] 


# Calcul de l'erreur 

# On peut utiliser la norme infinie : on cherche la plus grande différence entre les deux vecteurs yy et uu en valeur absolue
err_ep = np.linalg.norm(yy-uu_ep,np.inf) # max(np.abs(yy-uu_ep))
err_er = np.linalg.norm(yy-uu_er,np.inf) # max(np.abs(yy-uu_er))

# On peut aussi utiliser la norme 2 : on somme les carrés des différences, on prend la racine carrée et on divise par la norme de yy
# err_ep = np.linalg.norm(yy-uu_ep)/np.linalg.norm(yy)
# err_er = np.linalg.norm(yy-uu_er)/np.linalg.norm(yy)

On compare les graphes des solutions exacte (en bleu) et approchées (en rouge) et on affiche le maximum de l’erreur:

Code

plt.figure(figsize=(9, 7))

plt.plot(tt,yy   ,'b-*',label='Exacte')
plt.plot(tt,uu_ep,'r:+',label=f'Approchée par Euler Explicite, erreur max = {</span>err_ep<span class="sc">:g}')
plt.plot(tt,uu_er,'m:v',label=f'Approchée par Euler Implicite, erreur max = {</span>err_er<span class="sc">:g}')
plt.grid()
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');

5 Convergence empirique des schémas d’Euler

Définition: schéma convergent et ordre de convergence
Une méthode numérique est convergente si

\( |y_n-u_n|\le C(h)\xrightarrow[h\to0]{}0 \qquad\forall n=0,\dots,N \)

Si \(C(h) = \mathcal{O}(h^p)\) pour \(p > 0\), on dit que la convergence de la méthode est d’ordre \(p\).

Remarque : \(N\to+\infty\) lorsque \(h\to0\).

Exemple: considérons à nouveau le problème de Cauchy

\( \begin{cases} y'(t) = 2ty(t), &\forall t \in I=[0,1],\\ y(0) = 1. \end{cases} \)

On se propose d’estimer l’ordre de convergence des méthodes d’Euler.

Pour chaque schéma, on calcule la solution approchée avec différentes valeurs de \(h\) (ce qui correspond à différentes valeurs de \(N\)). On sauvegarde les valeurs de \(h\) dans le vecteur H.

Pour chaque valeur de \(h\), 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}|\) avec \(N_k\).

Code

H = []
err_ep = []
err_er = []

from numpy.linalg import norm

for k in range(1,7):
    tt_tmp = np.linspace(t0,tfinal,k*N+1) # on raffine la grille
    H.append(tt_tmp[1]-tt_tmp[0])
    yy = sol_exacte(tt_tmp)
    uu_ep_tmp = euler_progressif(phi,tt_tmp,y0)
    uu_er_tmp = euler_regressif(phi,tt_tmp,y0)    
    err_ep.append( norm(uu_ep_tmp-yy,np.inf) ) # idem que max(abs(uu_ep-yy))
    err_er.append( norm(uu_er_tmp-yy,np.inf) ) # idem que max(abs(uu_er-yy))

Pour afficher l’ordre de convergence on affiche les points (h[k],err_ep[k]) en échelle logarithmique: 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)\), autrement dit, en échelle logarithmique, \(p\) représente la pente de la ligne droite \(\ln(\text{err})\).

Code

plt.figure(figsize=(10, 7))
plt.loglog(H,err_ep, 'r-+', label='Euler Explicite')
plt.loglog(H,err_er, 'g-v', label='Euler Implicite')
plt.xlabel('$h$')
plt.ylabel('erreur')
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 échelle logarithmique. Pour estimer la pente globale de cette droite (par des moindres carrés) on peut utiliser la fonction polyfit (du module numpy).

Code

# ln(e) = a ln(h) + b
a_ep, b_ep = np.polyfit(np.log(H),np.log(err_ep), 1) # polyfit ( [liste des abscisses], [liste des ordonnées], degré du polynome)
print (f'Euler progressif: ordre {</span>a_ep <span class="sc">:1.2f}')
a_er, b_er = np.polyfit(np.log(H),np.log(err_er), 1)
print (f'Euler regressif: ordre {</span>a_er<span class="sc">:1.2f}')

Euler progressif: ordre 0.91
Euler regressif: ordre 1.12

On peut bien-sûr afficher la droite obtenue par régression linéaire en même temps que les points.

Soit on affiche les logarithmes des données et l’équation de la droite avec la commande plot:

Code

plt.figure(figsize=(18, 7))

plt.subplot(1,2,1)
plt.plot(np.log(H),np.log(err_ep), 'ro',label='Euler Explicite')
plt.plot(np.log(H),[a_ep*np.log(h)+b_ep for h in H])
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.legend(loc='upper left')
plt.grid(True);

plt.subplot(1,2,2)
plt.plot(np.log(H),np.log(err_er), 'ro',label='Euler Implicite')
plt.plot(np.log(H),[a_er*np.log(h)+b_er for h in H])
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.legend(loc='upper left')
plt.grid(True);

Soit on affiche les données en échelle logarithmique et l’équation de l’erreur \(Ch^p\) avec \(C=\ln(b)\) et \(p=a\) avec l’instruction loglog:

Code

plt.figure(figsize=(18, 7))

plt.subplot(1,2,1)
plt.loglog(H,err_ep, 'ro',label='Euler Explicite')
plt.loglog(H,[np.exp(b_ep)*(h**a_ep) for h in H])
plt.xlabel('$h$')
plt.ylabel('$e$')
plt.legend(loc='upper left')
plt.grid(True);

plt.subplot(1,2,2)
plt.loglog(H,err_er, 'ro',label='Euler Implicite')
plt.loglog(H,[np.exp(b_er)*(h**a_er) for h in H])
plt.xlabel('$h$')
plt.ylabel('$e$')
plt.legend(loc='upper left')
plt.grid(True);

6 Annexe: implémentation d’un schéma “moyen” des deux schémas précédents

On obtient la méthode de Crank-Nicolson (ou du trapèze)

Schéma CN:

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

Code

def CN(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.zeros_like(tt) # array de même taille que tt
    uu[0] = y0
    for i in range(len(tt)-1):
        temp = fsolve( lambda x: -x+uu[i]+h/2*(phi(tt[i],uu[i])+phi(tt[i+1],x)) , uu[i] )
        uu[i+1] = temp[0]
        # Rmk: pour ce phi linéaire en y, on peut calculer uu[i+1] directement
        # uu[i+1] = (1+h*tt[i])/(1-h*tt[i+1])*uu[i]
    return uu


N  = 8
tt = np.linspace(t0,tfinal,N+1)
uu_cn = CN(phi,tt,y0)
yy = sol_exacte(tt) 

plt.figure(figsize=(9, 7))

plt.plot(tt,yy   ,'b-',label='Exacte')
plt.plot(tt,uu_ep,'r:+',label=f'Approchée par Euler Explicite, erreur max = {</span><span class="bu">max</span>(np.<span class="bu">abs</span>(uu_ep<span class="op">-</span>yy))<span class="sc">:g}')
plt.plot(tt,uu_er,'m:v',label=f'Approchée par Euler Implicite, erreur max = {</span><span class="bu">max</span>(np.<span class="bu">abs</span>(uu_er<span class="op">-</span>yy))<span class="sc">:g}')
plt.plot(tt,uu_cn,'go',label=f'Approchée par Crank-Nicolson, erreur max = {</span><span class="bu">max</span>(np.<span class="bu">abs</span>(uu_cn<span class="op">-</span>yy))<span class="sc">:g}')
plt.grid()
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');

Code

H_cn = []
err_cn = []


for k in range(1,7):
    tt = np.linspace(t0,tfinal,k*N+1) # on raffine la grille
    h  = tt[1]-tt[0]
    yy = sol_exacte(tt)
    uu_cn = CN(phi,tt,y0)
    H_cn.append(h)
    err_cn.append( np.linalg.norm(uu_cn-yy) )

# ln(e) = a ln(h) + b
a_cn, b_cn = np.polyfit(np.log(H_cn),np.log(err_cn), 1) 
print (f'Crank Nicolson: ordre {</span>a_cn <span class="sc">:1.2f}') 


plt.figure(figsize=(10, 7))
plt.loglog(H, err_ep, 'o',label=f'Euler Explicite : pente = {</span>a_ep<span class="sc">:1.2f}')
plt.loglog(H,[np.exp(b_ep)*(h**a_ep) for h in H])
plt.loglog(H,err_er, 'd',label=f'Euler Implicite : pente = {</span>a_er<span class="sc">:1.2f}')
plt.loglog(H,[np.exp(b_er)*(h**a_er) for h in H])
plt.loglog(H_cn,err_cn, 'v',label=f'Cranc Nicolson : pente = {</span>a_cn<span class="sc">:1.2f}')
plt.loglog(H_cn,[np.exp(b_cn)*(h**a_cn) for h in H_cn])

plt.xlabel('$h$')
plt.ylabel('$e$')
plt.legend(loc='upper left',bbox_to_anchor=(1.04,1))
plt.grid(True);

Crank Nicolson: ordre 1.62

Retour au sommet