import numpy as np
import math

#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

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

#schéma du point-milieu
def runge(f, t0, y0, h, tfin):
    t=np.arange(t0,tfin+h,h) # python n'atteint pas la borne superieure :-()
    y = np.zeros(len(t))
    y[0] = y0
    for n in range(len(t)-1):
        u = y[n] + (h/2)*f(t[n], y[n])
        y[n+1]= y[n]+ h*( f(t[n]+h/2,u) )
    return y

#schéma de Heun
def heun(f, t0, y0, h, tfin):
    t=np.arange(t0,tfin+h,h) # python n'atteint pas la borne superieure :-()
    y = np.zeros(len(t))
    y[0] = y0
    for n in range(len(t)-1):
        u2 = y[n] + (h/3)*f(t[n], y[n])
        u3 = y[n] + (2*h/3)*f(t[n]+h/3, u2)
        y[n+1]= y[n]+ h*( (1/4)*f(t[n],y[n])+ (3/4)*f(t[n]+(2*h)/3,u3) )
    return y

#schéma de Runge-Kutta 4
def RK4(f, t0, y0, h, tfin):
    t=np.arange(t0,tfin+h,h) # python n'atteint pas la borne superieure :-()
    y = np.zeros(len(t))
    y[0] = y0
    for n in range(len(t)-1):
        u2 = y[n] + (h/2)*f(t[n], y[n])
        u3 = y[n] + (h/2)*f(t[n]+h/2, u2)
        u4 = y[n] + h*f(t[n]+h/2, u3)
        y[n+1]= y[n]+ h*((1/6)*f(t[n],y[n])+ (2/6)*f(t[n]+h/2,u2) + (2/6)*f(t[n]+h/2,u3) + (1/6)*f(t[n]+h,u4) )
    return y

t0, tfin = 0, 7
y0,h=1, 1
ye=euler(f, t0, y0, h, tfin)
ym=runge(f, t0, y0, h, tfin)
yh=heun(f, t0, y0, h, tfin)
yrk=RK4(f, t0, y0, h, tfin)
t=np.arange(t0,tfin+h,h)
tt=np.arange(t0,tfin+0.01,0.01)
## sol exacte
z=np.exp(tt)

import matplotlib.pyplot as plt
fig1 = plt.figure(figsize=(10, 8))
plt.plot(t,ye,'-or')
plt.plot(t,ym,'-og')
plt.plot(t,yh,'-ob')
plt.plot(t,yrk,'-om')
plt.plot(tt,z,'-k')
plt.legend(['Euler','Point milieu','Heun','RK4','solution exacte'])
plt.show()

pas = [1/4, 1/8, 1/16,  1/32, 1/64, 1/128]
err = np.zeros((len(pas), 4))
for i, h in enumerate(pas):
    yrk = RK4(f, t0, y0, h, tfin)
    yh = heun(f, t0, y0, h, tfin)
    yr = runge(f, t0, y0, h, tfin)
    ye = euler(f, t0, y0, h, tfin)
    t = np.arange(t0, tfin+h, h)
    z = np.exp(t)
    erreurrk = np.max(np.abs(yrk-z))
    erreurh = np.max(np.abs(yh-z))
    erreurr = np.max(np.abs(yr-z))
    erreure = np.max(np.abs(ye-z))
    err[i, :] = np.array([erreure, erreurr, erreurh, erreurrk])
# calcul des pentes avec les droites de régression
pe = np.polyfit(np.log(pas), np.log(err[:, 0]), 1)
pm = np.polyfit(np.log(pas), np.log(err[:, 1]), 1)
ph = np.polyfit(np.log(pas), np.log(err[:, 2]), 1)
prk = np.polyfit(np.log(pas), np.log(err[:, 3]), 1)
plt.loglog(pas, err[:, 0], '-or')
plt.loglog(pas, err[:, 1], '-og')
plt.loglog(pas, err[:, 2], '-ob')
plt.loglog(pas, err[:, 3], '-om')
#plt.loglog(pas, pas, '--or')
plt.loglog(pas, np.power(pas, 2), '--og')
plt.loglog(pas, np.power(pas, 3), '--ob')
plt.loglog(pas, np.power(pas, 4), '--om')
plt.legend([r'Euler pente=%.1f' % (pe[0], ), r'point milieu pente=%.1f' % (pm[0], ), r'Heun pente=%.1f' % (
    ph[0], ), r'RK4 pente=%.1f' % (prk[0], ),  'pente 2', '3', '4'], loc='best')
plt.xlabel('log pas')
plt.ylabel('log erreur')
plt.title('convergence comparée RK4, Heun, Point milieu, Euler')
plt.show()

Les principaux schémas à un pas¶

Tracé des solutions approchées¶

Tracé de l'erreur globale en log log¶