COURS ÉCRIT

1 Le modèle de Lorenz et les systèmes dissipatifs . . . . . . . . . 2

2 Le pendule simple et les systèmes conservatifs . . . . . . . . . . 10
3 Autres oscillateurs et phénomènes non linéaires . . . . . . . . . 16
4 Exploration numérique de système dynamiques . . . . . . . . . 21

Introduction

L’objet de ce chapitre est de donner quelques éléments utiles pour la compré-

ension des instabilités observées en mécanique des fluides.
L’exemple du modèle de Lorenz est instructif à plusieurs titres. Il permet tout
d’abord de montrer que l’étude de stabilité d’un équilibre se ramène au calcul
des valeurs propres de la matrice du système linéarisé. Ensuite, l’examen de
la première bifurcation de ce modèle est cohérent avec le fait que la “bifur-
cation fourche” est la bifurcation stationnaire générique lorsqu’une symétrie
est brisée. Enfin, le comportement chaotique observé pour certaines valeurs
des paramètres permet de mettre en évidence la sensibilité aux conditions ini-
tiales dans un système déterministe, ce que l’on appelle l’effet papillon dans
le contexte de la météorologie dont est d’ailleurs issu ce modèle.
L’exemple du pendule simple permet d’introduire la famille des système dy-
namique dépendant d’un potentiel qui est un cas particulier de la famille
des système dynamiques hamiltoniens. Comme l’énergie et le volume dans
l’espace des phases est invariant, ces systèmes sont appelés système conser-
vatifs. Ils présentent la particularité de ne pas posséder d’attracteurs, c’est-
à-dire de sous-ensembles de l’espace des phases vers lesquels convergent les
trajectoires de leurs voisinages. Au contraire, les systèmes dissipatifs ont la
propriété de pouvoir exhiber de tels attracteurs.
La présentation suivante étudie les autres exemples classiques d’oscillateurs :

1
2 APM-INPT thu-sysdyn (2003), O. Thual July 10, 2005

pendule amorti, pendule résonant, oscillateur de Van der Pol et oscillateur

paramétrique. Ces exemples permettent de mettre en évidence un certain
nombre de phénomènes que l’on retrouve dans de nombreuses système physi-
ques sujets à des oscillations : portraits de phases en forme d’ellipses conver-
gentes ou divergentes, phénomène de résonance, bifurcation de Hopf, phéno-
mène de résonance paramétrique, exposants de Floquet, etc. La détermination
de l’amplitude des oscillations bifurquées au voisinage du seuil d’instabilité
est effectuée à l’aide d’une méthode asymptotique aux échelles multiples pour
l’oscillateur de Van der Pol et l’oscillateur paramétrique.

1 Le modèle de Lorenz et les systèmes dissipatifs

Le modèle de Lorenz est l’un des systèmes dynamique les plus simples per-
mettant d’illustrer la compétition entre un forçage extérieur et la dissipi-
sation interne d’un système physique, tout en exhibant des comportements
dynamiques non triviaux.
En 1963, Lorenz a étudié numériquement un système de trois équations
différentielles censé représenter très grossièrement la convection thermique
dans l’atmosphère [E. N. Lorenz, J. Atmos. Sci. 20 (1963) 130]. La démarche
permettant de passer des équations fluides à ce système est exposée dans de
nombreux ouvrages, en plus de l’article original de Lorenz et n’est pas abordée
ici. Seul les régimes dynamiques de ce modèle sont exposés ici.

1.1 Les équations du modèle de Lorenz

On se contente ici d’énoncer que modèle de Lorenz est le système dynamique

qui s’écrit

ẋ = −σ x + σ y
ẏ = −x z + r x − y
ż = x y − b z (1)

avec X(t) = [x(t), y(t), z(t)] ∈ IR3 , r ≥ 0, σ ≥ 0 et b ≥ 0.

L’espace des phases, constitué des vecteurs X ∈ IR3 , est de dimension 3. On
dit aussi que le système dynamique possède trois degrés de liberté.
3 . Le choix
L’espace de contrôle est constitué des vecteurs µ = (r, σ, b) ∈ IR+
de Lorenz dans son article de 1963 est σ = 10 et b = 8/3 avec r variable.
Physiquement, r est proportionnel au gradient thermique vertical imposé au
fluide, σ au nombre de Prandtl et b à l’élongation de la boı̂te contenant le
1. LE MODÈLE DE LORENZ ET LES SYSTÈMES DISSIPATIFS 3

fluide. La solution d’équilibre trivial x = y = z = 0 du système correspond

physiquement à un régime où le fluide est au repos et où la chaleur se transmet
uniquement par diffusion moléculaire. C’est l’état de base conductif. On va
voir que pour r suffisamment grand, cet équilibre est instable et qu’il laisse la
place à des régimes où le transfert de chaleur est réalisé par diffusion et par
convection.

1.2 Équilibres du modèle

Le modèle de Lorenz est donc un système dynamique de la forme Ẋ =

F (µ; X). Il est autonome dans la mesure où la fonction F ne dépend pas du
temps. Les équilibres, points fixes, points critiques ou encore solutions sta-
tionnaires de ce système sont les points X e = (xe , ye , ze ) tels que F (X e ) = 0.
La résolution du système F (X e ) = 0 conduit au tableau suivant de solutions :

0 ≤ r ≤ 1 : Une solution 1 ≤ r : Trois solutions

Xe = 0

X e ∈ {X − , X + }

p
avec X − = (−a, −a, r − 1), X + = (a, a, r − 1) et a = b(r − 1).
La stabilité d’un équilibre X e s’obtient en étudiant les valeurs propres de la
jacobienne DF (X e ) obtenue en linéarisant le système autour du point X e . Le
calcul des composantes ∂Fi /∂Xj de cette jacobienne conduit à
−σ σ 0
 

DF (X e ) =  r − ze −1 −xe  . (2)
ye xe −b

1.3 Stabilité de l’équilibre nul

La stabilité de l’équilibre trivial X e = 0 s’obtient en étudiant les valeurs

propres s de la matrice
−σ σ 0
 

DF (0) =  r −1 0  (3)
0 0 −b
dont l’équation caractéristique est
h i
(s + b) s2 + (σ + 1)s + σ(1 − r) = 0 . (4)
4 APM-INPT thu-sysdyn (2003), O. Thual July 10, 2005

sigma=10, b=2.67, r=0.3 sigma=10, b=2.67, r=0.3

0.3 0.3

0.25 0.25

0.2 0.2

z
0.15 0.15

0.1 0.1

0.05 0.05

0.2 0.2
0 0
−0.3 0 −0.3 0
−0.2 −0.2
−0.1 −0.1
0 −0.2 0 −0.2
0.1 0.1
0.2 0.2
0.3 0.3

y y
x x

Figure 1: Tracé de particule pour r = .3. a) t=0, b) t=.3, c) t=1.2. L’équilibre

0 est stable.

Une première valeur propre est s = −b. Le discriminant de l’équation du

second degré permettant de déterminer les deux autres est ∆ = (σ −1)2 +4σr.
Comme il est positif, il y a deux racines réelles de somme −(σ+1) et de produit
σ(1 − r). On peut donc résumer le signe des trois valeurs propres réelles de
DF (0) sur le tableau suivant :

0 ≤ r ≤ 1 : 0 est stable 1 ≤ r : 0 est instable

signe de s ∈ {−, −, −} signe de s ∈ {−, −, +}

Lorsque r ≤ 0, les trois valeurs propres sont négatives : l’équilibre X e = 0 est

donc stable. Lorsque r ≥ 0, une des valeurs propres est positive : l’équilibre
est donc instable. Pour r = 1, la stabilité est marginale.
Notons sb = −b la valeur propre indépendante de r, sc (r) la valeur propre
toujours négative et sd (r) la valeur propre qui devient positive pour r > 1.
Les vecteurs propres associés sont notés φb , φc (r) et φd (r). Ces vecteurs sont
définis à une constante près et on peut par exemple choisir φb = (0, 0, 1)
pour le vecteur propres associé à la valeur propre sb = −b. L’expression
analytique des autres valeurs propres et vecteurs propres est triviale mais
n’est pas explicitée ici par souci de concision.

1.4 Portraits de phase au voisinage de l’équilibre nul

Si nous notons U (t) les trajectoires X(t) = 0 + U (t) proches de 0, l’approxi-

mation de linéarisation indique qu’elles sont solutions du système dynamique
linéaire U̇ (t) = DF (0)U ce que l’on peut expliciter de la manière suivante :

u̇ = −σ u + σ v
v̇ = r u − v
1. LE MODÈLE DE LORENZ ET LES SYSTÈMES DISSIPATIFS 5

sigma=10, b=2.67, r=1

0.9

0.8

0.7

0.6

z
0.5

0.4

0.3

0.2

0.1
1
0 0.5
-1 0
-0.5
0 -0.5
0.5
1 -1

y
x

Figure 2: Vecteurs propres de la statibilité de 0 pour r = 1 : φb , φc (1) et

φd (1).

ẇ = −b w , (5)

où l’on a noté U = (u, v, w). Par rapport au système initial, on voit que l’on
a négligé les termes non linéaires −x z et x y. On peut alors poser

U (t) = B(t) φb + C(t) φc + D(t) φd (6)

où (B, C, D) sont les composantes de U dans la base propre de la matrice

DF (0). En reportant cette décomposition dans le système linéaire ci-dessus
gouverné par cette matrice on obtient trivialement les équations

Ḃ = sb B , Ċ = sc C , Ḋ = sb D , (7)

ce qui se résout en

B(t) = B(0)esb t , C(t) = C(0)esc t , D(t) = D(0)esd t . (8)

Ce résultat démontre le fait que 0 est stable lorsque les trois valeurs propres
sont négatives et qu’il suffit d’une valeur propre positive pour obtenir une
divergence des trajectoires dans la direction instable.
Il permet aussi de tracer le portraits de phase, c’est-à-dire l’ensemble de tra-
jectoires, dans le voisinage de 0. Le cas le plus intéressant à examiner est
obtenu pour r proche de la valeur critique r = 1. Pour cette valeur, les
vecteurs propres associés aux valeurs propres (sb , sc , sd ) = [−b, −(σ + 1), 0]
sont
0 σ 1
     

φb = 0 , φc (1) = −1 , φd (1) = 1  .
     (9)
1 0 0
Dans un plan affine engendré par les vecteurs [φb (1), φc (1)], les trajectoires
convergent toutes vers l’intersection de ce plan de coordonnées (B, C) avec
6 APM-INPT thu-sysdyn (2003), O. Thual July 10, 2005

0
0 0

sd s(+)
d s(+)
sd d
s(+)
-2
-2 -2
s
s(+) b b
sb b
-4
-4 -4

s(+)

s(+)
s

-6 -6 -6

-8 -8 -8

sc s(+) sc
-10 -10 c -10 s(+)
c

-12 -12 -12

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

r r r

Figure 3: Valeurs propres en fonction de r ∈ [0, 2]. a) Stabilité de 0 : sd (r),

(+) (+) (+)
sb (r) et sc (r) b) Stabilité de X (+) : sd (r), sb (r) et sc (r)
c) Superposition des deux graphes.

σ+1
l’axe φd le long des courbes d’équations C/C(0) = [B/B(0)] b . Cette con-
vergence est exponentielle.
On voit donc que toutes les trajectoires convergent vers l’axe engendré par
φd (1), au moins pour celles qui sont proches de 0. On montre que c’est en fait
le cas pour toutes les trajectoires, même éloignées de 0, mais il faut utiliser
pour cela un raisonnement (fonctionnelle de Liapounov) qui dépasse le cadre
de cette présentation.
Pour r légèrement inférieur à un, l’orientation des vecteurs propres est légère-
ment modifiée. En revanche, les trajectoires sur l’axe φd convergent vers zéro,
alors qu’elles étaient réduites à un continuum de points (tout l’axe) pour r
exactement à un. Toutes les trajectoires du système linéaire convergent donc
vers 0 qui est donc stable.
Pour r légèrement supérieur à un, les trajectoires divergent, lentement, sur
l’axe φd . Il arrive alors un moment où l’approximation linéaire n’est plus
valable, si bien que ces trajectoires quittent l’axe φd . Il est naturel de penser
que ces trajectoires vont converger vers l’un des équilibres bifurqués à partir
de la valeur critique r = 1. Il convient maintenant d’étudier la stabilité de
ces nouveaux équilibres.

1.5 Stabilité des équilibres bifurqués

p
La stabilité de l’équilibre X e = X + = (a, a, r−1) avec a = b(r − 1) s’obtient
en étudiant les valeurs propres s de la matrice

−σ σ 0
 

DF (X + ) =  1 −1 −a  (10)
a a −b
1. LE MODÈLE DE LORENZ ET LES SYSTÈMES DISSIPATIFS 7

dont l’équation caractéristique est

s3 + s2 (σ + b + 1) + s b(σ + r) + 2b(r − 1)σ = 0 . (11)

Suivant les valeurs des paramètres de contrôles (r, σ, b), ce polynôme de degré
trois peut avoir trois racines réelles
h ou bien une racineiréelle et deux racines
(+) (+) (+)
complexes conjuguées. Notons sb (r), sc (r), sd (r) ces trois racines. Le
signe de la partie réelle des trois racines de X e détermine la stabilité du point
d’équilibre X e . Sans entrer dans une discussion paramétrique complète sur
cette stabilité, nous allons nous concentrer ici sur les cas où r est dans le
voisinage de la valeur critique r = 1.
Comme X + = 0 pour r = 1, on peut faire le lien avec les trois valeurs propres
déterminant la stabilité de 0 et ordonner les notations de telle sorte que l’on
ait :
(+) (+)
sb (1) = −b , s(+)
c (1) = −(σ + 1) et sd (1) = 0 . (12)

En posant r = 1 + ǫ où ǫ > 0 est un petit paramètre, on peut donc effectuer

les développements limités suivants :
(+)
sb (1 + ǫ) = −b + β ǫ + O(ǫ2 )
s(+) 2
c (1 + ǫ) = −(σ + 1) + γ ǫ + O(ǫ )
(+)
sd (1 + ǫ) = δ ǫ + O(ǫ2 ) . (13)
(+) (+)
Pour ǫ suffisamment petit, les valeurs propres sb et sc restent réelles et
(+)
négatives. La stabilité de X e est dépend donc du signe de sd et donc du signe
(+)
de δ. En reportant le développement de sd dans le polynôme caractéristique
et en identifiant l’ordre dominant en ǫ, on voit facilement que δ = −2σ/(σ+1).
On peut donc affirmer que X + est stable lorsque r est dans un voisinage
de la valeur critique r = 1. L’étude de stabilité de X − conduit au même
polynôme caractéristique et donc à la même conclusion. En fait, la symétrie
S qui à X = (x, y, z) associe S X = (−x, −y, z) est une symétrie du système
dynamique. En effet, on voit que F (S X) = S F (X), si bien que si X(t)
est un solution, alors Y (t) = S X(t) en est une autre. On remarque que
l’équilibre X e = 0 vérifie S X e = X e , tandis que les deux équilibres bifurqués
à partir de la valeur r = 1 vérifient S X + = X − et S X − = X + . On dit
que la bifurcation obtenue pour r = 1 brise la symétrie S. L’ensemble de
ces considérations montre que l’on est en présence d’une bifurcation fourche
supercritique.
Les trajectoires dans le voisinage de l’équilibre X + s’obtiennent en traçant
tout d’abord les vecteurs propres de la matrice DF (X + ). Pour r proche de
8 APM-INPT thu-sysdyn (2003), O. Thual July 10, 2005

0
(+)
sd
-2
(+)
sb
-4
s(+)

-6

-8

s(+)
-10 c

r*
-12
0 5 10 15 20 25 30

Figure 4: Valeurs propres en fonction de r ∈ [0, 30]. Stabilité de X (+) :

(+) (+) (+)
sd (r), sb (r) et sc (r)

sigma=10, b=2.67, r=1.5 sigma=10, b=2.67, r=1.5 sigma=10, b=2.67, r=1.5

1.5 1.5 1.5

1 1 1
z

0.5 0.5 0.5

1 1 1
0 0 0
-1.5 0 -1.5 0 -1.5 0
-1 -1 -1
-0.5 -0.5 -0.5
0 -1 0 -1 0 -1
0.5 0.5 0.5
1 1 1
1.5 1.5 1.5

y y y
x x x

Figure 5: Tracé de particule pour r = 1.5. a) t=0, b) t=3, c) t=6. L’équilibre

X (+) est stable.
1. LE MODÈLE DE LORENZ ET LES SYSTÈMES DISSIPATIFS 9

1, ces vecteurs propres sont proches des vecteurs propres φb , φc (1) et φd (1)
que l’on a déjà explicités pour étudier les trajectoires autour de 0. Dans
un plan affine engendré par φb et φc , les trajectoires sont identiques à celles
obtenues dans le voisinage de 0. Sur l’axe φd , les trajectoires convergent
vers l’équilibre X + . Le portrait de phase dans le voisinage de X − et pour
r toujours proche de 1, s’obtient en traçant les trajectoires par symétrie par
rapport aux trajectoires du voisinage de X + .
Pour r proche de 1, on peut réunir les trajectoires du voisinage de 0 et les
trajectoires des voisinages de X + et X − en traçant une trajectoire partant
de 0 dans la direction du vecteur propre de DF (0) proche de φd et rejoignant
X + le long de la direction propre de DF (X + ) proche de φd . La trajectoire
symétrique relie 0 et X + de la même manière.

1.6 Autres régimes lorsque r devient grand

Lorsque r augmente au-delà du voisinage de r = 1, on peut chercher s’il existe

une valeur critique r = r∗ pour laquelle les équilibres X + et X − deviennent
instables (simultanément à cause de la symétrie). On rappelle que la stabilité
de ces équilibres dépend du signe de la partie réelle des racines du polynôme
caractéristique :

s3 + s2 (σ + b + 1) + s b(σ + r) + 2b(r − 1)σ = 0 . (14)

La déstabilisation de ces équilibres par changement de signe d’une valeur

propre réelle est impossible car si s = 0 on a forcément r = 1. On peut donc
rechercher pour quelles valeurs des paramètres on peut obtenir deux racines
iω et −iω de partie réelle nulle. En reportant la valeur s = iω dans l’équation
caractéristique, on obtient les deux conditions −ω 2 (σ + b+ 1)+ 2b(r − 1)σ = 0
et −iω 3 +iωb(σ +r) = 0. En éliminant ω 2 entre les deux équations, on obtient
finalement
σ(σ + b + 3)
r∗ = . (15)
(σ − b − 1)
Pour les valeurs classiques σ = 10 et b = 8/3, on obtient r∗ = 470/19 ∼ 24, 73.
On peut remarquer qu’entre le voisinage de r = 1 où les trois racines sont
réelles, et le voisinage de r = r∗ où il y a une racine réelle et deux racines
complexes conjuguées, il existe une valeur r = rd telle que le polynôme ad-
mette un racine double négative. Pour r légèrement supérieur à rd , il y a deux
valeurs propres complexes conjuguées de partie réelle négative et les portraits
de phase dans le voisinages de X + et X − sont des spirales elliptiques con-
vergeant exponentiellement vers ces équilibres stables.
10 APM-INPT thu-sysdyn (2003), O. Thual July 10, 2005

sigma=10, b=2.67, r=5 sigma=10, b=2.67, r=5 sigma=10, b=2.67, r=5

5 5 5

4.5 4.5 4.5

4 4 4

3.5 3.5 3.5

3 3 3
z

z
2.5 2.5 2.5

2 2 2

1.5 1.5 1.5

1 1 1

0.5 0.5 0.5

5 5 5
0 0 0
-5 0 -5 0 -5 0
0 0 0
5 -5 5 -5 5 -5

y y y
x x x

Figure 6: Tracé de particule pour r = 5. a) t=0, b) t=1, c) t=3. L’équilibre

X (+) est stable.
sigma=10, b=2.67, r=30 sigma=10, b=2.67, r=30 sigma=10, b=2.67, r=30

30 30 30

25 25 25

20 20 20
z

z
15 15 15

10 10 10

5 5 5

20 20 20
0 0 0
-30 0 -30 0 -30 0
-20 -20 -20
-10 -10 -10
0 -20 0 -20 0 -20
10 10 10
20 20 20
30 30 30

y y y
x x x

Figure 7: Tracé de particule pour r = 30 > r∗ . a) t=0, b) t=.1, c) t=2. Le

régime est chaotique.

La déstabilisation de des équilibres X + et X − correspond donc à une insta-

bilité oscillatoire et donc une bifurcation de Hopf. Deux valeurs propres com-
plexes conjuguées traversent l’axe des imaginaires lorsque le paramètre r fran-
chit la valeur critique r∗ . Dans le quart de plan (σ, b), il existe deux domaines
connexes dans lesquels cette bifurcation est supercritique, le complémentaire
correspondant au cas où la bifurcation est sous-critique. La détermination
analytique de ces domaines dépasse le cadre de ce cours (calcul de formes
normales). Indiquons simplement que les valeurs classiques σ = 10 et b = 8/3
correspondent au cas où la bifurcation est sous-critique.
Dans ce cas, la simulation numérique du système dynamique montre que l’on
transite abruptement vers un régime chaotique où la trajectoire tourne autour
d’un des équilibres instables X + et X − comme si elles y convergeait avant
de basculer aléatoirement vers l’autre équilibre pour y répéter le même type
de comportement. On observe que toutes les trajectoires convergent vers
une trajectoire chaotique que l’on nomme “attracteur étrange”. On montre
que la dimension fractale de cette trajectoire est supérieure à un et que la
distance entre deux conditions initiales initialement très proches s’amplifie
très rapidement. Cette sensibilité aux conditions initiales, ainsi que le chaos
déterministe observé avec ce système dynamique simple, ont servi de base
2. LE PENDULE SIMPLE ET LES SYSTÈMES CONSERVATIFS 11

à ce que l’on appelle “l’effet papillon”. Si l’évolution de l’atmosphère est

comparable à l’évolution d’un système dynamique comme celui de Lorenz,
la prise en compte ou non de l’effet du battement d’une aile de papillon
dans les conditions initiales peut changer radicalement l’évolution du temps
(atmosphérique) au bout d’un temps (en nombre de jours) estimé à environ
7 jours.

2 Le pendule simple et les systèmes conservatifs

Le pendule simple est le système dynamique prototype des modèles conser-

vatifs, qui décrivent des systèmes pour lesquels le forçage extérieur et la dissi-
pation interne sont nulles. C’est un cas particulier des systèmes dynamiques
dépendant d’un potentiel qui sont eux-mêmes un cas particulier des systèmes
hamiltoniens.

2.1 Équation du pendule simple

l
u

Figure 8: Pendule simple

On considère un pendule simple formé d’une masse m assimilée à un point

et suspendue à un fil de masse négligeable et de longueur l. En l’absence de
frottement, l’angle u(t) que fait le pendule avec la verticale vérifie l’équation
différentielle ordinaire
ü + ω 2 sin u = 0 . (16)
p
où ω = g/l et g est la gravité. Cette équation ne fait pas intervenir la
masse.
Cette équation différentielle de degré deux peut se transformer en un système
dynamique de degré un Ẋ = F (X) en posant X = (q, p) et en écrivant

q̇ = p
12 APM-INPT thu-sysdyn (2003), O. Thual July 10, 2005

ṗ = −ω 2 sin q . (17)

L’espace des phases de ce système est l’ensemble des X ∈ IR2 et l’espace de

contrôle est constitué de l’ensemble des valeurs du paramètre ω 2 ∈ IR+ .
On remarque l’existence d’une invariance qui est la translation T2π X = X +
(2π, 0). Si X(t) est une trajectoire, alors Y (t) = T2π X(t) = X(t) + (2π, 0)
en est une autre. On dit aussi que l’espace des phases est l’ensemble des
X ∈ T 1 × IR où T 1 est le cercle de périmètre 2π.

2.2 Stabilité de l’équilibre haut

Les deux équilibres de ce système, obtenus en résolvant F (X e ) = 0, sont

X b = (0, 0) et X h = (π, 0). Ces équilibres correspondent à la position basse
(X b ) et à la position haute (X h ) du pendule.
La jacobienne de l’application F est la matrice

0 1
 
DF (X) = . (18)
−ω 2 cos X1 0

La stabilité de X h est donc obtenue en calculant les valeurs propres de la

matrice
0 1
 
DF (X h ) = , (19)
ω2 0
c’est-à-dire en résolvant l’équation s2 − ω 2 = 0. Les valeurs propres sont donc
s = ω, à laquelle est associé le vecteur propre φ+ = (1, ω) et s = −ω à laquelle
est associé le vecteur propre φ− = (1, −ω). L’équilibre est donc instable, la
direction d’instabilité étant le long de l’axe φ+ . Le portrait de phase autour
de X h est constitué d’hyperboles dont les asymptotes sont données par les
directions φ+ et φ− .

2.3 Stabilité de l’équilibre bas

La stabilité de X b = 0 est obtenue en calculant les valeurs propres de la

matrice
0 1
 
DF (X b ) = (20)
−ω 2 0
c’est-à-dire en résolvant l’équation s2 + ω 2 = 0. Les valeurs propres sont donc
s = iω, à laquelle est associé le vecteur propre φ = (1, iω) et s∗ = −iω à
laquelle est associé le vecteur propre φ∗ = (1, −iω).
2. LE PENDULE SIMPLE ET LES SYSTÈMES CONSERVATIFS 13

L’étude de la stabilité de X b = 0 revient à poser X(t) = X b + U (t) et à

négliger les termes d’ordre U 2 . On se ramène donc à l’étude de l’équation du
pendule linéaire ü + ω 2 u = 0 qui s’écrit aussi sous la forme
u̇ 0 1 u
    
= . (21)
v̇ −ω 2 0 v

On peut indiquer directement la solution de l’équation du second ordre en

posant u(t) = a cos(ω t + ϕ) ou bien poser Z(t) = u + i ω1 v pour combiner
astucieusement les équations du système d’ordre un.
On obtient alors Ż = i ωZ qui se résout en Z(t) = Z(0)eiωt = Reiωt+iϕ .
Comme u(t) = Z(t) + Z ∗ (t) = 2R cos(ωt + iϕ), ces deux approches se re-
joignent en posant a = 2R. On a alors v(t) = u̇(t) = −a ω sin(ωt + iϕ). Les
trajectoires sont donc des ellipses d’axes parallèles aux axes des coordonnées
(u, v), engendrées par les vecteurs a(1, 0) et a(0, ω). Le sens de rotation est
contraire au sens trigonométrique.
Une autre manière de présenter la deuxième approche consiste à poser
Z(t) + Z ∗ (t) u
   
∗ ∗
U (t) = Z(t) φ + Z (t) φ = = . (22)
iω Z(t) − iω Z ∗ (t) v
En reportant dans le système linéaire et en utilisant le fait que φ et φ∗ sont
vecteurs propres, on obtient le système

Ż = i ω Z et Z˙∗ = −i ω Z ∗ . (23)

On remarque que les ellipses formées par les trajectoires sont engendrées
par des couples de vecteurs proportionnels aux vecteurs V r = Re(φ) et
V i = Im(φ). On vérifie facilement que cette propriété s’applique de manière
très générale aux portraits de phase associés à des valeurs propres complexes
conjuguées obtenues en étudiant la stabilité d’un équilibre.

2.4 Systèmes dépendant d’un potentiel

Le pendule simple est un cas particulier de la famille des équations différentielles

s’écrivant sous la forme ü = −V ′ (u) où V (u) est une fonction dérivable quel-
conque appelée potentiel. On peut poser par exemple V (u) = ω 2 (1 − cos u)
pour le pendule non-linéaire et V (u) = 21 ω 2 u2 pour le pendule linéaire appelé
aussi pendule harmonique.
On peut ramener ces équations à des systèmes de degré un Ẋ = F (X) en
posant X = (q, p) et en écrivant

q̇ = p
14 APM-INPT thu-sysdyn (2003), O. Thual July 10, 2005

25 25

20 20

15 15

10 10

5 5
p

p
0 0

-5 -5

-10 -10

-15 -15

-20 -20

-25 -25
-6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6

q q

Figure 9: a) Portrait de phase du pendule simple. b) Portait de phase du

pendule amorti pour ω = 2π et λ = 1.

ṗ = −V ′ (q) . (24)

Les équilibres sont des points X e = (qe , 0) tels que Xe est un extremum de la
fonction V . La stabilité de ces équilibres est obtenue en calculant les valeurs
propres de la matrice

0 1
 
DF (X e ) = , (25)
−V ′′ (qe ) 0

c’est-à-dire en recherchant les solutions de s2 + V ′′ (qe ) = 0. Si V ′′ (qe ) < 0,

la valeur qe correspond à un maximum du potentiel V . Les valeurs propres
sont alors s = ± −V (qe ) et l’équilibre est instable. Si V ′′ (qe ) = ωe2 > 0, la
p
′′

valeur qe correspond à un minimum du potentiel V . Les valeurs propres sont

alors s = ±iωe et l’équilibre est marginal.
En restant dans le cadre d’un potentiel V (u) quelconque, on considère la
fonction E(t) = 12 p(t)2 + V [q(t)] associée à une trajectoire [q(t), p(t)] solution
du système. On montre facilement que Ė(t) = p(t)ṗ(t) + q̇(t) V ′ [q(t)] =
−p(t) V ′ [q(t)] + p(t) V ′ [q(t)] = 0. On appelle énergie la quantité E(t) = E0
qui est donc un invariant du mouvement. Cette énergie s’écrit aussi E(t) =
1 2 1 2
2 u̇ + V (u) et peut être vue comme la somme de l’énergie cinétique 2 u̇ et
de l’énergie potentielle V (u).
En utilisant l’invariance de l’énergie on peut affirmer que les trajectoires sont
les courbes H(q, p) = 21 p2 + V (q) = E0 . Grâce à cette interprétation, on peut
tracer le portrait de phase d’un système défini par une fonction V quelconque.
Il suffit pour cela d’imaginer une bille qui roule le long d’une colline dont
l’altitude est le potentiel V .
2. LE PENDULE SIMPLE ET LES SYSTÈMES CONSERVATIFS 15

Figure 10: Sytème dépendant d’un potentiel. a) Forme d’un potentiel V (q).
b) Portrait de phase dans le plan (q, p).

2.5 Systèmes dynamiques hamiltoniens

Les systèmes dépendant d’un potentiel sont des cas particuliers de la famille
des systèmes hamiltoniens qui s’écrivent
∂H
q̇ = (q, p)
∂p
∂H
ṗ = − (q, p) (26)
∂q
où H(q, p) est une fonction dérivable quelconque.
En posant H(q, p) = 12 p2 + V (q) on retrouve le cas des systèmes dépendant
d’un potentiel. Mais on peut considérer d’autres formes de l’Hamiltonien
H(q, p) comme par exemple H = q p.
Étant donnée une solution [q(t), p(t)], on montre que la quantité E(t) =
H[q(t), p(t)] est encore un invariant du mouvement. En effet,
∂H ∂H ∂H ∂H ∂H ∂H
Ė = q̇ + ṗ = − =0 (27)
∂q ∂p ∂p ∂q ∂q ∂p

On montre de plus que les aires sont conservées dans l’espace des phases. En
effet, si on note Ẋ = F (X) avec X = (q, p) le système dynamique hamiltonien,
on peut calculer la divergence
∂2H ∂2H
div F (X) = − =0. (28)
∂q ∂p ∂p ∂q
Cette propriété entraı̂ne qu’un volume (ici une aire) de conditions initiales
reste constant au cours du temps, comme dans le cas d’un écoulement isochore.
On parle alors de systèmes conservatifs. De tels systèmes n’admettent pas
d’attracteurs, c’est-à-dire d’ensembles vers lesquelles toutes les trajectoires
d’un bassin d’attraction pourraient converger.
16 APM-INPT thu-sysdyn (2003), O. Thual July 10, 2005

L’invariance de l’hamiltonien H et des aires dans l’espace des phases se

généralise au cas des systèmes hamiltoniens H(q, p) dans IR2n dont l’expression
générale est
∂H
q̇i = (q, p)
∂pi
∂H
ṗi = − (q, p) (29)
∂qi
pour i = 1, ..., n avec q ∈ IRn et p ∈ IRn .

3 Autres oscillateurs et phénomènes non linéaires

On considère maintenant des exemples d’oscillateurs qui prennent en compte

les phénomènes de frottement ainsi que des forçages externes. Contrairement
aux systèmes conservatifs, ces systèmes sont dissipatifs. Dans certains cas,
il s’établit un équilibre (au sens large) entre le forçage et la dissipation. Cet
équilibres sont parfois des attracteurs dans la mesure où un ensemble de con-
ditions initiales peuvent y converger. La particularité des exemples ci-dessous
est le fait que ces équilibres se traduisent par des régimes oscillant en temps.

3.1 Le pendule amorti

Lors de l’oscillation du pendule, le frottement de l’air peut se modéliser par

une force proportionnelle à la vitesse. On obtient donc l’expression

ü + 2λ u̇ + ω 2 sin u = 0 . (30)

Le paramètre λ > 0 est donc un facteur d’amortissement du mouvement du

pendule.
Les équilibres sont les mêmes que pour le pendule non amorti. Leur stabilité
et le portrait de phase de leur voisinage s’obtient en linéarisant les équations
autour de la position haute (instable) ou de la position basse (stable).
Une manière rapide d’étudier la stabilité de la position basse u = 0 consiste
à chercher des solutions complexes de la forme u(t) = um est avec um ∈ C I et
s∈C 2
I dans l’équation linéarisée autour de zéro qui s’écrit ü + 2λ u̇ + ω u = 0.
Comme l’équation est linéaire, la partie réelle d’une solution complexe est une
solution. On doit alors chercher les racines de l’équation s2 + 2λs + ω 2 = 0.
Lorsque λ < ω, il existe deux racines complexes conjuguées
q
s = −λ ± i ω 1 − λ2 /ω 2 . (31)
3. AUTRES OSCILLATEURS ET PHÉNOMÈNES NON LINÉAIRES 17

Le retour à l’équilibre s’effectue par des oscillations exponentiellement amor-

ties. Lorsque λ est petit, les oscillations sont faiblement amorties et la pulsa-
tion est proche de ω.
Lorsque λ > ω, il existe deux racines réelles
 q 
s = −λ 1 ± 1− ω 2 /λ2 . (32)

Le retour à l’équilibre exponentiellement sans oscillation. On parle alors de

“sur-relaxation”.
L’étude de stabilité autour de la position haute fait apparaı̂tre uniquement des
valeurs propres réelles dont l’une est positive. En connectant les portraits de
phase des voisinages des deux équilibres, on peut tracer toutes les trajectoires
de l’espace des phases IR2 de coordonnées (u, u̇).

Figure 11: portrait de phase du pendule amorti pour T = 2π/ω = 1 et λ = 1.

3.2 Le pendule résonant

On peut maintenant considérer le cas où les oscillations du pendule sont en-
tretenues en exerçant sur la masse ponctuelle une force oscillant à la pulsa-
tion ωe (en imposant par exemple un champ magnétique lorsque la masse est
aimantée). Les équations du mouvement peuvent alors s’écrire

ü + 2λ u̇ + ω 2 sin u = f cos(ωe t) . (33)

Cette équation différentielle ordinaire de degré deux peut se mettre sous la

forme d’un système dynamique de degré un Ẋ = F (X, t) avec X ∈ IR2 . Ce
système n’est pas autonome dans la mesure où F dépend du temps.
La position basse u = 0 n’est plus un équilibre du système dynamique si
f 6= 0 mais on peut quand même s’intéresser aux trajectoires situées dans le
voisinage de X = 0, surtout si le forçage f est un petit paramètre. Dans ce
cas, ces trajectoires sont solutions de l’équation linéarisée

ü + 2λ u̇ + ω 2 u = f cos(ωe t) . (34)

Sauf cas particulier discuté ci-dessous, les solutions de cette équation linéaire
(non homogène) sont de la forme u(t) = um eiωe t + u∗m e−iωe t où um ∈ C I et
s∈C I sont des coefficients complexes. En reportant dans l’équation linéaire
on obtient
(ω 2 − ωe2 + 2 i λ ωe )um = f /2 . (35)
18 APM-INPT thu-sysdyn (2003), O. Thual July 10, 2005

En notant ∆ = ωe2 /ω 2 et Λ = λ2 /ω 2 , l’amplitude de la réponse du pendule

au forçage est donc
f 1
|um | = . (36)
2ωn2 (1 − ∆)2 + 4 Λ ∆
p

Cette réponse est maximale pour ∆ = 1−2Λ et s’écrit alors |um | = 2ωf 2 4 Λ(1−Λ)
1
.
Lorsque Λ est très petit (faible dissipation) la réponse est très grande pour
pour pour ∆ ∼ 1 c’est-à-dire ωe ∼ ω. On est en présence du phénomène de
résonance.
Lorsque λ = 0 et ω = ωe = ω, c’est-à-dire pour le cas de la résonnance, on
iω t t
∗ −iω


cherche une solution sous la forme u(t) = t vm e + vm e avec vm ∈ C I.
f f
On trouve alors vm = −i 4ω ce qui conduit à la solution u(t) = 2ω t sin(ω t).
La réponse du pendule croı̂t donc linéairement avec le temps à la résonance.
Lorsque les trajectoires s’éloignent du voisinage de u = 0, l’approximation
linéaire n’est plus valable et il faut tenir compte de la nonlinéarité contenue
dans le terme sin u. La réponse du pendule lorsque ∆ varie et pour différentes
valeurs de f peut être obtenu numériquement ou par un développement
asymptotique lorsque f est petit. On peut alors observer le phénomène
d’hysteresis correspondant au fait que deux réponses sont possibles pour une
même valeur des paramètres de contrôle.

3.3 L’oscillateur de Van der Pol

Une autre façon de forcer le pendule peut être obtenue en imposant une
force proportionnelle à la vitesse et dans le même sens, au moins pour les
petites oscillations. On peut imaginer une telle force sur un pendule en util-
isant l’induction magnétique que provoque la masse ponctuelle lorsqu’elle est
électriquement chargée. Mais les équations du pendule avec une dissipation
négative pour les petites oscillations apparaissent naturellement dans un mon-
tage électronique portant le nom d’oscillateur de Van der Pol. Ces équations
s’écrivent
ü + (k u2 − 2 µ) u̇ + ω 2 sin u = 0 . (37)
La stabilité de l’équilibre u = 0 s’obtient en étudiant la partie réelle des
valeurs propres du système linéarisé qui s’écrivent s = µ ± iω. La position
basse u = 0 est donc stable pour µ < 0 (on retrouve la notion de dissipation
en posant µ = −λ) et instable pour µ > 0. La valeur critique de la bifurcation
est µ = 0. On s’attend alors à être en présence d’une bifurcation de Hopf.
Il est possible de calculer l’amplitude de saturation de l’instabilité par les
termes nonlinéaires au voisinage de la bifurcation µ = 0, grâce à des déve-
loppements asymptotiques où µ est le petit paramètre. On peut par exemple
3. AUTRES OSCILLATEURS ET PHÉNOMÈNES NON LINÉAIRES 19

effectuer un développement en échelles multiples en posant µ = M ǫ et

√ h i
u(t) = ǫ u0 (t, ǫ t) + ǫ u1 (t, ǫ t) + o(ǫ2 ) (38)

où les fonctions ui (t, T ) sont considérées comme des fonctions de deux vari-
ables indépendantes et doivent rester d’ordre 1. Ce choix de développement
est conforme à ce que l’on attend d’une bifurcation de Hopf où une oscil-
lation sur une échelle de temps d’ordre 1 coexiste avec une amplification
sur une échelle de temps inversement proportionnelle à l’écart au seuil et où
l’amplitude de l’oscillation bifurquée croı̂t comme la racine carrée de l’écart
au seuil.
Il convient alors de calculer les expressions suivantes :
√ ∂u0 ∂u0 ∂u1
  
  
u̇ = ǫ +ǫ + + o ǫ2
" ∂t ∂T ∂t ! #
√ ∂ 2 u0 2
∂ u0 ∂ 2 u1  
2
ü = ǫ +ǫ + +o ǫ
∂t2 ∂t ∂T ∂t2
1 3 √ 1 3
    
4
sin u = u − u + o(u ) = ǫ u0 + ǫ u1 − u0 + o ǫ2
h 6 i 6
2 2
u = ǫ u0 + o(ǫ) . (39)
√ 2
En reportant dans l’équation, l’ordre dominant en ǫ s’écrit ∂∂tu20 + ω 2 u0 = 0
ce qui conduit à la solution u0 (t, T ) = A(T )eiω t + A∗ (T )e−iω t où A(T ) est un
√
fonction complexe arbitraire de la variable T . À l’ordre suivant, en ǫ ǫ, on
obtient l’équation

∂ 2 u1 ∂ 2 u0 ∂u0 ∂u0 1 3
2
+ ω 2 u1 = −2 − k u20 + 2M + u0
∂t ∂t ∂T ∂t ∂t 6
n=3
Cn (T ) ei n ω t
X
= (40)
n=−3

où les fonctions Cn (T ) s’expriment en fonction de A(T ) et de ses dérivées.

Pour que u1 (t, T ) reste d’ordre 1, comme supposé dans le développement,
il est nécessaire d’imposer C1 (T ) = C−1 (T ) = 0. En effet, en l’absence de
cette condition dite de compatibilité, le terme u1 serait proportionnel à t et
donc non borné. On parle aussi de condition de résonance. En calculant
l’expression de C1 (T ), cette condition s’écrit donc
dA k 3
 
=M A− + i |A|2 A . (41)
dT 3 2ω
√
En posant Z(t) = ǫ A(ǫ t) eiωt , on obtient finalement la solution
 
u(t) = Z(t) + Z ∗ (t) + O |Z|3 (42)
20 APM-INPT thu-sysdyn (2003), O. Thual July 10, 2005

où l’amplitude complexe Z(t) est solution de l’équation

k 3
Ż = (µ + iω) Z − α|Z|2 Z avec α= + i. (43)
3 2ω
On voit donc que la bifurcation de Hopf est supercritique pour k > 0 et
sous-critique pour k < 0.

3.4 L’oscillateur paramétrique

Une autre façon d’exciter un pendule consiste à faire osciller sa longueur

ou bien le champ de gravité apparent en faisant osciller verticalement son
point d’attache. On note ωe la pulsation de ces oscillations de forçage dit
paramétrique dans la mesure où il modifie les paramètre de contrôle du pend-
ule simple. Lorsque ce forçage est petit, l’équation du pendule paramétrique
s’écrit :
ü + 2λ u̇ + ω 2 [1 + h cos(ωe t)] sin u = 0 . (44)
L’espace de contrôle de ce système est constitué des paramètres (λ, ωωe , h) ∈
3 . En effet, un changement d’unité de temps permet de voir que seul le
IR+
rapport ω/ωe est pertinent pour l’étude des solutions.
Cette équation différentielle ordinaire de degré deux se ramène à un système
dynamique d’ordre un de la forme Ẋ = F (X, t) avec X ∈ IR2 qui est
dépendant du temps. La stabilité de l’équilibre X = 0 s’obtient en linéarisant
les équations ce qui conduit ici au système linéaire

0 1
 
U̇ = L(t) U avec L(t) = 2 . (45)
−ω h cos(ωe t) −2 λ

L’étude des systèmes dynamiques linéaires U̇ = L(t)U avec U ∈ IRn et où

L(t) est une matrice réelle à coefficients périodiques de période Te est appelée
problème de Floquet. On montre que l’espace vectoriel des solutions com-
plexes est engendré par n solutions de la forme φj (t) esj t pour j = 1, ..., n où
les φj (t) sont des fonctions complexes périodiques de période Te et les sj des
nombres complexes appelés exposants de Floquet.
La démonstration de ce résultat est simple en considérant la matrice R telle
que U (Te ) = R U (0). On peut construire numériquement cette matrice en
intégrant le système dynamique linéaire de 0 à Te à partir d’un ensemble de
conditions initiales formant une base de IRn . Si on note rj pour j = 1, ..., n
les valeurs propres complexes de R, les coefficients de Floquet sont obtenus
en posant sj = T1e Ln rj pour j = 1, ..., n. La propriété U (t + Te ) = R U (t)
permet de conclure en construisant les fonctions φj (t) à partir des solutions
3. AUTRES OSCILLATEURS ET PHÉNOMÈNES NON LINÉAIRES 21

intégrées entre 0 et Te en prenant successivement les n vecteurs propres de R

comme conditions initiales.
Même pour le cas simple du pendule paramétrique (on parle d’équations de
Mathieu), la détermination analytique des exposants de Floquet fait l’objet
de développements mathématiques complexes (déterminants infinis, ...). Dans
l’espace de contrôle (λ, ωωe , h), on montre que l’équilibre se déstabilise prin-
cipalement pour ω/ωe ∼ 1/2, sur un intervalle d’autant plus large que h est
grand, à condition de dépasser une valeur une valeur critique d’autant plus
grande que la dissipation λ est grande. En plus de la résonnance paramétrique
ω/ωe ∼ 1/2, on observe les résonances ω/ωe ∼ p avec p ∈ IN ∗ , mais sur des
domaines de paramètres d’autant plus restreints que p est grand.
On peut effectuer un développement asymptotique lorsque l’intensité h du
forçage, le coefficient λ de la dissipation et le detuning δ = ω 2 /ωe2 − 1/4 sont
petits. On suppose que h = ǫ F , λ = ǫ Λ et δ = ∆ǫ
√ h i
u(t) = ǫ u0 (t, ǫ t) + ǫ u1 (t, ǫ t) + o(ǫ2 ) (46)

comme précédemment. À l’ordre dominant, on obtient u0 (t, T ) = A(T )eiω0 t +

A∗ (T )e−iω0 t avec ω0 = ωe /2. La réponse rapide du pendule s’effectue donc à
la pulsation moitié de celle du forçage, c’est-à-dire avec une période double.
L’ordre suivant nécessite d’imposer une relation de compatibilité qui conduit
à l’équation d’amplitude
dA
= (−Λ + i ν)A + µ A − iβ|A|2 A (47)
dT
avec ν = 2∆ ω0 , µ = ω0 F/4 et β = ω0 /4. L’étude de cette équation, en posant
par exemple A = R exp(iθ), permet de tracer le diagramme de stabilité dans
le plan (ν, F ) pour différentes valeurs de Λ. L’instabilité est supercritique
pour ν > 0 et sous-critique pour ν < 0.

6
F

instable
sous−
critique

stable

-
0 ν
Figure 1: stabilité de l’équilibre du pendule paramétrique
22 APM-INPT thu-sysdyn (2003), O. Thual July 10, 2005

4 Exploration numérique de système dynamiques

L’objectif de ce paragraphe est de motiver la programmation de systèmes

dynamiques pour observer numériquement la richesse et l’universalité des
comportements identifiés par l’approche théorique. Dans cet esprit, nous
rappelons quelques méthodes numériques simples pour l’approximation des
équations d’évolution. Nous énumérons ensuite une liste de systèmes dy-
namiques qui offre un large champ d’expérimentations numériques. Cette
liste est complétée par les exercices et problèmes présentés à la fin.

4.1 Méthodes numériques simples

On considère un système d’équations différentielles ordinaires couplées de la

forme
Ẋ = F (X) où X(t) ∈ IRN . (48)
En notations indicées ce système s’écrit

Ẋi = Fi (X1 , X2 , ..., XN ) pour i = 1, 2, ..., N . (49)

On suppose que F est différentiable. Le théorème de Cauchy indique alors

qu’une condition initiale X(0) donne naissance à une et une seule solution
X(t).
On présente ici quelques schémas d’approximation des solutions de ces systèmes
en discrétisant le temps en pas de temps égaux à ∆t. On considère alors une
succession d’états X 0 , X 1 , X 2 , ..., X n , ... qui doivent approximer une solution
X(t) aux temps tn = n ∆t.
Le schéma d’Euler consiste à effectuer l’itération
X n+1 − X n
= F (X n ) (50)
∆t
à partir de la condition initiale X 0 = X(0).
Le schéma Euler arrière s’écrit
X n+1 − X n
= F (X n+1 ) (51)
∆t
et nécessite donc de résoudre une équation implicite pour X n+1 à chaque pas
de temps.
Le schéma Leap-Frog (saute-mouton) est défini par l’équation explicite
X n+1 − X n−1
= F (X n ) (52)
2∆t
4. EXPLORATION NUMÉRIQUE DE SYSTÈME DYNAMIQUES 23

mais il fait intervenir les deux pas de temps précédant le nouveau pas de temps
calculé. Il doit donc être initialisé avec un schéma à un pas de temps (par
exemple Euler). Il est bon de savoir qu’il a tendance à découpler les pas de
temps pairs et impairs, ce qui se traduit par des oscillations numériques que
l’on peut éliminer en mélangeant, à intervalles réguliers ou systématiquement,
les pas de temps pairs ou impairs.
Le schéma Adams-Bashforth fait lui aussi intervenir deux pas de temps,
sans présenter le phénomène d’oscillation numériques du schéma précédent.
Il s’écrit
X n+1 − X n 3 1
= F (X n ) − F (X n−1 ) . (53)
∆t 2 2

Le schéma “prédicteur-correcteur” fait intervenir une étape intermédiaire

issue du schéma d’Euler. Il s’écrit
XI − Xn X n+1 − X n
= F (X n ) et = F (X I ) . (54)
∆t ∆t

Un schéma bien connu pour sa précision est le schéma de Runge-Kutta

quatrième ordre qui s’écrit

XI − Xn X II − X n
= F (X n ) = F (X I )
∆t ∆t/2
X III − X n X IV − X n
= F (X II ) = F (X III )
∆t/2 ∆t
∆t h i
X n+1 = Xn + F (X I ) + 2F (X II ) + 2F (X III ) + F (X IV ) (55)
6

Notons que l’on peut combiner plusieurs de ces schémas pour une même
équation. Si on décompose par exemple F (X) = A X + G(X) où A est
un opérateur linéaire, le schéma “Leap-Frog / Crank-Nicolson” s’écrit

X n+1 − X n−1 X n+1 + X n−1

=A + G(X n ) (56)
2∆t 2

Tous ces schémas peuvent se mettre sous la forme X n+1 = C(X n , X n−1 ; ∆t),
après avoir résolu les équations implicites lorsque c’est le cas. On peut alors
rappeler sommairement les notions de consistance, stabilité et convergence
d’un schéma.
Un schéma est consistant si

Lim kX(t + ∆t) − C [X(t), X(t − ∆t); ∆t]k = 0

∆t → 0 (57)
24 APM-INPT thu-sysdyn (2003), O. Thual July 10, 2005

pour une solution X(t). Il est précis à l’ordre p si cette quantité tend vers
zéro comme (∆t)p .
Un schéma est stable si pour tout temps T il existe une borne supérieure M
finie telle que kX n k < M pour tous les ∆t et tous les n tels que n∆t ≤ T .
Enfin, un schéma est convergeant si pour toutes les conditions initiales
l’approximation X n converge vers la solution X(n∆t) quand ∆t tend vers
zéro et n vers l’infini avec t = n∆t fixé.
Le théorème de Lax indique que pour un problème bien posé et pour un
schéma consistant, la stabilité est équivalente à la convergence.
On montre que les schémas d’Euler et Euler arrière sont précis au premier
ordre, que les schémas de Crank-Nicolson, Leap-Frog et prédicteur correcteur
sont précis au second ordre, et que le schéma de Runge-Kutta est précis au
quatrième ordre.
L’analyse de la stabilité des schémas est du ressort de l’étude des systèmes
dynamiques discrets, encore appelés “mappings”.

4.2 Le zoo des systèmes dynamiques continus

On présente ici plusieurs spécimens de systèmes dynamiques c’est-à-dire de

systèmes d’équations différentielles ordinaires couplées. Ces exemples peu-
vent servir de champ expérimental pour des simulations numériques ainsi que
d’illustration des concepts qui seront développés dans ce cours.
Cage des formes normales
On présente ici quatre systèmes dynamiques à un degré de liberté, c’est-à-
dire ne comprenant qu’une seule équation différentielle ordinaire d’ordre 1.
Ces équations sont étudiées en détail dans le chapitre intitulé “Bifurcations
génériques de l’équilibre”. Ce sont des formes normales dans la mesure où
elles sont représentatives de bifurcations de l’équilibre observées dans des
systèmes dynamiques plus complexes pouvant avoir un grand nombre de
degrés de liberté. La théorie des formes normales montrent qu’en effectu-
ant un développement asymptotique au voisinage de la bifurcation on ramène
la dynamique du système à celle de l’une ces formes normales.

1. Bifurcation noeud-col : x(t) ∈ IR

ẋ = µ + α x2 . (58)

2. Bifurcation fourche : x(t) ∈ IR

ẋ = µ x + α x3 . (59)
4. EXPLORATION NUMÉRIQUE DE SYSTÈME DYNAMIQUES 25

3. Bifurcation de Hopf : z(t) ∈ C

I

ż = (µ + i ω)z + (α + i β)|z|2 z . (60)

4. Catastrophe en fronce : x(t) ∈ IR

ẋ = −q + p x − x3 . (61)

Cage des vedettes du chaos

On présente maintenant deux systèmes dynamiques à trois degrés de liberté
souvent utilisés pour illustrer le concept de chaos déterministe. En simu-
lant numériquement ces systèmes pour certaines valeurs des paramètres de
contrôle, on montre que le comportement des solutions est erratique, sensible
aux conditions initiales et que les trajectoires simulées décrivent un ensemble
que l’on appelle souvent “attracteur étrange”. L’étude du premier exemple
fera l’objet du chapitre “Modèle de Lorenz”.

1. Modèle de Lorenz : X(t) = [x(t), y(t), z(t)] ∈ IR3

ẋ = −σ x + σ y avec σ ≥ 0
ẏ = −x z + r x − y avec r ≥ 0
ż = xy−bz avec b ≥ 0 . (62)

Le choix de Lorenz dans son article de 1963 est σ = 10 et β = 8/3. Pour

r suffisamment grand, ce système a un comportement chaotique.

2. Modèle de Rossler : X(t) = [x(t), y(t), z(t)] ∈ IR3

ẋ = −y − x
ẏ = x+ay
ż = b + x z − c z (63)

Un choix de paramètres conduisant à un comportement chaotique est

a = b = 0.2 et c = 5.

Cage des oscillateurs

On présente un certain nombre de système dynamique exprimés sous la forme
d’une équation différentielle d’ordre 2. En posant v = u̇, on se ramène à un
système dynamique d’ordre 1 à deux degrés de libertés. Ces systèmes ont
souvent un comportement oscillant, et décrivent donc des oscillateurs. Leur
étude fait l’objet du chapitre intitulé “Oscillateurs”.
26 APM-INPT thu-sysdyn (2003), O. Thual July 10, 2005

1. Pendule simple : u(t) ∈ IR

ü + ω 2 sin u = 0 . (64)

2. Oscillateur de Duffing : u(t) ∈ IR

ü + ω 2 u = ǫ u3 . (65)

3. Pendule linéaire amorti : u(t) ∈ IR

ü + 2λ u̇ + ω 2 u = 0 . (66)

4. Pendule tournant : u(t) ∈ IR

ü + 2λ u̇ + (ωn2 − ωe2 cos u) sin u = 0 . (67)

5. Pendule résonant : u(t) ∈ IR

ü + 2λ u̇ + ωn2 sin u = f cos(ωe t) . (68)

6. Oscillateur de Van der Pol : u(t) ∈ IR

ü + (k u2 − 2 µ) u̇ + ωn2 u = 0 . (69)

7. Oscillateur forcé : u(t) ∈ IR

ü + 2λ u̇ + ωn2 sin u = f cos(ωe t) . (70)

8. Oscillateur paramétrique : u(t) ∈ IR

ü + 2λ u̇ + ωn2 [1 + h cos(ωe t)] sin u = 0 . (71)

Conclusion

Plusieurs éléments de base dans l’étude des systèmes dynamiques ont été
abordé à travers l’examen du modèle de Lorenz, du pendule simple et autres
oscillateurs. La linéarisation de ces systèmes autour de leurs équilibre a été ef-
fectuée systématiquement. L’étude du système linéaire ainsi obtenu s’effectue
à l’aide des vecteurs propres et valeurs propres de la matrice jacobienne. On
obtient ainsi la stabilité de l’équilibre étudié ainsi que la topologie des trajec-
toires autour de cet équilibre.
Les changements de stabilité et de topologie des trajectoires obtenus en faisant
varier un paramètre de contrôle sont appelés des bifurcations. Nous avons
4. EXPLORATION NUMÉRIQUE DE SYSTÈME DYNAMIQUES 27

ainsi observé une “bifurcation fourche” lors de la déstabilisation de l’équilibre

trivial du modèle de Lorenz et une “bifurcation de Hopf” pour la déstabilisation
du repos dans l’oscillateur de Van der Pol. Ces bifurcations sont propres aux
systèmes dissipatifs qui modélisent la compétition entre un forçage extérieur
et une dissipation interne.
Les systèmes conservatifs modélisent des systèmes physiques isolés de tout
forçage et pour lesquels on néglige la dissipation d’énergie. Les systèmes
hamiltoniens forme une large classe de modèles conservatifs qui incluent les
sytèmes dépendant d’un potentiel et le pendule simple.
À travers les quelques exemples présentés, plusieurs concepts importants ont
été rapidement abordés : sensibilité aux conditions initiales, phénomène de
résonance, instabilité paramétrique, etc. Le calcul d’équations d’amplitude
par la méthode des échelles a été rapidement esquissé sur deux exemples.
Enfin, une initiation à l’exploration numérique des systèmes dynamiques a
été effectuée.
Cette introduction est une invitation à approfondir l’étude des systèmes dy-
namique dont la compréhension est essentielles pour comprendre la dynamique
de nombreux systèmes physiques, en particulier dans le cas des la mécanique
des fluides.