Méthode d’Euler

Dans la vidéo ci-dessous : Présentation des méthodes des équations aux différences et de la méthode d’Euler avec beaucoup d’exemples d’illustration. Rappels succincts sur les méthodes, puis application à 4 cas concrets à partir de 46’30 (codes à faire apparaître et à exécuter dans Basthon plus bas dans cette page) :

  • Une équation du premier ordre: Vitesse limite en chute libre avec frottements
  • Système d’équations du premier ordre: Système Proie-Prédateur
  • Une équation du second ordre: Position au cours du temps d’un pendule
  • Système d’équations du second ordre: Trajectoire d’un objet en balistique

On voit comment programmer sois même la méthode « Euler explicite » pour comprendre ce qu’il y a derrière, puis on utilise le module préprogrammé de Python: odeint

PDF de la présentation de la vidéo :

proxy.php?id=pdf-info-simu-euler-presentation

Mon cours sur Euler :

proxy.php?id=pdf-info-simu-euler-cours

Mon résumé sur Euler :

proxy.php?id=pdf-info-simu-euler-resume
Application 1 de la vidéo : Faire apparaître sujet-corrigé et Basthon
Application 2 de la vidéo : Faire apparaître sujet-corrigé et Basthon
Application 3 de la vidéo : Faire apparaître sujet-corrigé et Basthon
Application 4 de la vidéo : Faire apparaître sujet-corrigé et Basthon
Codes non exécutables indépendants de Basthon et Dropbox (peuvent ne pas être à jour)

Code 1

'''
Résolution d'une équation différentielle du premier ordre
Vitesse limite en chute libre
'''

## Imports

import matplotlib.pyplot as plt
plt.close('all')

## Affichage

def Affiche(fig,lx,ly,leg):
    plt.figure(fig)
    plt.plot(lx,ly,label=leg)
    plt.legend()
    plt.show()
    plt.pause(0.01)

## Fonction spécifique au problème traité

def F1(Y,t):
    m = 0.1
    g = 9.81
    k = 0.001
    Yp = (k/m)*Y**2-g
    return Yp

## Euler explicite

# Euler_Explicite

def Euler_Explicite(F,t0,tf,y0,dt):
    Lt = [t0]
    LY = [y0]
    t = t0
    Y = y0
    i = 0
    while t <= tf:
        Yp = F(Y,t)
        Y = Y + Yp*dt # Yp nombre
        i += 1
        t = t0 + i*dt
        Lt.append(t)
        LY.append(Y)
    return Lt,LY

# Résolution avec CI

t0 = 0
tf = 10
Y0 = 0

Ldt = [0.1,0.01,0.001]
for dt in Ldt:
    Lt,LY = Euler_Explicite(F1,t0,tf,Y0,dt)
    Affiche(1,Lt,LY,dt)

## odeint

N = 1000
dt = (tf-t0)/(N-1)
Lt = [t0+i*dt for i in range(N)]
from scipy.integrate import odeint
LY = odeint(F1,Y0,Lt)
Affiche(11,Lt,LY,dt)

Code 2

'''
Résolution d'un système d'équations différentielles du premier ordre
Proie-Prédateur
'''

## Imports

import matplotlib.pyplot as plt
plt.close('all')

## Affichage

def Affiche(fig,lx,ly,leg):
    plt.figure(fig)
    plt.plot(lx,ly,label=leg)
    plt.legend()
    plt.show()
    plt.pause(0.01)

## Fonction spécifique au problème traité

def F2(Y,t):
    H,A = Y
    a = 3
    b = 1
    c = 1
    d = 2
    dH = a*H-b*A*H
    dA = -c*A+d*A*H
    return [dH,dA]

## Euler explicite

# Euler_Explicite

def Euler_Explicite(F,t0,tf,y0,dt):
    Lt = [t0]
    LY = [y0]
    t = t0
    Y = y0
    i = 0
    while t <= tf:
        Yp = F(Y,t)
        Y = [Y[i] + Yp[i]*dt for i in range(len(Y))]
        i += 1
        t = t0 + i*dt
        Lt.append(t)
        LY.append(Y)
    return Lt,LY

# Résolution avec CI

t0 = 0
tf = 30
H0,A0 = 10,1
Y0 = [H0,A0]

Ldt = [0.0001]
for dt in Ldt:
    Lt,LY = Euler_Explicite(F2,t0,tf,Y0,dt)
    LH = [t[0] for t in LY]
    LA = [t[1] for t in LY]
    Affiche(1,Lt,LH,'H')
    Affiche(1,Lt,LA,'A')
    Affiche(2,LA,LH,'H=f(A)')

## odeint

N = 1000
dt = (tf-t0)/(N-1)
Lt = [t0+i*dt for i in range(N)]
from scipy.integrate import odeint
LY = odeint(F2,Y0,Lt)
LH = [t[0] for t in LY]
LA = [t[1] for t in LY]
Affiche(11,LA,LH,'H=f(A)')

Code 3

'''
Résolution d'une équation différentielle du second ordre
Pendule
'''

## Imports

import matplotlib.pyplot as plt
plt.close('all')

## Affichage

def Affiche(fig,lx,ly,leg):
    plt.figure(fig)
    plt.plot(lx,ly,label=leg)
    plt.legend()
    plt.show()
    plt.pause(0.01)

## Fonction spécifique au problème traité

from math import sin
def F3(Y,t):
    teta,tetap = Y
    m = 0.1
    R = 0.1
    k = 0.001
    g = 9.81
    J = m*R**2
    tetapp = (-k*tetap-R*m*g*sin(teta))/J
    return [tetap,tetapp]

## Euler explicite

# Euler_Explicite

def Euler_Explicite(F,t0,tf,y0,dt):
    Lt = [t0]
    LY = [y0]
    t = t0
    Y = y0
    i = 0
    while t <= tf:
        Yp = F(Y,t)
        Y = [Y[i] + Yp[i]*dt for i in range(len(Y))]
        i += 1
        t = t0 + i*dt
        Lt.append(t)
        LY.append(Y)
    return Lt,LY

# Résolution avec CI

t0 = 0
tf = 20
from math import pi
teta = 45*pi/180
tetap = 0
Y0 = [teta,tetap]

Ldt = [0.001]
for dt in Ldt:
    Lt,LY = Euler_Explicite(F3,t0,tf,Y0,dt)
    Lteta = [t[0] for t in LY]
    Ltetap = [t[1] for t in LY]
    Affiche(1,Lt,Lteta,dt)

## odeint

N = 1000
dt = (tf-t0)/(N-1)
Lt = [t0+i*dt for i in range(N)]
from scipy.integrate import odeint
LY = odeint(F3,Y0,Lt)
Lteta = [t[0] for t in LY]
Affiche(11,Lt,Lteta,dt)

Code 4

'''
Résolution d'un système d'équationz différentiellez du second ordre
Balistique
'''

## Imports

import matplotlib.pyplot as plt
plt.close('all')

## Affichage

def Affiche(fig,lx,ly,leg):
    plt.figure(fig)
    plt.plot(lx,ly,label=leg)
    plt.legend()
    plt.axis('scaled')
    plt.show()
    plt.pause(0.01)

## Fonction spécifique au problème traité

from math import sqrt
def F4(Y,t):
    x,xp,y,yp = Y
    beta = 1
    m = 1
    g = 9.81
    racine = sqrt(xp**2+yp**2)
    xpp =   -(beta/m)*xp*racine
    ypp = -g-(beta/m)*yp*racine
    return [xp,xpp,yp,ypp]

## Euler explicite

# Euler_Explicite

def Euler_Explicite(F,t0,tf,y0,dt):
    Lt = [t0]
    LY = [y0]
    t = t0
    Y = y0
    i = 0
    while t <= tf:
        Yp = F(Y,t)
        Y = [Y[i] + Yp[i]*dt for i in range(len(Y))]
        i += 1
        t = t0 + i*dt
        Lt.append(t)
        LY.append(Y)
    return Lt,LY

# Résolution avec CI

t0 = 0
tf = 0.84
x0 = 0
y0 = 0
xp0 = 10
yp0 = 10
Y0 = [x0,xp0,y0,yp0]

Ldt = [0.0001]
for dt in Ldt:
    Lt,LY = Euler_Explicite(F4,t0,tf,Y0,dt)
    Lx = [t[0] for t in LY]
    Ly = [t[2] for t in LY]
    Affiche(1,Lx,Ly,'y=f(x)')

## odeint

N = 1000
dt = (tf-t0)/(N-1)
Lt = [t0+i*dt for i in range(N)]
from scipy.integrate import odeint
LY = odeint(F4,Y0,Lt)
Lx = [t[0] for t in LY]
Ly = [t[2] for t in LY]
Affiche(11,Lx,Ly,'y=f(x)')

Les Sciences Industrielles de l'Ingénieur (SII) et l'Informatique du Tronc Commun (ITC) en CPGE par Denis DEFAUCHY