CPGE-AGADIR MP-PSI

RESOLUTION NUMERIQUE DES EQUATIONS DIFFERENTIELLES

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).
𝒅𝒕
𝒚(𝒕𝟎 ) = 𝒚𝟎

a. 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 [a,b].

avec une condition initiale : y(a)=y0
 On considère une subdivision régulière t = (t0, t1, …, tn) de [a,b] de pas h =(b-a)/n, on a donc :

ti =t0+ih =a+ih

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 as np
def euler(f,a,b,y0,n):
'''integration de y'(t)=f(t,y(t)) par la méthode d'Euler sur un intervalle [a,b] avec
n subdivisions et la condition initiale y(a)=y0 (a,y0)'''
h=(b-a)/n
t=a
y=y0
T=[a]
Y=[y0]
for k in range(n):
y=y+h*f(t,y)
[Link](y)
t=t+h
[Link](t)
return(T,Y)
Résolution d’équations différentielles -1- [Link]
Exemple 1 :
Résoudre l’équation : y′=2.x.y2 sur[0, 2] (on prendra une subdivision de 100 points) avec la condition
initiale y(0)=0.2 :

import [Link] as plt

f=lambda x,y: 2*x*y**2 #équivalent à def f(y,x) : return 2*x*y**2
(X,Y)= euler(f,0,2, 0.2,100)
[Link](X,Y)
[Link]()
[Link]()

Exemple 2 :
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.

import [Link] as plt

f = lambda x,y: -y #équivalent à def f(y,x) : return -y
(X,Y)= euler(f,0,5, 1,100)
[Link](X,Y)
[Link]()
[Link]()

b. Fonction odeint

 La bibliothèque [Link] contient la fonction odeint, qui résout numériquement des

équations différentielles. On commence donc par la charger via :

from [Link] import odeint

 L'utilisation se fait sous la forme :

odeint(f, y0,t)

#utilisation fonction odeint

Exemple 1 :
Résoudre l’équation : y′=2.x.y2 sur[0, 2] (on prendra une subdivision de 100 points) avec la condition
initiale y(0)=0.2 en utilisant la fonction odeint:

from [Link] import odeint

X = [Link](0,2,100)
f=lambda y,x: 2*x*y**2
Y=odeint(f,.2,x)
[Link](X,Y)
[Link]()
[Link]()

Résolution d’équations différentielles -2- [Link]

Exemple 2 :
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 en utilisant la fonction odeint:
from [Link] import odeint
f=lambda y,x: -y
a=0
b=5
n=100
T=[Link](a,b,n+1)
Y=odeint(f,1,T)
[Link](T,Y)
[Link]()
[Link]()

Exercice 1 :
A l'aide de la fonction odeint( ) du sous-module integrate de scipy, et du module [Link], donner le
code permettant de tracer la courbe représentative sur l'intervalle [0,4] (on prendra une subdivision de 100 points)
𝒚′ + √𝒕 . 𝒚 = 𝟎
de la solution approchée de l’équation suivante : {
𝒚(𝟎) = 𝟏
Correction
import [Link] as plt
from scipy import integrate
phi = lambda y,t: -y * t**0.5
T = [Link](0,4,100)
y0 = 1
Y = [Link](phi,y0,T
[Link](T,Y)
[Link]()

Exercice 2 : Evolution d’une température en chimie

On note y(t) la température en degrés Celsius d’une réaction chimique en fonction du temps t, t étant
exprimé en heures.
Après étude, on constate que la température est solution de l’équation différentielle :

y '  y  e 0, 25t

Avec la condition initiale y(0)=20°C.

- Ecrire un programme qui donne un tableau composé des températures de la réaction chimique chaque
demi-heure pendant la première journée de la réaction.
- Tracer la courbe représentative de cette solution
import numpy as np
from [Link] import odeint
import [Link] as plt

Résolution d’équations différentielles -3- [Link]

 T=[Link](0,24.5,0.5)
phi = lambda y,t:([Link](-0.25*t))-y
Y = odeint(phi,20,t)
[Link](T,Y)
[Link]()

2. Equations différentielles d’ordre 2

Si on cherche à intégrer l’équation différentielle d’ordre 2 : y′′= g(y, y′ , x) on la transforme en une équation

différentielle d’ordre 1 en posant z =(y, y′) donc z′ =(y′,y′′), l’équation différentielle peut alors s’écrire :

z′= f(z , x)
en posant f( (y, dy) , x) = (dy , g(y, dy ,x)).

La subdivision[x0, x1,..., xn] de l’intervalle sur lequel on intègre l’équation différentielle reste la même, par contre la

condition initiale doit être un couple (y(x0), y′(x0)) et la valeur de retour est une matrice numpy dont la première

colonne est constituée des valeurs y(xk) et la seconde des valeurs y′(xk).

Exemple 1: résolution de : y′′+ y′+ xy =0 sur le segment [1, 2] avec la condition initiale (y(0)=1, y′(0)=0).

a. Fonction odeint

from [Link] import odeint

import [Link] as plt
import numpy as np

def f(z,x):
(y,dy)=z
return [Link]([dy,-dy-x*y]) #retourne un tableau

>>>X=[Link](1.,2.,4)
>>>S=odeint(f,(1,0),X)
>>>print(S)
[[ 1. 0. ]
[ 0.94502133 -0.32624253]
[ 0.78778807 -0.60463152]
[ 0.55333316 -0.78078941]]

>>>print(S[:,0]) #première colonne: valeurs prises par S

[ 1. 0.94502133 0.78778807 0.55333316]

>>>[Link](X,S[:,0])
>>>[Link]() ; [Link]()

Résolution d’équations différentielles -4- [Link]

 b. Méthode d’Euler
La fonction euler utilisée pour résolution de l’équation différentielle d’ordre un reste la même
pour résoudre l’équation d’ordre 2 à condition que la fonction f(y,t) retourne un tableau array.

def euler(f,a,b,y0,n):
h=(b-a)/n
t=a
y=y0
T=[a]
Y=[y0]
for k in range(n):
y=y+h*f(t,y)
[Link](y)
t=t+h
[Link](t)
return(T,Y)

def f(x,z):
(y,dy)=z
[Link]([dy,-dy-x*y])

(X,Y)= euler(f,1,2,(1,0),10)
[Link](X,Y)
[Link]() ; [Link]()

Exemple 2 : résolution de l’équation y''=-y dans l’intervalle [0,2] avec y(0)=0 et y’(0)=1

Résolution d’équations différentielles -5- [Link]

 def f(x,z):
(y,dy)=z
[Link]([dy,-y])

(X,Y)=euler(f,0,2*[Link],(0,1),100)
[Link](X,Y)
[Link](X,[Link](X))
[Link]('x') #Etiquette sur l'axe x
[Link]('y(x)') #Etiquette sur l'axe y

[Link](('y','sin(x)'),'upper right', shadow=True)

[Link]('Résoluton de y\'\'=-y avec y(0)=0 et y\'(0)=1')
[Link]()
[Link]()

Exercice 3 :
Résoudre l’équation : 𝒚̈ + 𝟐 𝒚̇ + 𝝎𝟐𝟎 𝒚 = 𝑨 𝒄𝒐𝒔(𝝎𝒕) dans l’intervalle [0,4] avec une subdivision de 100
points ; avec : y(0)=1 ; y'(0)=0 ; A=1.5 ; 𝝎 = 𝟐 ;  = 𝟎. 𝟎𝟏 ; 𝝎𝟎 = 𝟏
Puis tracer la solution :

Résolution d’équations différentielles -6- [Link]

 #importations
from [Link] import *
from numpy import *
import [Link] as plt

#fonction f(z,t)
def f(z,t) :
y,dy=z
return [Link]([dy, -2*lamb*dy-w0**2*y+A*cos(w*t)])

#Appel
lamb,w0,A,w =0.01,1.0,1.5,2.
z=[1.,0.]
T=linspace(0,4*pi,100)
S=odeint(f,z,t)
[Link](t,S[:,0])
[Link]()
[Link]()

Exercice 4 : Le parachute

La trajectoire suivie par un objet relié à un parachute est un axe vertical noté (O, i ).
  
A un instant donné, le vecteur vitesse V de l’objet est défini par V(t)   (t ) i où v
est une fonction de la variable réelle positive t (temps).

Dans ces conditions de l’expérience, le vecteur R représentant la résistance de
 
l’air est défini par R  k V où k est un nombre réel strictement positif. (Voir
Figure ci-après)
On admet que la position x(t) de l’objet et du parachute vérifie l’équation
différentielle suivante :
𝒎𝒙̈ (𝒕) + 𝒌𝒙̇ (𝒕) = 𝒎𝒈
où m est la masse totale de l’objet et du parachute et g le coefficient de l’accélération de la pesanteur.
A partir de l’équation différentielle donner un programme de résolution numérique qui nous permet
de donner comme résultat la distance x de la chute (l’objet + parachute) après 2 minutes avec un pas de 2
seconds.
On prendra :
m=8 kg, g=10 m s-2, k=25 unités SI, position initiale x(0)=0m et vitesse initiale v(0)=v0=3,2 m s-1.

Résolution d’équations différentielles -7- [Link]

 #importations
from numpy import arange
from [Link] import odeint
import [Link] as plt
#fonction
def phi(z,t):
y, dy= z
return([dy , g-((k*dy)/m)])

#Résolution

t=arange(0.,121.,2)
x0=0.
v0=3.2
m=8
g=10
k=25
s=odeint(phi,(x0,v0),t)
x=s[:,0]
print(x)
print(x[-1])
[Link](t,x)
[Link]()
[Link]()

Résolution d’équations différentielles -8- [Link]

Exercice 5 : Etude d'un circuit RC
L'équation différentielle qui régit l'évolution de la tension aux bornes du condensateur d'un circuit RC en
régime libre s'écrit sous la forme:
𝒅𝑼𝒄 (𝒕) 𝟏
+ 𝑼 ( 𝒕) = 𝟎
𝒅𝒕 𝑹𝑪 𝒄
1. Ecrire une fonction qui implémente la méthode d'Euler-Cauchy d'intégration d'une équation-différentielle
du premier ordre.
2. On suppose qu'initialement Uc (t = 0) = 1V ; R = 1kΩ; C = 1nF. Réaliser l'intégration de l'équation-
différentielle précédente par la méthode d'Euler-Cauchy entre t=0 et 0.1ms, avec une subdivision de 1000
segments.
3. Donner le code permettant d'effectuer le tracé de la courbe représentative de Uc (t) pour t variant de 0 à
0.1ms.
Exercice 6 : Circuit RLC en série.
On considère un circuit électrique constitué d'une source de tension V, d'une
résistance R, d'une bobine L et d'une capacité C montés en série.
𝒅𝟐 𝒖𝑳 𝑹 𝒅𝒖𝑳 𝟏 𝒅𝟐 𝑽
UL satisfait l’équation différentielle: + + 𝒖𝑳 =
𝒅𝒕𝟐 𝑳 𝒅𝒕 𝑳𝑪 𝒅𝒕𝟐

C’est une équation différentielle du second ordre à coefficient constant avec second
membre :
𝑹 𝟏 𝟐 𝟐
𝒅𝟐 𝑽
𝒚̈ 𝒚̇ + 𝒚 = −𝟒 𝒇 𝒔𝒊𝒏(𝟐𝐟𝐭) = 𝟐
𝑳 𝑳𝑪 𝒅𝒕

Question : Tracer la courbe de la tension aux bornes de la bobine pour V(t) = sin(2ft) sur un intervalle de
temps de 1 seconde (on prendra 10000 points) pour f = 100Hz, R = 10, L = 870mH, C = 650mF. A l'instant
initial aux bornes de la bobine la tension est nulle ainsi que sa dérivée.

Exercice 7 : Pendule simple

Pour le pendule simple (non amorti) l’angle avec la verticale suit une équation différentielle de la forme :

𝜃̈ = − sin(𝜃) [1]
Q 1 : Résoudre numériquement et représenter la solution vérifiant les conditions initiales :
𝜽(𝟎) = 0.2 et 𝜽̇(𝟎) = 𝟎. 𝟐
Q 2 : Dans le cas des petits angles, on a sin(𝜃) ≈ 𝜃 , et l’on retrouve l’équation différentielle d’un oscillateur
harmonique :

𝜃̈ = − 𝜃 [2]
Résoudre numériquement et représenter la solution vérifiant les conditions initiales 𝜽(𝟎) = 𝟎. 𝟐 et 𝜽̇(𝟎) = 𝟎. 𝟐.
Y a-t-il une différence importante entre les solutions de [1] et [2] ?
Q 3 : En prenant toujours 𝜃 (0) = 0, représenter les solutions de [1] vérifiant les conditions initiales :
𝜽̇(𝟎) = {𝟏, 𝟏. 𝟗, 𝟐, 𝟐. 𝟏, 𝟑}.
Q 4 : Avec les différentes conditions initiales vues plus haut, représenter le portrait de phase c’est-à-dire l’ensemble
des points des couples (𝜽, 𝜽̇) .

Q 5 : Reprendre les questions précédentes avec un pendule amorti : 𝜃̈ = − sin(𝜃 ) − 𝑘𝜃

On pourra prendre k = 1 dans un premier temps, et éventuellement jouer sur la valeur de ce paramètre (qui
modélise l’amortissement).

Résolution d’équations différentielles -9- [Link]

 Système des équations différentielles

Exercice1 : Mise sous tension d’un moteur à courant continu :

On alimente un moteur initialement à l’arrêt par une tension u(t) dont la

représentation graphique est donnée ci-contre. On note w(t) la vitesse de
rotation du moteur et i(t) l’intensité du courant. On a les équations
différentielles:
𝑑𝑖
𝑢 (𝑡) = 𝑅𝑖 (𝑡) + 𝐿 + 𝑘𝜔 (𝑡)
𝑑𝑡
𝑑𝜔
𝑗 = 𝐾𝑖(𝑡) − 𝑓𝜔(𝑡)
𝑑𝑡
Où:
• J est le moment d’inertie J = 3.2.10-4 kg.m2 ;
• f est un coefficient de frottement f = 2.4 .10-3 [Link]-1 ;
• R est la résistance du bobinage R= 0.67 ;
• L est l’inductance du bobinage L=3.8 .10-3 H;
• K est la constante de force électromotrice du couple K=0.12N/A.

a. Définir la fonction u(t).

i(t) et w(t) sont solutions d’un système d’équations différentielles:
𝑑𝑖
=⋯
𝑑𝑡

𝑑𝜔
{ 𝑑𝑡 = ⋯
i
b. On note Y = ( ). Définir la fonction F(𝒀, 𝒕) associée au système.
ω
c. Résoudre le système avec la condition initiale 𝜔(0)= 0 et i (0)= 0 sur l’intervalle [0,0.4] et représenter les fonctions
obtenues.
d. Écrire la fonction associée au système sous la forme F(Y,t,K) afin que la force électromotrice K soit prise comme
un paramètre. On place une charge sur le moteur ce qui a pour effet d’augmenter son inertie à K1 =0.2N/A.
Résoudre alors le système pour cette valeur et représenter sur le même graphique les deux vitesses de rotation.
Quelle est l’influence de K sur la vitesse de rotation?

Exercice2 : Propagation d’une épidémie /Système d’équations différentielles :

On propose le modèle simple suivant représentant la propagation d’une épidémie (appelé modèle SIR):
— La population est divisée en trois groupes : sains(S), infectés(I) et résistants(R), dont les effectifs au temps t sont
notés respectivement s(t), i (t) et r (t);
— Les fonctions s, i et r sont solution d’un système d’équations différentielles du type:
𝑠 ′ (𝑡) = −𝑏𝑠(𝑡)𝑖 (𝑡)
{ 𝑖 ′ (𝑡) = 𝑏𝑠(𝑡)𝑖 (𝑡) − 𝑐𝑖 (𝑡)
𝑟 ′ (𝑡) = 𝑐𝑖 (𝑡)

Le coefficient b est un taux de contagion entre un individu sain et un individu infecté et le coefficient c est un taux de
guérison d’un individu infecté;
— On prendra comme conditions initiales s(0) =0.99, i (0) = 0.01 et r (0) = 0.
Représenter graphiquement les fonctions s et i pour t 𝝐 [0,50] et (b,c) = (0.9,0.1).
Faire de même lorsque (b,c) = (0.6,0.4).

Résolution d’équations différentielles -10- [Link]