Chapitre 4 Résolution numérique des équations différentielles ordinaire

Chapitre 4

Résolution Numérique des équations différentielles ordinaire

4.1 Introduction
Soit le problème de Cauchy
𝑑𝑦
𝑦′ = = 𝑓(𝑡, 𝑦)
(1) { 𝑑𝑡
𝑦(𝑡0 ) = 𝑦0
Où 𝑓 : [𝑡0 , 𝑡𝑁 ] × 𝐼𝑅 → 𝐼𝑅 . Le problème de Cauchy doit satisfaire les conditions qui assurent
l’existence et l’unicité de la solution (théorique) de ce problème :
L’objectif de ce chapitre est décrire un certain nombre de méthodes permettant de résoudre
numériquement le problème de Cauchy :
Etant donne une subdivision 𝑡0 < 𝑡1 < ⋯ < 𝑡𝑁 = 𝑡0 + 𝑇 de [𝑡0 , 𝑡0 + 𝑇], on cherche à
déterminer des valeurs approchées 𝑦0 , 𝑦1 , … , 𝑦𝑁 des valeurs 𝑦(𝑡𝑛 ) prises par la solution exacte
𝑦. On notera les pas successifs.
ℎ𝑛 = 𝑡𝑛+1 − 𝑡𝑛 , 0≤𝑛 ≤𝑁−1
Et ℎ𝑚𝑎𝑥 = 𝑚𝑎𝑥(ℎ𝑛 )le maximum des pas
On appelle méthode à un pas une méthode permettant de calculer 𝑦𝑛+1 à partir de la seule
valeur antérieure 𝑦𝑛 . Une méthode à 𝑟 pas est au contraire une méthode où le calcul de 𝑦𝑛+1
nécessite la mémoration des valeurs𝑦𝑛 , 𝑦𝑛−1 , … , 𝑦𝑛−𝑟+1.
Ce chapitre est limité aux méthodes numériques à un pas.
Les méthodes à un pas sont les méthodes de résolution numérique qui peuvent s’écrire sous la
forme : 𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑛 𝑓(𝑡𝑛 , 𝑦𝑛 , ℎ𝑛 ). 0≤𝑛≤𝑁
Où 𝑓 ∶ [𝑡0 , 𝑡0 + 𝑇] × 𝐼𝑅 × 𝐼𝑅𝐼𝑅 est une fonction que l’on supposera continue.
Dans toutes les méthodes numériques développées par la suite, on subdivise l’intervalle
(𝑡0 +𝑇)−𝑡0 𝑇
(𝑡0 , 𝑡0 + 𝑇) en 𝑁 intervalles de longueur ℎ = =𝑁, limités par les points
𝑁
𝑡𝑛 = 𝑡0 + 𝑛ℎ, 0 ≤ 𝑛 ≤ 𝑁.

4.2 Équations différentielles ordinaire d’ordre 1

4.2.1 Méthode d’Euler

En 𝑡0 on connait 𝑦0 , donc aussi 𝑦′(𝑡0 ) = 𝑓(𝑡0 , 𝑦0 ).Si 𝑦(𝑡) est la solution exacte de (1).

𝑦(𝑡)est approchée sur l’intervalle [𝑡0 , 𝑡1 ] par sa tangente en 𝑡0 .Et ainsi, on a :

𝑦1 = 𝑦0 + ℎ𝑓(𝑡0 , 𝑦0 ).

Sur l’intervalle[𝑡1 , 𝑦1 ].

Ceci conduit à l’Algorithme d’Euler :

𝑦 = 𝑦𝑖 + ℎ𝑓(𝑡𝑖 , 𝑦𝑖 ), 0≤𝑖 ≤𝑁−1
(E) { 𝑖+1
𝑡𝑖+1 = 𝑡𝑖 + ℎ
Précision de la méthode d’Euler

Dr. DOUAK MOHAMED Page 1

Chapitre 4 Résolution numérique des équations différentielles ordinaire

La méthode d’Euler est une méthode du premier ordre, c’est –à-dire que l’erreur au point 𝑡𝑛
s’exprime par l’inégalité
|𝑦𝑖 − 𝑦(𝑡𝑖 )| ≤ 𝑘ℎ
Où 𝑦𝑛 est la valeur approchée définie par l’algorithme d’Euler, 𝑦(𝑡0 ) est la valeur exacte de la
solution du problème de Cauchy au point 𝑡 = 𝑡𝑖 + 𝑖ℎ, et 𝑘 une constante indépendante de i et
de h.
Remarque : l’erreur dans la méthode d’Euler est relativement importante. Elle peut être
améliorée en choisissant le pas h le plus petit, ce qui augmente considérablement le volume
des calculs à effectuer, ou en approchant la solution du problème de Cauchy par des méthodes
permettant de réduire cette erreur.
4.2.2 Méthode d’Euler modifié

On a d’après la formule de Taylor :

ℎ2
𝑦(𝑡𝑖+1 ) = 𝑦(𝑡𝑖 + ℎ) = 𝑦(𝑡𝑖 ) + ℎ𝑦 ′ (𝑡𝑖 ) + 𝑦 ′′ (𝑡𝑖 ) + 𝑜(ℎ3 ) (1)
2!

𝑑𝑦 ′ 𝑦′(𝑡𝑖+1 )−𝑦′(𝑡𝑖 )
Or 𝑦 ′′ (𝑡𝑖 ) = = + 𝑜(ℎ)
𝑑𝑡 ℎ

Portons l’expression de 𝑦 ′′ (𝑡𝑖 ) dans l’équation (1) nous obtenons :

ℎ2 𝑦′(𝑡𝑖+1 ) − 𝑦′(𝑡𝑖 )
𝑦(𝑡𝑖+1 ) = 𝑦(𝑡𝑖 ) + ℎ𝑦 ′ (𝑡𝑖 ) + ( + 𝑜(ℎ)) + 𝑜(ℎ3 )
2! ℎ

ℎ
𝑦𝑖+1 = 𝑦𝑖 + (𝑦 ′ (𝑡𝑖+1 ) + 𝑦′(𝑡𝑖 ))
2
ℎ
𝑦(𝑡𝑖+1 ) = 𝑦(𝑡𝑖 ) + 2 (𝑓(𝑡𝑖 , 𝑦𝑖 ) + 𝑓(𝑡𝑖+1 , 𝑦𝑖+1 )) + 𝑜(ℎ3 ) (2)

Remarquons que cette formule ne donne pas yi+1 explicitement. Pour surmonter ce problème
on estime (prédit) par la méthode d’Euler. Appelons cette valeur estimée : yi+1,p. Cette valeur
peut être utilisée dans l’équation (2) pour obtenir une estimation meilleure de yi+1.
La méthode d’Euler modifié consiste donc à :
1/ prédire yi+1 par la méthode d’Euler : 𝑦𝑖+1,𝑝 = 𝑦𝑖 + ℎ𝑓(𝑡𝑖 , 𝑦𝑖 )
ℎ
2/ corriger yi+1,p par : 𝑦𝑖+1,𝑐 = 𝑦𝑖 + 2 (𝑓(𝑡𝑖 , 𝑦𝑖 ) + 𝑓(𝑡𝑖+1 , 𝑦𝑖+1,𝑝 ))

Ceci conduit à l’Algorithme d’Euler modifié :

ℎ
𝑦𝑖+1,𝑐 = 𝑦𝑖 + 2 (𝑓(𝑡𝑖 , 𝑦𝑖 ) + 𝑓(𝑡𝑖+1 , 𝑦𝑖+1,𝑝 )) , 0≤𝑖 ≤𝑁−1
(Em) { 𝑦𝑖+1,𝑝 = 𝑦𝑖 + ℎ𝑓(𝑡𝑖 , 𝑦𝑖 )
𝑡𝑖+1 = 𝑡𝑖 + ℎ
Dr. DOUAK MOHAMED Page 2
Chapitre 4 Résolution numérique des équations différentielles ordinaire

4.2.3 Méthode de Taylor (d’ordre 2)

Supposons que 𝑓soit de classe 𝐶 1 . Alors 𝑦(𝑡) est de class 𝐶 2 , et le développement de Taylor
d’ordre 2 implique :
ℎ2 ′′
𝑦(𝑡𝑖+1 ) = 𝑦(𝑡𝑖 + ℎ) = 𝑦(𝑡𝑖 ) + ℎ𝑦 ′ (𝑡𝑖 ) + 𝑦 (𝑡𝑖 ) + 𝑜(ℎ3 )
2!
Et comme :
D’une part, 𝑦 ′ (𝑡𝑖 ) = 𝑓(𝑡𝑖 , 𝑦(𝑡𝑖 )),
𝑑 𝑑 𝜕𝑓 𝜕𝑓 𝑑𝑦
D’autre part, 𝑦 ′′ (𝑡𝑖 ) = 𝑑𝑡 (𝑦′(𝑡))/𝑡=𝑡𝑖 = 𝑑𝑡 (𝑓(𝑡, 𝑦(𝑡))) = ( 𝜕𝑡 + 𝜕𝑦 . 𝑑𝑡 ) = 𝑓′(𝑡, 𝑦(𝑡))
𝑡=𝑡𝑖 𝑡=𝑡𝑖

ℎ
𝑦(𝑡𝑖+1 ) = 𝑦(𝑡𝑖 ) + ℎ [𝑓(𝑡𝑖 , 𝑦(𝑡𝑖 )) + 𝑓′(𝑡𝑖 , 𝑦𝑖 )] + 𝑜(ℎ3 )
2
Ceci conduit à l’Algorithme de Taylor (d’ordre 2):
ℎ
𝑦𝑖+1 = 𝑦𝑖 + ℎ (𝑓(𝑡𝑖 , 𝑦𝑖 ) + 2 𝑓′(𝑡𝑖 , 𝑦𝑖 )) , 0≤𝑖 ≤𝑁−1
(Taylor) {
𝑡𝑖+1 = 𝑡𝑖 + ℎ

4.2.4 Méthode de Runge-Kutta (d’ordre 2)

En pratique les méthodes de Taylor sont rarement utilisé car elles demandent le calcul de la
dérivée de la fonction f(t,y), ainsi que leurs valeurs aux points (ti, yi). Alors que, les méthodes
de Runge-kutta ont presque la même précision de la méthode de Taylor, et elles ne demandent
pas le calcul de la dérivée de f.
L’algorithme de Runge-Kutta (d’ordre 2) est donnée par :
1
𝑦𝑖+1 = 𝑦𝑖 + (𝑘1 + 𝑘2 )
2
𝑘1 = ℎ𝑓(𝑡𝑖 , 𝑦𝑖 )
𝑘2 = ℎ𝑓(𝑡𝑖+1 , 𝑦𝑖 + 𝑘1 )
4.2.4 Méthode de Runge-Kutta (d’ordre 4)
L’algorithme de Runge-Kutta (d’ordre 4) est donnée par :
1
𝑦𝑖+1 = 𝑦𝑖 + (𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4 )
6
𝑘1 = ℎ𝑓(𝑡𝑖 , 𝑦𝑖 )
ℎ 𝑘1
𝑘2 = ℎ𝑓 (𝑡𝑖 + , 𝑦𝑖 + )
2 2
ℎ 𝑘2
𝑘3 = ℎ𝑓 (𝑡𝑖 + , 𝑦𝑖 + )
2 2
𝑘4 = ℎ𝑓(𝑡𝑖 + ℎ, 𝑦𝑖 + 𝑘3 )

Dr. DOUAK MOHAMED Page 3

Chapitre 4 Résolution numérique des équations différentielles ordinaire

Exemples :
Soit le problème de Cauchy
𝑦 ′ = 𝑡 − 𝑦 + 1, 𝑡 ∈ [0, 1]
{ (1)
𝑦(0) = 1, ℎ = 0.1
Trouver les valeurs approchées de y(ti) pour : i=0,1,...,n, par:
1/ la méthode d’Euler
2/ la méthode d’Euler modifié
3/la méthode de Taylor
4/la méthode de Rung-Kutta-2
5/la méthode de Rung-Kutta-4
Solution
On a : 𝑦 ′ = 𝑓(𝑡, 𝑦) = 𝑡 − 𝑦 + 1, 𝑡 ∈ [0, 1]

1/ la méthode d’Euler
𝑦𝑖+1 = 𝑦𝑖 + ℎ𝑓(𝑡𝑖 , 𝑦𝑖 ), 𝑝𝑜𝑢𝑟 𝑖 = 0,1, … , 𝑛 − 1 = 9
𝑦(0) = 𝑦0 = 1
𝑦𝑖+1 = 𝑦𝑖 + ℎ(𝑡𝑖 − 𝑦𝑖 + 1)=0.9𝑦𝑖 + 0.1𝑡𝑖 + 0.1 𝑝𝑜𝑢𝑟 𝑖 = 0,1, … ,9
𝑦0 = 1
𝑦1 = 0.9𝑦0 + 0.1𝑡0 + 0.1 = 0.9(1) + 0.1(0) + 0.1 = 1.
𝑦2 = 0.9𝑦1 + 0.1𝑡1 + 0.1 = 0.9(1) + 0.1(0.1) + 0.1 = 1.010
𝑦3 = 0.9𝑦2 + 0.1𝑡2 + 0.1 = 0.9(1.01) + 0.1(0.2) + 0.1 = 1.029
𝑦4 = 0.9𝑦3 + 0.1𝑡3 + 0.1 = 0.9(1.029) + 0.1(0.3) + 0.1 = 1.056
𝑦5 = 0.9𝑦4 + 0.1𝑡4 + 0.1 = 0.9(1.056) + 0.1(0.4) + 0.1 = 1.090
𝑦6 = 0.9𝑦5 + 0.1𝑡5 + 0.1 = 0.9(1.090) + 0.1(0.5) + 0.1 = 1.131
𝑦7 = 0.9𝑦6 + 0.1𝑡6 + 0.1 = 0.9(1.131) + 0.1(0.6) + 0.1 = 1.178
𝑦8 = 0.9𝑦7 + 0.1𝑡7 + 0.1 = 0.9(1.178) + 0.1(0.7) + 0.1 = 1.230
𝑦9 = 0.9𝑦8 + 0.1𝑡8 + 0.1 = 0.9(1.230) + 0.1(0.8) + 0.1 = 1.287
𝑦10 = 0.9𝑦9 + 0.1𝑡9 + 0.1 = 0.9(1.287) + 0.1(0.9) + 0.1 = 1.349
Programme de la méthode d’Euler en FORTRAN
! Euler's method for solving an ordinary differential equation
program main
real :: a= 0.0, b= 1.0, y=1., h, t, f
integer :: n= 10, k
open(100,file='resultats.dat')
h = (b - a)/ real(n)
t=a
k=0

Dr. DOUAK MOHAMED Page 4

Chapitre 4 Résolution numérique des équations différentielles ordinaire

write(100,*)'k t y f(t,y)'
write(*,*)' k t y f(t,y)'
do k = 0,n
write(*,10)k,t,y,f(t,y)
write(100,10)k,t,y,f(t,y)
y = y + h*f(t,y)
t=t+h
end do
10 format(i3,2x,f5.2,2f11.4)
read(*,*)red
end program main

function f(t,y)
real, intent(in)::y, t
f = t-y+1
end function f

Résultats du programme de FORTRAN

k t y f(t,y)
0 0.00 1.0000 0.0000
1 0.10 1.0000 0.1000
2 0.20 1.0100 0.1900
3 0.30 1.0290 0.2710
4 0.40 1.0561 0.3439
5 0.50 1.0905 0.4095
6 0.60 1.1314 0.4686
7 0.70 1.1783 0.5217
8 0.80 1.2305 0.5695
9 0.90 1.2874 0.6126
10 1.00 1.3487 0.6513

Programme de la méthode d’Euler en MATLAB

clc;
clear all;
h=0.1;
x=0:h:1;
y=zeros(1,length(x)-1);
y(1)=1;
f_xy=@(x,y) x-y+1;
for i=1:(length(x)-1)
y(i+1)=y(i)+h*f_xy(x(i),y(i))
end
sol=[y]'

2/ la méthode d’Euler modifié

ℎ
𝑦𝑖+1,𝑐 = 𝑦𝑖 + (𝑓(𝑡𝑖 , 𝑦𝑖 ) + 𝑓(𝑡𝑖+1 , 𝑦𝑖+1,𝑝 ))
2
𝑦𝑖+1,𝑝 = 𝑦𝑖 + ℎ𝑓(𝑡𝑖 , 𝑦𝑖 ) = 0.9𝑦𝑖 + 0.1𝑡𝑖 + 0.1
ℎ
𝑦𝑖+1,𝑐 = 𝑦𝑖 + (𝑡𝑖 − 𝑦𝑖 + 1 + 𝑡𝑖+1 − 𝑦𝑖+1,𝑝 + 1)
2
ℎ
𝑦𝑖+1,𝑐 = 𝑦𝑖 + (𝑡𝑖 − 𝑦𝑖 + 1 + 𝑡𝑖+1 − (0.9𝑦𝑖 + 0.1𝑡𝑖 + 0.1) + 1)
2
𝑦𝑖+1,𝑐 = 0.905𝑦𝑖 + 0.095𝑡𝑖 + 0.1 𝑝𝑜𝑢𝑟 𝑖 = 0,1, … , 𝑛 − 1 = 9

Dr. DOUAK MOHAMED Page 5

Chapitre 4 Résolution numérique des équations différentielles ordinaire

Programme de la méthode d’Euler modifier en FORTRAN

! Euler's modifiè method for solving an ordinary differential equation
program main
real(4) :: a= 0.0, b= 1., y,yp, h, t, f
integer :: n= 10, k
open(100,file='resultats.dat')
h = (b - a)/ real(n)
t=a
k=0
y=1.;
write(*,*) ' k, t, y, f(t,y), yp, f(t+h,yp)'
write(100,*) ' k, t, y, f(t,y), yp, f(t+h,yp)'
do k = 0,n
yp = y + h*f(t,y)
write(*,10) k,t,y,f(t,y),yp,f(t+h,yp)
write(100,10) k,t,y,f(t,y),yp,f(t+h,yp)
y=y+(h/2)*(f(t,y)+f(t+h,yp))
t=t+h
end do
10 format(i3,2x,5f11.4)
read(*,*)red
end program main

function f(t,y)
real, intent(in)::y, t
f = t-y+1
end function f
Résultats du programme de FORTRAN
k, t, y, f(t,y), yp, f(t+h,yp)
0 0.0000 1.0000 0.0000 1.0000 0.1000
1 0.1000 1.0050 0.0950 1.0145 0.1855
2 0.2000 1.0190 0.1810 1.0371 0.2629
3 0.3000 1.0412 0.2588 1.0671 0.3329
4 0.4000 1.0708 0.3292 1.1037 0.3963
5 0.5000 1.1071 0.3929 1.1464 0.4536
6 0.6000 1.1494 0.4506 1.1945 0.5055
7 0.7000 1.1972 0.5028 1.2475 0.5525
8 0.8000 1.2500 0.5500 1.3050 0.5950
9 0.9000 1.3072 0.5928 1.3665 0.6335
10 1.0000 1.3685 0.6315 1.4317 0.6683
Programme de la méthode d’Euler modifier en MATLAB

clc;
clear all;
h=0.1;
x=0:h:1;
y=zeros(1,length(x)-1);
yp=zeros(1,length(x)-1);
yp(1)=1;y(1)=1;
f_xy=@(x,y) x-y+1;
for i=1:(length(x)-1)
yp(i)=y(i)+h*f_xy(x(i),y(i));
y(i+1)=y(i)+0.5*h*(f_xy(x(i),y(i))+f_xy(x(i)+h,yp(i)));
end
sol=[y]'

3/la méthode de Taylor

Dr. DOUAK MOHAMED Page 6

Chapitre 4 Résolution numérique des équations différentielles ordinaire

ℎ
𝑦𝑖+1 = 𝑦𝑖 + ℎ (𝑓(𝑡𝑖 , 𝑦𝑖 ) + 𝑓′(𝑡𝑖 , 𝑦𝑖 ))
2

Or 𝑓 ′ (𝑡, 𝑦) = 1 − 𝑦 ′(𝑡) = 1 − (𝑡 − 𝑦 + 1) = −𝑡 + 𝑦
ℎ
𝑦𝑖+1 = 𝑦𝑖 + ℎ (𝑡𝑖 − 𝑦𝑖 + 1 + (−𝑡𝑖 + 𝑦𝑖 ))
2

𝑦𝑖+1 = 0.905𝑦𝑖 + 0.095𝑡𝑖 + 0.1 𝑝𝑜𝑢𝑟 𝑖 = 0,1, … , 𝑛 − 1 = 9

Programme de la méthode de Taylor en FORTRAN
! Taylor series method (order 2) for solving an ordinary differential equation
program main
real :: a=0.0,b=1., y, h, t,f,df,red
integer :: n=10, k
open(100,file='resultats.dat')
h=(b-a)/real(n)
t=a; y=1.
write(*,*)' k, t, y, f(t,y), df(t,y)'
write(100,*)' k, t, y, f(t,y), df(t,y)'
do k=0,n
write(100,10) k,t,y,f(t,y),df(t,y)
write(*,10) k,t,y,f(t,y),df(t,y)
y = y + h*(f(t,y) + (h/2.0)*df(t,y))
t = t+h
end do
10 format(i3,4f11.4)
read(*,*) red
end program main
function f(t,y)
real, intent(in)::y, t
f = t-y+1
end function f
function df(t,y)
real, intent(in)::y, t
df = y-t
end function df

Résultats du programme de FORTRAN

k, t, y, f(t,y), df(t,y)
0 0.0000 1.0000 0.0000 1.0000
1 0.1000 1.0050 0.0950 0.9050
2 0.2000 1.0190 0.1810 0.8190
3 0.3000 1.0412 0.2588 0.7412
4 0.4000 1.0708 0.3292 0.6708
5 0.5000 1.1071 0.3929 0.6071
6 0.6000 1.1494 0.4506 0.5494
7 0.7000 1.1972 0.5028 0.4972
8 0.8000 1.2500 0.5500 0.4500
9 0.9000 1.3072 0.5928 0.4072
10 1.0000 1.3685 0.6315 0.3685
Programme de la méthode de Taylor en MATLAB

clc;
clear all;
h=0.1;

Dr. DOUAK MOHAMED Page 7

Chapitre 4 Résolution numérique des équations différentielles ordinaire

x=0:h:1;
y=zeros(1,length(x)-1);
y(1)=1;
f_xy=@(x,y) x-y+1;
df_xy=@(x,y) y-x; % dérivé de f_xy
for i=1:(length(x)-1)
y(i+1)=y(i)+h*(f_xy(x(i),y(i))+0.5*h*df_xy(x(i),y(i)));
end
sol=[y]'

4/la méthode de Rung-Kutta-2

𝑘1 = ℎ𝑓(𝑡𝑖 , 𝑦𝑖 ) = 0.1(𝑡𝑖 − 𝑦𝑖 + 1) = 0.1𝑡𝑖 − 0.1𝑦𝑖 + 0.1
𝑘2 = ℎ𝑓(𝑡𝑖+1 , 𝑦𝑖 + 𝑘1 ) = 0.1(𝑡𝑖+1 − (𝑦𝑖 + 𝑘1 ) + 1)
𝑘2 = 0.1(𝑡𝑖+1 − (𝑦𝑖 + 0.1𝑡𝑖 − 0.1𝑦𝑖 + 0.1) + 1) = 0.09𝑡𝑖 − 0.09𝑦𝑖 + 0.1
1
𝑦𝑖+1 = 𝑦𝑖 + (𝑘1 + 𝑘2 )
2
𝑦𝑖+1 = 0.905𝑦𝑖 + 0.095𝑡𝑖 + 0.1 , 𝑝𝑜𝑢𝑟 𝑖 = 0,1, … , 𝑛 − 1 = 9
Programme de la méthode de Rung-Kutta-2 en FORTRAN
program main
real::a=0.,b=1.,t,y=1.,f,h,k1,k2
integer::i,n=10
open(100,file='resultats.dat')
h=(b-a)/n
t=a
write(*,*)' i, t, y, f(t,y), k1, f(t+h,y+k1), k2'
write(100,*)' i, t, y, f(t,y), k1, f(t+h,y+k1), k2'
do i=0,n
k1=h*f(t,y)
k2=h*f(t+h,y+k1)
write(100,10)i,t,y,f(t,y),k1,f(t+h,y+k1),k2
write(*,10)i,t,y,f(t,y),k1,f(t+h,y+k1),k2
y=y+0.5*(k1+k2)
t=t+h
enddo
10 format(i3,6f11.4)
read(*,*)red
end program main
function f(t,y)
real,intent(in)::t,y
f=t-y+1
end function

Résultats du programme :
i, t, y, f(t,y), k1, f(t+h,y+k1), k2
0 0.0000 1.0000 0.0000 0.0000 0.1000 0.0100
1 0.1000 1.0050 0.0950 0.0095 0.1855 0.0185
2 0.2000 1.0190 0.1810 0.0181 0.2629 0.0263
3 0.3000 1.0412 0.2588 0.0259 0.3329 0.0333
4 0.4000 1.0708 0.3292 0.0329 0.3963 0.0396
5 0.5000 1.1071 0.3929 0.0393 0.4536 0.0454
6 0.6000 1.1494 0.4506 0.0451 0.5055 0.0506
7 0.7000 1.1972 0.5028 0.0503 0.5525 0.0553
8 0.8000 1.2500 0.5500 0.0550 0.5950 0.0595
9 0.9000 1.3072 0.5928 0.0593 0.6335 0.0633

Dr. DOUAK MOHAMED Page 8

Chapitre 4 Résolution numérique des équations différentielles ordinaire

10 1.0000 1.3685 0.6315 0.0631 0.6683 0.0668

Programme de la méthode de Rung-Kutta-2 en MATLAB

clc;
clear all;
h=0.1;
x=0:h:1;
y=zeros(1,length(x)-1);
y(1)=1;
f_xy=@(x,y) x-y+1;
for i=1:(length(x)-1)
k1=h*f_xy(x(i),y(i));
k2=h*f_xy(x(i)+h,y(i)+k1);
y(i+1)=y(i)+(1/2)*(k1+k2);
end
sol=[y]'

5/la méthode de Rung-Kutta-4

𝑘1 = ℎ𝑓(𝑡𝑖 , 𝑦𝑖 ) = 0.1𝑡𝑖 − 0.1𝑦𝑖 + 0.1
ℎ 𝑘1
𝑘2 = ℎ𝑓 (𝑡𝑖 + , 𝑦𝑖 + ) = 0.095𝑡𝑖 − 0.095𝑦𝑖 + 0.1
2 2
ℎ 𝑘2
𝑘3 = ℎ𝑓 (𝑡𝑖 + , 𝑦𝑖 + ) = 0.095𝑡𝑖 − 0.095𝑦𝑖 + 0.1
2 2
𝑘4 = ℎ𝑓(𝑡𝑖 + ℎ, 𝑦𝑖 + 𝑘3 ) = 0.091𝑡𝑖 − 0.091𝑦𝑖 + 0.1
1
𝑦𝑖+1 = 𝑦𝑖 + (𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4 )
6
𝑦𝑖+1 = 0.905𝑦𝑖 + 0.095𝑡𝑖 + 0.1 , 𝑝𝑜𝑢𝑟 𝑖 = 0,1, … , 𝑛 − 1 = 9

Programme de la méthode de Rung-Kutta-4 en FORTRAN

Runge-Kutta method of order 4 for solving an initial value problem (rk4,f)
program main
real :: K1, K2, K3, K4,a = 0.0, b = 1., y=1., h, t
integer :: n =10,i
open(100,file='resultats.dat')
h=(b-a)/n
t=a
write(*,*)'i, t,y,f(t, y),K1,f(t + h/2.0, y + 0.5*K1),K2,f(t + h/2.0, y + 0.5*K2),K3,f(t + h, y + K3),K4'
do i=0,n
K1 = h*f(t, y)
K2 = h*f(t + h/2.0, y + 0.5*K1)
K3 = h*f(t + h/2.0, y + 0.5*K2)
K4 = h*f(t + h, y + K3)
write(100,10)i,t,y,f(t, y),K1,f(t + h/2.0, y + 0.5*K1),K2,f(t + h/2.0, y + 0.5*K2),K3,f(t + h, y + K3),K4
write(*,10)i,t,y,f(t, y),K1,f(t + h/2.0, y + 0.5*K1),K2,f(t + h/2.0, y + 0.5*K2),K3,f(t + h, y + K3),K4
y = y + (K1 + 2*K2 +2*K3 + K4)/6.0
t=t+h
end do
10 format(i2,10f9.4)
read(*,*)red

Dr. DOUAK MOHAMED Page 9

Chapitre 4 Résolution numérique des équations différentielles ordinaire

end program main

function f (t,y)
real, intent(in) :: t, y
f = t-y+1
end function f

Résultats du programme de FORTRAN

i, t, y, f(t, y), K1, f(t + h/2.0, y + 0.5*K1),K2,f(t + h/2.0, y + 0.5*K2),K3, f(t + h, y + K3), K4
0 0.0000 1.0000 0.0000 0.0000 0.0500 0.0050 0.0475 0.0047 0.0952 0.0095
1 0.1000 1.0048 0.0952 0.0095 0.1404 0.0140 0.1381 0.0138 0.1813 0.0181
2 0.2000 1.0187 0.1813 0.0181 0.2222 0.0222 0.2202 0.0220 0.2593 0.0259
3 0.3000 1.0408 0.2592 0.0259 0.2962 0.0296 0.2944 0.0294 0.3297 0.0330
4 0.4000 1.0703 0.3297 0.0330 0.3632 0.0363 0.3615 0.0362 0.3935 0.0394
5 0.5000 1.1065 0.3935 0.0393 0.4238 0.0424 0.4223 0.0422 0.4512 0.0451
6 0.6000 1.1488 0.4512 0.0451 0.4786 0.0479 0.4773 0.0477 0.5035 0.0503
7 0.7000 1.1966 0.5034 0.0503 0.5282 0.0528 0.5270 0.0527 0.5507 0.0551
8 0.8000 1.2493 0.5507 0.0551 0.5731 0.0573 0.5720 0.0572 0.5935 0.0593
9 0.9000 1.3066 0.5934 0.0593 0.6138 0.0614 0.6127 0.0613 0.6322 0.0632
10 1.0000 1.3679 0.6321 0.0632 0.6505 0.0651 0.6496 0.0650 0.6672 0.0667

Programme de la méthode de Rung-Kutta-4 en MATLAB

clc;
clear all;
h=0.1;
x=0:h:1;
y=zeros(1,length(x)-1);
y(1)=1;
f_xy=@(x,y) x-y+1;
for i=1:(length(x)-1)
k1=h*f_xy(x(i),y(i));
k2=h*f_xy(x(i)+0.5*h,y(i)+0.5*k1);
k3=h*f_xy(x(i)+0.5*h,y(i)+0.5*k2);
k4=h*f_xy(x(i)+h,y(i)+k3);
y(i+1)=y(i)+(1/6)*(k1+2*k2+2*k3+k4);
end
sol=[y]'

Résolution des EDO (1) avec Matlab

Les principales fonctions de MATLAB qui permettent la résolution des équations
différentielles sont :
- ode23 : Méthode de Runge-Kutta dite d’ordre 2,3
- ode45 : Méthode de Runge-Kutta dite d’ordre 4,5 (la plus utilisée)
Exemple 1 :
𝑦 ′ = 𝑡 − 𝑦 + 1, 𝑡 ∈ [0, 1]
Pour résoudre cette équation différentielle d’ordre 1 à l’aïd du solveur ode, les différentes
étapes sont :
1. On écrit un fichier .m appelé equadiff.m dans lequel
(a) On définit le domaine d’étude : ti=tinitial, tf=tfinal
deltat=[ti tf]
(b) On définit la condition initiale : y0=y(0)

Dr. DOUAK MOHAMED Page 10

Chapitre 4 Résolution numérique des équations différentielles ordinaire

yinit=y0
(c) On appelle ode45 on écrivant :
[t,y]=ode45(f,deltat,yinit);
- t: vecteur contenant les instants auxquels la solution a été calculée,
- y : vecteur contenant la solution de l’équation différentielle
(d) on fait la représentation graphique de la solution en fonction de t à l’aide de
plot(t,y).
Le programme equadiff permettant de réaliser les différentes étapes est :
clc; clear all; close all;
f=@(t,y)t-y+1;
ti=0;
tf=1;
y0=1;
[t y]=ode45(f,[0 tf],y0);
figure;
plot(t,y)
xlabel('t');
ylabel('y');
grid on

4.3 Équations différentielles ordinaire d’ordre supérieur à 1 et les systèmes d’équations

différentielles

Considérons le problème suivant :

Dr. DOUAK MOHAMED Page 11

Chapitre 4 Résolution numérique des équations différentielles ordinaire

𝑦 (𝑚) (𝑡) = 𝑓(𝑡, 𝑦, 𝑦 ′ , … . , 𝑦 (𝑚−1) ), 𝑡 ∈ [𝑎, 𝑏]

(𝐼) {
𝑦(𝑎) = 𝛼1 ; 𝑦 ′(𝑎) = 𝛼2 ; … … … 𝑦 (𝑚−1) (𝑎) = 𝛼𝑚 ;

L’équation différentielle d’ordre m du problème (I) peut être transformer en un système de m

équations différentielles d’ordre 1, en effet posons :
𝑢1 (𝑡) = 𝑦(𝑡), 𝑢2 (𝑡) = 𝑦 ′ (𝑡); … … … , 𝑢𝑚 (𝑡) = 𝑦 (𝑚−1) (𝑡)
Nous obtenons le système suivant :
𝑑𝑢1
= 𝑢2
𝑑𝑡
𝑑𝑢2
= 𝑢3
𝑑𝑡
(𝐼) .
.
.
𝑑𝑢𝑚
{ = 𝑓(𝑡, 𝑢1 , 𝑢2 , … , 𝑢𝑚−1 )
𝑑𝑡

Les méthodes de résolution des EDO (1) : Euler, Euler modifier, Runge-kutta (2) et (4)
peuvent être généralisées pour résoudre le problème (II).

Résolution des EDO (2)

Soit l’équation différentielle suivante :
𝑑2 𝑦 𝑑𝑦
(𝐼) { 𝑑𝑡 2 = 𝑓 (𝑡, 𝑦, 𝑑𝑡 ) , 𝑡 ∈ [𝑎, 𝑏]
𝑦(𝑎) = 𝛼1 ; 𝑦 ′ (𝑎) = 𝛼2
Posons : 𝑢1 (𝑡) = 𝑦(𝑡), 𝑢2 (𝑡) = 𝑦 ′ (𝑡) on obtient le système d’équation suivant :
𝑑𝑢1
= 𝑢2 = 𝑓1 (𝑡, 𝑢1 , 𝑢2 )
𝑑𝑡
(𝐼𝐼) 𝑑𝑢2
= 𝑓2 (𝑡, 𝑢1 , 𝑢2 )
𝑑𝑡
{ 𝑢1 (𝑎) = 𝛼1 ; 𝑢2 (𝑎) = 𝛼2
Résolution du système d’équation (II)

1/ Méthode d’Euler :
𝑢1,𝑖+1 = 𝑢1𝑖 + ℎ𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) ; 𝑢10 = 𝛼1
𝑢2,𝑖+1 = 𝑢2𝑖 + ℎ𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) ; 𝑢20 = 𝛼2

2/ la méthode d’Euler modifié

𝑢1,𝑖+1,𝑝 = 𝑢1𝑖 + ℎ𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 )
𝑢2,𝑖+1,𝑝 = 𝑢2𝑖 + ℎ𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 )

Dr. DOUAK MOHAMED Page 12

Chapitre 4 Résolution numérique des équations différentielles ordinaire

ℎ
𝑢1,𝑖+1,𝑐 = 𝑢1𝑖 + (𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) + 𝑓1 (𝑡𝑖+1 , 𝑢1,𝑖+1,𝑝 , 𝑢2,𝑖+1,𝑝 ))
2
ℎ
𝑢2,𝑖+1,𝑐 = 𝑢2𝑖 + (𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) + 𝑓2 (𝑡𝑖+1 , 𝑢1,𝑖+1,𝑝 , 𝑢2,𝑖+1,𝑝 ))
2

3/la méthode de Taylor

ℎ
𝑢1,𝑖+1 = 𝑢1𝑖 + ℎ (𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) + 𝑓1′ (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ))
2

ℎ
𝑢2,𝑖+1 = 𝑢2𝑖 + ℎ (𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) + 𝑓2′ (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ))
2

4/la méthode de Rung-Kutta-2

𝑘11 = ℎ𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 )
𝑘12 = ℎ𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 )
𝑘21 = ℎ𝑓1 (𝑡𝑖 + ℎ, 𝑢1𝑖 + 𝑘11 , 𝑢2𝑖 + 𝑘12 )
𝑘22 = ℎ𝑓2 (𝑡𝑖 + ℎ, 𝑢1𝑖 + 𝑘11 , 𝑢2𝑖 + 𝑘12 )
1 1
𝑢1,𝑖+1 = 𝑢1𝑖 + 2 (𝑘11 + 𝑘21 ) ; 𝑢2,𝑖+1 = 𝑢2𝑖 + 2 (𝑘12 + 𝑘22 )

5/la méthode de Rung-Kutta-4

𝑘11 = ℎ𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 )
𝑘12 = ℎ𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 )
ℎ 𝑘11 𝑘12
𝑘21 = ℎ𝑓1 (𝑡𝑖 + , 𝑢1𝑖 + , 𝑢2𝑖 + )
2 2 2
ℎ 𝑘11 𝑘12
𝑘22 = ℎ𝑓2 (𝑡𝑖 + , 𝑢1𝑖 + , 𝑢2𝑖 + )
2 2 2
ℎ 𝑘21 𝑘22
𝑘31 = ℎ𝑓1 (𝑡𝑖 + , 𝑢1𝑖 + , 𝑢2𝑖 + )
2 2 2
ℎ 𝑘21 𝑘22
𝑘32 = ℎ𝑓2 (𝑡𝑖 + , 𝑢1𝑖 + , 𝑢2𝑖 + )
2 2 2
𝑘41 = ℎ𝑓1 (𝑡𝑖 + ℎ, 𝑢1𝑖 + 𝑘31 , 𝑢2𝑖 + 𝑘32 )
𝑘42 = ℎ𝑓2 (𝑡𝑖 + ℎ, 𝑢1𝑖 + 𝑘31 , 𝑢2𝑖 + 𝑘32 )
1
𝑢1,𝑖+1 = 𝑢1𝑖 + (𝑘11 + 2𝑘21 + 2𝑘31 + 𝑘41 )
6
1
𝑢2,𝑖+1 = 𝑢2𝑖 + (𝑘12 + 2𝑘22 + 2𝑘32 + 𝑘42 )
6

Exemple d’application :

Dr. DOUAK MOHAMED Page 13

Chapitre 4 Résolution numérique des équations différentielles ordinaire

𝑑𝑢1
= 𝑢1 𝑢2 + 𝑡 = 𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 )
𝑑𝑡
(𝐼𝐼) 𝑑𝑢2
= 𝑡𝑢2 + 𝑢1 = 𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 )
𝑑𝑡
{ 𝑢1 (0) = 𝑢1,0 = 1; 𝑢2 (0) = 𝑢2,0 = −1; ℎ = 0.1
Calculer u1(0.1) et u2(0.1) ?

1/ Méthode d’Euler :
𝑢1,𝑖+1 = 𝑢1𝑖 + ℎ𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) = 𝑢1𝑖 + 0.1(𝑢1𝑖 𝑢2𝑖 + 𝑡𝑖 )
𝑢1,1 = 𝑢1,0 + 0.1(𝑢1,00 𝑢2,0 + 𝑡0 ) = 1 + 0.1(−1 + 0) = 0.9
𝑢2,𝑖+1 = 𝑢2,𝑖 + ℎ𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) = 𝑢2𝑖 + 0.1(𝑡𝑖 𝑢2𝑖 + 𝑢1𝑖 )
𝑢2,1 = 𝑢20 + 0.1(𝑡0 𝑢2,0 + 𝑢1,0 ) = −1 + 0.1(0 + 1) = −0.9

Résultats du programme
k, t, u1, u 2, f1, f2
0 0.0000 1.0000 -1.0000 -1.0000 1.0000
1 0.1000 0.9000 -0.9100 -0.7190 0.8090
2 0.2000 0.8281 -0.8363 -0.4925 0.6608
3 0.3000 0.7788 -0.7751 -0.3037 0.5463
4 0.4000 0.7485 -0.7235 -0.1416 0.4591
5 0.5000 0.7343 -0.6790 0.0014 0.3948

2/ la méthode d’Euler modifié

𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) = 𝑢1,𝑖 𝑢2,𝑖 + 𝑡𝑖 ; 𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) = 𝑡𝑖 𝑢2,𝑖 + 𝑢1,𝑖
𝑢1,0 = 1; 𝑢2,0 = −1
𝑢1,𝑖+1,𝑝 = 𝑢1𝑖 + ℎ𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) = 𝑢1𝑖 + 0.1(𝑢1𝑖 𝑢2𝑖 + 𝑡𝑖 )
𝑢2,𝑖+1,𝑝 = 𝑢2𝑖 + ℎ𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) = 𝑢2𝑖 + 0.1(𝑡𝑖 𝑢2,𝑖 + 𝑢1,𝑖 )
ℎ
𝑢1,𝑖+1,𝑐 = 𝑢1𝑖 + (𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) + 𝑓1 (𝑡𝑖+1 , 𝑢1,𝑖+1,𝑝 , 𝑢2,𝑖+1,𝑝 ))
2
ℎ
𝑢2,𝑖+1,𝑐 = 𝑢2𝑖 + (𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) + 𝑓2 (𝑡𝑖+1 , 𝑢1,𝑖+1,𝑝 , 𝑢2,𝑖+1,𝑝 ))
2
Résultats du programme
k, t, u1, u2, f1, f2, up1, up2, f1p, f2p
0 0.0000 1.0000 -1.0000 -1.0000 1.0000 0.9000 -0.9000 -0.7100 0.8100
1 0.1000 0.9145 -0.9138 -0.7356 0.8231 0.8409 -0.8315 -0.4992 0.6746
2 0.2000 0.8528 -0.8420 -0.5180 0.6844 0.8010 -0.7735 -0.3196 0.5689
3 0.3000 0.8109 -0.7814 -0.3336 0.5765 0.7775 -0.7238 -0.1627 0.4880

Dr. DOUAK MOHAMED Page 14

Chapitre 4 Résolution numérique des équations différentielles ordinaire

4 0.4000 0.7861 -0.7294 -0.1734 0.4943 0.7687 -0.6800 -0.0227 0.4287

5 0.5000 0.7763 -0.6838 -0.0308 0.4344 0.7732 -0.6403 0.1049 0.3890

3/la méthode de Taylor

ℎ
𝑢1,𝑖+1 = 𝑢1𝑖 + ℎ (𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) + 𝑓1′ (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ))
2

ℎ
𝑢2,𝑖+1 = 𝑢2𝑖 + ℎ (𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) + 𝑓2′ (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ))
2
𝑑𝑢1 𝑑𝑢2
𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) = 𝑢1,𝑖 𝑢2,𝑖 + 𝑡𝑖 ; 𝑓1′ (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) = . 𝑢2,𝑖 + . 𝑢1,𝑖 = 𝑓1 . 𝑢2,𝑖 + 𝑓2 . 𝑢1,𝑖
𝑑𝑡 𝑑𝑡

𝑓1′ (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) = (𝑢1,𝑖 𝑢2,𝑖 + 𝑡𝑖 ). 𝑢2,𝑖 + (𝑡𝑖 𝑢2,𝑖 + 𝑢1,𝑖 ). 𝑢1,𝑖
𝑑𝑢2 𝑑𝑢1
𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) = 𝑡𝑖 𝑢2,𝑖 + 𝑢1,𝑖 ; 𝑓2′ (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) = 𝑢2,𝑖 + 𝑡𝑖 + = 𝑢2,𝑖 + 𝑡𝑖 𝑓2 + 𝑓1
𝑑𝑡 𝑑𝑡

𝑓2′ (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) = 𝑢2,𝑖 + 𝑡𝑖 (𝑡𝑖 𝑢2,𝑖 + 𝑢1,𝑖 ) + 𝑢1,𝑖 𝑢2,𝑖 + 𝑡𝑖

Résultats du programme :
k, t, u1, u2
0 0.0000 1.0000 -1.0000
1 0.1000 0.9150 -0.9100
2 0.2000 0.8538 -0.8354
3 0.3000 0.8126 -0.7728
4 0.4000 0.7884 -0.7193
5 0.5000 0.7793 -0.6727

4/la méthode de Rung-Kutta-2

𝑘11 = ℎ𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) = 0.1(𝑢1,𝑖 𝑢2,𝑖 + 𝑡𝑖 )
𝑘12 = ℎ𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 ) = 0.1(𝑡𝑖 𝑢2,𝑖 + 𝑢1,𝑖 )
𝑘21 = ℎ𝑓1 (𝑡𝑖 + ℎ, 𝑢1𝑖 + 𝑘11 , 𝑢2𝑖 + 𝑘12 ) = 0.1[(𝑢1𝑖 + 𝑘11 )(𝑢2𝑖 + 𝑘12 ) + 𝑡𝑖 + ℎ]
𝑘22 = ℎ𝑓2 (𝑡𝑖 + ℎ, 𝑢1𝑖 + 𝑘11 , 𝑢2𝑖 + 𝑘12 ) = 0.1[(𝑡𝑖 + ℎ)(𝑢2𝑖 + 𝑘12 ) + 𝑢1𝑖 + 𝑘11 ]
1 1
𝑢1,𝑖+1 = 𝑢1𝑖 + 2 (𝑘11 + 𝑘21 ) ; 𝑢2,𝑖+1 = 𝑢2𝑖 + 2 (𝑘12 + 𝑘22 )

Résultats du programme
i, t, u1, u2, k11, k12, k21, k22
0 0.0000 1.0000 -1.0000 -0.1000 0.1000 -0.0710 0.0810
1 0.1000 0.9145 -0.9095 -0.0732 0.0824 -0.0496 0.0676
2 0.2000 0.8531 -0.8345 -0.0512 0.0686 -0.0314 0.0572
3 0.3000 0.8118 -0.7716 -0.0326 0.0580 -0.0156 0.0494

Dr. DOUAK MOHAMED Page 15

Chapitre 4 Résolution numérique des équations différentielles ordinaire

4 0.4000 0.7877 -0.7179 -0.0165 0.0501 -0.0015 0.0437

5 0.5000 0.7787 -0.6710 -0.0022 0.0443 0.0113 0.0400

5/la méthode de Rung-Kutta-4

𝑘11 = ℎ𝑓1 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 )
𝑘12 = ℎ𝑓2 (𝑡𝑖 , 𝑢1𝑖 , 𝑢2𝑖 )
ℎ 𝑘11 𝑘12
𝑘21 = ℎ𝑓1 (𝑡𝑖 + , 𝑢1𝑖 + , 𝑢2𝑖 + )
2 2 2
ℎ 𝑘11 𝑘12
𝑘22 = ℎ𝑓2 (𝑡𝑖 + , 𝑢1𝑖 + , 𝑢2𝑖 + )
2 2 2
ℎ 𝑘21 𝑘22
𝑘31 = ℎ𝑓1 (𝑡𝑖 + , 𝑢1𝑖 + , 𝑢2𝑖 + )
2 2 2
ℎ 𝑘21 𝑘22
𝑘32 = ℎ𝑓2 (𝑡𝑖 + , 𝑢1𝑖 + , 𝑢2𝑖 + )
2 2 2
𝑘41 = ℎ𝑓1 (𝑡𝑖 + ℎ, 𝑢1𝑖 + 𝑘31 , 𝑢2𝑖 + 𝑘32 )
𝑘42 = ℎ𝑓2 (𝑡𝑖 + ℎ, 𝑢1𝑖 + 𝑘31 , 𝑢2𝑖 + 𝑘32 )
1
𝑢1,𝑖+1 = 𝑢1𝑖 + (𝑘11 + 2𝑘21 + 2𝑘31 + 𝑘41 )
6
1
𝑢2,𝑖+1 = 𝑢2𝑖 + (𝑘12 + 2𝑘22 + 2𝑘32 + 𝑘42 )
6
Résultats du programme
i, t, u1, u2, k11, k12, k21, k22, k31, k32, k41, k42
0 0.0000 1.0000 -1.0000 -0.1000 0.1000 -0.0852 0.0903 -0.0864 0.0910 -0.0730 0.0823
1 0.1000 0.9139 -0.9092 -0.0731 0.0823 -0.0612 0.0747 -0.0620 0.0753 -0.0510 0.0685
2 0.2000 0.8522 -0.8341 -0.0511 0.0685 -0.0411 0.0627 -0.0418 0.0631 -0.0325 0.0579
3 0.3000 0.8106 -0.7711 -0.0325 0.0579 -0.0240 0.0535 -0.0244 0.0538 -0.0164 0.0499
4 0.4000 0.7863 -0.7174 -0.0164 0.0499 -0.0089 0.0467 -0.0093 0.0470 -0.0021 0.0442
5 0.5000 0.7772 -0.6705 -0.0021 0.0442 0.0047 0.0420 0.0044 0.0422 0.0109 0.0405

Résolution des Equation différentielle du second ordre avec Matlab

Soit L’équation différentielle du second ordre connu sous le nom d’équation de l’oscillateur
harmonique :
𝑑2𝑦
+ 𝜔02 𝑦 = 0
𝑑𝑡 2
Pour résoudre cette équation différentielle d’ordre 2 à l’aïd du solveur ode, les différentes
étapes sont :
2. Transformation de l’équation différentielle d’ordre 2 en un système de deux équations
différentielles du premier ordre, pour cela :

Dr. DOUAK MOHAMED Page 16

Chapitre 4 Résolution numérique des équations différentielles ordinaire

(a) On crée 2 vecteurs de dimension 2 correspondant à l’ordre de l’équation

différentielle :
- Vecteur y
- y(1)=y contient la solution y(t)
- y(2)=dy/dt contient la dérivée de la solution dy/dt
- vecteur dydt
- dydt(1)=dy/dt=y(2)
- dydt(2)=d2y/dt2
(b) on transforme l’équation différentielle d’ordre 2 en un système de deux équations
différentielles du premier ordre :
L’équation différentielle
𝑑2𝑦
= −𝜔02 𝑦
𝑑𝑡 2
Est mise sous la forme vectorielle équivalente :
𝑑𝑦𝑑𝑡(1) = 𝑦(2)
{
𝑑𝑦𝑑𝑡(2) = −𝜔02 𝑦(1)
𝑦1 = 𝑦
𝑑𝑦
= 𝑦2
𝑑𝑡
𝑑𝑦2
= −𝜔02 𝑦
𝑑𝑡
{𝑦1 (0) = 1, 𝑦2 (0) = 0, 𝑡 ∈ [0, 1]

f=@(t,y)[y(2);-omega0.^2*y(1)];

3. On écrit un fichier .m appelé equadiff.m dans lequel

(d) On définit le domaine d’étude : ti=tinitial, tf=tfinal
deltat=[ti tf]
(e) On définit les conditions initiales : y0=y(0) et dydt0=dy/dt(0)
yinit=[y0 dydt0]
(f) On appelle ode45 on écrivant :
[t,y]=ode45(f,deltat,yinit);
- t: vecteur contenant les instants auxquels la solution a été calculée,
- y : matrice deux colonnes avec :
-la première colonne y( :,1) contient y(t),
- la seconde colonne y( :,2) contient la dérivée dy/dt.

Dr. DOUAK MOHAMED Page 17

Chapitre 4 Résolution numérique des équations différentielles ordinaire

(d) on fait la représentation graphique de la solution en fonction de t à l’aide de

plot(t,y).
Le programme equadiff permettant de réaliser les différentes étapes est :
clc
clear
close
omega0=2;
f=@(t,y)[y(2);-omega0.^2*y(1)];
ti=0;
tf=1;
yin=1;
dydtin=0;
deltat=[ti tf];
yinit=[yin dydtin];
[t,y]=ode45(f,deltat,yinit);
subplot(211)
plot(t,y(:,1))
xlabel('temps')
ylabel('y')
title(['\omega_0=',num2str(omega0)])
grid on
subplot(212)
plot(t,y(:,2))
xlabel('temps')
ylabel('dy/dt')
grid on

Dr. DOUAK MOHAMED Page 18