Systèmes Dynamiques

Francine Diener et Marc Diener

3 septembre 2017
2
Table des matières

1 Le modèle de croissance logistique 5

1.1 Croissance d’une population d’éléphants . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Le modèle logistique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Calcul des solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 Tracé de trajectoires avec Scilab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2 Equations différentielles en dimension 1 9

2.1 Définitions et premiers exemples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Le cas linéaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Etude qualitative dans le cas autonome . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

3 Solutions numériques approchées 13

3.1 Champs de vecteurs associé à une équation différentielle . . . . . . . . . . . . . . . . . . . 13
3.2 Solutions numériques approchées . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2.1 Algorithme de la méthode d’Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2.2 Evaluation de l’erreur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2.3 Algorithme de la méthode RK2 ou méthode du point milieu . . . . . . . . . . . . . 15
3.2.4 Méthode Runge-Kutta d’ordre 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

4 Existence et unicité des solutions 17

4.1 Enoncé du théorème fondamental . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.1.1 Unicité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.1.2 Cas non autonome . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.1.3 Durée de vie et explosion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.2 Dépendance par rapport à la condition initiale et à un paramètre . . . . . . . . . . . . . . 18
4.2.1 Déterminisme et chaos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

5 Le système de Lotka-Volterra 19
5.1 Un exemple en écologie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
5.1.1 Le modèle proie-prédateur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
5.1.2 Présence de solutions périodiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
5.1.3 Comment lire ce comportement directement sur l’équation ? . . . . . . . . . . . . . 20
5.1.4 Que s’est-il passé dans l’adriatique ? . . . . . . . . . . . . . . . . . . . . . . . . . . 20
5.2 Une théorie du cycle économique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.2.1 Introduction des variables pertinentes . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.2.2 Construction du modèle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.2.3 Les cycles économiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

6 Noeuds, cols, foyers et centres 23

6.1 Systèmes de deux équations différentielles . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
6.1.1 Trajectoires et champs de vecteurs . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
6.1.2 Etude qualitative, isoclines, équilibres . . . . . . . . . . . . . . . . . . . . . . . . . 24
6.1.3 Calcul de solutions approchées . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
6.2 La classification de Poincaré . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
6.2.1 Calcul des solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3
4 TABLE DES MATIÈRES

7 Linéarisé d’un système différentiel non linéaire 27

7.1 Linéarisé au voisinage d’un équilibre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
7.2 Exemples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
7.3 Le cas où le linéarisé est un centre : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

8 Dynamiques de populations structurées en âges 31

8.1 Le modèle de Leslie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
8.2 Le théorème de Perron Frobenius . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

9 Systèmes différentiels en dimension n 35

9.1 Systèmes dynamiques à n dimensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
9.2 Le cas linéaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
9.3 Stabilité locale dans le cas non linéaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
9.4 Loi de conservation ou intégrale première : . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

10 L’équation de Van der Pol 39

10.1 L’oscillateur linéaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
10.2 Le cycle unique de l’équation de Van der Pol . . . . . . . . . . . . . . . . . . . . . . . . . 39
10.3 Une bifurcation de Hopf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
Chapitre 1

Le modèle de croissance logistique

Dans ce chapitre on se propose de découvrir sur un exemple très simple qui modélise une situation
concrête, les notions d’équation différentielle, de solutions et de conditions initiales. On verra aussi que,
même s’il est possible dans ce cas simple de calculer explicitement l’ensemble des solutions, on pourrait
néanmoins étudier l’équation et ses solutions même si l’on ne pouvait pas les calculer.

1.1 Croissance d’une population d’éléphants

L’éléphant africain de la savane (loxodonta africana) se comptait par millions dans la savane africaine
avant qu’il ne soit décimé, durant des siècles, par des chasseurs, notamment pour exploiter l’ivoire de ses
défenses et prendre possession de ses territoires à des fins agricoles. A la fin du 19e siècle, cette population
étant pratiquement arrivée à extinction en Afrique du Sud, il fut décidé la création d’un parc naturel, le
parc Kruger à la frontière entre l’Afrique du Sud et le Mozambique. Le premier responsable du parc en
1903 ne trouva aucun éléphant à son arrivée mais un petit groupe de 10 éléphants furent repérés en 1905,
vraissemblablement venu du Mozambique. Des mesures de protection strictes, à la fois des animaux et de
leur habitat furent décidées dans ce parc et maintenues tout au long du 20e siècle. Elles permirent une
croissance naturelle de cette population, qui fut d’abord lente jusque dans les années 30, puis très rapide
jusque dans les années 60. C’est alors qu’on observa à la fois un ralentissement du taux de croissance et,
en même temps, un début de dégradation par les éléphants d’autres espèces de l’écosystème comme les
baobabs par exemple.

Figure 1.1 – Effectifs observés et effectifs théoriques de la population d’éléphants dans le parc Kruger.

Pour décider de l’attitude à adopter pour gérer au mieux les populations de ce parc, les responsables
eurent recours à un modèle mathématique appelé modèle logistique.
Le tableau suivant indique les effectifs observés y(t) ainsi que les effectifs théoriques y(t) calculés en
suivant ce modèle (et arrondis à l’entier le plus proche).

5
6 CHAPITRE 1. LE MODÈLE DE CROISSANCE LOGISTIQUE

t 1905 1923 1930 1939 1945 1950 1960 1970 1980 1990 2000
y(t) 10 13 29 450 980 3010 5800 6500 7400 7200 7310
y(t) 10 146 402 1346 2623 3994 6271 7186 7428 7484 7496
Cela permit de déterminer la valeur d’une taille limite, ici K = 7500, qui représente la taille de la
population en deça de laquelle il convient de rester si l’on veut préserver la cohabitation harmonieuse de
la population avec son écosystème. Le parc mit alors en place un programme d’abattage controlé destiné
à limiter la surpopulation en maintenant le nombre d’éléphants approximativement égal à cette valeur.

1.2 Le modèle logistique

La courbe représentant les effectifs théoriques est appelée une courbe logistique en raison de sa forme
particulière qui présente cette caractéristique de croissance amortie. Elle est en réalité le graphe de la
solution de l’équation différentielle  
′ y(t)
y (t) = ry(t) 1 − (1.1)
K
dite équation logistique de condition initiale y(t0 ) = 10, avec t0 = 1905. Cette équation comporte deux
paramètres, r = 0, 15 qui s’appelle le taux de croissance intrinsèque et K = 7500 qu’on appelle la capacité
biotique. L’idée du modèle logistique, introduit par Verhulst en 1836, est la suivante. Si la population
concernée pouvait croı̂tre indéfiniment, sans rencontrer aucune limitation de ressource ou d’espace, elle
aurait une croissance exponentielle. En effet elle serait solution de l’équation différentielle suivante
y ′ (t) = ry(t) (1.2)
que l’on appelle, comme dans le cas des modèles à temps discrets, un modèle de Malthus (en prenant
pour taux de croissance r = 0, 15). Or la solution de cette équation différentielle, de condition initiale
y(t0 ) = 10, s’écrit y(t) = y(t0 )er(t−t0 ) = 10e0,15(t−1905) . On peut le vérifier en remplaçant y(t) par cette
fonction dans l’équation différentielle (1.2). On obtiendrait alors pour y(t) les valeurs suivantes :

t 1905 1923 1930 1939 1945 1950 1960 1970 1980 1990 2000
10e0,15(t−1905) 10 149 425 1640 4034 8541 38276 171542 768799 3445519 15441745
La taille de la population dépasserait alors les 15 millions dès la fin de la période considérée et cela
ne s’arrêterait pas là. Une telle croissance exponentielle n’est donc pas adaptée aux effectifs observés
ici, à l’exception peut-être de la période située avant 1940. En effet, pour l’équation (1.2), le taux de
′
croissance yy(t)(t)
reste constant (égal à r) au cours du temps et ceci ne tient pas compte des limitations
environnementales qui, de fait, ralentissent la croissance lorsqu’on s’approche de la capacité biotique K.
D’où l’idée de remplacer ce taux constant r par un taux variable r(1 − y(t) K ) qui dépend de la taille de la
y
population. On voit en effet que le coefficient 1 − K reste proche de 1 lorsque la taille de la population
est très petite, ce qui explique le début de croissance exonentielle, puis il diminue jusqu’à tendre vers
0 lorsque la taille de la population augmente et tend vers la capacité biotique. En fait ce coefficient
représente la part de la capacité biotique encore disponible à chaque instant t. Plus cette part s’amenuise
et plus la croissance se ralentit.

1.3 Calcul des solutions

y
On peut calculer explicitement les solutions de l’équation différentielle logistique y ′ = ry(1 − K) de
la façon suivante :  
dy(t) y(t)
= ry(t) 1 −
dt K
peut se réécrire
1
  dy(t) = rdt
y(t)
y(t) 1 − K
1
1 1
Mais comme on a l’égalité y
y(1− K )
= y + K
y
1− K
, l’équation devient
1
dy(t) dy(t)
+ K y(t) = rdt
y(t) 1− K
1.3. CALCUL DES SOLUTIONS 7

d’où, en intégrant,
y(t)
ln y(t) − ln(1 − ) = rt + C ste
K
soit encore en prenant l’exponentielle
y(t) ste

y(t)
= ert eC
1− K
 
y(0)K
Il est facile de vérifier que la constante d’intégration vaut ici C ste = ln K−y(0) . D’où, après simplifica-
tions, la solution
y(0)K
y(t) = . (1.3)
y(0) + e−rt (K − y(0))
Pour savoir qu’elle est l’allure du graphe de cette solution (en réalité il y a autant de solutions que
de choix de valeurs initiale y(0)), on pourrait calculer sa dérivée (ce qui serait légèrement fastidieux...)
mais il est bien plus simple d’utiliser l’équation différentielle : en effet, l’équation différentielle donne la
dérivée y ′ comme une fonction de y, puisque y ′ = ry(1 − y/K). On voit donc, sans calcul, que
— y ′ s’annule lorsque y = 0 et y = K
— y ′ > 0 lorsque y est compris entre 0 et K
— y ′ < 0 sinon.

Figure 1.2 – La parabole qui s’annule en y = 0 et y = K, graphe de la fonction y ′ = f (y).

Il en résulte que, aussi longtemps que la population reste inférieure à sa capacité biotique K, elle
ne cesse de croı̂tre (puisque y ′ > 0). Et on calcule facilement la limite, quand t tend vers l’infini, de
l’expression trouvée pour la solution y(t), qui vaut précisément K. A l’inverse, si y(t) est supérieure à
cette capacité, y(t) décroı̂t (puisque y ′ < 0) et on vérifie facilement qu’elle tend également vers K.
L’examen du graphe de la fonction f (y) = ry(1 − y/K), qui représente la dérivée de y, renseigne
aussi sur le taux de croissance maximal d’une population soumise à une croissance logistique. En effet,
la maximum de f est atteint pour y = 12 K, ce qui signifie que c’est lorsque la taille de la population est
égale à la moitié de sa capacité biotique que sa croissance est la plus forte.

ModŁle de logistique

2.5

1.5

0.5

0 5 10 15 20
x

Figure 1.3 – Quatres solutions d’une dynamique logistique : elles sont croissantes si la taille de la
population est inférieure à sa capacité biotique et décroissantes si elle est supérieure. La capacité biotique
est un équilibre de la dynamique (tout comme la solution nulle). Chaque vecteur représenté au point de
coordonnées (t, y) la direction tangente à la dynamique : sa pente est donnée par la valeur de f(y) en ce
point.
8 CHAPITRE 1. LE MODÈLE DE CROISSANCE LOGISTIQUE

1.4 Tracé de trajectoires avec Scilab

Si l’on se donne une ”condition initiale” (t0 , y0 ), l’instruction Scilab
ode(y0,t0,t,f)
donne la valeur (approchée avec 8 chiffres significatifs) de la valeur de la solution à l’instant t de l’équation
y ′ = f (t, y) telle que y(t0 ) = y0 , pour t0 =t0 et y0 =y0. Si on remplace t par une liste de valeurs telle
tt=0:0.1:20, tous les instants compris entre 0 et 20 par pas de 0.1, c’est-à-dire 0, 0.1, 0.2, 0.3, . . . , 19.8, 19.9, 20,
la commande ode(y0,t0,tt,f) retourne la liste des valeurs de la solution pour toutes les valeur de t
dans la liste tt.
Notons que, pour ode, la fonction f qui définit l’équation différentielle doit être une fonction de deux
variables (ici t et y), par exemple f (t, y) = −4 ∗ y + 8.5 ∗ cos(t) et il faut veiller à la définir comme telle,
même si la dépendance en t n’existe pas.
Voici comment calculer, puis tracer le graphe de la solution de cette équation telle que y(0) = 5 dans
une fenêtre ”numéro 1”, pour t ∈ [0..2] :
t0=0 ;y0=5 ;tt=0:0.1:2; yy=ode(y0,t0,tt,fonct); xset("window",1); plot(tt,yy);
Pour faire un autre tracé dans une nouvelle fenêtre, on appelle d’abord l’ouverture d’une deuxième
fenêtre graphique avec xset("window",2) ;). Si on ne le fait pas, le tracé suivant se superpose au
précédent dans la dernière fenêtre ouverte.
On verra au chapitre suivant comment tracer le champs des vecteurs tangents à la dynamique comme
sur la figure précédente avec la commande Scilab fchamp.
Chapitre 2

Equations différentielles en
dimension 1

Pour modéliser une quantité qui évolue au cours du temps il est naturel de postuler une relation entre
cette quantité et sa dérivée : on obtient alors une équation différentielle. C’est l’exemple le plus simple
de système dynamique.

2.1 Définitions et premiers exemples

Considérons une quantité y(t) (taille d’une population, concentration d’une substance, ...) qui évolue
au cours du temps et sa dérivée y ′ (t) (lorsqu’il est raisonnable de supposer que cette dérivée existe).
Supposons qu’on soit conduit à postuler une relation entre cette quantité et sa dérivée de la forme

dy(t)
= f (t, y(t))
dt
pour une fonction f particulière. Cette relation est une équation différentielle du premier ordre 1 et la
résolution d’une telle équation consiste à trouver toutes les fonctions y(t) inconnues qui satisfont cette
équation. Voici deux exemples, le modèle exponentiel et le modèle logistique(déja rencontré au chapitre
précédent.
Exemple : Le modèle exponentiel, très rudimentaire, a été proposé pour représenter la croissance d’une
population par Thomas Malthus en 1798. Il suppose que la population possède un taux de reproduction r
constant, simple différence du taux de natalité et du taux de mortalité (la population est supposée isolée
c’est-à-dire qu’aucune migration n’est envisagée). Si y(t) désigne la taille de la population à l’instant t et
y ′ (t) sa dérivée, la formule y ′ (t) = ry(t) signifie que le taux de croissance y(t+δt)−y(t)
δt entre les instants t
et t + δt est proportionnel à y(t) et surtout que le coefficient de proportionnalité r ne varie pas au cours
du temps. On peut résoudre cette équation : sa solution est donnée par y(t) = y(0)ert où y(0) désigne
la taille de la population à l’instant t = 0 qu’on appelle condition initiale. Ce modèle correspond donc à
une croissance exponentielle de la population lorsque r > 0 d’où son nom de modèle exponentiel souvent
utilisé à la place de modèle malthusien. Notons qu’il peut s’agir aussi d’une décroissance exponentielle si
r est négatif.
A noter que l’équation différentielle y ′ (t) = ry(t) est définie par la fonction f (y) = ry qui est une
fonction linéaire. L’ensemble de ses solutions s’écrit {y(t) = y(0)ert , y(0) ∈ R} et comporte donc une
infinité de solutions différentes, autant que de valeurs possibles pour la condition initiale y(0).

Exemple : L’idée du modèle logistique, introduit par Verhulst en 1836, est la suivante. Si la population
pouvait croı̂tre indéfiniment, sans rencontrer aucune limitation de ressource ou d’espace, elle aurait une
croissance exponentielle. Mais une croissance exponentielle n’est pas adaptée aux populations que l’on
observe le plus souvent à l’exception peut-être d’une période initiale où la taille de la population est
encore petite, car elle ne tient pas compte des limitations environnementales qui, de fait, ralentissent la
croissance lorsqu’on s’approche de la taille normale de la population qu’on appelle sa capacité biotique
K. D’où l’idée de remplacer le taux constant r par un taux variable r(1 − y(t)
K ) qui dépend de la taille de

1. Les équations différentielles du 2e ordre font intervenir non seulement la fonction y et sa dérivée y ′ mais aussi sa
dérivée seconde y ′′ et les équations d’ordre n font intervenir les dérivées de la fonction jusqu’à l’ordre n.

9
10 CHAPITRE 2. EQUATIONS DIFFÉRENTIELLES EN DIMENSION 1

y
la population. Ce coefficient 1 − K reste proche de 1 lorsque la taille de la population est très petite, ce
qui explique le début de croissance exonentielle, puis il diminue jusqu’à tendre vers 0 lorsque la taille de
la population augmente et tend vers la capacité biotique. L’équation logistique est y ′ (t) = ry(t)(1 − y(t) K ).
En fait le coefficient r(1 − y/K)représente la part de la capacité biotique encore disponible à chaque
instant t. Plus cette part s’amenuise et plus la croissance se ralentit.
L’équation différentielle logistique, y ′ (t) = ry(t)(1 − y(t) y
K ) est définie par la fonction f (y) = ry(1 − K )
y(0)K
qui est un polynôme de degré deux. L’ensemble de ses solutions s’écrit {y(t) = y(0)+(K−y(0))e−rt , y(0) ∈
R}. Il y en a aussi une infinité.

2.2 Le cas linéaire

Lorsque l’équation différentielle est de la forme dy(t)
dt = a(t)y(t)+b(t), où a(t) et b(t) sont des fonctions
données, on dit que l’équation différentielle est linéaire. C’est le cas du modèle exponentiel pour lequel
a(t) = r et b(t) = 0 (a et b sont des fonctions constantes) mais ce n’est pas le cas du modèle logistique.
Les équations linéaires peuvent être résolues explicitement, c’est-à-dire qu’on peut écrire de façon
explicite l’ensemble de
R leurs solutions (que l’on appelle encore la solution générale). Tout d’abord, si
t
a(t)dt
b(t) = 0, y(t) = y(0)e 0 . Plus généralement,lorsque b(t) 6= 0, on doit d’abord rechercher une solution
particulière de l’équation, que l’on notera y ∗ (t), la solution générale s’écrivant alors
Rt
a(t)dt
y(t) = y ∗ (t) + (y(0) − y ∗ (0)) e 0 .

Par exemple, on peut vérifier que y ∗ (t) = te2t est une solution de l’équation y ′ = 2y + e2t et en déduire
que la solution générale de cette équation s’écrit y(t) = y(0)e2t + te2t .
Au dela des équations linéaires, il n’y a qu’un petit nombre d’autres équations différentielles qui
peuvent être résolues explicitement. Le plus souvent, les équations différentielles que l’on est amené à
utiliser ne peuvent pas être résolues explicitement. On a alors recours au calcul approché, nous verrons
comment plus loin, ou bien à l’étude qualitative des solutions, principalement centrée sur l’étude des
équilibres de l’équation et de leur stabilité.

2.3 Etude qualitative dans le cas autonome

Dans le modèle logistique, la fonction f qui définit la dynamique ne dépend que d’une seule variable
(f (y) et non f (t, y)), ce qui entraine que la dynamique ne dépend que de la taille de la population et non
du temps. Ce cas particulier important s’appelle le cas des équations autonomes.
La principale caractéristique du modèle logistique est qu’il présente un équilibre attractif vers lequel
tendent toutes les solutions du modèle, quelque soit leur condition initiale (sauf si y(0)=0 !). Or l’existence
d’équilibres et leurs propriétés (par exemple le fait que les autres solutions tendent vers lui) sont des
éléments que l’on peut souvent déduire directement de l’équation différentielle, même si l’on ne sait pas
calculer explicitement ses solutions. C’est l’étude du graphe de la fonction f (une droite dans le cas
exponentiel et une parabole du cas logistique) qui permet de le faire.
Définition : Pour une équation différentielle de la forme

dy(t)
= f (y(t)), (2.1)
dt
on appelle équilibre ou état stationnaire une valeur constante y ∗ de la quantité y telle que si y(0) = y ∗
alors y(t) = y ∗ pour tout t (la quantité reste à l’équilibre). Un équilibre est donc une solution constante
de l’équation différentielle. Une telle solution a nécessairement une dérivée nulle, c’est-à-dire que l’on a
f (y ∗ ) = 0 ; en d’autres termes y ∗ est aussi un zéro de la fonction f .
Ainsi dans le modèle exponentiel où f (y) = ry, il y a un seul équilibre y ∗ = 0 et dans le modèle
y
logistique où f (y) = ry(1 − K ), il y en a deux, y ∗ = 0 et y ∗ = K.
Dans un modèle de type (2.1), il y a autant d’équilibres différents qu’il y a de zéros différents de
la fonction f . On peut donc visualiser les différents équilibres de l’équation en traçant le graphe de la
fonction f . Les équilibres sont les abscisses des points d’intersection du graphe avec l’axe horizontal (qui
est ici l’axe des y). Et ce graphe permet en outre de visualiser, sur son axe horizontal, un schéma de la
dynamique : il suffit de mettre une flèche dans le sens des y croissants sur les segments de l’axe où f > 0
2.3. ETUDE QUALITATIVE DANS LE CAS AUTONOME 11

(c’est-à-dire où le graphe de f est au dessus de l’axe) et une flèche dans le sens des y décroissants sur les
segments de l’axe où f < 0. Parfois ce schéma de la dynamique est suffisant et peut remplacer à lui seul
une résolution de l’équation (qui, de toute façon, est bien souvent impossible).
Définition : On dit qu’un équilibre y ∗ pour lequelle on a f ′ (y ∗ ) < 0 est un équilibre stable car dans ce cas
l’évolution de toute solution dont la condition initiale est proche de l’équilibre y ∗ est de s’en rapprocher.
De façon analogue, on dit qu’un équilibre y ∗ pour lequelle on a f ′ (y ∗ ) > 0 est un équilibre instable car
dans ce cas l’évolution de toute solution dont la condition initiale est proche de l’équilibre y ∗ est de s’en
éloigner.
On peut vérifier en appliquant ce critère que l’unique équilibre du modèle exponentiel est stable
lorsque r < 0 (extinction) et instable lorque r > 0 (explosion) et de même, si l’on suppose r > 0, on peut
vérifier que l’équilibre y ∗ = K(> 0) du modèle logistique est un équilibre stable (capacité biotique) alors
que y ∗ = 0 est un équilibre instable.
Lorsque f ′ (y ∗ ) = 0, on ne peut pas savoir à partir de f ′ si l’équilibre est stable, instable ou ni l’un ni
l’autre.
La condition f ′ (y ∗ ) < 0 (resp. f ′ (y ∗ ) > 0) est donc un critère de stabilité (resp. d’instabilité) qui se
révèle très opérationnel puisqu’il se calcule facilement. Pour rendre ce critère intuitif, on se reportera à
nouveau au schéma de la dynamique obtenu à partir du graphe de f . On y voit que lorsque f ′ (y ∗ ) < 0
le graphe de f passe au point y ∗ de valeurs positives à des valeurs négatives et donc que la population
croı̂t tant qu’elle est plus petite que y ∗ (puisque f ′ (y) > 0) et décroit tant qu’elle est plus grande. Elle
tend donc dans tous les cas à se rapprocher de l’équilibre. On fait le même raisonnement, inversé cette
fois, dans le cas où f ′ (y ∗ ) > 0.

y
Figure 2.1 – Graphe de la fonction f (y) = ry(1 − K ) dans le plan (y, y ′ ) et esquisse des solutions de
l’équation différentielle y ′ (t) = ry(t)(1− y(t)
K ) dans le plan (t, y). Sur l’axe des y, les points représentent les
équilibres et les flèches indiquent le sens de variation des solutions (croissantes si y ′ > 0 et décroissantes
si y ′ < 0).

La figure ci dessus montre que pour une équation différentielle telle que (2.1), la détermination des
équilibres et du sens de variation des solutions suffit bien souvent pour tracer l’esquisse des solutions
de l’équation. C’est ce qu’on appelle l’étude qualitative. Notons que cette esquisse en dit souvent plus
sur le comportement des solutions que l’expression explicite de la solution générale (lorsqu’elle peut être
calculée) car son l’expression, éventuellement compliquée, se révèle souvent bien peu parlante.
12 CHAPITRE 2. EQUATIONS DIFFÉRENTIELLES EN DIMENSION 1
Chapitre 3

Solutions numériques approchées

Avant de présenter quelques méthodes numériques permettant de calculer des solutions approchées des
équations différentielles, précisons l’importante notion de champs de vecteurs associé à une telle équation.

3.1 Champs de vecteurs associé à une équation différentielle

Bien qu’il soit très rare que l’on puisse résoudre explicitement une équation différentielle donnée,
on peut souvent avoir une idée de l’allure des graphes des solutions en observant le champs de vecteurs
associé. En effet le graphe d’une solution de l’équation y ′ = f (t, y) est par définition tangent à son vecteur
vitesse (1, y ′ (t)) et donc au vecteur (1, f (t, y)). La connaissance de la fonction f en chaque point (t, y)
permet donc de représenter facilement ces vecteurs tangents même si l’on ne connait pas les solutions. Et
si l’on en trace un grand nombre, uniformément répartis dans le plan (t, y), on obtient une représentation
du champ de vecteurs associé à l’équation différentielle qui permet souvent de deviner les graphes des
solutions puisqu’il s’agit des courbes qui sont tangentes en tous leurs points aux vecteurs de ce champs
de vecteurs.
Les deux graphiques suivants sont des exemples de champs de vecteurs associés respectivement à
l’équation linéaire autonôme y ′ = y et à l’équation linéaire non autonôme y ′ = 2y − 5 cos t.

Champ de vecteurs de y’=y

–1 0 1 2 3 4 5

Nous allons voir maintenant que l’idée de la plus simple des méthodes numériques de calcul de solutions
approchées, appelée la méthode d’Euler, est construite sur la notion de champs de vecteurs.

3.2 Solutions numériques approchées

La toute première méthode numérique permettant de calculer une solution approchée d’une équation
différentielle y ′ = f (t, y), étant donnée une condition initiale (t0 , y0 ), a été trouvée par Leonard Euler
en 1768. Elle a été améliorée de façon décisive par deux mathématiciens allemands, Runge et Kutta en

13
14 CHAPITRE 3. SOLUTIONS NUMÉRIQUES APPROCHÉES

1901. C’est la méthode de Runge-Kutta du 4e ordre qui est la plus utilisée aujourd’hui mais il y en a bien
d’autres comme la méthode prédicteur-correcteur d’Adams qui est le choix par défaut de Scilab.
Toutes ces méthodes fournissent une suite de points (tn , yn ) issue du point donné (t0 , y0 ), qu’on peut
aussi voir comme une fonction linéaire par morceaux en reliant les points entre eux par des segments et
qui fournit une approximation de la solution exacte y(t) (solution dont on s’est assuré de l’existence et
de l’unicité). Cette approximation ne fournit pas la solution exacte et elle comporte donc une erreur ;
l’analyse numérique est la branche des mathématiques qui étudie les méthodes de calcul approché comme
celles d’Euler ou de Runge-Kutta, pour tenter de les améliorer mais surtout pour en prévoir et en maitriser
les erreurs.

3.2.1 Algorithme de la méthode d’Euler

Soit y(t) la valeur à l’instant t > t0 de la solution telle que y(t0 ) = y0 . On choisit un pas d’intégration
h, que l’on prendra ici égal à (t−t0 )/n pour un entier n assez grand, et on pose, pour tout i = 0, . . . , n−1

ti+1 = ti + h
(3.1)
yi+1 = yi + hf (ti , yi )

Cette formule est ce que l’on appelle un algorithme : elle permet de calculer, étant donné la condition
initiale (t0 , y0 ), la suite des points (ti , yi ) de proche en proche jusqu’au point (t = tn = t0 + nh, yn ), la
quantité yn étant la valeur approchée cherchée de la solution exacte y(t).
L’idée géométrique de la méthode d’Euler est très simple. Partons du point initial (t0 , y0 ) dans le plan
(t, y). La solution issue de ce point n’est pas connue mais le champs de vecteurs se déduit immédiatement
de l’équation différentielle ; en d’autres termes, on connait la tangente à la solution en ce point. On prend
alors simplement comme approximation de la solution jusqu’à l’instant suivant t1 = t0 + h sa tangente en
ce point. Si h n’est pas trop grand, l’erreur ne le sera pas trop non plus. Cette tangente a pour équation

y = f (t0 , y0 )t + C ste

puisque y ′ (t0 ) = f (t0 , y0 ) et on vérifie aisément que sa valeur en t1 = t0 + h vaut f (t0 , y0 )h + y0 .

3.2.2 Evaluation de l’erreur

Pour évaluer l’erreur commise lorsqu’on remplace la solution exacte y(t) par sa valeur obtenue par la
méthode d’Euler, il convient d’estimer tout d’abord l’erreur faite après un pas, erreur qu’on appelle erreur
locale. On verra ensuite que, comme les erreurs commises à chaque pas peuvent s’accumuler lorsqu’on
répète ce calcul un grand nombre de fois, elles forment ce que l’on appelle l’erreur globale.
En utilisant la valeur en t = t1 du développement de Taylor de la solution exacte y(t) au point initial
t0 (si l’on suppose que la solution est une fonction deux fois continûment dérivable)

h2 ′′
y(t1 ) = y(t0 ) + hy ′ (t0 ) + y (t0 ) + o(h2 ),
2

on voit que les deux premiers terme de ce développement sont précisément la valeur choisie pour y1 . On
a donc
h2
y(t1 ) = y1 + y ′′ (t0 ) + o(h2 ).
2
On voit donc que l’erreur commise sur un seul pas, y(t1 ) − y1 , est une quantité qui tend vers 0 lorsque h
tend vers 0 comme h2 . En d’autres termes, si l’on divise par dix le pas d’intégration h l’erreur est alors
divisée par cent (si h est assez petit).
Mais comme on fait des erreurs du même ordre à chaque pas, on peut montrer par un calcul semblable
que l’erreur globale, celle qui sera commise un fois parvenue au point t, est une erreur d’ordre h et non
plus d’ordre h2 : c’est intuitif si l’on observe que s’il y a n pas et si t − t0 = 1 par exemple, on accumulera
n = 1/h erreurs de taille de l’ordre de grandeur de h2 .
Finalement, on montre que l’erreur à l’instant t vaut yn − y(t0 + nh) = o(h). C’est pourquoi on dit
que la méthode d’Euler est une méthode d’ordre 1.
3.2. SOLUTIONS NUMÉRIQUES APPROCHÉES 15

3.2.3 Algorithme de la méthode RK2 ou méthode du point milieu

C’est une astuce très simple qui permet d’améliorer la méthode d’Euler pour obtenir une méthode,
presque aussi simple, mais qui soit cette fois d’ordre (global) 2 (c’est-à-dire dont l’erreur tend vers zéro
comme le carré du pas h). Voici l’algorithme de la méthode du point milieu :

ti+1 = ti + h
(3.2)
yi+1 = yi + hf (ti + h2 , yi + k21 )
où, à chaque pas, k1 désigne k1 = hf (ti , yi ).
Cette méthode porte le nom de méthode du point milieu car elle consiste à prendre comme pente de
la droite entre ti et ti+1 non pas la pente du champs de vecteurs au point (ti , yi ) comme dans le cas de
la méthode d’Euler mais la pente du champs de vecteurs au point situé au milieu du segment joignant
(ti , yi ) et (ti+1 , yi + hf (ti , yi )) d’Euler. Cela conduit à prendre une valeur intermédiaire entre les pentes
aux deux extrémités, raprochant ainsi le segment de la solution exacte.
On peut montrer que cette astuce permet de diminuer l’erreur commise car elle élimine le terme en
h2 de l’erreur locale, celle-ci devenant alors une erreur locale en h3 . Par contre cette méthode augmente
sensiblement le temps de calcul puiqu’il faut, à chaque pas, évaluer le champs non seulement à l’origine
(t, y) du mais au point (t + h/2, y + k1 /2) situé au milieu du segment utilisé par la méthode d’Euler.

3.2.4 Méthode Runge-Kutta d’ordre 4

La méthode précédente s’appelle aussi la méthode Runge-Kutta d’ordre 2 (RK2) car elle a été inventée
par ces deux auteurs qui ont proposé une famille de méthodes de complexité de plus en plus grande qui
permettent d’obtenir une approximation des solutions d’ordre d pour tout entier d ≥ 1, c’est-à-dire
avec une erreur globale en hd . Les méthodes RK1 et RK2 ne sont autres que la méthode d’Euler et la
méthode du point milieu. En général, dans les applications on n’utilise que rarement les deux méthodes
précédentes car elles produisent des erreurs trop grandes. La méthode RK4 est par contre très populaire
car elle correspond à un bon compromis entre complications de programmation et augmentation du temps
de calcul d’une part et taille de l’erreur commise d’autre part.
Voici les formules pour la méthode RK4 :

ti+1 = ti + h
(3.3)
yi+1 = yi + 16 (k1 + 2k2 + 2k3 + k4 )

où, à chaque pas (ti , yi ), 

 k1 = hf (ti , yi )
hf (ti + h2 , yi + k21 )

k2 =

(3.4)
k = hf (ti + h2 , yi + k22 )
 3


k4 = hf (ti + h, yi + k3 )
A noter que malgré la qualité de cet algorithme, l’utilisateur ne doit jamais perdre de vue qu’il ne
s’agit que d’une approximation de la solution et qu’il existe bien de cas où elle pourrait se révéler très
éloignée de la solution exacte. L’idéal est de compléter chaque fois qu’on le peut, le calcul numérique
approché par une étude qualitative qui permet de contrôler la validité de l’approximation fournie par la
méthode numérique.
16 CHAPITRE 3. SOLUTIONS NUMÉRIQUES APPROCHÉES
Chapitre 4

Existence et unicité des solutions

Le théorème fondamental de Cauchy-Lipschitz garantit l’existence et l’unicité des solutions d’équations

différentielles sous des conditions très peu contraignantes. On ne démontrera pas ici ce théorème car la
preuve peut être facilement trouvée dans tous les ouvrages consacrés aux Equations Différentielles ou,
par exemple, dans Wikipedia. Comme son nom l’indique, ce théorème a d’abord été montré par Cauchy
(en 1844) puis généralisé par Lipschitz peu après.

4.1 Enoncé du théorème fondamental

On rapelle qu’une fonction (t, y) 7→ f (t, y) est dite continuement dérivable ou de classe C 1 dans un
domaine D ⊆ R2 si elle admet des dérivées partielles sur D et si ces dérivées sont continues. En fait, toutes
les fonctions usuelles, polynôme, fraction rationnelle, exponentielle, logarithme, fonction trigonométrique,
etc...sont continuement dérivables sur leur domaine de définition.

Théorème 1 (Théorème d’existence et d’unicité) Soit f : D → R une fonction définie sur un

domaine ouvert D ⊆ R2 et y ′ = f (t, y) l’équation différentielle associée. Si la fonction f est continuement
dérivable alors, pour tout (t0 , y0 ) ∈ D, il existe un intervalle onvert I contenant t0 et une unique solution
y(t) de l’équation différentielle définie sur I et vérifiant la condition initiale y(t0 ) = y0 .

Si l’on attache son attention aux graphes des solutions, représentés dans un plan (t, y), le théorème
précédent affirme que par tout point (t0 , y0 ) passe une unique courbe (t, y(t)) qui est le graphe d’une
solution de l’équation différentielle et qui est définie sur un intervalle I qui contient l’instant intial t0 . Le
graphe d’une solution s’appelle aussi une trajectoire de l’équation différentielle.

4.1.1 Unicité
Ce théorème garantissant l’unicité de la solution correspondant à une condition initiale y0 donnée,
on peut donc de parler de la solution de condition initiale y0 . Autrement dit, à deux conditions ini-
tiales différentes correspondent forcément deux solutions différentes pour toutes les valeurs antérieures et
postérieures de t. Géométriquement, cela signifie que deux trajectoires partant de deux points initiaux
différents ne peuvent ni se couper ni même se toucher.

4.1.2 Cas non autonome

Le théorème précédent se généralise au cas des équations différentielles dépendant du temps qu’on
appelle équation différentielle non autonome, y ′ = f (t, y). Dans ce cas, on demande seulement à f d’être
une fonction continue sur son domaine comme fonction des deux variables t et y et d’être continuement
dérivable comme fonction de y uniquement. Le point initial (t0 , y0 ) doit en outre appartenir au domaine
de définition de f .
On verra plus loin qu’en fait ce théorème reste vrai aussi lorsque y ∈ Rn , c’est-à-dire dans le cas de
systèmes de n équations différentielles.

17
18 CHAPITRE 4. EXISTENCE ET UNICITÉ DES SOLUTIONS

4.1.3 Durée de vie et explosion

On remarquera que le théorème précédent n’affirme rien sur la taille de l’intervalle d’existence d’une
solution qui peut donc être petit : c’est la raison pour laquelle on parle de solution locale et de théorème
d’existence et d’unicité local.
Le fait que les solutions pourraient cesser d’exister au dela d’un certain intervalle de temps I est lié au
phénomène d’explosion qui est bien illustré par l’exemple suivant. Considérons la solution de l’équation
1
différentielle y ′ = y 2 de condition initiale y(0) = 10. Il est facile de voir qu’elle vaut y(t) = 0.1−t et qu’elle
n’est donc pas définie au dela de t = 0.1 car elle présente une asymptote verticale en ce point. Et ce sera
vrai aussi pour n’importe quelle autre solution de condition initiale y0 qui ne sera pas définie au dela de
l’instant t = 1/y0 .
Ce phénomène d’explosion est la signature d’un comportement non linéaire car on peut montrer
facilement (en examinant la forme de la solution explicite par exemple) que pour une équation différentielle
linéaire y ′ = a(t)y + b(t), si les fonctions a(t) et b(t) sont définies pour tout t, alors l’intervalle I peut être
choisi égal à R tout entier. En général, on choisit pour I le plus grand intervalle possible (avant explosion).
Par exemple l’intervalle ] − ∞, 1/y0[ dans le cas des solutions de l’équation y ′ = y 2 de condition initiale
y(0) = y0 . La solution s’appelle alors une solution maximale.

4.2 Dépendance par rapport à la condition initiale et à un pa-

ramètre
Comme chaque solution y(t) est caractérisée par sa condition initiale (t0 , y0 ), on peut considérer
la solution comme une fonction de la condition initiale également y(t, t0 , y0 ). A t0 et t fixés dans I,
l’application y0 7→ y(t, t0 , y0 ) est bien définie d’après le théorème précédent, et on peut montrer qu’elle
est continuement dérivable. En d’autres termes, la valeur à l’instant t d’une solution y(t) est donc une
fonction très régulière de sa valeur à l’instant t0 .
Lorsque l’équation différentielle dépend d’un paramètre a, c’est-à-dire lorsque f (y, a) est une fonction
d’un paramètre, comme c’est le cas de l’équation linéaire y ′ = ry ou de l’équation logistique y ′ =
ry(1 − y/K) qui dépendent des paramètres r et K, on peut aussi voir de quelle façon la solution y(t),
correspondant à une condition initiale fixée (t0 , y0 ), va dépendre de ce (ou ces) paramètres. On peut
s’assurer que, lorsque la fonction a 7→ f (y, a) est continuement dérivable alors la solution a 7→ y(t, a) le
sera également.

4.2.1 Déterminisme et chaos

Le théorème de Cauchy-Lipschitz montre qu’un phénomène régit par une équation différentielle est
parfaitement (et uniquement) déterminé dès que l’on connait sa condition initiale. Ce déterminisme,
qui est en fait la capacité du modèle à prévoir, prévalait à l’époque de Cauchy et de ses successeurs
qui pensaient que la seule limite que l’on avait à pouvoir tout prédire était l’éventuelle incapacité dans
laquelle on était temporairement de connaitre les équations différentielles régissant les lois de la nature.
C’est Poincaré qui, un peu plus tard, a remis en cause ce paradigme à l’occasion de son étude du système
solaire. Il écrivait :“Une cause très petite, qui nous échape, détermine parfois un effet considérable que
nous ne pouvons pas ne pas voir, et alors nous disons que cet effet est dû au hasard. Si nous connaissions
exactement les lois de la nature et la situation de l’univers à l’instant initial, nous pourrions prédire
exactement les lois de la nature et la situation de l’univers à tout instant ultérieur. Mais, lors même
que les lois naturelles n’auraient plus de secret pour nous, nous ne pourrions connaitre la situation
qu’approximativement. Si cela nous permet de prévoir la situation ultérieure avec la même approximation,
c’est tout ce qu’il nous faut, nous disons que le phénomène à été prévu, qu’il est régit par des lois ; mais
il n’en est pas toujours ainsi, il peut arriver que de petites différences dans les conditions initiales en
engendrent de très grandes dans les phénomènes finals ; une petite erreur sur les premières produirait
une erreur énorme sur les derniers. La prédiction devient impossible et nous avons le phénomène fortuit”.
Le théorème de Cauchy-Lipschitz indique bien qu’une prévision parfaite est possible, mais uniquement
sous réserve de connaitre parfaitement la condition initiale. Mais certains systèmes dynamiques sont
rapidement imprévisibles pour la raison qu’indique Poincaré. On utilise alors le terme de chaos pour
décrire cette situation.
Chapitre 5

Le système de Lotka-Volterra

Jusqu’à présent nous avons étudié les équations différentielles en dimention 1 qui fournissent des
modèles dynamiques pour une quantité y(t) évoluant au cours du temps. Ce chapitre aborde le cas de
la dimension 2, c’est-à-dire des modèles dynamiques pour deux quantités (x(t), y(t)) qui évoluent en
interagissant l’une avec l’autre au cours du temps. Pour cela nous allons étudier successivement deux
exemples, le premier qui vient de l’écologie, le modèle de Lotka Volterra, a donné son nom à ce type de
système différentiel et un second qui vient de l’économie et s’appelle le modèle de Goodwin.

5.1 Un exemple en écologie

5.1.1 Le modèle proie-prédateur
Le modèle que nous étudions ici a été proposé par Volterra (et indépendemment par Lotka) en 1926
dans un ouvrage intitulé “Théorie mathématique de la lutte pour la vie” qui est probablement le premier
traité d’écologie mathématique. Volterra avait été consulté par le responsable de la pêche italienne à
Trieste qui avait remarqué que, juste après la première guerre mondiale (période durant laquelle la pêche
avait été nettement réduite) le nombre de requins et autres prédateurs impropres à la consommation que
l’on relevait involontairement dans les filets parmi les poissons consommables était nettement supérieur à
ce qu’il avait été avant guerre alors que la population des sardines que l’on avait beaucoup moins pêché,
semblait avoir diminué. Ceci apparaissait comme un paradoxe que Volterra parvint à expliquer (voir
ci-dessous) avec le modèle qu’il proposa et qui porte aujourd’hui son nom.
Notons respectivement x(t) et y(t) la taille des deux populations à l’instant t, la seconde (ici des
requins, appelés les prédateurs) se nourissant de la première (ici des sardines, appelés les proies). On fait,
sur la dynamique de ces deux populations plusieurs hypothèses, inévitablement simplificatrices, qui vont
nous permettre d’écrire le modèle : on suppose d’une part que les proies disposent de nouriture en quantité
illimitée et que seuls les prédateurs s’opposent à leur croissance. On suppose aussi que le nombre de
prédateurs est limité par la quantité de proies dont ils disposent pour se nourir et qu’en l’absence de proies,
la population des prédateurs ne peut survivre. Enfin en ce qui concerne les interactions entre ces deux
populations, on suppose que le nombre de rencontres entre proies et prédateurs est à la fois proportionnel
à x(t) et y(t), donc proportionnel au produit x(t)y(t) et que le taux de disparition des proies ainsi que le
taux de croissance des prédateurs dus à ces rencontres sont l’un et l’autre proportionnels au nombre de
rencontres entre les deux populations. Ceci conduit aux deux équations différentielles suivantes :
 ′
x = α1 x − β1 xy
(5.1)
y ′ = −α2 y + β2 xy

où les quatre constantes α1 , α2 , β1 et β2 sont respectivement le taux de croissance des proies, le taux de
(dé)croissance (naturelle) des prédateurs et les coefficients d’interaction entre les deux populations.

5.1.2 Présence de solutions périodiques

Pour étudier plus facilement le système, donnons aux quatre paramètres α1 , α2 , β1 et β2 des valeurs
particulières (qui ne cherchent pas à être réalistes),

α1 = 0.8, α2 = 0.6, β1 = 1.6, β2 = 1.8. (5.2)

19
20 CHAPITRE 5. LE SYSTÈME DE LOTKA-VOLTERRA

ce qui permet de calculer des solutions avec Scilab. Si l’on trace les trajectoires dans le plan et aussi les
graphes des solutions en fonction du temps, on observe comme sur la figure ci-dessous que les trajectoires
sont fermées et que leurs deux composantes x(t) et y(t) sont périodiques.

Figure 5.1 – A gauche le graphe des effectifs des proies x(t) (en trait plein) et celui des prédateurs y(t)
(en pointillés), dans le cas où x(0) = 0.5 et y(0) = 0.3. A droite plusieurs trajectoires du système (5.1)
représentées dans le plan de phase (la trajectoire en gras est issue du point de coordonnées (0.5, 0.3) et
correspond aux deux graphes de gauche).
Il est intéressant de noter que les oscillations de la taille des deux populations ne sont pas dues à des
variations de leur environnement mais à leur interaction. En effet, en suivant ces évolutions sur la partie
gauche de la figure à partir de l’instant t = 0, on observe que la diminution du nombre de proies entraine,
avec un petit décalage dans le temps, une diminution du nombre de prédateurs qui en viennent à manquer
de nourriture, diminution qui, à son tour, rendra possible une nouvelle augmentation du nombre de proies
profitant de l’absence de prédateurs. Mais cette augmentation va permettre à son tour un redémarrage de
la croissance des prédateurs et ainsi de suite.

5.1.3 Comment lire ce comportement directement sur l’équation ?

Le commentaire ci-dessus du tracé des solutions peut se lire directement sur l’équation. En effet, il
porte sur la croissance et décroissance de la population x(t) de proies et celles de la population y(t)
de prédateurs. Or la croissance de la solution x se lit sur le signe de x′ , donc sur le signe de f (x, y) =
α1 x − β1 xy = x(α1 − β1 y), qui ne dépends que du signe de x et α1 − β1 y ; si ces deux quantités sont de
même signe x augmente, si elles sont de signe opposé x diminue. De même, la croissance de la population
y(t) se lit sur le signe de y ′ , donc sur le signe de g(x, y) = −α2 y + β2 xy = y(−α2 + β2 x). Cette remarque
est tout-à-fait générale : les courbes f (x, y) = 0 et g(x, y) = 0 découpent des régions du plan sur lesquelles
les solutions x et y sont monotones.
On peut voir aussi que, parmi les trajectoires, l’une d’elles joue un rôle particulier :  c’est l’équilibre

α2 /β2
(0.3; 0.5). Dans le cas général du système (5.1), l’équilibre (non nul) a pour coordonnées .
α1 /β1
Définition : On appelle équilibre toute solutions constante du système d’équations différentiel.
La recherche des équilibres se ramène donc simplement à celle des solutions du système d’équations
f (x, y) = 0 et g(x, y) = 0.

5.1.4 Que s’est-il passé dans l’adriatique ?

Que s’est-il passé dans l’adriatique durant la première guerre mondiale qui a vu diminuer la population
de sardines et augmenter celle des requins, alors que, justement la pêche avait été presque arrêtée ? L’arrêt
de la pêche a eu pour effet d’augmenter α1 , puisque la mortalité des sardines (par pêche) a diminué, et
5.2. UNE THÉORIE DU CYCLE ÉCONOMIQUE 21

de diminuer α2 pour les mêmes raisons. Par ailleurs, les coefficients d’interactions β1 et β2 n’ont pas été
modifiés par l’arrêt de la pêche. L’équilibre (qui est aussi le point autour duquel les oscillations se font)
a donc été déplacé durant la période de guerre, la valeur d’équilibre des proies ayant diminué (avec α2 )
et celle des prédateurs ayant au contraire augmenté.

5.2 Une théorie du cycle économique

Le modèle de Goodwin 1 a été l’un des premier modèle de cycle économique. Avant Goodwin on
considérait volontiers que les oscillations parfois présentes dans les évolutions de quantités économiques
que l’on étudie s’expliquaient avant tout par des chocs externes intervenus comme perturbations des
dynamiques de ces quantités. On pensait en effet que ces dynamiques, elles, n’avaient aucune raison de
présenter des oscillations par elles-mêmes. C’est la notion de cycles exogènes, par opposition à celle de
cylces endogènes qui, eux, sont des trajectoires oscillantes, naturellement présentes dans la dynamique. Le
modèle de Goodwin fut l’un des premiers exemples de cycles endogènes ; il s’agit de cycle de répartition
salaire/profit. Il montre une dynamique qui ne tend pas vers un équilibre mais au contraire est constituée
de cycles dont la période et l’amplitude dépendent des conditions initiales.

5.2.1 Introduction des variables pertinentes

Dans ce modèle, les variables Y , K et L représentent, comme c’est le cas habituellement, la production
nationale, le capital ou richesse du pays, et la force de travail ou nombre d’employés dans les entreprises. Il
s’agit d’un modèle dans lequel on ne suppose pas le plein emploi, c’est-à-dire que l’on aura à distinguer la
quantité L d’une autre quantité, notée N qui représente la demande d’emploi et qui n’a pas de raison d’être
identique. Les deux variables économiques dont on se propose d’étudier la dynamique sont respectivement
le taux d’emploi de la main d’oeuvre, noté x,

L
x :=
N
et la part des salaires dans le revenu national, noté z,

wL
z=
Y
où w représente le salaire moyen et donc wL le salaire total.

5.2.2 Construction du modèle

Cinq équations gouvernent l’évolution dans le temps de ces diverses quantités. Elles correspondent
chacunes à des hypothèses, bien entendu extrêmement simplificatrices, faites sur le modèle :
— La production nationale est proportionnelle au capital investi, ce qui permet d’écrire K comme
une fonction linéaire de Y , le coefficient v, représentant la productivité du capital, étant supposé
constant :
Y
=v (5.3)
K
— Ce revenu nécessite une main d’oeuvre dont la productivité évolue au cours du temps selon la
relation :
Y
= ueµt , µ > 0 et u > 0 constantes. (5.4)
L
— La demande de travail augmente au taux constant n (=natalité − mortalité) :

N ′ = nN (5.5)

— L’accroissement du capital (épargne) est égal au solde “revenus moins salaires” qui est ainsi supposé
entièrement réinvesti :
K ′ = Y − wL. (5.6)
1. Goodwin R.M. A growth cycle, in C.F. Feinstein ed. Socialism, capitalism and economic growth, Cambridge University
Press, 1967
22 CHAPITRE 5. LE SYSTÈME DE LOTKA-VOLTERRA

— Le taux d’accroissement du salaire moyen est une fonction croissante du taux d’emploi : c’est la loi
de Phillips 2 . Pour simplifier, on suppose ici que cette fonction croissante est simplement affine :

w′
= ax − b, a > 0 et b > 0 constantes. (5.7)
w
On déduit facilement 3 de ces cinq relations les deux équations différentielles qui gouvernent la dyna-
′
mique de z et de x. Tout d’abord on calcule, à partir de la définition de z, le rapport zz :

z′ w′ L′ Y′
= + − . (5.8)
z w L Y
Donc, de (5.7) et (5.4) on déduit une première équation

z′
= ax − b − µ. (5.9)
z
x′
Puis on calcule, à partir de la définition de x, le rapport x :

x′ L′ N′
= − . (5.10)
x L N
Donc, de (5.5), (5.3) et (5.4) on déduit :

x′ K′
= −µ + − n.
x K
K′ Y (1−z)
Mais comme, en appliquant (5.6) et (5.3), K = Y /v , on obtient une deuxième équation :

x′
= v(1 − z) − (µ + n). (5.11)
x
Finalement la dynamique du vecteur (x, z) est donc régie par le système différentiel suivant :
 ′
x = (v − (µ + n))x − vxz
(5.12)
z ′ = −(b + µ)z + axz

5.2.3 Les cycles économiques

Le sytème (5.12) s’appelle en économie le modèle de Goodwin mais il appartient clairement à la
même famille que le système différentiel de Lotka-Volterra (5.1)car il en a exactement la forme : ici
α1 = v − (µ + n), β1 = v, α2 = b + µ, et β2 = a. On connait donc sa dynamique : il présente en particulier
un équilibre (non nul)  ∗   
x (b + µ)/a
=
y∗ 1 − (µ + n)/v
et toutes les autres solutions (positives) sont périodiques et oscillent autour de cet équilibre. Ce sont des
cycles.

2. Phillips A.W. The relationship between unemployment and the rate of change of money wage rates in the United
Kingdom 1861-1957, Economica, vol 25, 1958
′ ′ ′ ′ ′ ′
3. Si γ = αβ alors γγ = αα + ββ , et de manière similaire si γ = α
β
alors γγ = αα − ββ . Ceci se démontre très facilement
γ′
en calculant la dérivée γ
de log(γ).
Chapitre 6

Noeuds, cols, foyers et centres

L’étude de l’exemple du système de Lotka Voterra montre qu’une bonne connaissance des équilibres
fournit des renseignements précieux sur le comportements des trajectoires proches. C’est ce que nous
allons voir, d’abord dans le cas des systèmes différentiels linéaires puis, dans le prochain chapitre, dans
le cas général des systèmes différentiels non nécessairement linéaires de dimension 2.

6.1 Systèmes de deux équations différentielles

On considère le système de deux équations différentielles suivant :
 ′
x = f (x, y)
(6.1)
y ′ = g(x, y)

où f et g sont deux fonctions que l’on supposera lisses (c’est-à-dire continûment dérivables).
On appelle solution du système (6.1) un vecteur (x(t), y(t)) dont les deux coordonnées sont des
fonctions du temps qui vérifient le système différentiel, c’est-à-dire telles que l’on a x′ (t) = f (x(t), y(t))
et aussi y ′ (t) = g(x(t), y(t)). On appelle condition initiale la valeur de la solution à l’instant initial (que
l’on choisit souvent égal à 0), c’est-à-dire le vecteur (x(0), y(0)).
Par exemple pour le système différentiel suivant, appelé oscillateur harmonique,
 ′
x = −y
(6.2)
y′ = x

on peut vérifier facilement que, pour toutes les valeurs de r ≥ 0 et θ ∈ [0, 2π[, le vecteur (x(t), y(t)) =
(r cos(t − θ), r sin(t − θ)) est une solution du système et aussi que, parmi toutes ces solutions, celle qui
correspond à r = 2 et θ = 0, (x(t), y(t)) = (2 cos t, 2 sin t), est la solution du système de condition initiale
(2, 0).
Comme pour les équations différentielles, on peut rarement calculer les solutions exactes d’un système
différentiel. Mais, comme pour les équations différentielles, on peut montrer que pour assurer l’existence
et l’unicité des solutions du système de deux équations différentielles (6.1), étant donnée une condition
initiale (x(0), y(0)), il suffit que les fonctions f et g qui définissent l’équation soient de classe C ∞ . On peut
donc, à défaut de savoir calculer des solutions exactes, chercher à décrire le comportement des solutions
soit par une étude qualitative, soit en calculant des solutions approchées (ou, mieux encore, lorsque c’est
possible, en combinant les deux approches).

6.1.1 Trajectoires et champs de vecteurs

On peut représenter géométriquement les solutions du système différentiel de deux façons différentes :
soit on trace les graphes de chacunes des deux composantes de la solution comme des fonctions du temps,
soit on trace la courbe image de t → (x(t), y(t)) qui est une courbe paramétrée dans le plan (x, y) qu’on
appelle une trajectoire du système.
Dans le cas de l’oscillateur harmonique, les trajectoires sont des cercles concentriques (pourquoi ?).
On sait que la “vitesse de déplacement” sur la courbe solution est donnée par le vecteur vitesse que l’on
peut calculer simplement à l’aide des dérivées des deux composantes de la solution
 ′ 
x (t)
V = .
y ′ (t)

23
24 CHAPITRE 6. NOEUDS, COLS, FOYERS ET CENTRES
p
A noter que plus sa longueur kV k = x′ (t)2 + y ′ (t)2 est grande et plus la courbe est parcourue rapidement
par la dynamique associée au système.
Bien qu’on ne connaisse pas en général les trajectoires, on connait néanmoins leurs vecteurs tangent
V (x, y) en tout point puisqu’il est donné par le système différentiel V (x, y) = (f (x, y), g(x, y)). Au
système différentiel correspond donc un champs de vecteurs dans le plan. Et les courbes paramétrées
t 7→ (x(t), y(t)) qui sont solutions du système différentiel sont les courbes tangentes en chacun de leurs
points au vecteur de coordonnées (f (x, y), g(x, y)).

6.1.2 Etude qualitative, isoclines, équilibres

L’étude qualitative du système consiste à déterminer, à partir d’un examen du système, un aperçu du
champs de vecteurs afin d’en déduire l’allure des trajectoires.
Pour cela on remarque que si f (x, y) = 0 en un point, le vecteur du champs de vecteur sera vertical
en ce point, et de même si g(x, y) = 0, il sera horizontal. On en déduit que la courbe d’équation g(x, y) =
0, appelée isocline y ′ = 0, est une courbe sur laquelle les solutions t 7→ (x(t), y(t)) ont une tangente
horizontale. De même la courbe d’équation f (x, y) = 0, appelée isocline x′ = 0, est une courbe sur
laquelle les solutions t 7→ (x(t), y(t)) ont une tangente verticale. Les points d’intersections des isoclines
sont les équilibres (x∗, y∗) du système c’est-à-dire les points tels que la trajectoire issue d’un tel point
reste en ce point pour tout t (car si x′ (t) et y ′ (t) sont nuls, x(t) et y(t) sont constants).
Dans chacune des régions du plan délimitées par les isoclines, les quantités f (x, y) et g(x, y) sont de
signe constant et on peut donc schématiser la direction du champs de vecteurs par une flèche de l’un des
quatres types suivants : vers la droite et vers le haut (si f > 0 et g > 0), vers la droite et vers le bas (si
f > 0 et g < 0), vers la gauche et vers le haut (si f < 0 et g > 0), ou vers la gauche et vers le bas (si
f < 0 et g < 0).
La représentation des deux isoclines y ′ = 0 et x′ = 0, des équilibres (à l’intersection des isoclines) et
du schéma des flèches du champs de vecteurs constituent l’étude qualitative du système d’où l’on pourra
le plus souvent déduire, en utilisant la propriété qu’ont les trajectoires de ne jamais se croiser, l’allure
des solutions en fonction de leur condition initiale (x(0), y(0)). A noter qu’en raison de l’unicité aucune
trajectoires ne traverse un équilibre.

6.1.3 Calcul de solutions approchées

Comme dans le cas d’une équation différentielle unique, la méthode d’Euler permet de calculer une
solution approchée (x̃(t), ỹ(t)) passant par un point donné (x0 , y0 ) pour un système différentiel. L’idée
est ici encore de suivre le champ de vecteur durant un pas de temps h, supposé assez petit, puis de
recommencer durant un second pas de temps h, en suivant à nouveau le champ de vecteur au point que
l’on a atteint. L’algorithme est le suivant : on définit une suite de points du plan de coordonnées (x0 , y0 ),
(x1 , y1 ) , ...(xn , yn ), ... par

xn = xn−1 + hf (xn−1 , yn−1 )
(6.3)
yn = yn−1 + hg(xn−1 , yn−1 )
et on relie ces points entre eux de façon à former une ligne brisée. Si h est suffisamment petit, cette ligne
brisée fournit une approximation de la solution exacte. En effet on peut montrer que si h = (t − t0 )/n,
alors la suite (xn , yn ) converge, quand n tend vers l’infini, vers la valeur exacte de la solution (x(t), y(t))
à l’instant t ≥ t0 .
On peut aussi vérifier que cet algorithme, comme dans le cas d’une équation différentielle unique,
est un algorithme d’ordre 1 (c’est-à-dire que l’erreur tend vers zéro comme le pas h) et qui peut être
remplacé par des algorithmes plus performants, notamment des algorithmes d’ordre plus élévé. Parmi les
plus utilisé on trouve l’algorithme de Runge-Kutta d’ordre 4.

6.2 La classification de Poincaré

Henri Poincaré a introduit une classification des champs de vecteurs linéaires du plan qui regroupe ces
champs en un nombre fini de classes selon l’aspect géométrique de la famille des trajectoires du champs.
Cette classification est très importante car on l’utilise non seulement pour l’étude de champs linéaires
mais surtout pour l’étude de champs non linéaires dont on approxime l’allure des trajectoires au voisinage
des points d’équilibre par celle des trajectoires du linéarisé du champs comme nous allons le voir plus
loin.
6.2. LA CLASSIFICATION DE POINCARÉ 25

Valeurs propres réelles λ et µ

0<λ<µ Noeud instable

0 < λ = µ, A diagonalisable Noeud dégénéré instable

0 < λ = µ, A non diagonalisable Noeud instable

λ<0<µ Col

λ = µ < 0, A diagonalisable Noeud dégénéré stable

λ = µ < 0, A non diagonalisable Noeud stable

µ<λ<0 Noeud stable

Valeurs propres complexes (α ± iω, ω 6= 0)

α>0 ou Foyer instable

α=0 Centre

α<0 ou Foyer stable

Figure 6.1 – Classification de Poincaré des systèmes linéaires

26 CHAPITRE 6. NOEUDS, COLS, FOYERS ET CENTRES

Un système différentiel (6.1) est dit linéaire lorsque les deux fonctions f et g sont des fonctions linéaires
de x et de y. Il s’écrit alors sous forme matricielle de la façon suivante :
 ′      
x a b x x
= = A
y′ c d y y

où A est une matrice réelle 2 × 2. On supposera que A est non dégénérée, c’est-à-dire que 0 n’est pas
une valeur propre. On notera λ et µ les deux valeurs propres de A lorsqu’elles sont réelles et on notera
α ± iω ces deux valeurs propres de A lorsqu’elles sont complexes. On sait qu’il existe une base de R2 dans
laquelle l’application linéaire associée à A a pour matrice l’une des suivantes :
     
λ 0 λ 1 α ω
0 µ 0 λ −ω α

En notant U et V les coordonnées dans cette base, il est aisé de résoudre le système en U et V : dans
le premier cas on a (U, V ) = (eλt U0 , eµt V0 ), dans le second (U, V ) = eλt (U0 + tV0 , V0 ) et enfin dans le
troisième     
U cos ωt sin ωt U0
= eλt .
V − sin ωt cos ωt V0
Il est dès lors facile d’en déduire le comportement des solutions. La figure 6.1 indique ce comportement,
selon les valeurs respectives de λ, µ et α, ainsi que les noms donnés par Poincaré à ces divers cas.

6.2.1 Calcul des solutions

Que peut-on dire des solutions d’un système linéaire lorsqu’on connait les valeurs propres et vecteurs
propres du système ?    
′ x(t) a b
Dans le cas d’un système linéaire V = AV , où V = et A = , on a les résultat
y(t) c d
suivants :
— L’ensemble des solutions du système est un espace vectoriel de dimension 2 (car toute combi-
naison linéaire de solutions est encore une solution) et il suffit donc de connaitre deux solutions
linéairement indépendantes pour connaitre toutes les solutions (qui seront des combinaison linéaires
de ces deux solutions linéairement indépendantes).  
xλ
— Si λ est une valeur propre de la matrice A, de vecteur propre Vλ = , alors le vecteur V (t) =
yλ
e Vλ est une solution puisqu’on a V (t) = λe Vλ et aussi AV = Ae Vλ = eλt AVλ = eλt λVλ (car
λt ′ λt λt

Vλ est un vecteur propre de A). On a donc bien Vλ′ = AVλ .

— Les deux observations précédentes permettent de résoudre explicitement le système linéaire lorsque
la matrice A possède deux valeurs propres réelles λ1 et λ2 et deux vecteurs propres associés V1 et
V2 linéairement indépendants. Toutes les solutions du système sont alors des combinaisons linéaires
des deux solutions particulières eλ1 t V1 et eλ2 t V2 .
— Dans le cas où les deux valeurs propres sont complexes conjuguées, si λ = α + iω est l’une de ces
deux valeurs propres, et Vλ un vecteur propre associé, on peut vérifier facilement que la partie
réelle et la partie imaginaire de la solution complexe eλt Vλ sont des solutions (réelles) du système
et elles sont linéairement indépendantes. Donc ces deux solutions Re(eλt Vλ ) et Im(eλt Vλ ) forment
donc une base de l’ensemble des solutions.
Chapitre 7

Linéarisé d’un système différentiel

non linéaire

La classificiation de Poincaré permet de distinguer, pour les systèmes différentiels linéaires, leur nature,
selon qu’ils présente un noeud, un col, un foyer ou un centre ainsi que leur stabilité. On va voir à présent
qu’elle va aussi permettre de prédire le comportement local des trajectoires (et donc des solutions) de
tout système différentiel non linéaire, c’est-à-dire le comportement au voisinage des équilibres. En effet,
au voisinage de chaque équilibre, le système est localement linéaire. La nature d’un équilibre d’un système
non linéaire sera simplement la nature de son linéarisé au voisinage de cet équilibre.

7.1 Linéarisé au voisinage d’un équilibre

Supposons que (x∗ , y ∗ ) soit un équilibre du système différentiel
 ′
x = f (x, y)
(7.1)
y ′ = g(x, y)

c’est-à-dire un zéro commun de∗ f et g. Soit ε > 0 un très petit paramètre. Effectuer le changement de
∗ y−y
variables X := x−x ε , Y := ε revient à regarder à la loupe au voisinage de l’équilibre (x∗ , y ∗ ). En effet,
∗ ∗
lorsque x − x et y − y sont très petits, de l’ordre de ε, X et Y sont alors des grandeurs appréciables
et donc les dessins obtenus dans le plan (X, Y ) correspondent à l’image de points (x, y) très proches de
l’équilibre. Après calculs, on constate que le sytème obtenu sous la loupe peut s’écrire sous la forme
 ′
X = aX + bY + o1 (ε)
(7.2)
Y ′ = cX + dY + o2 (ε)

où o1 (ε) et o2 (ε) sont des expressions qui contiennent ε en facteur et qui donc tendent vers 0 avec ε. Si l’on
néglige ces termes, le système différentiel devient linéaire (cela signifie que quand on regarde à la loupe
un système différentiel au voisinage d’un de ses équilibres, on voit un système différentiel pratiquement
linéaire), c’est-à-dire qu’il peut s’écrire sous la forme
 ′      
X a b X X
= = A
Y′ c d Y Y

La matrice A s’appelle la matrice jacobienne du système initial. On peut calculer facilement cette matrice
A à partir des dérivées partielles de f et g calculées au point d’équilibre (x∗ , y ∗ ). En effet on a
!
∂f ∗ ∗ ∂f ∗ ∗
∗ ∗ ∂x (x , y ) ∂y (x , y )
A = A(x , y ) = ∂g ∗ ∗ ∂g ∗ ∗
∂x (x , y ) ∂y (x , y )

On a vu que les équilibres d’une dynamique linéaire sont principalement de 4 types, noeuds, cols,
foyers et centres, les noeuds et les foyers se divisant eux-même en deux catégories selon qu’ils sont stables
ou instables. Par extension, on dit qu’un équilibre d’un système non linéaire est un noeud, un col, un foyer
ou un centre selon le type de son linéarisé. Le type de l’équilibre s’appelle sa nature. La connaissance de la
nature des équilibres d’une dynamique apporte souvent des renseignement précieux sur le comportement
des trajectoires au voisinage des points d’équilibre.

27
28 CHAPITRE 7. LINÉARISÉ D’UN SYSTÈME DIFFÉRENTIEL NON LINÉAIRE

coexistence des deux especes coexistence improbable: extinction de l’une des deux especes

2 1

0.8

1.5

0.6

y y
1

0.4

0.5

0.2

0 0.5 1 1.5 2 0 0.2 0.4 0.6 0.8 1

x x

7.2 Exemples
A titre d’exemples, considérons les deux systèmes différentiels suivants :
 ′
x = (2 − x − 2y/3)x
(7.3)
y ′ = (2 − 2x/3 − y)y
 ′
x = (1 − x − 2y)x
(7.4)
y ′ = (1 − 2x − y)y
dont la figure ci-dessus représente les champs de vecteurs associé. Ces deux exemples présentent chacun 3
équilibres situés sur les axes de coordonnées, (0, 0), (0, 2), (2, 0) et un quatrième équilibre de coordonnées
( 56 , 65 ) dans le premier exemple et de coordonnées ( 13 , 13 ) dans le second. On détermine la nature de ces
équilibres en linéarisant le système au voisinage de chacun de ces équilibres et en calculant, dans chaque
cas, les valeurs propres de la matrice jacobienne A trouvée. On peut ensuite controler en vérifiant si la
nature trouvée est compatible avec ce que l’on peut voir sur la figure. On établit ainsi que le quatrième
équilibre est un noeud stable dans le premier système et un col dans le second.

7.3 Le cas où le linéarisé est un centre :

Il est important de comprendre que si la classification de Poincarré des champs linéaires du plan
permet en général d’en déduire la nature des équilibres des champs non linéaires par simple étude de leur
linéarisé, cela n’est pas le cas lorsqu’il s’agit d’un centre. En effet si le linéarisé d’un champ au voisinage
d’un de ses équilibre est un noeud stable ou un noeud instable, un col, ou bien encore un foyer stable ou
un foyer instable, alors il en sera de même de l’équilibre du champs non linéaire.
Par contre si le linéarisé est un centre, alors il y a trois possibilités pour l’équilibre du champs non
linéaire : soit il est un centre, soit il est un foiyer stable, soit il est un foyer instable. En effet, lorsqu’on
perturbe légèrement un champs au voisinage d’un centre en lui ajoutant des termes non linéaires, cela
suffit le plus souvent à casser la périodicité des solutions qui continuent à tourner autour de l’équilibre
mais en suivant une spirale qui, soit tend vers l’équilibre (foyer stable), soit s’en éloigne (foyer instable).
On n’est donc jamais assuré que l’équilibre d’un champ non linéaire dont le linéarisé est un centre, est
lui-même un centre. C’est rarement le cas.
L’exemple du système de Lotka Volterra est pourtant un exemple où le champ non linéaire présente,
comme son linéarisé, un équilibre qui est un centre. Rappelons son expression :
 ′
x = α1 x − β1 xy
(7.5)
y ′ = −α2 y + β2 xy

Ce système (7.5) présente deux équilibres, l’origine (qui n’a guère d’intérêt dans ce modèle) et l’équilibre
(α1 /β1 , α2 /β2 ) autour duquel tournent les autres trajectoires. Un examen du linéarisé du système montre
que le premier est un col (deux valeurs propres réelles de signe opposé) et le second un centre (deux
valeurs propres imaginaires pures).
Mais bien que le système (5.1) ait un équilibre dont le linéarisé est un centre, rien ne permet, à priori,
d’affirmer que ses solutions sont périodiques, c’est-à-dire que l’équilibre (α1 /β1 , α2 /β2 ) est effectivement
7.3. LE CAS OÙ LE LINÉARISÉ EST UN CENTRE : 29

un centre et non un foyer. Pour s’en assurer, on a recours à l’existence d’une intégrale première pour ce
système.
Définition : Soient ((x(t), y(t)) la dynamique de deux espèces comme par exemple les proies et les
prédateurs du système de Lotka Volterra, on dit que la fonction H(x, y) est une intégrale première ou loi
de conservation de cette dynamique lorsque la quantité H(x(t), y(t)) reste constante au cours du temps.
Par exemple si l’on considère la dynamique de l’oscillateur harmonique, dont les solutions sont de
la forme (x(t), y(t)) = (r cos t, r sin t), la fonction H(x, y) = x2 + y 2 est une loi de conservation puisque
H(x(t), y(t)) = r2 .
Dans le cas du système de Lotka-Volterra (5.1), il n’est pas difficile de vérifier que la fonction suivante
est une loi de conservation :

H(x, y) = α1 ln y − β1 y + α2 ln x − β2 x.

Il suffit en effet de calculer le produit scalaire du gradient de H par le champs de vecteur (x′ (t), y ′ (t))
qui est donné par le système différentiel.
L’importance des lois de conservation pour l’étude des systèmes différentiels comme ceux de Lotka-
Volterra est facile à comprendre. Dès que la fonction H est connue, on peut, en utilisant ses dérivées
partielles, tracer ses courbes de niveau et en déduire les trajectoires de la dynamique. Alors qu’une
étude qualitative permet de prévoir l’oscillation des deux populations (car les trajectoires tournent dans
le plan (x, y)), elle ne permet pas de s’assurer que la dynamique est réellement périodique, c’est-à-
dire que les trajectoire se referment effectivement après un tour. Au contraire cette information découle
immédiatement de l’existence de la loi de concervation H.
30 CHAPITRE 7. LINÉARISÉ D’UN SYSTÈME DIFFÉRENTIEL NON LINÉAIRE
Chapitre 8

Dynamiques de populations
structurées en âges

8.1 Le modèle de Leslie

Le modèle présenté ici est dû à Sir Paul Leslie (1945) et il est l’un des plus utilisé en dynamique des
populations et en démographie. Il suppose que la population étudiée est constituée de plusieurs groupes
d’individus à des stades différents ou classes d’ages différentes (oeufs, oisillons, oiseaux, par exemple ou
bien graines, rosettes, plantes en fleurs, etc...). Les effectifs de chacune des classes évoluent de façons
différentes mais pas indépendemment les unes des autres. On va étudier la dynamique de ce type de
modèle et notamment chercher à répondre aux deux questions suivantes :
1. l’effectif total, somme des effectifs des différentes classes, a-t-il une croissance exponentielle avec
un taux de croissance constant, et dans ce cas, comment calculer ce taux ?
2. La répartition des individus dans les différentes classes, la distribution initiale, se maintient-elle
au cours du temps ou bien se modifie-t-elle et de quelle façon ?

Exemple : Pour commencer examinons un exemple. Il s’agit d’une population de rongeurs ayant un
cycle de reproduction de 3 ans. On ne considère ici que la sous-population formée des individus femelles.
On suppose que chaque femelle donne en moyenne naissance à 6 femelles durant sa deuxième année et à
10 femelles durant sa troisième année. Cependant, seul un rongeur sur deux survit au dela de sa première
année et seul 40% de ceux qui survivent la deuxième année survivront jusqu’à la troisième année.
Si l’on désigne respectivement par j(t), p(t) et a(t) les effectifs à l’instant t des femelles juvéniles,
des femelles préadultes (rongeurs de 1 an) et des femelles adultes (rongeurs de 2 ans), les informations
précédentes peuvent s’écrire : 
 j(t + 1) = 6p(t) + 10a(t)
p(t + 1) = 0, 5j(t) (8.1)
a(t + 1) = 0, 4p(t)


Ces formules (8.1) permettent, à partir des effectifs initiaux des trois classes, (j0 , p0 , a0 ), de calculer les
effectifs (j1 , p1 , a1 ) à l’instant suivant t = 1, puis, (j2 , p2 , a2 ) à l’instant t = 2 et ainsi de suite. Si l’on
désigne par N (t) = j(t) + p(t) + a(t) l’effectif total de la population à l’instant t (et donc N0 l’effectif
initial), on peut également calculer à partir de (8.1) les termes successifs de la suite (N (t)), ce qui permet
d’apréhender aussi la dynamique de cette population dans son ensemble. Pour avoir une idée du taux de
croissance de chacune des classes, on peut calculer les quotients j(t+1) p(t+1) a(t+1)
j(t) , p(t) et a(t) pour t = 0, 1, 2, ...
mais le résultat est très irrégulier et on voit mal sur ces premiers termes quel taux de croissance on
pourrait retenir pour rendre compte de la dynamique de ces différentes classes d’age. Et si l’on considère
la population dans son ensemble, les quotients NN(t+1) (t) ne sont pas plus réguliers.
Par contre si on laisse le temps augmenter, on constate que ces taux tendent tous vers la même valeur
λ, ici λ = 2, c’est-à-dire qu’après un certain temps, la dynamique considérée consiste simplement en
une multiplication par un facteur 2 des effectifs de chaque classe d’une période à la suivante. Ce facteur
multiplicatif, qui correspond à un taux de croissance asymptotique peut être calculé facilement comme
nous allons le voir.

31
32 CHAPITRE 8. DYNAMIQUES DE POPULATIONS STRUCTURÉES EN ÂGES

Si l’on s’intéresse maintenant non plus à la dynamique des effectifs mais à l’évolution au cours
du temps de la répartition des individus entre les diverses classes, on peut aussi calculer, à partir
de la répartition initiale des individus selon ces trois classes v0 = N10 (j0 , p0 , a0 ) l’évolution de cette
répartition au cours du temps v(t) = N1(t) (j(t), p(t), a(t)). On constate que cette répartition tend vers
1
une répartition asymptotique qui est celle du vecteur v = 130 (100, 25, 5), c’est-à-dire la répartition
100 25 5
( 130 , 130 , 130 ) ≃ (0.77, 0.192, 0.038). Cette répartition particulière a en outre la propriété remarquable
que, sur une population initiale répartie de cette façon, la dynamique est exactement le comportement
asymptotique indiqué plus haut, à savoir une multiplication des effectifs par 2.
On peut écrire le modèle précédent en utilisant une notation matricielle de la façon suivante :
     
j(t + 1) 0 6 10 j(t)
 p(t + 1)  =  0, 5 0 0  ·  p(t) 
a(t + 1) 0 0, 4 0 a(t)
Si l’on introduit une notation vectorielle X(t) pour le vecteur colonne des effectifs des trois classes à
l’instant t, et un nom L pour cette matrice, la dynamique peut donc se réécrire d’une façon qui est très
semblable aux dynamiques linéaires d’une population à une seule classe :
X(t + 1) = L · X(t). (8.2)
La matrice L est un exemple de matrice de Leslie.
On appelle plus généralement matrice de Leslie toute matrice de la forme
 
f1 f2 f3 ... fn
 p1 0 0 ... 0 
 
 0 p2 0 . . . 0 
 
 ... ... ... ... ... 
0 . . . . . . pn−1 0
Elle permet de modéliser par la dynamique (8.2) une population structurée en n classes d’age : la première
ligne contient les coefficients de fertilité de chaque classe d’age f1 , f2 , ...fn et la sous diagonale les
probabilités (ou taux) de survie p1 , p2 , ...,pn−1 d’une classe d’age à la suivante. Les matrices de Leslie
ont tous leurs coefficients positifs ou nuls (mais elles ne sont pas pour autant des matrices stochastiques
car elles n’ont pas généralement la somme des coefficients de leurs lignes égale à 1).

8.2 Le théorème de Perron Frobenius

Théorème de Perron Frobenius
C’est le théorème de Perron Frobenius qui va nous permettre dans la plupart des cas de décrire les
dynamiques de Leslie. On dit qu’une matrice de Leslie est primitive lorsque l’une de ses puissances L, L2 ,
L3 , L4 , etc...a tous ses coefficients strictement positifs. C’est le cas de la matrice de l’exemple puisque sa
puissance L5 est à coefficients strictement positifs comme on peut le vérifier facilement.
Le théorème de Perron Frobenius affirme qu’une matrice primitive possède une valeur propre positive
strictement plus grande que toutes les autres valeurs propres que l’on appelle valeur propre dominante
λ∗ à laquelle est associé un vecteur propre X ∗ dit vecteur propre dominant dont tous les coefficients sont
positifs et de somme 1. De plus si X(0) est un vecteur initial dont tous les coefficients sont strictement
positifs, si X(t) = (x1 (t), x2 (t), . . . , xn (t) est sa dynamique, et si N (t) = x1 (t) + x2 (t) + . . . + xn (t)) la
somme de ses coefficients, on a les propriétés suivantes :
xi (t+1)
1. pour tout i = 1..n, xi (t) →t→∞ λ∗
X(t)
2. N (t) →t→∞ X ∗ .
Ce résultat important permet d’affirmer que si la matrice de Leslie d’un modèle dynamique (8.2)
est primitive, alors cette dynamique présente lorsque t augmente, un comportement asymptotique de
croissance exponentielle et la population se répartit selon une répartition particulière qui ensuite sera
invariante au cours du temps. De plus le calcul de ce taux de croissance et de cette répartition asymp-
totique se fait simplement en recherchant la valeur propre dominante λ∗ de la matrice de Leslie et un
vecteur propre X ∗ associé de somme 1.
8.2. LE THÉORÈME DE PERRON FROBENIUS 33

Exercice 1 : Une scientifique étudie une colonie de souris. Elle note qu’elles produisent en moyenne
une fille par femelle pendant leur première année de vie et 8 pendant leur seconde année. Elle note
aussi qu’elles sont seulement 25% à survivre une seconde année et aucune ne survivra au delà.
1. Ecrire le système dynamique modélisant cette population de souris en indiquant quelle est la
matrice de Leslie L du système. Cette matrice est-elle primitive ?
2. Pour une population initiale de 10 souris, toutes de la première classe d’âge, que pouvez-vous
dire de l’évolution du système ?
Exercice 2 : Pour la matrice suivante, calculer L2 et L3 et en déduire qu’elle n’est pas primitive.
Calculer les images successives Ln V d’un vecteur V quelconque par cette matrice.
 
0 1 0
L= 0 0 1 
1 0 0

Exercice 3 : On reprend le même exercice sur une population de saumons femelles que dans le TP
6, mais avec cette fois pour matrice
 
0 0 2000
L =  0, 005 0 0 
0 0, 1 0
1. Préciser quel est le sens des 3 coefficients non nuls de cette matrice.
2. Etant donnée une population initiale de 1000 individus pour chacune des 3 classes d’age,
calculer combien y aura-t-il d’ individus pour chaque classes d’age au temps t = 1, puis aux
temps t = 2 et t = 3. Qu’en déduisez-vous pour sa dynamique pour t ≥ 3 ?
3. Pensez-vous qu’on puisse appliquer le théorème de Perron-Frobenius à cette dynamique ?
Exercice 4 : A bird species has a maximum life span of 3 years. On average, each pair of birds in
their first year will produce 2 offsprings. A typical sample of 8 birds in their 2nd year will produce
a total of 15 more offsprings. only 40% of birds in their first year will survive to their 2nd year
and only 30% in their 2nd year will survive to their third year. Survival rates do not depend on
gender.
Describe how this bird population evolves over time. In particular, describe whether the population
remains stable, decreases or increases and at what rate. Furthermore, if possible, describe the
relative proportions of each age group after a number of years have passed.
34 CHAPITRE 8. DYNAMIQUES DE POPULATIONS STRUCTURÉES EN ÂGES
Chapitre 9

Systèmes différentiels en dimension n

Nous avons étudié jusqu’ici des équations différentielles et des systèmes de deux équations différentielles
pour lesquels nous avons présenté notamment la classification de Poincaré. Cette classification, et plus
généralement l’ensemble des méthodes d’étude que nous avons présentées, se généralisent à des systèmes
à n équations. Bien entendu en dimension n les dynamiques sont souvent beaucoup plus complexes ce
qui rend par exemple une étude qualitative souvent trop difficile. Nous allons voir que cela n’empêche
pas néanmoins l’étude de ces systèmes.

9.1 Systèmes dynamiques à n dimensions

La forme générale d’un système de n équations différentielles à n inconnues, notées x1 (t), x2 (t), ....,
xn (t), est la suivante  ′
x = f1 (x1 , x2 , . . . , xn )
 1′


x2 = f2 (x1 , x2 , . . . , xn )
(9.1)

 ... = ...................
 ′
xn = fn (x1 , x2 , . . . , xn )
Si l’on désigne par x le vecteur de composantes x1 , x2 , . . . , xn et par f : Rn → Rn la fonction de
composantes f1 , f2 , . . . , fn , on peut écrire vectoriellement le système x′ = f (x). A noter qu’une telle
fonction f donne, comme dans le cas n = 2, un champs de vecteurs dans Rn , puisqu’en chaque point
M = (x1 , x2 , . . . , xn ) de Rn , f (M ) est un vecteur de coordonnées (f1 , f2 , . . . , fn ).
Tout comme dans le cas n ≤ 2, un théorème d’existence et d’unicité des solutions permet d’assurer
que, si f est continuement dérivable, il existe, pour toute condition initiale M0 ∈ Rn , une unique solution
isuue de ce point.
Les équilibres (ou solutions stationnaires) sont des solutions constantes, c’est-à-dire telles que xi (t) =
x∗i = C ste pour tout i = 1, . . . , n. Ces équilibres sont donc les points M ∗ = (x∗1 , x∗2 , . . . , x∗n ) solutions du
système de n équations à n inconnues


 f1 (x1 , x2 , . . . , xn ) = 0
f2 (x1 , x2 , . . . , xn ) = 0

(9.2)

 ........ = ...
fn (x1 , x2 , . . . , xn ) = 0

.

9.2 Le cas linéaire

Lorsque la fonction f qui permet de définir le système différentiel est une fonction linéaire, de matrice
A = (aij )1≤i≤n,1≤j≤n , le système s’écrit
 ′
x = a11 x1 + a12 x2 + . . . + a1n xn
 1′


x2 = a21 x1 + a22 x2 + . . . + a2n xn
(9.3)

 .... = ...................
 ′
xn = an1 x1 + an2 x2 + . . . + ann xn

ou, plus simplement, x′ = Ax et on dit qu’il est linéaire.

35
36 CHAPITRE 9. SYSTÈMES DIFFÉRENTIELS EN DIMENSION N

Les systèmes linéaires sont importants puisqu’ils fournissent des approximations locales de la dyna-
mique d’un système non linéaire au voisinage de ses équilibres.
On peut vérifier facilement que l’ensemble des solutions d’un système (9.3) forment un espace vectoriel
de dimension n et donc, si l’on connait n solutions linéairement indépendantes, on peut calculer n’importe
quelle solution (comme combinaison linéaire de ces n solutions. Il est le plus souvent facile de calculer une
base de l’espace des solutions dès que l’on connait le spectre de A, c’est-à-dire l’ensemble de ses valeurs
propres. On a le résultat important suivant :

Proposition 2 Pour un système linéaire (9.3), on a les propriétés suivantes :

— Si toutes les valeurs propres de A sont de partie réelles strictement négatives, alors l’équilibre
(0, 0, . . . , 0) du système (9.3) est asymptotiquement stable, c’est-à-dire que toutes les solutions
convergent vers cet équilibre lorsque t tend vers l’infini. On dit qu’il s’agit d’un noeud stable
lorsque toutes les valeurs propres sont réelles. Dans ce cas, la convergence vers l’équilibre est
sans oscillations. Si certaines valeurs propres sont complexes la convergence peut comporter des
oscillations.
— Si au moins une valeur propre de A est de partie réelle strictement positive, la plupart des solutions
tendent vers l’infini (et non pas vers l’équilibre) quand t tend vers l’infini.
— Les valeurs propres imaginaires pures (de partie réelle nulles) conduisent à des comportements de
type centre, c’est-à-dire des oscillations périodiques autour de l’équilibre.

La preuve est facile puisqu’il est possible de calculer les solutions explicitement : en effet, si λ1 , λ2 , . . . , λn
sont les n valeurs propres de A (supposées distinctes) de vecteurs propres V1 , V2 , . . . , Vn , alors les n vec-
teurs eλ1 t V1 , eλ2 t V2 , . . . , eλn t Vn forment une base de l’ensemble des solutions.

9.3 Stabilité locale dans le cas non linéaire

Dans le cas d’un système non linéaire (9.1), il suffit en général, comme en dimension 1 ou 2, de linéariser
le système en ses points d’équilibre pour connaitre la nature locale de la dynamique au voisinage de ses
points d’équilibre. Rappelons que le linéarisé d’un système (9.1) au voisinage d’un équilibre M ∗ est donné
par le système x′ = Ax, où A est la matrice jacobienne de f calculée au point M ∗ . On a le résultat suivant :

Théorème 3 Si un système (9.1) possède un équilibre M ∗ , cet équilibre est asymptotiquement stable
(c’est-à-dire que toute solution issue d’un point suffisamment proche de l’équilibre tend vers l’équilibre
quand t tend vers l’infini) si et seulement si toutes les valeurs propres de son linéarisé en ce point M ∗
ont leurs parties réelles strictement négatives.

A noter que ce théorème donne des informations sur la dynamique du système non linéaire mais
elle n’est valable qu’au voisinage de ses équilibres. C’est pour cela qu’on appelle ce critère un critère de
stabilité locale.
Pour étendre ces informations locales au dela d’un petit voisinage de l’équilibre et montrer par exemple
la stabilité d’un équilibre dans une région D contenant l’équilibre, l’une des rares méthodes à notre
disposition est l’utilisation d’une fonction de Lyapounov.

Fonctions de Lyapounov
On appelle fonction de Liapounov du système (9.1) dans un domaine D ⊆ Rn contenant un équilibre
M , une fonction Φ : D → R, continument dérivable sur D qui présente un minimum au point M ∗
∗

et qui soit décroissante sur les trajectoires du système, c’est-à-dire telle que Φ(x(t)) soit une fonction
décroissante de t.
Comme on a
∂Φ ′ ∂Φ ′ ∂Φ ∂Φ
(Φ(x(t)))′ = x (t) + . . . + x (t) = f1 (x(t)) + . . . + fn (x(t)),
∂x1 1 ∂xn n ∂x1 ∂xn
cette décroissance s’éxprime par le fait que l’inégalité Grad(Φ) · f (x) ≤ 0 doit être satisfaite pour tout
x ∈ D. Géométriquement cela signifie que les trajectoires du système descendent les surfaces de niveau de
la fonction de Liapounov Φ, ce qui les conduit donc inéluctablement, lorsque t augmente, vers le minimum
de Φ dans le domaine D qui est l’équilibre M ∗ du système différentiel.

Exercice Le modèle suivant est appelé le modèle d’exclusion compétitive. Il modélise la compétition
entre n entreprises qui se partagent un marché. On désigne par x1 (t), x2 (t), .... , xn (t) la part de marché
détenue à l’instant t par chacunes d’elles. Le taux de croissance de chaque entreprise comporte une part
9.4. LOI DE CONSERVATION OU INTÉGRALE PREMIÈRE : 37

de croissance propre x′i (t) = βi xi (t) (βi représente le taux naturel de croissance propre de l’entreprise en
l’absence de compétiteurs) mais sa croissance est limitée (un peu comme dans un modèle logistique) par
la concurrence des autres à travers un terme x′i (t) = −γi Fi (x(t))xi (t) qu’on supposera ici linéaire pour
simplifier F (x) = α1 x1 + . . . + αn xn . Ceci conduit au système suivant :
 ′
 x1′ = (β1 − γ1 F (x)) x1

x2 = (β2 − γ2 F (x)) x2

(9.4)
 .... = ...................

 ′
xn = (βn − γn F (x)) xn
βi
Montrer que, si l’on suppose que les γi sont tous différents, une seule entreprise, celle qui a le plus
grand coefficient βγii , supposons que ce soit celle qui porte le numéro i = 1, va survivre, les n − 1 autres
étant conduites à la disparition. Cela se traduit par la présence d’un unique équilibre stable dont toutes
les coordonnées sont nulles sauf la première. On pourra vérifier pour cela que, sous les hypothèses faites,
la matrice jacobienne du système est une matrice triangulaire avec des coefficients nuls sous la diagonale
et des coefficients strictement négatifs sur la diagonale.

9.4 Loi de conservation ou intégrale première :

Le système (5.1) ayant un équilibre dont le linéarisé est un centre, rien ne permet, à priori, d’affirmer
que ses solutions sont périodiques, c’est-à-dire que l’équilibre (α1 /β1 , α2 /β2 ) est effectivement un centre
et non un foyer. Pour s’en assurer, on a recours à l’existence d’une intégrale première pour ce système.
Définition : Soient ((x(t), y(t)) la dynamique de deux espèces comme par exemple les proies et les
prédateurs du système de Lotka Volterra, on dit que la fonction H(x, y) est une intégrale première ou loi
de conservation de cette dynamique lorsque la quantité H(x(t), y(t)) reste constante au cours du temps.
Par exemple si l’on considère la dynamique de l’oscillateur harmonique, dont les solutions sont de
la forme (x(t), y(t)) = (r cos t, r sin t), la fonction H(x, y) = x2 + y 2 est une loi de conservation puisque
H(x(t), y(t)) = r2 .
Dans le cas du système de Lotka-Volterra (5.1), il n’est pas difficile de vérifier que la fonction suivante
est une loi de conservation :

H(x, y) = α1 ln y − β1 y + α2 ln x − β2 x.

Il suffit en effet de calculer le produit scalaire du gradient de H par le champs de vecteur (x′ (t), y ′ (t))
qui est donné par le système différentiel.
L’importance des lois de conservation pour l’étude des systèmes différentiels comme ceux de Lotka-
Volterra est facile à comprendre. Dès que la fonction H est connue, on peut, en utilisant ses dérivées
partielles, tracer ses courbes de niveau et en déduire les trajectoires de la dynamique. Alors qu’une
étude qualitative permet de prévoir l’oscillation des deux populations (car les trajectoires tournent dans
le plan (x, y)), elle ne permet pas de s’assurer que la dynamique est réellement périodique, c’est-à-
dire que les trajectoire se referment effectivement après un tour. Au contraire cette information découle
immédiatement de l’étude de H.

Exercices :
1. En prenant les valeurs particulières des paramètres (5.2), faire l’étude qualitataive du système,
calculer le linéarisé au voisinage de ses deux équilibres, en déduire la nature de l’un d’eux puis
vérifier que H est bien une intégrale première et en déduire la nature du deuxième.

Exercices :
1. Montrer que le système de Goodwin (5.12)possède effectivement des solutions périodiques et expli-
quez les oscillations du taux d’emploi de la main d’oeuvre et de la part des salaires dans le revenu
national durant un cycle économique.
2. Pour le système de Lotka Voltera perturbé (par un paramètre α) étudié en TP, calculer le signe
de la dérivée par rapport à t de la quantité H(x(t), y(t)) en fonction du paramètre α et en déduire
une nouvelle façon de prouver la stabilité de l’équilibre (la fonction H s’appelle une fonction de
Liapounov du système perturbé).
38 CHAPITRE 9. SYSTÈMES DIFFÉRENTIELS EN DIMENSION N
Chapitre 10

L’équation de Van der Pol

10.1 L’oscillateur linéaire

On appelle oscillateur linéaire l’équation différentielle du second ordre
ẍ + k ẋ + x = 0, où k est un nombre fixé. (10.1)
Cette équation est dite du second (ou deuxième) ordre car elle fait intervenir les deux premières dérivées
de la fonction inconnue x ; elle est dite linéaire car l’inconnue et ses dérivées interviennent linéairement
dans l’équation. Nous n’avons pour le moment que des équations et des systèmes du premier ordre mais,
en posant y = ẋ, on voit que l’équation (10.1) est équivalent au système différentiel linéaire

ẋ = y
(10.2)
ẏ = −x − ky
 
0 1
puisque ẍ = ẏ. La matrice de ce système différentiel linéaire est et ses valeurs propres sont
  −1 −k
donc 12 −k ± (k − 2)(k + 2) . Il en découle que le système (10.2)est une noeud ou un foyer, selon que
p

la valeur absolue de k est plus grande ou plus petite que 2, et surtout que l’équilibre (0, 0) est stable ou
instable selon que k est positif ou négatif. Le cas k = 0 est une centre : toute les solutions sont périodiques
et correspondent aux solutions x = r cos(t + φ), où r et φ sont des constantes arbitraires. Le cas k = 0
s’appelle de ce fait l’oscillateur linéaire.
Exercice : Démontrez ce qui vient d’être affirmé sur la nature du système (10.2) dans la classification
de Poincaré ainsi que sur les solutions lorsque k = 0.

10.2 Le cycle unique de l’équation de Van der Pol

L’inconvénient de l’oscillateur harmonique lorsque qu’on souhaite modéliser un phénomène périodique
est que son amplitude (r ou 2r, selon les auteurs) est fixée par la condition initiale, or souvent on observe
des phénomènes qui tendent vers un cycle-limite, c’est à dire un régime périodique vers lequel tendent
toutes les solutions, quelques que soient leur condition initiale. Notons que le système de Lotka-Volterra
de la relation proies-prédateurs présente également cette caractéristique. C’est Balthasar van der Pol qui
a introduit l’exemple suivant d’équation différentielle (non linéaire) comportant un cycle-limite :
ẍ + (x2 − 1)ẋ + x = 0 (10.3)
auquel on associe généralement le système de Van der Pol
x3

ẋ = y − F (x) , avec F (x) = −x
3 (10.4)
ẏ = −x
.
Exercice : Vérifier que le système (10.6) se déduit de l’équation de Van des Pol (10.3) en posant
y = ẋ + F (x).
L’idée heuristique de Van der Pol est que pour x voisin de 0 le facteur de ẋ est négatif ce qui rend
l’équilibre (0, 0) instable, et pour x grand il est positif ce qui force x à diminuer : ceci crée un cycle-limite.

39
40 CHAPITRE 10. L’ÉQUATION DE VAN DER POL

Figure 10.1 – Le cycle du système de Van der Pol. Le champ est vertical sur la cubique d’équation
y = F (x).

10.3 Une bifurcation de Hopf

Si l’idée heuristique qu’il y a un cycle-limite car le facteur de ẋ est négatif lorsque x est petit, on
devrait pouvoir moduler la taille du cycle en augmentant ou diminuant la région où ce facteur est négatif.
Ceci conduit à considérer l’équation de Van der Pol à paramètre suivante :

ẍ + (x2 − a)ẋ + x = 0 (10.5)

auquel on associe le système de Van der Pol à paramètre

x3

ẋ = y − Fa (x) , avec Fa (x) = − ax.
3 (10.6)
ẏ = −x
√ √
On voit que le facteur (x2 − 1) de ẋ est négatif sur l’intervalle ] − a, + a[ et donc à condition que
a > 0. Voici les cycles qu’on observe pour diverses valeurs de a ; on représente à chaque fois la cubique
d’équation y = Fa (x)

Figure 10.2 – Le cycle du système de Van der Pol. Le champ est vertical sur la cubique d’équation
y = Fa (x). Les valeurs de a sont, de gauche à droite et de haut en bas, 1.0, 0.5, 0.2, 0.1, 0, et −0.5. On
observe que, tant que a > 0, il subsiste un cycle, dont la taille diminue progressivement avec a. Ce cycle
est de moins en moins attractif et les autres trajectoires s’en distinguent de plus en plus. Lorsque a = 0
le cycle disparaı̂t est devient un équilibre stable, mais faible, au point que la portion de trajectoire tracée
(pour une durée de 50) reste clairement à l’écart de l’équilibre (0, 0). Lorsque a < 0 l’équilibre devient
plus attractif. Cette disparition d’un cycle-limite par diminution de sa taille et son remplacement par un
équilibre stable s’appelle une bifurcation de Hopf.

Exercice : Étudier, en fonction de la valeur du paramètre a, la nature dans la classification de Poincaré,

du linéarisé du champ (10.6) au voisinage de l’équilibre (0, 0).