CPGE OUJDA SPE

Résolution approchée d'une équation différentielle

1. Equations différentielles du premier ordre

De nombreux problèmes conduisent à la résolution d'équations différentielles du premier ordre du type :

Ecrite aussi  :(1).

Remarque :

Une équation différentielle du premier ordre normalisée peut toujours s'écrire sous la forme (1).

Par exemple : s'écrira : avec .

Parfois la fonction fne dépend que de y ; n’empêche, la fonction fdevra toujours être déclarée avecles deux variables t et y (dans cet ordre : variable puis fonction).

Par exemple : s'écrira :

Ø Fonction odeint

La bibliothèque scipy.integrate contient la fonction odeint, qui résout numériquement des équations différentielles. On commence donc par la charger via :

fromscipy.integrateimportodeint

L'utilisation se fait sous la forme :

odeint(f, y0,t)

Exemple1 :

importmatplotlib.pyplotasplt

fromscipy.integrateimportodeint

import numpy asnp

f = lambdax,y: -y #équivalent à def f(x,y) : return -y

a=0

b=5

t=np.linspace(a,b,101)

y=odeint(f,1,t)

print(t,y)

plt.plot(t,y)

plt.grid()

plt.show()

Exemple 2 :

Le code ci-dessous permet de résoudre l’équation différentielle du premier ordre suivante :

fromscipyimport*

fromscipy.integrateimportodeint #Pour résoudre l’équation différentielle

importmatplotlib.pyplotasplt #Permet de tracer des graphes

defderiv (y , x ) :

return y*sin(y) + sin (x)

#Choix de l’intervalle d’intégration

x0 = 0

xmax = 20

npoints = 1000

x = linspace( x0 , xmax , npoints )

#Condition initiale

y0 = 2 # y(x0) = y0

solution = odeint( deriv , y0 , x)

y = solution [ : , 0 ]

plt.plot (x , y)

plt.show ()

I. Méthode d'Euler

Ø Lorsqu'il n'est pas possible d'obtenir une solution explicite sous forme de fonctions usuelles, on détermine une solution approchée (numérique). Il existe plusieurs méthodes permettant de réaliser ceci. L'une d'entre elles est la méthode d'Euler qui sera étudiée ici.

Ø On veut déterminer une solution approchée de l’équation différentielle :

Où : est une fonction continue sur l’intervalle.

Avec une condition initiale :

· On considère une subdivision régulière de de, on a donc :

Si h est suffisamment petit on a :

· Les approximations sont alors calculées de proche en proche par :

Voici une implémentation de la méthode Euler en python

import numpy asnp

defeuler(a,b,f,y0,n):

'''integration de y'(t)=f(t,y(t)) par la méthode d'Euler sur un intervalle [a,b] avecn points et la condition initiale y(a)=y0 (a,y0)'''

x=np.linspace(a,b,n+1)

y=np.zeros(n+1) # tableau de n+1 éléments

y[0]=y0

h=(b-a)/n

for k inrange(n):

y[k+1]=y[k]+h*f(x[k],y[k])

return(x,y)

Illustrons l'utilisation de la fonction Euler pour trouver une solution approchée de l'équation :

y'=-y intégrée sur [0,5] avec la condition initiale y(0)=1.

importmatplotlib.pyplotasplt

f = lambdax,y:-y #équivalent à def f(x,y) : return -y

(x,y)=euler(0,5,f,1,100)

print(x,y)

plt.plot(x,y)

plt.grid()

plt.show()

2. Résolution d’équations différentielles du deuxième ordre

Exemple :

On attache une bille au pendule. En l’absence de frottements, les oscillations du pendule vérifient l’équation différentielle : (2)

On prendra dans toute la suite ω₀ = 3 rad · s⁻¹.

Dans le cas des petits angles , on a , et l’on retrouve l’équation différentielle d’un oscillateur harmonique : (3)

L’équation 2 n’a pas de solutions analytiques, on se propose donc de l’intégrer avec SciPy pour comparer ses solutions à celles de l’équation 3.

On prendra comme conditions initiales (t = 0) = 0,1 rad et(t = 0) = 0.

fromscipyimport*

fromscipy.integrateimportodeint

omega0 = 3#rad/s

defderiv (TH, t ) :

# prend comme arguments :

# TH : une liste contenant theta et thetap

# t : un instant donné

# retourne la liste : [thetap ,-omega**2*theta]

theta= TH[ 0 ]

thetap= TH[ 1 ]

return [thetap,-omega0**2*theta]

#définition de l’intervalle de temps

t0 = 0

tmax = 15

npoints = 1000

t = linspace( t0 , tmax , npoints )

#Définition des conditions initiales

theta0 = .1 #rad

thetap0 = 0#rad/s

ci = [ theta0 , thetap0 ]

solution = odeint( deriv , ci , t )

#la résolution proprement dite

theta = solution [ : , 0 ]

thetap = solution [ : , 1 ]

plt.plot(t ,theta)

plt.show ()