Introduction à l’approximation numérique d’EDO

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’existence et l’unicité d’une solution \(y\) sur \(I\) aient été établies.

Nous disposons donc de :

• une valeur initiale \(y_0\) pour \(y(t)\) à un certain instant \(t_0\),
• l’EDO elle-même \(y'(t) = \varphi(t, y(t))\) qui relie la variation de \(y\) par rapport à \(t\) à la fonction \(\varphi\).

Exemple : rechercher la fonction \(t\mapsto y(t)\) pour \(t\in[0;2]\) telle que \(y'(t) = -2y(t)\) et \(y_0 = y(t_0) = 17\). La solution existe, est unique, et est donnée par \(y(t) = 17 \exp(-2t)\).

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.

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.

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 14})

def euler_explicit(f, y0, t_values):
    h = t_values[1] - t_values[0]
    y_values = np.empty_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')

    plt.scatter(t_exact_discrete, 0*t_exact_discrete, label='Discrétisation', color='black', marker='x')
    for i in range(len(t_exact_discrete)):
        plt.annotate('$t_{'+str(i)+'}$', (t_exact_discrete[i]-0.04, - 0.75), color='black')

    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_{'+str(i)+'}$', (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_{'+str(i)+'}$', (t_exact_discrete[i]-0.04, y_approx[i] - 0.75), color='green')

    plt.xticks([0])
    plt.yticks([0])
    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);

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 “à droite” comme suit :

\( \begin{cases} y'(t_n) = \varphi(t_n,y(t_n))\\ y'(t_n) = \displaystyle \lim_{k\to 0}\frac{y(t_n+k)-y(t_n)}{k} = \frac{y(t_{n}+k)-y(t_n)}{k}+\mathscr{O}(k) \end{cases} \) Pour \(k=h>0\), on obtient \( \varphi(t_n,y(t_n)) = \frac{y(t_{n+1})-y(t_n)}{h}+\mathscr{O}(h) \quad\leadsto\quad \varphi(t_n,y(t_n)) \approx \frac{y(t_{n+1})-y(t_n)}{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).

Une autre méthode classique est la méthode d’Euler implicite (ou rétrograde, de l’anglais backward), qui s’obtient en considérant l’équation différentielle en chaque nœud \(t_{n+1}\) et en remplaçant la dérivée exacte \(y'(t_{n+1})\) par le taux d’accroissement “à gauche” :

\( \begin{cases} y'(t_{n+1}) = \varphi(t_{n+1},y(t_{n+1}))\\ y'(t_{n+1}) = \displaystyle \lim_{k\to 0}\frac{y(t_{n+1}+k)-y(t_{n+1})}{k} = \frac{y(t_{n+1}+k)-y(t_{n+1})}{k}+\mathscr{O}(k) \end{cases} \)

Pour \(k = -h <0\), on obtient \( \varphi(t_{n+1},y(t_{n+1})) = \frac{y(t_{n})-y(t_{n+1})}{-h}+\mathscr{O}(h) \quad\leadsto\quad \varphi(t_{n+1},y(t_{n+1})) \approx \frac{y(t_{n+1})-y(t_n)}{h}. \) ce qui donne encore une suite récurrente pour construire une solution numérique :

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}\). Pour cela on peut utiliser des méthodes itératives comme la méthode de la dichotomie ou une méthode de point fixe comme la méthode de Newton. Néanmoins, nous verrons que les méthodes implicites jouissent de meilleures propriétés de stabilité que les méthodes explicites.

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

1. la fonction fsolve du module scipy.optimize pour résoudre les équations implicites présentes dans le schéma implicite ;
2. le module matplotlib et le module numpy:
   • soit en important séparément matplotlib.pyplot et numpy (avec les alias classiques plt et np respectivement) ;
   • soit (déconseillé) en important 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 redéfinit toutes les fonctions mathématiques du module math et ces fonctions sont vectorisées (e.g. on pourra écrire directement yy = np.sin(xx) avec xx une liste/tuple/array au lieu d’écrire yy = [np.sin(x) for x in xx]).

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({'font.size': 14})

On initialise le problème de Cauchy

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

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

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

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.empty_like(tt) # array vide de même taille que tt
    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)\).

def euler_regressif(phi,tt,y0):
    h  = tt[1]-tt[0]
    N = len(tt)-1 # len(tt) = N+1 car tt contient N+1 points
    uu = np.empty_like(tt) # array vide de même taille que tt
    uu[0] = y0
    for n in range(N): # range(N) = 0,1,...,N-1
        fct_implicit = lambda x: -x + uu[n] + h*phi(tt[n+1], x)
        temp = fsolve( fct_implicit , 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:

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. Pour cela il faut décider si on veut une erreur absolue ou relative et quelle norme utiliser. Ici on choisit l’erreur absolue maximale sur tous les nœuds.

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:

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 = {err_ep:g}')
plt.plot(tt,uu_er,'m:v',label=f'Approchée par Euler Implicite, erreur max = {err_er:g}')
plt.grid()
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');

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 : pour un intervalle fixé, \(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\).

from numpy.linalg import norm

Ntest = 6

H_ep_er = np.empty(Ntest)
err_ep = np.empty(Ntest)
err_er = np.empty(Ntest)


for k in range(1,Ntest+1):

    tt_tmp = np.linspace(t0,tfinal,k*N+1) # on raffine la grille
    yy = sol_exacte(tt_tmp)
    uu_ep_tmp = euler_progressif(phi,tt_tmp,y0)
    uu_er_tmp = euler_regressif(phi,tt_tmp,y0)    

    H_ep_er[k-1]= tt_tmp[1]-tt_tmp[0]
    err_ep[k-1] = norm(uu_ep_tmp-yy,np.inf) # idem que max(abs(uu_ep-yy))
    err_er[k-1] = 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})\).

plt.figure(figsize=(10, 7))
plt.loglog(H_ep_er,err_ep, 'r-+', label='Euler Explicite')
plt.loglog(H_ep_er,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).

# ln(e) = a ln(h) + b
a_ep, b_ep = np.polyfit(np.log(H_ep_er),np.log(err_ep), 1) # polyfit ( [liste des abscisses], [liste des ordonnées], degré du polynome)
print (f'Euler progressif: ordre {a_ep :1.2f}')
a_er, b_er = np.polyfit(np.log(H_ep_er),np.log(err_er), 1)
print (f'Euler regressif: ordre {a_er: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 :

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

plt.subplot(1,2,1)
plt.plot( np.log(H_ep_er), np.log(err_ep), 'ro',label='Euler Explicite')
plt.plot( np.log(H_ep_er), a_ep*np.log(H_ep_er)+b_ep ) # la droite de régression
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_ep_er), np.log(err_er), 'ro',label='Euler Implicite')
plt.plot( np.log(H_ep_er), a_er*np.log(H_ep_er)+b_er ) # la droite de régression
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 :

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

plt.subplot(1,2,1)
plt.loglog(H_ep_er,err_ep, 'ro',label='Euler Explicite')
plt.loglog( H_ep_er, np.exp(b_ep)*(H_ep_er**a_ep) ) # la droite de régression
plt.xlabel('$h$')
plt.ylabel('$e$')
plt.legend(loc='upper left')
plt.grid(True);

plt.subplot(1,2,2)
plt.loglog(H_ep_er,err_er, 'ro',label='Euler Implicite')
plt.loglog( H_ep_er, np.exp(b_er)*(H_ep_er**a_er) ) # la droite de régression
plt.xlabel('$h$')
plt.ylabel('$e$')
plt.legend(loc='upper left')
plt.grid(True);

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

def CN(phi,tt,y0):
    h  = tt[1]-tt[0]
    N = len(tt)-1 # len(tt) = N+1 car tt contient N+1 points
    uu = np.empty_like(tt) # array de même taille que tt
    uu[0] = y0
    for i in range(N): # range(N) = 0,1,...,N-1
        fc_implicit = lambda x: -x + uu[i] + h/2*( phi(tt[i],uu[i]) + phi(tt[i+1], x) )
        temp = fsolve( fc_implicit , 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 = {max(np.abs(uu_ep-yy)):g}')
plt.plot(tt,uu_er,'m:v',label=f'Approchée par Euler Implicite, erreur max = {max(np.abs(uu_er-yy)):g}')
plt.plot(tt,uu_cn,'go',label=f'Approchée par Crank-Nicolson, erreur max = {max(np.abs(uu_cn-yy)):g}')
plt.grid()
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');

H_cn = np.empty(Ntest)
err_cn = np.empty(Ntest)


for k in range(1,Ntest+1):

    tt = np.linspace(t0,tfinal,k*N+1) # on raffine la grille
    yy = sol_exacte(tt) # on calcule la solution exacte sur cette grille
    uu_cn = CN(phi,tt,y0) # on calcule la solution approchée par Crank-Nicolson
    
    H_cn[k-1] = tt[1]-tt[0] # on stocke le pas de temps
    err_cn[k-1] = np.linalg.norm(uu_cn-yy,np.inf) # on stocke l'erreur (norme infinie)
    
# 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 {a_cn :1.2f}') 


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

plt.plot(np.log(H_ep_er), np.log(err_ep), 'o',label=f'Euler Explicite : pente = {a_ep:1.2f}')
plt.plot( np.log(H_ep_er), a_ep*np.log(H_ep_er)+b_ep ) # la droite de régression

plt.plot(np.log(H_ep_er),np.log(err_er), 'd',label=f'Euler Implicite : pente = {a_er:1.2f}')
plt.plot(np.log(H_ep_er), a_er*np.log(H_ep_er)+b_er )

plt.plot(np.log(H_cn),np.log(err_cn), 'v',label=f'Crank Nicolson : pente = {a_cn:1.2f}')
plt.plot(np.log(H_cn), a_cn*np.log(H_cn)+b_cn )

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

Crank Nicolson: ordre 2.01