Le schéma d'Euler

$$y'(t)= \frac{dy}{dt} \approx \frac{y(t+h) - y(t)}{h}$$

$$y(t+h) \approx y(t) + h\,f(t,y(t))$$

On génère une suite de valeurs $y_n$ aux instants $t_n = t_0 + n\cdot h$

le schéma s'écrit $$y_{n+1} = y_n + h\cdot f(t_n,y_n)$$ où $y_0$ est la valeur initiale $y(t_0)$

import numpy as np
#schéma d'Euler
def euler(f, t0, y0, h, tfin):
    t=np.arange(t0,tfin+h,h) # la subdivision, python n'atteint pas la borne superieure :-(
    y = np.zeros(len(t)) # initialisation
    y[0] = y0
    for n in range(len(t)-1):
        y[n+1] = y[n] + h*f(t[n], y[n])
    return y

import math
# definition de la fonction dans l'EDO
def f(t,y):
    return y

import matplotlib.pyplot as plt
# tinitial, donnee initiale, tfinal
t0,y0=0.0,1.0 
def quadeuler(  h=1.0, tfin = 4.0):
# appel de la fonction
    y=euler(f, t0, y0, h, tfin)
    t=np.arange(t0,tfin+h,h)
    plt.plot(t,y,'-ob')
    z=y0*np.exp(t)
    plt.plot(t,z,'or')
    tt=np.linspace(t[0],t[-1],100)
    #pour avoir exactement le meme intervalle que y
    plt.plot(tt,y0*np.exp(tt),'-r')
    plt.legend(['schéma euler','valeurs exactes','solution exacte'])
quadeuler()   

On peut faire varier le pas de la subdivision

from ipywidgets import interact
import ipywidgets as widgets

interact(quadeuler,  h=(1/64, 1.0), tfin=(1.0,7.0));

Widget Javascript not detected.  It may not be installed or enabled properly. Reconnecting the current kernel may help.

Convergence quand le pas de temps $h\to 0$

On peut essayer avec des pas de plus en plus petits.

tfin=4.0
pas = [1, 1/2, 1/4,  1/8, 1/16, 1/32,1/64]
for  h in pas:
    y=euler(f, t0, y0, h, tfin)
    t=np.arange(t0,tfin+h,h)
    plt.plot(t,y,'-ob',markersize=2)
    tt=np.arange(t0,tfin+1e-1,1e-1)
    plt.plot(tt,y0*np.exp(tt),'-r')
    plt.legend(['schéma euler','solution exacte'])

Calcul des erreurs globales quand le pas de temps $h\to 0$

tfin=2.0
pas = [1, 1/2, 1/4,  1/8, 1/16, 1/32, 1/64]
err=np.zeros(len(pas))
for i, h in enumerate(pas):
    y=euler(f,t0,y0,h,tfin)
    t=np.arange(t0,tfin+h,h)
    z=np.exp(t)
    erreur= np.max(np.abs(y-z))
    err[i]=erreur
for i, h in enumerate(pas):
    print('pas', pas[i],'erreur',round(err[i],4))   

pas 1 erreur 3.3891
pas 0.5 erreur 2.3266
pas 0.25 erreur 1.4286
pas 0.125 erreur 0.8058
pas 0.0625 erreur 0.4304
pas 0.03125 erreur 0.2228
pas 0.015625 erreur 0.1134

On trace en log log

plt.loglog(pas, err,'-or')
plt.loglog(pas,pas,'-ob')
plt.legend(['erreur','droite de pente 1'])
plt.xlabel('log pas')
plt.ylabel('log erreur')
plt.title('convergence Euler')
plt.show()

Le fait que la pente vaut approximativement 1 traduit que le schéma d'Euler est d'ordre 1:

si $err(h) = C h + o(h)$ alors $\ln err(h) = \ln h + \ln C + o(1)$

En coordonnée log-log, on a bien une droite de pente 1.