Approximation de solutions d’équations différentielles,

schémas numériques.
C. Dossal
Février 2013

1 Le cadre général
On va chercher à approcher numériquement les solutions d’équations différentielles de la
forme :
y0 (t) = f (t, y(t)) avec y(t0 ) = y0 . (1)
où f est une fonction continue Lipschitz par rapport à la deuxième variable. Le théorème de
Cauchy Lipschitz assure l’existence et l’unicité de la solution. Il arrive que l’on sache résoudre ce
problème de manière analytique mais dans un grand nombre de cas on ne connait pas de forme
explicite de la solution. On peut alors essayer d’approcher la solution par un schéma numérique.
Le principe de toutes les méthodes que nous allons voir est le même : On cherche à estimer la
solution y de (EqRef1) en des points régulièrement espacés (ti )i∈I en partant de y(t0 ) qui est connu
et en estimant de proche en proche les valeurs de y(ti ) en utilisant le fait que
Z ti+1
y(ti+1 ) = y(ti ) + f (t, y(t))dt. (2)
ti

De fait, ne connaissant pas y sur l’intervalle [ti ,ti+1 ] on doit souvent approcher cette intégrale
avec les moyens du bord. L’étude des schémas numériques consiste à étudier la convergence et à
borner l’erreur de ces schémas. Dans la suite on notera h la distance entre deux points ti et ti+1 ;
∀i ∈ N, h = ti+1 −ti . On peut faire des schémas avec des pas variables mais nous concentrerons
Rt
sur
les schémas à pas constant. Chaque méthode va essayer d’approcher cette intégrale tii+1 f (t, y(t))dt
par une fonction h × φ (ti , y(ti ), h).

2 Préliminaire
Les 4 équations de référence que nous utiliserons pour tester les différentes métodes sont les
suivantes :
y0 (t) = −y(t) sous la contrainte y(0) = 1. (EqRef1)
y0 (t) = y(t) sous la contrainte y(0) = 1. (EqRef2)
0 2
y (t) = −y (t) sous la contrainte y(0) = 1. (EqRef3)
y0 (t) = − tan(t)y(t) sous la contrainte y(0) = 1. (EqRef4)
1. Résoudre ces 4 équations différentielles.
2. Créer les fonctions f1.m et solf1.m suivantes

1
 function z=f1(t,y)
z=-y;
et
function z=solf1(t);
z=-exp(t);
3. Créer les fonctions analogues pour les autres équations de référence.

3 Schéma d’Euler explicite

Euler a été le premier
Rt
à proposer un schéma d’approximation. Ce schéma consiste à appro-
cher l’intégrale tii+1 f (t, y(t))dt par φ (ti , y(ti ), h) = f (ti , y(ti )). Le schéma dit d’Euler explicite
s’écrit alors
yn+1 = yn + h × f (tn , y(tn )). (3)
4. On considère la solution approchée par la méthode d’Euler de l’équation (EqRef1). Si on
pose h = Tn montrer que pour tout T , la solution obtenue par la métode d’Euler converge
vers la solution analytique au point T quand n tend vers +∞.
5. Montrer que si n est fixé et que T tend vers +∞, yn s’éloigne de 0. On dit que la méthode
d’Euler est un schéma instable. On ne l’utilise qu’en temps fini.
En Matlab la méthode d’Euler peut se coder de la manière suivante :
function y=MethodeEuler(name,t0,T,y0,h)
ff=str2func(name);
N=floor(T/h)+1;
t=[t0:h:t0+(N-1)*h];
y=zeros(N,1);
y(1)=y0;
for k=1:N-1
y(k+1)=y(k)+h*ff(t(k),y(k));
end
On pourra utiliser la fonction suivante pour tester la méthode :
function er=TestMethode(Methode,ff,solff,t0,T,y0,h)
solf=str2func(solff);
N=floor(T/h)+1;
t=[t0:h:t0+(N-1)*h];
if Methode(1)==’A’;
Meth=str2func(Methode(1:end-1));
r=str2num(Methode(end));
[yn,tn,fn]=Meth(r,ff,t0,y0,h,N);
else
Meth=str2func(Methode);
y=Meth(ff,t0,T,y0,h);
end
plot(t,y);
z=zeros(N,1);
for k=1:N
z(k)=solf(t(k));

2
 end
hold on;plot(t,z,’r’);title(Methode)
hold off;
er=max(abs(z-y));
6. Tester la méthode d’Euler sur les équations de référence en faisant varier la valeur du pas
h.
7. Pour ces différentes équations comment évolue l’erreur maximale en fonction du pas h ?

4 Méthode du point milieu

Dans la méthode
Rt
du point milieu on essaie d’être un peu plus précis dans l’estimation de
l’intégrale tii+1 f (t, y(t))dt pour cela on utilise un point intermédiaire entre yn et yn+1 que
nous noterons yn+ 1 :
2
h
yn+ 1 = yn + f (tn , yn ).
2 2
On a alors
h
yn+1 = yn + h f (tn + , yn+ 1 ) et toujours tn+1 = tn + h.
2 2

8. A partir du programme précédent créer une nouvelle fonction MethodePointMilieu qui

résout numériquement une EDO par cette méthode.
9. Tester cette nouvelle méthode sur les équations de référence en faisant varier le pas h.
Comparer les erreurs de cette méthode et de la méthode d’Euler.
10. Comment varie l’erreur en fonction de h ?
11. Un bon indicateur de la complexité d’une méthode est le nombre de fois par étape où l’on
doit évaluer la fonction f . Pour un nombre égal d’évaluations de f quelle méthode assure
l’erreur la plus faible ?

5 Consistance, stabilité et convergence

L’objectif d’un schéma numérique est de fournir une solution approchée qui converge vers
la solution du problème de Cauchy (1) quand h tend vers 0.
Pour estimer cette convergence on définit la notion d’erreur de consistance. Si une méthode
de résolution numérique est de la forme yn+1 = yn + hφ (tn , yn , hn ),
Définition 1. L’erreur de consitance en relative à la solution exacte z de (1) est définie
par
en = z(tn+1 ) − z(tn ) − hφ (tn , z(tn ), h)
C’est l’erreur que l’on commetrait dans le calcul de z(tn+1 ) si on connaissait exactement
z(tn ) par la méthode numérique utilisée. En pratique on connait rarement la valeur exacte
de z(tn ) ...
Définition 2. On dit que la méthode numérique est consistante si pour toute solution
exacte z, limN→∞ ∑06n6N |en | = 0.

3
 Définition 3. On dit que la méthode est stable s’il existe une constante S > 0, appelée
constante de stabilité, telle que pour toutes suites (yn ) et (ỹn ) définie par

yn+1 = yn + hφ (tn , yn , h), 06n6N

ỹn+1 = ỹn + hφ (tn , ỹn , h) + εn , 06n6N

on ait !
max |ỹn − yn | 6 S |ỹ0 − y0 | + ∑ |εn |
06n6N 06n6N

Ainsi si la méthode est stable, même si on commet une petite erreur à chaque étape par
rapport au schéma souhaité, et en pratique c’est toujours le cas à cause des erreurs d’ar-
rondis, on peut contrôler l’erreur globale pourvu que la somme de chacune des erreurs soit
contrôlée.
Définition 4. On dit que la méthode est convergente si pour toute solution exacte z, la
suite yn+1 = yn + hφ (tn , yn , h) vérifie

max |yn − z(tn )| → 0

06n6N

quand y0 → z(t0 ) et quand h → 0.

La quantité max |yn − z(tn )| s’appelle l’erreur globale, c’est cette erreur qui importe dans
06n6N
la pratique.
12. En prenant ỹn = z(tn ) montrer que si la méthode est stable et consistante alors elle est
convergente.
T
13. Soit z une solution exacte de (1), soit T > 0 fixé et N ∈ N, on note h = N et

en = z(tn+1 ) − z(tn ) − hφ (tn , z(tn ), h).

(a) Comment appelle t-on en ?

(b) Montrer que pour tout n < N, il existe cn ∈]tn ,tn+1 [ tel que en = h( f (cn , z(cn ))+φ (tn , z(tn ), h)).
(c) On pose

αn = f (cn , z(cn )) − φ (cn , z(cn ), 0) et βn = φ (cn , z(cn ), 0) − φ (tn , z(tn ), h)

Montrer que pour tout ε > 0 il existe h0 > 0 tel que si h 6 h0 , ∀n 6 N, |βn | 6 ε.
(d) En déduire que si h < h0

| ∑ |en | − ∑ h|αn || 6 T ε.
06n<N 06n<N

(e) Justifier que

Z T
lim ∑ h|αn | = | f (t, z(t)) − φ (t, z(t), 0)|dt.
N→∞ t0
06n<N

(f) En déduire que la méthode est consistante si et seulement si

∀t ∈ [t0 , T ], φ (t, z(t), 0) = f (t, z(t)).

4
 Lemme de Gronval discret

14. Soit (θn )n∈N ∈ R+ et (εn )n∈N ∈ R deux suites telles que θn+1 6 (1 + Λh)θn + |εn |.
Justifier que 1 + Λh 6 eΛh .
15. Montrer par récurrence que pour tout n ∈ N
n−1
θn 6 eΛ(tn −t0 ) θ0 + ∑ eΛ(tn −ti+1 ) |εi |. (4)
i=0

16. Dans cette question on suppose que φ est Lipschitz de constante Λ par rapport à la deuxième
˜ n ) définies par
variable. On considère deux suites (yn ) et (y

yn+1 = yn + hφ (tn , yn , h)
ỹn+1 = ỹn + hφ (tn , ỹn , h) + εn

(a) Montrer que la suite θn = |ỹn − yn | vérifie

θn+1 6 (1 + Λh)θn + |εn |.

(b) En utilisant le lemme de Gronval discret que

!
N
ΛT
max |ỹn − yn | 6 e |ỹ0 − y0 | + ∑ |εn | .
06n6N n=0

(c) En déduire que la méthode numérique associée à la fonction φ est stable de constante
de stabilité eΛT .
17. Que peut-on dire d’une méthode de résolution associée à une fonction φ Lipshitz par rapport
à la deuxième variable et telle que φ (t, y, 0) = f (t, y) ?
18. En déduire que si f est Lipschitz de constante k par rapport à la deuxième variable, le
schéma d’Euler explicite est convergent et donner sa constante de stabilité.
19. On rappelle que pour la méthode du point milieu
 
h h
φ (t, y, h) = f t + , y + f (t, y) .
2 2

Si f est k-Lipschitz montrer que

 
h
|φ (t, y1 , h) − φ (t, y2 , h)| 6 k |y1 − y2 | + | f (t, y1 ) − f (t, y2 )| .
2

20. En déduire que si f est k−Lipschitz la méthode du point milieu est stable et donner sa
constante de stabilité.

5
 6 Ordre d’une méthode numérique de résolution d’une
EDO
Définition 5. On dit qu’une méthode à 1 pas est d’ordre p si pour toute solution exacte
z de (1) où f est de classe C p , il existe une constante C telle que l’erreur de consistance
relative à z vérifie
|en | 6 Ch p+1 , ∀n < N.
21. Montrer que si une méthode est d’ordre p alors
N−1
∑ |en | 6 CT h p .
n=0

22. Ordre de la méthode d’Euler.

On suppose que f est de classe C1 et on note z une solution exacte de (1) avec f ∈ C1 et
on note en l’erreur de consistance de la méthode d’Euler.
(a) Montrer que
1
en = h2 z(2) (tn ) + o(h2 ).
2
En déduire une expression de en en fonction des dérivées partielles de f au point (tn , yn )
et que la méthode d’Euler est d’ordre 1.
23. Ordre de la méthode du point milieu.
On suppose maintenant que f est de classe C2 et on considère la méthode du point milieu.
Si z est une solution de (1), on a z0 (t) = f (t, z(t)) et z est C3 , on en déduit que la dérivée
troisième z(3) de z solution exacte de (1) s’exprime comme combinaison linéaire de produits
de dérivées partielles de f d’ordre 1 et 2 et est donc bornée [t0 ,t0 + T ].
(a) Montrer que l’erreur de consistance en s’écrit en = εn + εn0 avec

h
εn = z(tn+1 ) − z(tn ) − hz0 (tn + )
2
h
εn0 = hz0 (tn + ) − (yn+1 − z(tn )).
2

(b) Montrer qu’il existe C1 tel que εn 6 C1 h3 .

(c) Montrer que     
h h h
εn0 = h f tn + , z(tn + ) − f tn + , yn+ 1
2 2 2 2

(d) En utilisant l’égalité des accroissements finis montrer qu’il existe C2 tel que εn0 6 C2 h3 .
(e) En déduire que la méthode du point milieu est une méthode d’ordre 2.
La fonction suivante permet d’estimer numériquement l’ordre d’une méthode en uti-
lisant plusieurs valeurs du pas h.
function r=OrdreMethode(Methode,fonction,solution,H,t0,y0,T)
f=str2func(fonction);
sol=str2func(solution);
K=length(H);
for k=1:K
h=H(k);

6
 N=floor(T/h)+1;
t=[t0:h:t0+(N-1)*h];
z=zeros(N,1);
for j=1:N
z(j)=sol(t(j));
end
if Methode(1)==’A’;
Meth=str2func(Methode(1:end-1));
r=str2num(Methode(end));
y=Meth(r,fonction,t0,T,y0,h);
else
Meth=str2func(Methode);
y=Meth(fonction,t0,T,y0,h);
end
er(k)=max(abs(y-z));
end
plot(log(er),log(H),’ro’);title(’Diagramme Log log de l Erreur en fonction du pas’);
hold on;
plot(log(er),log(H));
hold off;
coef=polyfit(log(H),log(er),1);
r=coef(1);
On peut l’utiliser en lançant par exemple :
>> r=OrdreMethode(’MethodeEuler’,’f1’,’solf1’,[0.0001,0.001,0.01,0.1],0,1,1)
(f) Utiliser cette fonction pour estimer numériquement l’ordre des méthodes d’Euler et
du point milieu sur les 4 équations de référence.
(g) Cela correspond il aux calculs théoriques ?
(h) Que se passe t-il si on utilise de très petites valeurs du pas pour estimer l’ordre ?
Expliquer le phénomène.
24. La méthode de Heun consite à utiliser la formule de récurrence suivante
1
yn+1 = yn + h ( f (tn , yn ) + f (tn + h, yn + h f (tn , yn ))) .
2
Ecrire un programme qui met en place cette méthode
25. Tester cette nouvelle méthode sur les équations de référence et évaluer numériquement son
ordre.

7
 7 Modèle de Volterra pour les systèmes proie-prédateur
En 1926, Volterra propose un modèle simple de système proie-prédateur afin d’expliquer
les oscillations des campagnes de pêche dans l’adriatique. Il s’écrit
 0
S (t) = S(t)(a − bR(t))
(5)
R0 (t) = R(t)(−c + dS(t))
où t représente le temps, S(t) et R(t) représente respectivement le nombre de proies et de
prédateur à l’instant t et a, b, c et d sont des paramètres positifs.
26. Montrer qu’en posant u(τ) = dc S(t), v(τ) = ab R(t), τ = at et α = ac on obtient le système
suivant.  0
u (τ) = u(τ)(1 − v(τ))
(6)
v0 (τ) = αv(τ)(u(τ) − 1)
27. Soit H(t) = αu(t) + v(t) − α ln(u(t)) − ln(v(t)).
28. Calculer H 0 (t). On dit que H est une intégrale première du système.
29. Soit f la fonction définie de ]0, +∞[ dans R par f (x) = x − ln x. Etudier les variations de f .
30. En déduire que quelques soient les valeurs initiales de u(0) et v(0) choisies, les fonctions u
et v sont bornées supérieurement par un réel M > 0 et inférieurement par un réel m > 0.
31. Montrer qu’à une valeur de u(t) correspondent au plus 2 valeurs de v(t).
32. Expliquer ce que réalise le programme suivant :
function [u,v,H]=Voltera(a,u0,v0,T,h)
N=ceil(T/h);
S=zeros(N,1);
R=zeros(N,1);
H=zeros(N,1);
u(1)=u0;
v(1)=v0;
for k=2:N
u(k)=u(k-1)+h*u(k-1)*(1-v(k-1));
v(k)=v(k-1)+a*h*v(k-1)*(u(k-1)-1);
H(k)=a*v(k)+u(k)-a*log(u(k))-log(v(k));
end
%close all;
plot(u,v);
figure;plot(H);
A quel schéma d’approximation numérique correspond t-il ?
33. Faire un test avec a = 1, u0 = 0.4, v0 = 0.4, T = 100 et h = 0.01 et interpréter graphiquement
les courbes.
34. Ces courbes sont-elles cohérentes avec l’analyse théorique ? Pourquoi ? Tracer plusieurs
courbes avec les mêmes paramètres avec T = 20 en faisant varier h. Qu’observe t-on ?
35. Reprendre les valeurs précédentes en prenant v0 = 2.
36. Que se passe t-il si u0 = 1 et v0 = 1 ?
37. Tester également u0 = 0.95 et différentes valeurs de a ?
38. Modifier le programme de manière à approcher la solution réelle en utilisant la méthode
du point milieu.
39. Comparer les trajectoires et l’évolution de la valeur de H pour différentes valeurs de h avec
cette nouvelle méthode et la méthode précédente. Qu’en concluez vous ?

8
 8 Equations différentielles d’ordre 2.
L’exemple précédent sur le modèle de Volterra a permis de voir que les méthodes numé-
riques pouvaient être utilisés pour les systèmes différentiels d’ordre 1.
Cette remarque permet d’utiliser ces méthodes numériques pour résoudre certaines équa-
tions d’ordre 2. En effet si on doit résoudre une équation sclalaire du type

y(2) = f (y0 , y,t) avec y0 (t0 ) = z0 et y(t0 ) = y0 (7)

on peut écrire le système différentiel d’ordre 1 équivalent :

 0
z = f (z, y,t) avec z(t0 ) = z0
(8)
y0 = z avec y(t0 ) = 0.

40. Ecrire le système différentiel équivalent à l’équation suivante :

ty0
y(2) = −y+3 avec y(0) = 1 et y0 (0) = 0.
2

41. Après avoir remarqué que z(t) = t 2 + 1 est solution, résoudre numériquement ce système en
utilisant la méthode du point milieu.
Equation de la gravitation
Si on note Y le vecteur tridimensionnel de la position d’une planète dans un référentiel
héliocentrique (le soleil est au point (0, 0, 0)) la force de gravitation du soleil sur ce corps
Y
est égal à −GmMS , où m désigne la masse du corps, MS celle du soleil et G est la
kY k3
constante de gravitation. Nous rappelons que la relation fondamentale de la dynamique
assure que la somme des forces qui s’exercent sur un corps est égale à la masse multipliée
par l’accéleration de ce corps.
42. Montrer que le vecteur Y est solution d’une équation différentielle de la forme
−KY
Y (2) = .
kY k3

Pour simplifier les calculs on va se placer dans le plan, c’est à dire que Y = (x, y).
43. Ecrire un système différentiel d’ordre 1 dont x et y sont les solutions.
44. Compléter le programme suivant en utilisant la méthode d’Euler pour résoudre ce système
numériquement.

function UnCorpsEuler(x01,x02,v01,v02,h,T)
K=2;
z0=0;y0=0;
N=floor(T/h);
x1=zeros(N,1);
x2=x1;v1=x1;v2=x1;
x1(1)=x01;
x2(1)=x02;
v1(1)=v01;
v2(1)=v02;
for k=1:N

9
 x1(k+1)=x1(k)+
x2(k+1)=x2(k)+
v1(k+1)=v1(k)-
v2(k+1)=v2(k)-
end
V=400;
a=6;
c1=0.5*(z0+x01);
c2=0.5*(y0+x02);
for k=1:V:N
plot(x1(k),x2(k),’ro’);title(’Methode d Euler’);
hold on;
plot(x1(1:k),x2(1:k),’k’);
plot(y0,z0,’ko’);

axis([c1-a c1+a c2-a a+c2])

hold off;
pause(0.001);
end

45. A quoi correspondent les 4 données intitiales ?

46. Tester différentes vitesses initiales en utilisant des paramètres assurant que la planète fasse
plusieurs tours du soleil.
47. Proposer un programme analogue pour la méthode du point milieu et comparer les deux
méthodes.
48. Que se passe t-il si la vitesse initiale est trop importante ?

10