Equ Diff - Matl
Transféré par
Samir2020Equ Diff - Matl
Transféré par
Samir2020Partie 5 :
Équations différentielles
1 Exemples et motivations
Exemple 1 (Biologie) Considérons une population d’animaux dans un milieu ambiant où au plus
animaux peuvent coexister. On suppose que initialement la population soit et que le facteur de
croissance des animaux soit égal à une constante . Dans ce cas, l’évolution de la population au cours du
temps sera proportionnelle au nombre d’animaux existants, sans toutefois que ce nombre ne dépasse la
limite . Cela peut s’exprimer au moyen de l’équation
(1.1)
La résolution de cette équation permet de trouver l’évolution de la population au cours du temps.
On suppose maintenant que deux populations, "! et # , soient en compétition entre elles. L’évolution des
deux populations est alors décrite par le système d’équations différentielles
!
$! "!%
'& (*)!$"!%
,+-##
.
#
/# #
'& (*)0#%#
'*+1!$"!%
. (1.2)
! # 1
+ ! -
+ #
où et sont les facteurs de croissance des deux populations, et tiennent compte de l’interaction
entre les deux populations, tandis que ) ! et ) # sont liés à la quantité de nourriture disponible pour chaque
population. Ce système d’équations différentielles est connu comme modèle de Lotka-Volterra.
2 Introduction
On considère une fonction continue 2,3 465 7 . Pour 8 :9 donné, on cherche 8*3; 9
4<7=8?>@;-A 9 qui satisfait le problème suivant, appelé problème de Cauchy :
8 >@;-A?BC2?>D;FEG8?>@;-A%A si ;IHKJ
(2.1)
8>LJABM8
N
8?>@;-A
où 8 @> ;-AB N .
;
1
Exemple 2
O Un problème de Cauchy peut être linéaire, comme par exemple :
P% '*P% si (Q
(2.2)
pour lequel R
P *P% et dont la solution est
T (* U P
VXWXY[Z\ ]Z* U P .
O On peut avoir aussi desS problèmes non-linéaires, comme :
0 S
]_^
si
` (2.3)
avec R a b ^ . Ce problème admet les trois solutions suivantes : ac , a d% W U e f ,
] 0S %d W U e S f .
O Le problème suivant :
#
](Z\ si
` (2.4)
admet comme solution la fonction gFh%ij avec lk ?knm # , c’est-à-dire une solution locale.
Théorème 2.1 (de Cauchy-Lipschitz) Si 2 est continue sur 4o5 et s’il existe une constante
p
HqJ telle que r r r r
w p y 9 y
2?>D;FE/s'A]tc2?>D;FE/uvA s:txu s[E/u E ;zHqJ{E (2.5)
alors le problème de Cauchy (2.1) admet une solution globale (c-à-d. pour tout ;zH|J ) et elle est
r r r r r r
unique.
Dans l’exemple (2.2), 2>D;FE/s'At,2?>@;FE/u:A B >L}~st*}~;-AtM>L}~ut,};-A Bq} stu . Donc, la relation
p
(2.5) est satisfaite pour B} .
3 Dérivée numérique (Sec. 9.1 du livre)
!
Soit 23 E7 , de classe et vBc E/ ! E~]E/'vB , nœuds équirépartis dans E~- .
On note B>$atA% la distance entre deux nœuds consécutifs.
La dérivée 2 >@'$A est donnée par
2?>D'cAt2>D'$A
2 >D'$AB ¢¤£ ¡ -¥ E
2>D'$At2>D't§A
B ¢¤£ ¡ ~¦ E
2?>D'cA`tc2?>D't§A
B ¢¤£ ¨
Soit, maintenant, >L©\2jA0 une approximation de 2 >D'$A . On appelle
2
ª différence finie progressive l’approximation
2?>@'4 ! At§2?>@'0A
>$©\2[A¬« B E ®BqJ{E~¯E/<t° (3.1)
ª différence finie rétrograde l’approximation
2?>@'$At2?>@'² ! A
>$©\2jA0± B E ®B³~E~¯E/ (3.2)
ª différence finie centrée l’approximation
2?>@'4 ! At2?>@'² ! A
>$©\2[A¬´ B ¨ E ®B~E~]E/<tC (3.3)
#
Si 2 9 > A , pour tout 9 , il existe un µ entre ' et tel que l’on a le développement de
Taylor
2 D> µ[A #
2>DA?BC2?>@'$A]c2 >@'0A >@\t¶'LA ¨ >@·tx'0A (3.4)
ª Pour Bc'4 ! dans (3.4) on obtient pour la différence finie progressive
>$©o2jA « B2 >D'$A] ¨ 2 D> µ[AE (3.5)
r type
ce qui conduit à une estimation du r
w
r r 2 >D'LAt°>$©\2[A¬« ¸E
!
où B # º¹»¼½¤¾ ¼¿¡À ¼¿ ¥ÁD 2 >@A .
ª Pour Bc'² ! dans (3.4) on obtient pour la différence finie rétrograde
>L©\2jA ± B2 >D'0At ¨ 2 >Dµ[AE
r r
ce qui conduit à
r r 2 >D'$At°>$©\2jA0± w ¸ E
!
où B # º¹» ¼ ½¤¾ ¼¿ ¦ Á À ¼¿  2 >DA .
ª Pour Bc'4 ! et ÃBc'² ! avec un développement d’ordre ¨ (si 2 9 W )
2 D> 'LA # 2 @> µ ! A W
A B2?>D'LA]§2 >D'0AÄ
2?>D'4 ! ¨ Å E
2 >D'0A # 2 >Dµ # A W
2?>D'² ! AB2?>D'LAt2 >@'$A ¨ t Å E
on obtient
2 D> µ ! A§2 Æ @> µ # A #
>L©\2jA ´ B2 >D'$A] ¨
et donc on a l’estimation suivante r r
w #
r 2 r >D'LAt°>$©\2[A ´ ¸ E
!
où B Ç º¹» ¼½¤¾ ¼¿ ¦ Á À ¼¿ ¥Á@ 2 Æ >DA .
3
r r
Définition 3.1 La différence 2 >D'0AtC>L©\2jA « (et celles correspondantes aux autres différences
finies) est appelée erreur de troncature au point ' .
¨
L’erreur de troncature est d’ordre pour les formules progressive et rétrograde et d’ordre pour
la formule centrée.
4 Méthodes de différences finies (Sec. 10.2 du livre)
Soient JBT; ÉÈ ; !ÊÈ ~ È ;G È ;G4 !(È une suite de nombres réels. On note :BM;G4 ! tË;G .
On notera par :
Ì une approximation de 8>D;G{A
Dans le problème de Cauchy (2.1), pour ;ÍBM;G on a
8 >D;GA?BC2?>D;GE/8>D;G{A%A%
On veut alors également approcher la dérivée 8 D> ;GA au point ;G . Cela se fait en utilisant des
schémas de dérivation numérique.
Schéma d’Euler progressif :
Ì 4 ! t Ì ¨
B2>D;GE Ì A pour ,BJ{E~E
(4.1)
Ì BT8
Schéma d’Euler rétrograde :
Ì 4 ! t Ì ¨
BC2>D;G4 ! E Ì 4 ! A pour *BJ{E~E
(4.2)
Ì Bc8
Remarque 4.1
ª Le schéma d’Euler progressif est un schéma explicite, car il permet de calculer Ì Î4 ! en
fonction de Ì :
Ì 4 ! B Ì Ïc"2?>@;G"E Ì {A%
ª Le schéma d’Euler rétrograde est un schéma implicite, car il ne permet pas de calculer
Ì 4 ! en fonction de Ì :
Ì 4 ! B Ì Ïc"2?>@;G4 ! E Ì Î4 ! A%E
(sauf pour des cas particuliers de fonctions 2 , par exemple lorsque 2 est linéaire par
rapport à sa deuxième variable).
4
Exemple 3 On considère l’équation différentielle suivante
#
Q 0
?
e (4.3)
On résout cette équation avec les méthodes d’Euler progressive et Euler rétrograde, dans l’intervalle & .
0Ð
avec un pas de temps Ñ Ò @e (ce qui équivaut à approcher la solution exacte
aux instants ·
@e @Ó ).
Ð
Euler progressive
>> f=’-t*y.ˆ2’;
>> h=0.2; t0=0; tf=4; y0=2;
>> [t_EP, y_EP] = euler_prog(t0,y0,tf,h,f)
La variable t EP contient la suite des instants et la variable y EP les valeurs Ô correspondantes
calculées par la méthode.
Euler rétrograde
La fonction euler retr utilise la méthode de Newton pour résoudre à chaque pas de temps
l’équation non-linéaire qu’on obtient en appliquant le schéma d’Euler rétrograde. C’est pourquoi il
faut donner aussi la fonction Õ-Ö .
Õ-×
>> df=’-2*t*y’;
>> [t_ER, y_ER] = euler_retr(t0,y0,tf,h,f,df)
#
La solution exacte du problème (4.3) est
e U ÊZ\
. La figure 1 obtenue grâce aux commandes
>> t=[0:0.05:4]; sol_ex=’2./(1+t.ˆ2)’;
>> plot(t,eval(sol_ex),’b’,t_EP,y_EP,’ro-’,t_ER, y_ER’,’go-’)
>> legend(’sol-ex’,’Euler prog.’,’Euler retr’)
montre la solution exacte ainsi que les deux solutions calculées.
5 La propriété de stabilité (Sec. 10.3.3 du livre)
Pour Ø·H|J donné, on considère le problème :
8 D> ;-A?BCt¸Ø{8?>@;-A si ;zHKJ
(5.1)
8?>ÙJA?BM8
dont la solution est
²Û Y
8?>@;-A?BÚ 8 E ¡£ Ü 8>D;-A?BKJ{
Yl
donc
y
Posons JÝBn; ÄÈ ; !zÈ È ;G È ;GÎ4 !zÈ ~ tels que ;GBn? , c-à-d. oBÞ , où <HÞJ est
donné et s’appelle pas de temps.
5
y ’ (t) = − t y2
2
1.8 sol−ex
Euler prog.
1.6 Euler retr
1.4
1.2
y(t)
0.8
0.6
0.4
0.2
0
0 0.5 1 1.5 2 2.5 3 3.5 4
t
F IG . 1: Comparaison entre la solution exacte et celles obtenues par les méthodes d’Euler pro-
gressive et rétrograde.
6
Définition 5.1 On dira qu’une méthode de différences finies est stable si, pour suffisamment
petit, on a que
Ì :7 J lorsque ;GÝ7 ß
ª Pour le schéma d’Euler progressif :
Ì 4 ! ³ Ì vB>/(tØ'A 8 E y
B r >/rtØ[A Ì E d’où xàJ{ (5.2)
Si ?tØ' È t alors Ì 7 ß quand x7 ß , donc le schéma d’Euler progressif est instable.
Pour assurer la stabilité, on a besoin de limiter le pas de temps , en imposant la condition de
r r
stabilité :
¨
?tØ' È d’où È Ø[
ª Pour le schéma d’Euler rétrograde :
Ì Î4 ! B Ì Ì vB y
et donc 8 E xàJ{
§Ø' §Ø'
Et comme ¡ ¸ £lÜ Ì BJ , on a la stabilité sans condition sur .
Exemple 4 On résout le problème (5.1) pour á c e et § dans l’intervalle & %. avec les méthodes
@â
d’Euler progressive et rétrograde et Ñ et Ñ @ . Voilà les commandes Matlab pour le cas Ñ @â :
>> f=’-2*y’; h=0.9; t0=0; tf=10; y0=1;
>> [t_ep, y_ep] = euler_prog(t0,y0,tf,h,f);
>> df=’-2’; [t_er, y_er] = euler_retr(t0,y0,tf,h,f,df);
>> t=[0:0.1:tf]; sol_ex=’exp(-2*t)’;
>> plot(t,eval(sol_ex),’b’,t_ep,y_ep,’ro-’,t_er, y_er’,’go-’)
La figure 2 montre les solutions obtenues pour Ñ ã @â (à gauche) et Ñ ã @ (à droite) ainsi que la
solution exacte.
Exemple 3 (suite) Utilisons un pas de temps beaucoup plus grand : Ñ ³ @d . La figure 3 montre les
solutions obtenues par les deux méthodes d’Euler progressive et rétrograde. On remarque que la solution
obtenue par la méthode d’Euler progressive est devenue instable.
6 Convergence de la méthode d’Euler progressive (Sec. 10.3.2 du
livre)
Dans (2.1), pour ;ÍBM;G on a
8 >D;GA?BC2?>D;GE/8>D;G{A%A% (6.1)
7
y ’ (t) = −2 y(t)
y ’ (t) = −2 y(t) 6
1
sol−ex sol−ex
Euler prog. Euler prog.
0.8 4 Euler retr
Euler retr
0.6
2
0.4
0.2
y(t)
0
y(t)
−2
−0.2
−0.4
−4
−0.6
−0.8 −6
0 2 4 6 8 10 0 2 4 6 8 10
t t
F IG . 2: Comparaison entre les solutions obtenues par les méthodes d’Euler progressive et
rétrograde pour oBqJ{Dä (à gauche) et ËB³~$ (à droite).
y ’ (t) = − t y2
5
0
y(t)
−5
sol−ex
Euler prog.
Euler retr
−10
0 0.5 1 1.5 2 2.5 3 3.5 4
t
F IG . 3: Comparaison entre les deux solutions obtenues par les méthodes d’Euler progressive et
rétrograde pour oBqJ{Då .
8
On sait qu’il existe µ entre ;G et ;G4 ! BM;Gɧ (voir formule 3.5) tel que
#
8?>@;GA?Bc8>D;G{A]§8 >@;GA] ¨ 8 @> µ{A
On en déduit que
8>D;G4 ! A tx8?>@;G{A
8 >D;GA(B t ¨ 8 @> µ{A
En conséquence, (6.1) devient :
8>D;GÎ4 ! At¶8>D;G{A
B2?>@;G"E/8>D;G{A/A] ¨ 8 D> µ{A (6.2)
Pour le problème (2.1), la méthode d’Euler progressive avec un pas de temps constant , s’écrit :
Ì 4 ! t Ì ¨
B2>D;GE Ì A pour ,BJ{E~E
(6.3)
Ì BT8
En imposant à la solution exacte de vérifier la méthode numérique, on introduit donc une erreur
de troncature. Plus précisément, on appelle erreur locale de troncature la fonction æÎ4 ! >$A telle
que
8?>@;G4 ! At¶8>D;G{A
BC2?>@;G"EG8?>@;GA/A]xæ~4 ! >LA (6.4)
Donc :
æ~4 ! >$A(B ¨ 8 @> µ{A
On va démontrer le résultat de convergence suivant :
# p
Théorème 6.1 Si 8 9èç > 4[A et 2 est lipschitzienne,
r r étant sa constante de Lipschitz, alors
y Ì w
xàJ{E 8?>@;GAt >D;GAE (6.5)
où r r
Ú é Y
ê ° t
o>@;GA?B ¨ p º¹» 8 @> ;-A
Y
On dit qu’une méthode est convergente dès que r r
y Ì
xàCJ{E £l
Ü 8?>@;GAt BJ{E
Démonstration. On définit l’erreur au pas par :
ÚvBc8>D;G{At Ì
9
De (6.3) et (6.4) on obtient :
ÚÎ4 ! t Ú
Bq2?>D;GE/8>D;G{A%At2>D;GE Ì A],æÎ4 ! >$AE
(6.6)
Ú BJ{
r
Comme la fonction 2 est lipschitzienne, on a r r r
w p
2?>@;GE/8?>@;G{A%At2?>@;G"E Ì {A Ú
Compte tenu de cette inégalité,r (6.6)r permet d’écrire
r r quer : r
w p
ÚÎ4 ! >% A Ú § æ~4 ! >LA
On définit l’erreur globale de troncature æ?>$A par :r r
r r æ?>$A(B º ¹» æ[>$A
Soit ëÊìÊB Úì . Alors on a les inégalités suivantes :
w p
ëÊÎ4 ! >/§ Aëɧæ(>LA
w p p
>/§ A ¡>/?c A%ëÊβ ! cæ(>$A/{§æ(>LA
# p Î4 !
w p p p >/a§ A t°
C>/§ A >%§ A §~1>%§ A æ(>LA(B p æ?>$A
(6.7)
p w ¢ y p w Yê ¢
Mais ?c Ú é E ·HKJ , d’où /> ?c A ° t Ú é t°I
r BCÚ é r t° . Donc
Y
ê t°
Y ê tC
w Ú é w Ú é
ëÊ p æ >$A
? p ¨ º ¹» Ì >@;-A
Y
Remarque 6.1 Le même type de résultat peut être établi pour la méthode d’Euler rétrograde.
7 Méthodes de Runge-Kutta d’ordre 2
Si on intègre l’équation 8 >D;-A(BC2?>@;FE/8?>@;-A%A entre ;G et ;GÎ4 ! on obtient :
Yê
¥ Á
N
8>D;G4 !
A t¶8?>@;G{AB 2?>D;FE/8>D;-A/A F; (7.1)
Yê
En utilisant la formule du trapèze, on trouve le schéma implicite suivant, appelé schéma de
Crank-Nicolson ou du trapèze :
Ì Î4 ! t Ì B y
¨ 2?>D;GE Ì {A]§2?>@;G4 ! E Ì Î4 ! A%~E xàJ{ (7.2)
En modifiant le schéma (7.2) afin de le rendre explicite, on identifie la méthode de Heun :
Ì Î4 ! t Ì vB
¨ 2>D;GE Ì A]c2?>D;GÎ4 ! E Ì Ïc2?>@;G"E Ì {A%A/ (7.3)
10
Les deux méthodes (de Crank-Nicolson et de Heun) sont d’ordre 2 (par rapport à oB º¹»íI ).
Si on utilise dans (7.1) la méthode du rectangle on trouve
Ì Î4 ! t Ì vBqI2?>D; Á Ì Á
Î4 î E 4 î A (7.4)
Si maintenant on approche Ì 4 !
ï¬# par
Ì Ì
Á
4 î B ¨ 2?>@;G"E Ì {AE
on trouve la méthode d’Euler modifiée :
Ì 4 ! t Ì vBCI2
; 4 îÁ E Ì ¨ 2?>@;G"E Ì {A (7.5)
Les méthodes de Heun et d’Euler modifiée sont des cas particuliers dans la famille de méthodes
de Runge-Kutta d’ordre 2.
Il existe d’autres méthodes plus compliquées, comme par exemple la méthode de Runge-Kutta
d’ordre 4 suivante :
ð !
BC2?>D;GE Ì {A
ð # Ì ð !
BC2?>@;GÏ ¨ E ¨ EA
Ì :7 ð ð #
W CB 2?>@;GÏ ¨ E Ì ¨ E A
ð¸ñ ð
BC2>D;G4 ! E Ì § W A
Ì 4 ! B Ì Ï ð ! ¨ ð # ¨ ð ð¸ñ
Å > W A
Exemple 5 On considère le problème de Cauchy
T @\òó i[ ?
(7.6)
On a résolu ce problème par les méthodes d’Euler progressive et de Heun sur l’intervalle & e%. avec un
pas de temps Ñ .
Ð
>> f=’(cos(t) - 0.1)*y’;
>> h=0.4; t0=0; tf=12; y0=1;
>> % euler progressif
>> [t_ep, y_ep] = euler_prog(t0,y0,tf,h,f)
>> % heun
>> [t_heun, y_heun] = heun(t0,y0,tf,h,f)
Lafigure
o
ô ! 4Îõ öÆ÷ø 4 montre les solutions obtenues par les deux méthodes ainsi que la solution exacte
V ² Y YÙù . On remarque que la solution obtenue par la méthode de Heun est beaucoup plus précise que
celle d’Euler progressive.
D’ailleurs, si on réduit le pas de temps, la solution obtenue par la méthode d’Euler progressive s’approche
de la solution exacte. La figure 5 montre les solutions obtenues avec Ñ @e @ @ ú .
Ð
11
y ’ (t) = − y ( 0.1 − sin (t) )
2.5
sol−ex
Euler prog.
Heun
2
1.5
y(t)
0.5
0
0 2 4 6 8 10 12
t
F IG . 4: Comparaison entre les solutions obtenues par les méthodes d’Euler progressive et de
Heun pour BJ{Dû .
>> plot(t,eval(sol_ex),’b--’); hold on;
>> h=0.4;
>> for i=1:4
[t_ep, y_ep] = euler_prog(t0,y0,tf,h,f);
plot(t_ep,y_ep)
h=h/2;
end
On veut, maintenant, estimer l’ordre de convergence de ces deux méthodes sur cet exemple. Pour cela, on
résout le problème avec différents pas de temps et on compare les résultats obtenus à l’instant (cÓ avec
la solution exacte.
>> h=0.4; t=6; y6=eval(sol_ex)
>> for i=1:5
% euler progressif
[t_ep, y_ep] = euler_prog(t0,y0,tf,h,f);
err_ep(i) = abs(y6 - y_ep(fix(6/h)+1));
% Heun
[t_heun, y_heun] = heun(t0,y0,tf,h,f);
err_heun(i) = abs(y6 - y_heun(fix(6/h)+1));
h=h/2;
end
>> h=[0.4, 0.2, 0.1, 0.05, 0.025];
12
Euler progressive
y ’ (t) = − y ( 0.1 t − sin (t) )
2.5
sol. ex.
1.5
y(t)
h=0.05
h=0.1
1
h=0.2
h=0.4
0.5
0
0 2 4 6 8 10 12
t
F IG . 5: Solutions obtenues par la méthode d’Euler progressive pour différents pas de temps.
0
10
−1
10
−1
−2
| y(6) − un (6) |
10
−3
10
−4
10
−2 euler prog.
heun
−5
10 −2 −1 0
10 10 10
h
F IG . 6: Erreurs en échelle logarithmique commises par les méthodes d’Euler progressive et de
Å
Heun dans le calcul de 8?> A .
13
>> loglog(h,err_ep,’b’,h,err_heun,’r’)
La figure 6 montre en échelle logarithmique les erreurs commises par les deux méthodes en fonction de Ñ .
On voit bien que la méthode d’Euler progressive converge à l’ordre 1 tandis que celle de Heun à l’ordre 2.
8 Applications
On revient à l’exemple proposé au début du chapitre.
Exemple 1 (suite) On considère d’abord l’équation scalaire (1.1). Prenons une population initiale de
Ð
lapins dont le facteur de croissance est c @ d (l’unité de temps est 1 mois) et la population maximale
de
f lapins. On résout l’équation par la méthode de Heun avec Ñ T
mois sur une période de trois
ans :
>> f=’0.08*y.*(1-(y/70))’; h=1; t0=0; tf=36; y0=40;
>> [t, y] = heun(t0,y0,tf,h,f); plot(t,y)
La solution est montrée en figure 7.
70
65
60
population de lapins
55
50
45
40
0 6 12 18 24 30 36
mois
F IG . 7: Solution de l’équation scalaire de Lotka-Volterra.
Dans Matlab il y a plusieurs algorithmes déjà implémentés pour résoudre des équations différentielles
(voir par exemple ode45, ode23, ode23s, ode15s). En particulier, ode45 résout l’équation par une
méthode d’ordre 4 similaire à la méthode de Runge-Kutta d’ordre 4 ; en plus, il utilise un pas de temps
variable et il choisit à chaque instant le pas de temps le plus convenable de façon à satisfaire une tolérance
fixée. Pour utiliser cette commande, il faut d’abord écrire une fonction Matlab qui contient la définition
de la fonction R
:
function dy=fun(t,y)
dy = 0.08*y.*(1-(y/70));
return
14
Ensuite, on tape les commandes
>> options=odeset(’RelTol’,1e-4);
>> [t,y] = ode45(’fun’,[t0,tf], y0, options);
ñ
La première commande permet de fixer des options ; dans notre cas on a défini une tolérance relative de
² (taper help odeset pour savoir toutes les options qu’on peut définir). La deuxième commande
permet de résoudre l’équation différentielle définie dans la fonction fun, dans l’intervalle [t0, tf] et
avec comme condition initiale y0.
On considère, maintenant, le système (1.2). Prenons une population initiale !
de 40 lapins et une
population #
de 20 renards et le système de Lotka-Volterra
! @ dí"! * @ "!- #
# l @ Óí#1 Z* @ Ð eí"!- # (8.1)
On souhaite étoudier l’évolution des deux populations sur une période de 10 ans. Si on introduit les
vecteurs
"!%
ý @ dí"!
* @ "!-
#
ü]
ü
#
l @ Óí#
Z* @ Ð eí"!
#
on peut écrire le système (8.1) sous la forme très générale :
ü ý ]ü ü] `K& "!% # .þ (8.2)
Toutes les méthodes qu’on a vues jusqu’à maintenant sont applicables aussi au système (8.2). Par exemple,
la méthode d’Euler progressive s’écrit
ÿ 4 ! ÿ
ý ÿ
Ñ
ce qui équivaut au schéma
Î4 !
Ô ! Ô ! Q @ d Ô ! * @ Ô !
Ô #
Ñ! Ð
# Î4 #
Ô
Ñ
Ô @ Ó Ô # *
Q
Z @ e
Ô ! Ô #
Ô ! !
Ô # #
La commande Matlab ode45 permet de résoudre aussi des systèmes d’équations différentielles. Il faut
d’abord écrire une fonction qui définisse le système :
function dy=fun2(t,y)
dy=zeros(2,1);
dy(1) = 0.08*y(1) - 0.004*y(1)*y(2);
dy(2) = -0.06*y(2) + 0.002*y(1)*y(2);
return
Ensuite on peut résoudre le système :
15
>> y0=[40;20]; t0=0; tf=120;
>> options=odeset(’RelTol’,1e-4);
>> [t,y] = ode45(’fun2’,[t0,tf], y0, options);
>> plot(t,y(:,1),’b’, t,y(:,2),’r’)
La première colonne de y contient la solution "! tandis que la deuxième colonne contient # . La figure
8 montre l’évolution des deux populations. On remarque que la distribution des instants de discrétisation
n’est pas uniforme.
40
lapins
renard
35
30
25
20
15
0 20 40 60 80 100 120
F IG . 8: Évolution des populations de lapins et de renards en 10 ans.
16