Cours GNI
Cours GNI
Cours GNI
P. Chartier
24 février 2016
2
Table des matières
3
4 TABLE DES MATIÈRES
5
1. CADRE GÉNÉRAL ET RAPPEL SUR LES ÉQUATIONS DIFFÉRENTIELLES ORDINAIRES7
où f est une fonction définie sur D et (t0 , y0 ) un point de R × D. On suppose en outre que f
est continue et localement lipschitzienne, de sorte que pour tout (t0 , y0 ) ∈ R × D, le système
(1.1) admet une solution maximale unique sur un intervalle ouvert J(t0 , y0) ⊂ R (théorème de
Cauchy-Lipschitz). Alors, l’application
qui associe la valeur en t de la solution de (1.1) est bien définie sur l’ouvert
On note en outre
Ω0 = {(t, y0 ) ∈ R × D; t ∈ J(0, y0)}
Ω0 → Rn
(t, y0 ) 7→ ϕt (y0 ) = y(t; 0, y0)
et qu’on a en outre, pour tout (t, t0 , y0 ) ∈ Ω, (t − t0 , y0) ∈ Ω0 et y(t; t0 , y0 ) = ϕt−t0 (y0 ), ce qui
justifie la définition de ϕt comme une application indépendante de t0 .
Remarque 1.2 On peut aussi définir le flot d’un système non-autonome, mais on rappelle que
tout système non-autonome peut se réécrire comme un système autonome par l’adjonction de la
variable t. L’hypothèse d’autonomie ne constitue donc pas une restriction.
Proposition 1.3 L’application (t, y) 7→ ϕt (y) de Ω0 dans Rn est continue. En particulier, pour
tout (t, y) ∈ Ω0 , il existe un voisinage V de y tel que l’application ϕt (·) soit définie et continue sur
V.
8
R → Diff(Rn )
t 7→ ϕt (·)
est un homomorphisme du groupe (R, +) dans le groupe (Diff(Rn ), ◦), où Diff(Rn ) est l’ensemble
des difféomorphismes de Rn dans lui-même.
où A(t, y0 ) = ∂f
∂y
(ϕt (y0 )) ne dépend pas seulement du temps t mais aussi d’un “paramètre” y0 . La
résolvante exhibe donc elle-même une dépendance en y0 , mais aucune en t0 (seul t − t0 compte),
et on la note pour cette raison S(t; y0 ). D’après la discussion précédente, on a alors
∂ϕt
S(t; y0 ) = (y0 ).
∂y0
Cette égalité est en fait justifiée rigoureusement dans le théorème suivant (donné sans démonstration,
voir le cours de L3) :
2. EXEMPLES DE SYSTÈMES HAMILTONIENS 9
Théorème 1.5 Soit f une fonction de D dans Rn continue et localement lipschitzienne, telle
que ∂f∂y
(y) existe et soit continue sur D. Alors, le flot de (1.1) est une application continûment
différentiable par rapport à y et sa dérivée Ψt (y0 ) = ∂ϕ t
∂y0
(y0 ) vérifie l’équation variationelle as-
sociée à (1.1) :
(
Ψ̇t (y0) = ∂f ∂y
(ϕt (y0 ))Ψt (y0 )
(1.4)
Ψ̇0 (y0) = In
Théorème 1.6 Si f est de classe C k sur D, alors (t, y) 7→ ϕt (y) est également de classe C k sur
Ω0 .
ẏ = J −1 ∇H(y)
Dans les applications à la physique, le vecteur y ∈ R2n est souvent partitionné en ses composantes
position q ∈ Rn et quantité de mouvement p ∈ Rn . On peut alors écrire y = (p, q) et le système
sous la forme :
q̇ = ∇p H(q, p),
ṗ = −∇q H(q, p).
Lorsque l’hamiltonien est séparable, i.e. qu’il peut s’écrire sous la forme
q̈ = −∇q U(q).
10
Le pendule simple est constitué d’une masse fixée à l’extrémité d’un fil sans masse, inexten-
sible et sans raideur. Cette masse supposée ponctuelle oscille sous l’effet de l’attraction terrestre.
Pour mettre en équation le mouvement, on repère la position du pendule simple par l’angle θ qu’il
fait avec la verticale orientée vers le bas. descendante. L’accélération est notée g et vaut approxi-
mativement g ≃ 9, 81 m.s−2 en unités S.I. Deux forces s’exercent sur la masse : le poids mg et la
tension T du fil, perpendiculaire au mouvement de la masse. Dans ce modèle simplifié, les frotte-
ments sont négligés. L’énergie du pendule (cinétique + potentielle) s’écrit alors en fonction de de
la longueur du fil l, de la vitesse angulaire θ̇ comme :
1 1
H(mlθ̇, lθ) = ml2 θ̇2 + mgl(1 − cos θ) = (mlθ̇)2 + mgl(1 − cos((lθ)/l)) (2.5)
2 2m
Ici, la variable position q coı̈ncide avec lθ et la variable quantité de mouvement p avec mlθ̇. Les
équations hamiltoniennes s’écrivent ainsi :
∂H
q̇ = = p/m
∂p
∂H
ṗ = − = −mg sin(q/l)
∂q
θ̈ + ω 2 sin θ = 0
Pendule double
où r = kx1 −x2 k2 . Le corps 1 exerce sur le corps 2 une force F1→2 portée par le vecteur u = x1 −xr
2
,
m1 m2
de direction la direction de u, et de norme G r2 . Le corps 2 exerce sur le corps 1 une force F2→1
portée par le vecteur −u, de direction la direction de −u, et de norme G m1rm 2
2
. On a ainsi la figure
suivante :
m1 F2→1 F1→2 m2
✲ ✛
m2
mẌ(t) = −G 3 X(t),
r
où r(t) = kX(t)k2 . Cette composante décrit un mouvement à force centrale, que nous allons
maintenant analyser.
12
Remarque 2.3 On peut également travailler avec l’énergie massique Hm (t) = 21 kẊ(t)k22 −
m
G kX(t)k = E(t)/m. La constante G est la constante gravitationelle. Dans le système S.I., elle
vaut G = 6, 6742867 × 10−11 m3 kg−1 s−2 .
u2 v3 − u3 v2
1. On rappelle que le produit vectoriel est défini par u∧v = u3 v1 − u1 v3 . Il est bilinéaire et antisymétrique.
u1 v2 − u2 v1
2. EXEMPLES DE SYSTÈMES HAMILTONIENS 13
Lemme 2.4 Soit X(t) ∈ R2 le vecteur coordonnées du centre de gravité du système, satisfaisant
m2
mẌ(t) = −G X(t),
r3
où r(t) = kX(t)k2 . Alors, il existe des constante E0 et L0 de R telles que pour tout t ≥ 0, on a
H(t) = E0 ,
L(t) = L0 .
et d’autre part
L(t) = m det(X(t), Ẋ(t))
r(t) cos(θ(t)) ṙ(t) cos(θ(t)) − r(t) sin(θ(t))θ̇(t)
= m
r(t) sin(θ(t)) ṙ(t) sin(θ(t)) + r(t) cos(θ(t))θ̇(t))
= mr 2 (t)θ̇(t)
L0 L0
autrement dit, θ̇(t) = mr 2 (t)
. En remplaçant θ̇(t) par mr 2 (t)
dans l’expression de H(t), il vient alors
2 !
1 L0 m2
m ṙ 2 (t) + r 2 (t) −G = H0
2 mr 2 (t) r
soit
L20 2m 2H0
ṙ 2 (t) + 2 2
−G =
m r (t) r m
On observe alors que θ̇(t) conserve un signe constant, celui de L0 , que l’on suppose non nul (si L0
est nul, alors Ẋ(t) et X(t) sont colinéaires pour tout t, et donc le mouvement de X est rectiligne).
On peut alors faire le changement de changement de variable t = t(θ) et exprimer l’équation en
fonction de θ. Il vient :
dr dr dθ dr L0 L0 d(1/r)
ṙ = = = 2
=−
dt dθ dt dθ mr m dθ
de sorte qu’en posant z = 1/r, on obtient finalement l’équation :
2
dz 2Gm3 2H0 m
(θ) + z 2 (θ) − 2
z(θ) =
dθ L0 L20
Gm3
Un dernier calcul avec α = L20
β = 2HL02m et u = z − α montre que
,
0
2
du
(θ) + u2 (θ) = α2 + β
dθ
et en dérivant par rapport à θ :
d2 u
(θ) + u(θ) = 0
dθ2
p la solution générale peut s’écrire : u(θ) = uθ0 cos(θ − θ0 ), avec la condition initiale uθ0 =
dont
± β + α2 . Ainsi, r s’écrit :
p
r(θ) =
1 ± e cos(θ − θ0 )
q
L2 2H L2
avec p = 1/α = Gm0 3 et e = 1 + m50G20 . Suivant la valeur de e, la trajectoire de X(t) est
donc une parabole, une hyperbole, ou une ellipse d’excentricité e (on peut vérifier en effet que
2H L2
1 + m50G20 > 0).
2. EXEMPLES DE SYSTÈMES HAMILTONIENS 15
Système gravitationnel
On considère donc un système à N particules de positions qi ∈ R3 et moment cinétiques
pi ∈ R3 , i = 1, . . . , N, soumises à l’attraction gravitationnelle. L’hamiltonien s’écrit alors
N
1X 1 T X
H(p, q) = pi pi + Gmi mj U(kqi − qj k)
2 i=1 mi 1≤i<j≤N
où mi désigne la masse de la particule i et U(r) = − 1r . Le système peut être aisément si-
mulé (voir http ://w3.bretagne.ens-cachan.fr/math/people/gilles.vilmart/sunjupitersaturn.html pour
le problème à trois corps Soleil-Jupiter-Saturne). Un tel système comporte trois invariants :
1. l’hamiltonien (énergie totale du système)
P
2. la quantité de mouvement totale P (p) = N i=1 pi
PN
3. le moment angulaire total L(p, q) = i=1 qi ∧ pi
La conservation de l’énergie est une propriété générale des systèmes hamiltoniens (voir ci-après).
On vérifie par ailleurs aisément la conservation des deux autres quantités en dérivant par rapport à
t. D’une part
XN XN X
d
P (p) = ṗk = −G ∇q k mi mj U(kqi − qj k)
dt k=1 k=1 1≤i<j≤N
X
′
= −G mi mk U (kqi − qk k)∇qk kqi − qk k
1≤i<k≤N
X
−G mi mk U ′ (kqj − qk k)∇qk kqj − qk k
1≤k<j≤N
X
= −G mi mk U ′ (kqi − qk k)∇qk kqi − qk k = 0
1≤i6=k≤N
car
1
∇qk kqi − qk k = (qk − qi ).
kqi − qk k
D’autre part
XN XN
d 1
L(p, q) = (q̇i ∧ pi + qi ∧ ṗi ) = pi ∧ pi + qi ∧ ṗi
dt i=1 i=1
mi
16
Dynamique moléculaire
N
1X 1 T X
H(p, q) = pi pi + U(kqi − qj k)
2 i=1 mi 1≤i<j≤N
avec
U(r) = 4E0 (r0 /r)12 − (r0 /r)6 .
Le potentiel de Lennard-Jones est utilisé pour décrire les interactions entre deux atomes au sein
d’un gaz monoatomique de type gaz rare. Le terme (r0 /r)6 est attractif et dominant à grande dis-
tance : c’est l’interaction de Van der Waals. L’exposant 12 du terme répulsif, dominant à courte
distance, est empirique.
Cette définition implique que la quantité I est conservée le long de toute solution de ẏ = f (y),
quelle que soit la condition initiale associée.
Exemple 3.2 L’énergie (l’hamiltonien) d’un système hamiltonien est une intégrale première (voir
paragraphe suivant).
Exemple 3.3 Quantité de mouvement et moment angulaire dans les systèmes à N particules.
3. PROPRIÉTÉS GÉOMÉTRIQUES DU FLOT 17
n
X
∂f ∂fi
∀y ∈ D, div (f )(y) = Tr (y) = (y) = 0
∂y i=1
∂yi
son volume. Le flot ϕt (·) considéré comme application de Rn dans Rn envoie chaque point y de A
sur un point ϕt (y) de ϕt (A) et il est naturel de considérer le volume de l’ensemble image ϕt (A), à
savoir :
Z
dy
ϕt (A)
On a alors :
Théorème 3.4 Pour un système différentiel de la forme ẏ = f (y), avec f de classe C 1 sur D telle
que divf ≡ 0, alors
∂ϕt
Preuve. Soit Ψt (y) = ∂y
(y). Ψt est solution de l’équation variationnelle
Ψ̇t (y) = ∂f ∂y
(ϕt (y))Ψt(y),
Ψ0 (y) = IRn
Pour une matrice M ∈ GLn (R), on a :
∀ H ∈ Mn (R), (dM det)H = (det M)T r(M −1 H)
En effet :
det(M + tH) − det(M) −1 n 1 −1
= t det(M) t det( IRn + M H) − 1
t t
= t det(M) (1 + t T r(M −1 H) + . . . + tn det(M −1 H) − 1)
−1
Il vient donc 2 :
d d d
det(Ψt (y)) = (dΨt det) Ψt = det(Ψt ) T r(Ψ−1 t Ψt ),
dt dt dt
= det(Ψt ) T r(Ψ−1
t (∂y f (ϕt (y))) Ψt ),
= det(Ψt ) T r(∂y f (ϕt (y))) = 0.
Ainsi, det(Ψt (y)) = det(Ψ0 (y)) = det(IRn ) = 1, et
Z Z ∂ϕ Z
t
dz = det (y) dy = dy.
ϕt (A) A ∂y A
ẏ = J −1 ∇y H(y),
où H(·) est une fonction de R2n dans R et où J est la matrice de M2n (R)
0 In
J=
−In 0
Notons que pour un système hamiltonien, on a
∂f
divf = Tr( ) = Tr(J −1 ∇2 H) = Tr(∇2 HJ −T ) = −Tr(J −1 ∇2 H) = −divf = 0
∂y
de sorte que le flot d’un système hamiltonien préserve le volume. On a en outre :
Théorème 3.5 Soit ϕt le flot associé à un système hamiltonien. Alors, pour tout (t, y) ∈ Ω0 ,
H(ϕt (y)) = H(y).
2. On peut aussi faire un calcul direct. Soient Ψ1 , . . . , Ψn les vecteurs colonnes de Ψt et αi,j , 1 ≤ i, j ≤ n, les
coefficients de la matrice Θ = Ψ−1
t ∂y f (ϕt (y))Ψt . On a bien sûr pour tout j = 1, . . . , n
X n
d
Ψj = (∂y f (ϕt (y)))Ψj = αi,j Ψi .
dt i=1
P = {tξ + sη | 0 ≤ t, s ≤ 1}
En dimension n > 1, on remplace cette expression par la somme ω(ξ, η) des aires orientées des
projections sur les plans (pi , qi ) de P :
n
X n
X
ξip ηip
ω(ξ, η) = = (ξip ηiq − ξiq ηip ).
ξiq ηiq
i=1 i=1
ω définit ainsi une forme bilinéaire antisymétrique, que l’on peut encore écrire :
ω(ξ, η) = ξ T Jη,
Transformations symplectiques
Définition 3.6 Une application linéaire A : R2n → R2n (confondue une fois encore avec sa
représentation matricielle de GL2n (R)) est dite symplectique si :
AT JA = J,
R2n−2
ξ
η
qI
pI
(ξip ηiq − ξiq ηip )
F IGURE 3 – Application ω
Définition 3.7 Une application g de U ouvert de R2n dans R2n , de classe C 1 sur U, est dite sym-
plectique si sa matrice jacobienne g ′ (y) est symplectique pour tout y de U, c’est-à-dire si :
T
∀y ∈ U, (g ′ (y)) Jg ′(y) = J,
ou de manière équivalente si :
Soit M une variété bidimensionnelle de U, telle qu’il existe une “carte globale” :
M = ψ(K),
où K est un compact de R2 et ψ(s, t) un difféomorphisme de K dans M. M peut être vue comme
la limite de l’union de “petits parallélogrammes” engendrés par les vecteurs
∂ψ ∂ψ
ds et dt.
∂s ∂t
Alors, la somme des aires orientées des projections sur les plans (pi , qi ) de tous ces parallélogrammes
s’écrit :
ZZ
∂ψ ∂ψ
Ω(M) = ω( (s, t), (s, t))dsdt.
K ∂s ∂t
3. PROPRIÉTÉS GÉOMÉTRIQUES DU FLOT 21
∂ψ ∂ψ
ds dt ∂g◦ψ
∂s ∂t g ∂s
ds
R2d−2
∂g◦ψ
∂t
dt
(p0 , q0 )
g(p0 , q0 )
qI
I
q
pI pI
Théorème 3.8 Soit U un ouvert de R2n et g une application de U dans R2n , de classe C 1 . Si g est
symplectique sur U, alors elle préserve Ω(M), c’est-à-dire :
Ω(g(M)) = Ω(M)
Théorème 3.9 (Poincaré 1899) Soit H(·) une fonction de classe C 2 d’un ouvert D de R2n dans R
(telle que ∇y H(·) soit localement lipschitzienne). Alors pour tout (t, y) ∈ Ω0 , ϕt est une transfor-
mation symplectique.
∂ϕt
Preuve. La matrice Ψt (y) = ∂y
est solution de l’équation variationnelle :
Ψ̇t (y) = J −1 ∇2 H(ϕt (y))Ψt (y)
Ψ0 (y) = I2n
22
T
Or ∇2 H(ϕt (y)) est symétrique, i.e. ((∇2 H(ϕt (y))) = ∇2 H(ϕt (y)). D’où :
d
ΨTt (y)JΨt (y) = Ψ̇Tt (y)JΨt (y) + ΨTt (y)J Ψ̇t (y),
dt
= ΨTt (y)(∇2 H)T |J −T T −1 2
| {zJ}(∇ H)Ψt (y)
{z J} Ψt (y) + Ψt (y) J
=−J −1 J=−I =I
= 0
En outre, pour t = 0 on a
ΨT0 (y)JΨ0(y) = I T JI = J
ce qui permet de conclure.
Théorème 3.10 Soit D un ouvert connexe de R2n et f une fonction C 1 de D dans R2n . On suppose
qu’il existe un t > 0 et un ouvert U simplement connexe 3 ou étoilé 4 tels que [0, t] × U ⊂ Ω0 et
que pour tout 0 ≤ s ≤ t et tout y ∈ U, ϕs (y) est symplectique. Alors, ẏ = f (y) est un système
hamiltonien sur U, c’est-à-dire qu’il existe une fonction H de classe C 2 sur U telle que
∀y ∈ U, f (y) = J −1 ∇y H(y).
Lemme 3.11 (Lemme d’intégrabilité) Soit U un ouvert de R2n et f : U → R2n une fonction de
classe C 1 , telle que f ′ (y) soit symétrique pour tout y ∈ U. Alors, pour tout y0 de D, il existe un
voisinage V(y0 ) et une fonction H définie sur V(y0 ) telle que :
Preuve. Soit B0 ⊂ U une boule de centre y0 contenue dans U. On définit H sur B0 par :
Z 1
H(y) = (y − y0 )T f (y0 + t(y − y0 ))dt.
0
Remarque 3.12 Lorsque l’ouvert considéré est étoilé par rapport à l’un de ses points y0 ∈ U,
alors on peut définir H sur tout U.
Remarque 3.13 Dans le cas général d’un ouvert quelconque, le résultat du théorème 3.10 est
faux. Considérons par exemple le système :
p
ṗ = p2 +q 2
q (3.7)
q̇ = p2 +q2
défini sur U = {(p, q) ∈ R2 ; (p, q) 6= (0, 0)}. Pour (p0 , q0 ) ∈ U, le flot s’écrit :
p0
ϕt (p0 , q0 ) = α(t, r0 ) ,
q0
p p
avec α(t, r0 ) = 1 + 2t/r02 et r0 = p20 + q02 . Sa dérivée ∂(p∂ϕ t
0 ,q0 )
est de la forme :
" #
p2
∂ϕt (∂p0 α)p0 + α (∂q0 α)p0 (∂r0 α) r00 + α (∂r0 α) pr00q0
= = q2 .
∂(p0 , q0 ) (∂p0 α)q0 (∂q0 α)q0 + α (∂r0 α) pr00q0 (∂r0 α) r00 + α
C’est une matrice symplectique si α2 + αr0 (∂r0 α) = 1, ce qui se vérifie par un simple calcul.
Localement, on peut écrire le système comme un système hamiltonien. Par exemple, dans le demi-
plan p > 0, on peut prendre H(p, q) = − arctan pq et vérifier que :
1 1 p
−∇q H(p, q) = q 2 = = ṗ (3.8)
p 1+ 2 p + q2
2
p
q 1 q
∇p H(p, q) = 2 2 = = q̇ (3.9)
p 1+ 2 q p + q2
2
p
24
δ(t)
δ0
A contrario, supposons qu’il existe un H de classe C 2 sur U, tel que le champ de vecteur f (p, q) =
(p2 + q 2 )−1 [p, q]T s’écrive f (p, q) = J −1 ∇p,q H(p, q). Considérons alors la forme différentielle
et calculons son intégrale le long du chemin Γ paramétré par (p, q) = (cos(θ), sin(θ)) :
Z Z 2π Z 2π
− sin(θ)
ω= ωcos(θ),sin(θ) dθ = (− sin2 (θ) − cos2 (θ))dθ = −2π
Γ 0 cos(θ) 0
Or on aurait par ailleurs ωp,q = (∂p H)(p, q)dp + (∂q H)(p, q)dq = dH dont l’intégrale sur Γ est
nulle. H ne peut donc pas être défini sur tout U. L’hypothèse de simple connexité est essentielle
pour cela.
L’une des propriétés fondamentales des transformations symplectiques est qu’elles préservent la
structure hamiltonienne. Plus encore, c’est une propriété caractéristique. On a en effet le théorème
suivant :
Théorème 3.14 Soit Ψ un difféomorphisme d’un ouvert U de R2n dans V = Ψ(U). Alors, si Ψ
est symplectique, le système hamiltonien
ẏ = J −1 ∇y H(y)
ż = J −1 ∇z K(z) (3.11)
Or on a, en omettant les arguments des fonctions, ∇y H = (∂y (K ◦ Ψ))T = ((∂z K)(∂y Ψ))T =
(Ψ′ )T (∇z K). D’où :
ż = Ψ′ (Ψ−1 (z))J −1 (Ψ′ (Ψ−1 (z)))T ∇z K(z). (3.12)
Si Ψ est symplectique, (Ψ′ )T J(Ψ′ ) = J, donc Ψ′ J −1 (Ψ′ )T = J −1 et
ż = J −1 ∇z K(z).
Réciproquement, si le système en z (3.12) est hamiltonien pour tout K, alors Ψ′ J −1 (Ψ′ )T = J −1
donc (Ψ′ )T J(Ψ′ ) = J et Ψ est symplectique.
Si f et g sont supposées de classe C ∞ , alors on peut itérer l’opérateur Lf et considérer ses puis-
sances Lkf en définissant
Lk+1 k
f [g] = Lf [Lf [g]], k = 1, . . . , ∞
où l’argument ϕt (y) de g ′′, f ′ et f a été omis dans la seconde ligne pour alléger les notations, et
plus généralement
dk
k
g(ϕt(y)) = (Lkf [g])(ϕt (y))
dt
de sorte que le développement en série de Taylor de g(ϕt (y)) en t = 0 s’écrit :
X tk
g(ϕt (y)) = (Lkf [g])(y) = exp(tLf )[g](y),
k≥0
k!
et pour g(y) = y :
Remarque 3.16 L’opérateur Lf fait apparaı̂tre la dérivée première de son argument g et est dit
d’ordre 1 pour cette raison. Le k-ième itéré de Lf fait apparaı̂tre la dérivée k-ième de son argument
g et est dit d’ordre k.
d’où il apparaı̂t clairement que les opérateurs Lf1 et Lf2 ne commutent pas en général. Les deux
opérateurs Lf1 Lf2 et Lf2 Lf1 sont d’ordre 2, mais de manière remarquable,
Les opérateurs apparaissent dans l’ordre inverse car (3.13) s’obtient en considérant g(y) = ϕ2t (y) et
en développant g(ϕ1s (y)). Les séries considérées sont des séries formelles, en général non conver-
gentes, car les opérateurs Lf1 et Lf2 sont non bornés. Par ailleurs, si Lf1 et Lf2 ne commutent pas,
alors
exp(sLf1 ) exp(tLf2 ) 6= exp(sLf1 + tLf2 ).
3. PROPRIÉTÉS GÉOMÉTRIQUES DU FLOT 27
s2 t st2
L(s, t) = sLf1 + tLf2 + st[Lf1 , Lf2 ] + [Lf1 , [Lf1 , Lf2 ]] + [Lf , [Lf2 , Lf1 ]] + . . . (3.14)
12 12 2
où [Lf1 , Lf2 ] = Lf1 Lf2 − Lf2 Lf1 . Cette identification devient rigoureuse si elle est interprétée
comme l’identification de développements de Taylor au voisinage de s = t = 0. Les termes
suivants du développement sont donnés par la formule Baker-Campbell-Hausdorff (BCH). La for-
mule (3.14) conduit à penser que si [Lf1 , Lf2 ] = 0, alors les flots ϕ1s et ϕ2t commutent. Le résultat
est vrai mais ne peut être démontré à partir de (3.14) puisque les séries considérées ne sont pas
convergentes.
Théorème 3.17 Soient f1 et f2 deux champs de vecteurs de C 1 (D, Rn ). La composition des flots
ϕ1s et ϕ2t est commutative si et seulement si
∂· ∂f1 ∂f2
[Lf1 , Lf2 ] = f2 − f1 = 0.
∂y ∂y ∂y
Preuve. Quitte à reparamétriser s en multipliant f1 par une constante (ce qui ne change rien à la
nullité du crochet de Lie), on peut considérer s = t. La loi de groupe ϕ1t ◦ ϕ1s = ϕ1t+s = ϕ1s ◦ ϕ1t
donne par dérivation par rapport à t
ϕ̇1t ◦ ϕ1s = ∂y (ϕ1s ) ◦ ϕ1t · ϕ̇1t
i.e.
f1 ◦ ϕ1t ◦ ϕ1s = ∂y (ϕ1s ) ◦ ϕ1t · (f1 ◦ ϕ1t )
et donc pour t = 0
f1 ◦ ϕ1s = ∂y (ϕ1s ) · f1 .
De même,
d2 2
ϕt ◦ ϕ1s = f2′ f1 .
dsdt s=t=0
Preuve. [Preuve alternative] Quitte à reparamétriser s en multipliant f1 par une constante (ce qui
ne change rien à la nullité du crochet de Lie), on peut considérer s = t. On remarque tout d’abord
que, d’après la formule (3.14), pour tout y de Rn , et pour h > 0 suffisament petit :
(ϕ2h ◦ ϕ1h − ϕ2h ◦ ϕ1h )(y) = O(h3 ) (3.15)
où la constante contenue dans O est obtenue par majoration des dérivées sur un compact contenant
(ϕ2u ◦ ϕ1v )(y) et (ϕ1v ◦ ϕ2u )(y) pour tout 0 ≤ u, v ≤ t. Prenons h = t/N. Il vient :
ϕ2t ◦ ϕ1t = ϕ2t−h ◦ ϕ2h ◦ ϕ1h ◦ ϕ1t−h
= ϕ2t−h ◦ ϕ1h ◦ ϕ2h ◦ ϕ1t−h + O(h3 )
= ϕ2t−h ◦ ϕ12h ◦ ϕ2h ◦ ϕ1t−2h + 2O(h3 )
= ···
= ϕ2t−h ◦ ϕ1N h ◦ ϕ2h + NO(h3 ).
En répétant l’opération N fois, on obtient :
(ϕ2t ◦ ϕ1t − ϕ1t ◦ ϕ2t ) = N 2 O(h3 ) = O(1/N) (3.16)
4. MÉTHODES NUMÉRIQUES GÉOMÉTRIQUES 29
La trajectoire exacte est une ellipse dans l’espace des phases (p, q) qui dépend des conditions
initiales. On teste maintenant le comportement qualitatif de trois schémas numériques :
1. Méthode d’Euler explicite :
pi+1 = pi + h(−ω 2 qi ) = pi + (−hω 2 )qi
qi+1 = qi + h(pi ) = hpi + qi
Tous ces schémas prennent la forme (particulièrement simple ici en raison du caractère linéaire du
problème) suivante :
pi+1 pi
= M(hω)
ωqi+1 ωqi
0
v
−1
−2
−3
−3 −2 −1 0 1 2 3
u
1 1 2 1
H(p1 , p2 , q 1 , q 2 ) = [(p ) + (p2 )2 ] − p ,
2 (q 1 )2 + (q 2 )2
= T (p) + V (q).
Le système différentiel correspondant est donc de la forme :
ṗ = − ∂H = −V ′ (q)
∂q
⇐⇒ q̈ = −V ′ (q)
q̇ = ∂H∂p
= p
La trajectoire exacte est une ellipse pour H < −0.5 dans l’espace des phases qui dépend des
conditions initiales. La figure 7 démontre clairement que la méthode du point milieu (qui est sym-
plectique comme on le verra) se comporte de bien meilleure façon que les deux autres.
Elle est donc déterminée par la donnée des coefficients A = (aij ) ∈ Rs×s et b = (bj ) ∈ Rs . Pour
un système différentiel de la forme
ṗ = f (p, q),
(4.17)
q̇ = g(p, q),
4. MÉTHODES NUMÉRIQUES GÉOMÉTRIQUES 31
1.5 −1
2
1 −2
1
0.5 −3
q2
q2
H
0 0 −4
−0.5 −5
−1
−1 −6
−2
−1.5 −7
−2 −3 −8
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 0 50 100 150 200 250 300 350 400 450
q1 q1 Time
F IGURE 7 – Solutions du problème de Kepler (pas h = 0.01 pour Euler, 0.1 pour le point milieu).
Définition 4.1 Soit une méthode de Runge-Kutta partitionnée de coefficients A, Â, b, b̂. Pour |h|
suffisamment petit, les équations précédentes définissent une fonction “flot numérique” Φh , telle
que yi+1 = Φh (yi ) définit la solution approchée. Si cette fonction est symplectique dès lors que le
système est hamiltonien, alors la méthode est dite symplectique.
Remarque 4.2 La méthode est explicite lorsque A et  sont triangulaires inférieures strictes.
Q(p, q) = pT Cq,
32
où B = diag(b) et B̂ = diag(b̂), sont satisfaites. Si en outre dans (4.17), f ne dépend que de q et
g ne dépend que de p, alors la première condition est suffisante.
Or
X
p0 = P j − h ajk f (Pk , Qk ), j = 1, . . . , s,
k
X
q0 = Qj − h âjk g(Pk , Qk ), j = 1, . . . , s,
k
donc
X
Q(p1 , q1 ) = pT0 Cq0 + h bj f T (Pj , Qj ) C Qj + b̂j PjT C g(Pj , Qj )
j | {z }
terme nul car f T (p,q)Cq+pT Cg(p,q)=0
X
−h2 bj âj k f T (Pj , Qj ) C g(Pk , Qk ) + b̂j aj k f T (Pk , Qk ) C g(Pj , Qj )
j, k
X
+h2 bj b̂k f T (Pj , Qj ) C g(Pk , Qk )
j, k
X
= pT0 Cq0 −h 2
bj âj k + b̂k ak j − bj b̂k f T (Pj , Qj ) C g(Pk , Qk )
j, k | {z }
=0 par hypothèse
D’après le théorème 3.4, si Tr(M) = 0, alors det(Y (t)) = det(Y (0)) = 1 pour tout t ∈ R. Dès
que n ≥ 3, det Y est donc un invariant polynomial de degré au moins 3. Nous allons montrer
qu’aucune méthode de Runge-Kutta ne peut conserver cet invariant génériquement, c’est-à-dire
pour toute matrice M.
Lemme 4.4 Soit z 7→ R(z) une fonction analytique au voisinage de z = 0, satisfaisant R(0) = 1
et R′ (0) = 1. Soit n ≥ 3. Alors les deux assertions suivantes sont équivalentes :
1. ∀M ∈ M(Rn ), Tr(M) = 0 ⇒ det(R(M)) = 1
2. ∀z ∈ C, R(z) = exp(z)
Preuve. Il est clair que 2. implique 1., car la solution de l’équation différentielle en t = 1 est
donnée par
Y (1) = etM t=1 = R(M)
et d’après le théorème 3.4, 1 = det(Y (0)) = det(Y (1)) = det(R(M)). Montrons maintenant la
réciproque : pour M = Diag(ν, µ, −ν − µ, 0, . . . , 0) on a
R(M) = Diag(R(ν), R(µ), R(−ν − µ), 1, . . . , 1),
de sorte que 1 = det(R(M)) = R(ν)R(µ)R(−ν − µ). Pour ν = 0, on obtient R(−µ) = 1/R(µ),
puis
R(ν)R(µ) = R(ν + µ).
La fonction z 7→ R(z) est donc un morphisme de groupe de (C, +) dans (C, ×) satisfaisant R(0) =
R′ (0) = 1, donc R(z) = exp(z).
Théorème 4.5 Une méthode de Runge-Kutta consistante (d’ordre au moins 1) ne peut pas conser-
ver les invariants polynomiaux de degrés d ≥ 3.
Preuve. Soit une méthode de Runge-Kutta de coefficients (A, b). Appliquée au système matriciel
Ẏ = MY , la méthode s’écrit simplement :
Xs
Yk,i = Yk + hM aij Yk,j , i = 1, . . . , s
j=1
Xs
Yk+1 = Yk + hM bj Yk,j ,
j=1
de sorte que
Yk+1 = R(hM)Yk ,
34
avec R(z) = 1 + zbT (In − zA)−1 1. On a R(0) = 1 et la méthode étant consistante, R′ (0) = bT 1 =
1. Si la méthode de Runge-Kutta conserve les invariants de degrés d ≥ 3, alors elle satisfait la
proposition 1. du lemme précédent, et donc R(z) = exp(z), ce qui entre en contradiction avec la
forme rationnelle de R(z) obtenue ci-dessus.
avec B = diag(b), B̂ = diag(b̂). Si l’hamiltonien est séparable, i.e. H(p, q) = T (p) + V (q), la
première condition est nécessaire et suffisante.
Preuve. Pour démontrer que la condition est nécessaire, il suffit d’interpréter la condition de
symplecticité comme un invariant quadratique. Pour un système hamiltonien
ẏ(t) = J −1 ∇H(y(t))
(4.25)
y(0) = y0
on a en effet
d T
(Ψ JΨt ) = 0 (4.26)
dt t
∂ϕt (y0 )
où Ψt = ∂y0
. Maintenant, si Φh est le flot numérique d’une méthode de Runge-Kutta,
(Φh (y0 ), ∂y0 Φh (y0 ))
est, d’après le lemme 4.6, la solution numérique à l’instant h du système
ẏ = J −1 ∇H(y), y(0) = y0
−1 2
Ψ̇t = J (∇ H(y))Ψt, Ψ(0) = I
pour lequel ΨTt JΨt est un invariant quadratique. Si la méthode satisfait les conditions de sym-
plecticité, cet invariant est donc conservé. Pour montrer que la condition énoncée est en outre
nécessaire, il faut introduire une notion de méthode de Runge-Kutta “réduite”, traduisant le fait
que certaines étapes internes (autrement dit les Yi ) peuvent être redondantes.
Exemples
Il existe plusieurs classes de méthodes symplectiques. Citons en particulier les méthodes de
Gauss et les méthodes de Lobatto IIIA-IIIB (voir []), qui sont des méthodes partitionnées. La
méthode de Gauss la plus simple est la méthode du point milieu, qui est d’ordre 2 :
h
y1 = y0 + hf (Y ), Y = y0 + f (Y ).
2
Il est trivial de vérifier qu’elle est symplectique car ici A = 1/2 et b = 1. La méthode de Lobatto
IIIA-IIIB la plus simple est la méthode d’Euler symplectique que l’on peut écrire :
p1 = p0 + hf (p1 , q1 ), q1 = q0 + hg(p0, q0 ),
36
qui devient explicite dès lors que f de dépend que de q et g que de p. On a ici A = 1, = 0 et
b = b̂ = 1.
Théorème 4.9 Sout ϕt le flot exact de l’équation différentielle. On considère une méthode numérique
Φh d’ordre p telle que pour tout y dans Rn ,
Φh (y) = ϕh (y) + hp+1 C(y) + O(hp+2 ).
Alors la méthode adjointe Φ∗h , définie par Φ∗h = Φ−1
−h est du même ordre p et on a
Preuve. On a clairement
Φ−h (y) = ϕ−h (y) + (−h)p+1 C(y) + O(hp+2 ),
donc
(ϕh ◦ Φ−h )(y) = y + (−h)p+1 C(y) + O(hp+2 ),
puis
(ϕh ◦ Φ−h )−1 (y) = y − (−h)p+1 C(y) + O(hp+2 ),
ce qui implique enfin que
(Φ∗h ◦ ϕ−1
h )(y) = y − (−h)
p+1
C(y) + O(hp+2 ).
En prenant y = ϕh (z), il vient finalement
Φ∗h (z) = ϕh (z) + (−1)p hp+1 C(z) + O(hp+2 ).
Si la méthode Φh est symétrique, alors Φ∗h = Φh et on a nécessairement (−1)p = 1, i.e. p pair.
Notons que ces équations n’ont pas de solutions réelles pour les valeurs impaires de p. En revanche,
si p est pair et s = 3 par exemple, on obtient la solution
1 21/(p+1)
γ1 = γ3 = et γ 2 = − .
2 − 21/(p+1) 2 − 21/(p+1)
Il est donc possible à partir d’une méthode de base symétrique de construire une méthode d’ordre
au moins p + 1, donc p + 2 en raison de la symétrie de la méthode de composition obtenue. En
répétant la procédure, il est possible de construire des méthodes d’ordre arbitrairement élevé.
5 Analyse rétrograde
5.1 Idée générale
L’idée originale de l’analyse rétrograde est d’interpréter la solution numérique comme la solu-
tion exacte d’un problème approché.
A NALYSE R ÉTROGRADE
M
ẏ = f (y) numéthode
ériqu
e
q y0 , y1 , y2 , y3 , . . .
=
✶ ỹ(0), ỹ(h), ỹ(2h), . . .
ỹ˙ = f˜h (ỹ) t i on
Solxuacte
e
38
Le champ de vecteurs f˜h (ỹ) est déterminé sous la forme d’un développement en puissances
de h en général non convergent ; l’équation différentielle modifiée (satisfaite par l’approximation
numérique y0 , y1 , y2 , ...) est de la forme :
ỹ˙ = f (ỹ) + hf2 (ỹ) + h2 f3 (ỹ) + · · · , ỹ(0) = y0 .
Conséquence : L’erreur numérique sur l’hamiltonien est en Chp si p est l’ordre de la méthode
et ce pour des intervalles de temps exponentiellement longs.
5. ANALYSE RÉTROGRADE 39
2
2 2
1 1
−2
0 0
−4
v
v
−6
−1 −1
−8
−2 −2
−10
−3 −3
−12
−14 −4 −4
−20 −15 −10 −5 0 5 −7 −6 −5 −4 −3 −2 −1 0 1 2 3 −8 −6 −4 −2 0 2 4
u u u
F IGURE 8 – Solutions exactes des équations modifiées (en rouge) et solutions numériques (en
vert).
Si le champ de vecteur f est supposé (réel-)analytique dans une boule complexe B2R (y0 ) et est
majoré par M sur B2R (y0 ), alors pour h suffisamment petit, le flot numérique Φh (y) associé à la
méthode de Runge-Kutta est analytique et les coefficients dj (y) de son développement
R (n) (0)
Preuve. Soit h0 = γ 2κM pour γ < 1. On définit les fonctions Θi (h, y) par Θi (h, ·) ≡ idCm et
pour n ∈ N par :
(n+1)
Θi : {|h| ≤ h0 } × BR (y0 ) −→ Cm
X (n)
(h, y) 7−→ y + h ai,j f (Θj (h, y)),
j
m étant la dimension du système d’équations différentielles. Si pour tout j, tout |h| ≤ h0 et tout
(n) (n+1)
y ∈ BR (y0 ), Θj (h, y) ∈ B3R/2 (y0 ), alors il en est de même pour Θi . En effet :
(n+1)
X (n) R
max kΘi (h, y) − yk ≤ max kh ai j f (Θj (h, y))k ≤ h0 κ M ≤ .
|h|≤h0 , y∈BR (y0 ) |h|≤h0 , y∈BR (y0 )
j
2
40
(n)
Les fonctions Θi (h, ·) prennent donc leurs valeurs dans B3R/2 (y0 ). En dérivant par rapport à z la
fonction α(z) = f (y + z∆y) pour y ∈ B3R/2 (y0 ), k∆yk ≤ 1 et z ∈ BR/2 (0), une majoration de
Cauchy donne
M 2M
kα′(0)k = kf ′ (y)∆yk ≤ = . (5.30)
R/2 R
L’opérateur f ′ est donc borné par (2M)/R sur B3R/2 et il vient alors pour tous m ≥ n :
(m) (n)
X (m−1) (n−1)
max kΘi (h, y) − Θi (h, y)k ≤ max kh ai,j (f (Θj (h, y)) − f (Θj (h, y))k
i, |h|≤h0 , y∈BR (y0 ) i, |h|≤h0 , y∈BR (y0 )
j
2M (m−1) (n−1)
≤ h0 κ max kΘ (h, y) − Θj (h, y)k
R j, |h|≤h0, y∈BR (y0 ) j
...
(m−n) (0)
≤ γn max kΘi (h, y) − Θi (h, y)k
i, |h|≤h0 , y∈BR (y0 )
R
≤ γn
2
(n)
La suite des Θi converge donc uniformément vers une fonction Θ∞
i analytique (en h et y),
(n)
puisque chacune des Θi est analytique. Finalement,
kΦh (h, y) − yk ≤ µM|h| (5.31)
et
j−1
1 dj µ M h0 2κM
kdj (y)k = (Φh (y) − y) ≤ ≤ µM . (5.32)
j! dhj h=0 hj0 γR
Le résultat s’obtient alors en faisant tendre γ vers 1.
où les Li ≡ Lfi sont les opérateurs dérivée de Lie associés aux champs de vecteurs fi .
Or Lif˜h 2
= L1 + hL2 + h L3 + . . ., et il vient donc :
!i
X ti X
ϕ̃ht (y) = hk−1 Lk [id](y)
i≥0
i! k≥1
X ti X
= y+ i i!
hk1 +...+ki Lk1 · · · Lki [id](y)
i≥1
h k ≥1,...,k ≥1
1 i
j
!
X X ti X
j
= y+ h Lk1 · · · Lki [id](y) .
j≥1 i=1
hi i! k
1 +...+ki =j,k1 ≥1,...,ki ≥1
En identifiant les puissances de h, on a maintenant que ϕ̃hh (y) = Φh (y) si on a pour tout j ≥ 1 :
j
X 1 X
dj (y) = Lk1 · · · Lki [id](y). (5.39)
i=1
i! k
1 +...+ki =j,k1 ≥1,...,ki ≥1
Le terme de la somme correspondant à i = 1 n’est autre que Lj [id](y) = fj (y) et le résultat est
alors immédiat.
42
Théorème 5.5 On suppose que le champ de vecteur f est analytique sur B2R (y0) et borné par
M sur cette boule. On suppose par ailleurs, que la majoration (5.29) des coefficients dj du
développement en série de Taylor de la méthode numérique est valable sur la boule BR (y0 ). Alors
les coefficients fj de l’équation modifiée tels que calculés dans le lemme 5.3 sont majorés par :
j−1
ηMj
∀ y ∈ BR/2 (y0 ), kfj (y)k ≤ (log 2)ηM , (5.42)
R
où η = 2 max(κ, µ/(2 log 2 − 1)).
R
Preuve : Soit J ≥ 2 un entier fixé et δ = 2(J−1) . On définit, pour j = 1, . . . , J, les rayons
j−1 R J −j
rj = + R (5.43)
J −1 2 J −1
5. ANALYSE RÉTROGRADE 43
i
1 Y
kLk1 · · · Lki [id]krj ≤ kfkl krkl . (5.44)
δ i−1
l=1
En insérant cette majoration dans l’expression de fj calculée dans le lemme 5.3, on obtient alors,
pour 1 ≤ j ≤ J :
j i
X 1 X Y
kfj krj ≤ kdj kR + kfkl krkl . (5.45)
i=2
δ i−1 i! k
1 +...+ki =j l=1
on peut montrer par une récurrence facile, que 5 kfj krj ≤ δβj . Afin d’estimer les βj , on définit
maintenant la fonction génératrice
X
b(ζ) = βj ζ j . (5.47)
j≥1
De (5.46), on tire :
j−1 Xj
X Yi
µM
j 2κM 1
βj ζ = ζ ζ + βkl ζ kl ,
δ R i=2
i! k +...+k =j l=1
1 i
puis
X j−1 j i
µM X
j 2κM XX 1 X Y
βj ζ = ζ ζ + βkl ζ kl ,
j≥1
δ j≥1
R j≥2 i=2
i!
k1 +...+ki =j l=1
µM
5. Notons au passage que β1 = δ et que kf1 kr1 = kf kR ≤ kf k2R ≤ M .
44
K(Blog 2 )
B2 log 2−1
L
c’est-à-dire :
X1X X Yi
µM 1
b(ζ) = ζ + βkl ζ kl ,
δ 1 − 2κM
R
ζ i≥2
i! j≥i k +...+k =j l=1 1 i
X1 X i
Y
µM 1
= ζ 2κM
+ βkl ζ kl ,
δ 1− R ζ i≥2
i! k ,...,k l=1
1 i
µM 1 X 1
= ζ 2κM
+ b(ζ)i .
δ 1− R ζ i≥2
i!
où ϕt (y) désigne le flot exact de l’équation différentielle ẏ = f (y). Alors l’équation modifiée est
de la forme :
Preuve. C’est une conséquence immédiate de la définition des coefficients de l’équation modifiée.
On a en effet :
ϕ̃hh (y) = Φh (y) = ϕh (y) + hp+1 δp+1 (y) + O(hp+2 )
et donc :
j
!
X X 1 X X hj
hj
Lk1 · · · Lki [id](y) = Ljf [id](y) + hp+1 δp+1 (y) + O(hp+2 )
j≥1 i=1
i! k j≥1
j!
1 +...+ki =j,k1 ≥1,...,ki ≥1
46
Théorème 5.7 Soit Φh le flot numérique associé à une méthode de Runge-Kutta symplectique
appliquée à un système hamiltonien d’hamiltonien H de classe C ∞ sur R2n . Alors l’équation
modifiée est elle-aussi hamiltonienne, c’est-à-dire que pour tout j ∈ N∗ , il existe un hamiltonien
Hj de classe C ∞ sur R2n tel que
fj (y) = J −1 ∇Hj (y). (5.54)
Preuve. La preuve se fait par récurrence. Pour j = 1, f1 ≡ f et le résultat est vrai. Supposons
maintenant que pour r ≥ 1, la relation (5.54) soit établie pour tout 1 ≤ j ≤ r. Le flot ϕ̃r,h (y)
correspondant à l’équation modifiée
ỹ˙ = f (ỹ) + hf1 (y) + . . . + hr−1 fr (y) (5.55)
est le flot d’un système hamiltonien, donc est symplectique. Or :
Φh (y) = ϕ̃r,h (y) + hr+1 fr+1 (y) + O(hr+2 ) (5.56)
donc
Φ′h (y) = (ϕ̃r,h )′ (y) + hr+1 fr+1
′
(y) + O(hr+2 ). (5.57)
Les flots Φh (y) et ϕ̃r,h (y) étant symplectiques, on a alors :
T T
J = (Φ′h (y)) JΦ′h (y) = (ϕ̃r,h )′ (y)) J ϕ̃r,h )′ (y)
h T i
r+1 ′ T ′ ′ ′
+h (ϕ̃r,h ) (y)) Jfr+1 (y) + fr+1 (y) J(ϕ̃r,h ) (y)
+O(hr+2 )
(ϕ̃r,h )′ (y)=I+O(h)
h T i
r+1 ′ ′
= J +h Jfr+1 (y) + fr+1 (y) J + O(hr+2 ).
′
En conséquence, la matrice Jfr+1 (y) est symétrique et d’après le lemme d’intégrabilité, on a :
′
Jfr+1 (y) = ∇Hr+1 (y), (5.58)
où H est défini sur R2n et de classe C ∞ .
Théorème 5.8 Soit Φh le flot numérique associé à une méthode de Runge-Kutta préservant le
volume appliquée à un système dont le champ de vecteur est à divergence nulle. Alors l’équation
modifiée est elle-aussi à divergence nulle, c’est-à-dire que pour tout j ∈ N∗ , on a
divfj (y) = 0. (5.59)
5. ANALYSE RÉTROGRADE 47
Preuve. A nouveau, la preuve se fait par récurrence. Pour j = 1, f1 ≡ f et le résultat est vrai.
Supposons maintenant que pour r ≥ 1, la relation (5.59) soit établie pour tout 1 ≤ j ≤ r. Le flot
ϕ̃r,h (y) correspondant à l’équation modifiée
ỹ˙ = f (ỹ) + hf1 (y) + . . . + hr−1 fr (y) (5.61)
est le flot d’un système à divergence nulle, donc préserve le volume. De l”égalité
Φ′h (y) = ϕ̃′r,h (y) + hr+1 fr+1
′
(y) + O(hr+2 ),
on tire alors
det(Φ′h (y)) = det(ϕ̃′r,h (y) + hr+1 fr+1
′
(y) + O(hr+2 )).
Or det(M + H) = det(M) + (dM det)(H) + O(kHk2), soit ici avec M = ϕ̃′r,h (y) et H =
′
hr+1 fr+1 (y) + O(hr+2 ) :
′
1 = det(Φ′h (y)) = det(ϕ̃′r,h (y)) + hr+1 (dϕ̃′r,h (y) det)fr+1 (y) + O(hr+2 ),
= 1 + hr+1 Tr (ϕ̃′r,h (y))−1fr+1 ′
(y) + O(hr+2 ),
′
= 1 + hr+1 Tr fr+1 (y) + O(hr+2 ),
car ϕ̃′r,h (y) = I + O(h). Donc Tr fr+1
′
(y) = 0. En procédant de manière similaire avec I
I(y) = I (Φh (y)) = I (ϕ̃r,h (y)) + hr+1 (∇y I)T (ϕ̃r,h (y)) fr+1 (y) + O(hr+2 ),
on prouve la seconde partie du théorème.
et soit ϕ̃N,h (y) le flot exact correspondant. On souhaite ici déterminer s’il existe un N minimisant
l’erreur entre le flot numérique Φh et le flot ϕ̃h,N (y). Notons tout d’abord que si Φh admet un
développement de la forme (5.28) et que les termes du développement sont majorés comme en
48
R
(5.29), alors le dit développement est convergent pour |h| ≤ 2κM . Symétriquement, le flot ϕ̃N,h (y)
se développe en puissance de h et de l’opérateur Lf˜h = hL1 + . . . hN −1 LN :
N
!i
X hi X X
ϕ̃h,N = id + hk−1 Lk [id] = id + hj Gj [id]
i≥1
i! k=1 j≥1
Afin de majorer les Gj [id] sur la boule BR/4 (y0 ), on procède comme dans le théorème 5.5, avec
δ = R/(4(i − 1)). Pour i ≥ 2 indices kl , 1 ≤ kl ≤ N de somme j, on a :
i
1 Y
kLk1 · · · Lki [id]kR/4 ≤ kfkl kR/2
δ i−1 l=1
i−1 j−1
4(i − 1) log 2 ηMN
≤ (ηM log 2) .
N R
Il vient alors :
j−1 X
j i−1
ηMN 1 4(i − 1) log 2
kGj [id]kR/4 ≤ (ηM log 2) Card S(i, j, N) (5.64)
R i=1
i! N
où S(i, j, N) = {(k1 , . . . , ki ) ∈ [1, N]i , k1 + . . . + ki = j}, de sorte que Card S(i, j, N) ≤ N i−1
et donc :
j−1 X j j−1
ηMN (4(j − 1) log 2)i−1 16ηMN
kGj [id])kR/4 ≤ (ηM log 2) ≤ (ηM log 2) (5.65)
R i=1
i! R
R
En particulier, le développement de ϕ̃N,h est convergent pour y ∈ BR/4 (y0 ) et |hN| ≤ 16ηM
.
Preuve. Les champs de vecteurs fj , 1 ≤ j ≤ N sont construits de manière à ce que les termes des
développements de Φh et ϕ̃N,h coincident jusqu’à l’ordre N. On a donc :
X
Φh − ϕ̃N,h = hj (dj − Gj [id]) (5.67)
j≥N +1
Ainsi :
X
kΦh − ϕ̃N,h kR/4 ≤ hj kdj − Gj [id]kR/4
j≥N +1
X
j−1 j−1!
2κM 16ηMN
≤ hj µM + (ηM log 2)
j≥N +1
R R
X 16ηMNh j
≤ (2ηM log 2) h
j≥N
R
Pour Nh ≤ h0 , on a alors :
N
hN e (2eηM log 2) −N
kΦh − ϕ̃N,h kR/4 ≤ (2ηM log 2) h ≤ he (5.68)
eh0 e−1 e−1
car la fonction x 7→ (ǫx)x admet pour ǫ < 1 un minimum en x = (ǫe)−1 , et pour 0 < x ≤ (ǫe)−1 ,
(ǫx)x ≤ e−x (ici, ǫ = hh0 e ). Pour N = E(h0 /h) on obtient finalement :
(2e2 ηM log 2) − h0
kΦh − ϕ̃N,h kR/4 ≤ he h .
e−1
h0
sur des intervalles de temps nh ≤ e h .
Preuve : On peut majorer uniformément f = J −1 ∇H, par M, sur le compact KR . Toutes les
majorations obtenues successivement dans les théorèmes 5.29, 5.5 et 5.9 sont donc applicables
uniformément sur K avec un h0 et un N = E(h0 /h) uniformes eux-aussi. Maintenant,
n
X Xn
H̃(yn ) − H̃(y0 ) = H̃(yj ) − H̃(yj−1 ) = H̃(yj ) − H̃(ϕ̃N,h (yj−1 )) (5.72)
j=1 j=1
Xn
= H̃(Φh (yj−1 )) − H̃(ϕ̃N,h (yj−1)) (5.73)
j=1
On rappelle que l’exponentielle exp(A) d’une matrice A ∈ Mn (K) est définie par la série conver-
gente
X∞
1 k
exp(A) = A .
k!
k=0
Lorsque A et B sont deux matrices de Mn (K) qui commutent, alors exp(A+B) = exp(A) exp(B) =
exp(B) exp(A). L’objectif de cette section est d’établir une formule, dite formule de Baker-Campbell-
Hausdorff (BCH en abrégé), dans le cas où A et B ne commutent pas, puis d’approcher exp(t(A +
B)) par des produits de matrices de la forme exp(taA) et exp(tbB).
Pour Ω ∈ Mn (K) fixée, on définit l’opérateur adΩ comme l’application linéaire suivante :
et on note adiΩ l’opérateur composé adΩ ◦ . . . ◦ adΩ . La norme considérée sur l’espace vectoriel
| {z }
i fois
L(Mn (K)) est la norme subordonnée, qu’on notera encore k.k : en particulier, on a
kadΩ (A)k
kadΩ k = sup
A∈Mn (K),kAk6=0 kAk
dΩk
Exemple 6.1 Pour H ∈ Mn (K), on peut calculer dΩ
H pour k = 1, 2, 3. Il est clair que
dΩ
dΩ
H = H. Pour k = 2, 3, le calcul donne :
dΩ2
H = HΩ + ΩH,
dΩ
3
dΩ
H = HΩ2 + ΩHΩ + Ω2 H.
dΩ
En effet, on a HΩ + ΩH = 2HΩ + (ΩH − HΩ) = 2HΩ + adΩ (H) puis on vérifie que
Preuve. On a :
(Ω + tH)k+1 − Ωk+1 (Ω + tH)(Ω + tH)k − ΩΩk
=
t t
(Ω + tH)k − Ωk
= Ω + H(Ω + tH)k
t
En passant à la limite on obtient le résultat.
Preuve. La preuve se fait par récurrence. Pour k = 1, 2, 3, le résultat est établi (voir calculs
précédents). Supposons la formule établie jusqu’à un rang k ≥ 3. Alors
k+1 dΩ dΩk
dΩ k
H = H Ω +Ω H
dΩ dΩ dΩ
Xk−1
k
= Ω adiΩ (H)Ωk−1−i + HΩk
i+1
i=0
Lemme 6.4 Soit d expΩ l’application linéaire de Mn (K) dans Mn (K) définie par
X 1
∀H ∈ Mn (K), d expΩ H = adkΩ (H)
k≥0
(k + 1)!
Preuve. Notons tout d’abord que l’opérateur adiΩ est borné par 2i kΩki . En effet,
∀H ∈ Mn (K), kadΩ (H)k ≤ kΩkkHk + kHkkΩk = 2kΩkkHk
kadΩ (H)k
Donc, supkHk6=0 kHk
≤ 2kΩk et
kadiΩ (H)k i−1
kadΩ (H)k
sup ≤ 2kΩk sup ≤ . . . ≤ 2i kΩki .
kHk6=0 kHk kHk6=0 kHk
Maintenant, pour Ω non nulle et N quelconque, il vient
N
X N
X
1 k 1 e2kΩk − 1
k adΩ (H)k ≤ 2k kΩkk kHk ≤ kHk ≤ kHke2kΩk
(k + 1)! (k + 1)! 2kΩk
k=0 k=0
La borne reste valable pour Ω = 0, donc la série converge et l’opérateur d expΩ est borné par e2kΩk .
Preuve. Notons tout d’abord que toutes les séries considérées sont absolument
0 convergentes. On
P−1
peut donc écrire (par convention i=0 = 0 ce qui correspond au fait que dΩ dΩ
H = 0) :
X 1 dΩk
d exp(Ω)
H = H
dΩ k! dΩ
k≥0
X 1 X k−1
k
= adiΩ (H)Ωk−1−i
k! i+1
k≥0 i=0
X X 1
= adiΩ (H)Ωk−1−i
i≥0 k≥i+1
(i + 1)!(k − i − 1)!
XX 1
= adiΩ (H)Ωj
i≥0 j≥0
(i + 1)!j!
Cette dernière expression n’est autre que le produit de Cauchy des deux séries d expΩ (H) et
exp(Ω).
inverse. Notons tout d’abord que 0 est valeur propre de adΩ . Les vecteurs propres associés sont les
matrices M qui commutent avec Ω. Par exemple, pour n = 2 et
0 1
Ω= ,
−1 0
Preuve. La valeur propre λ = 0 de adΩ est “envoyée” sur la valeur P propre 1 de d expΩ . En
1
effet, pour λ = 0 et H non nulle commutant avec Ω, d expΩ H = k≥0 (k+1)! λk H = H. Soit
Q
P (λ) = ri=1 (λ − λi )mi le polynôme caractéristique de adΩ . Notons alors
r
Y
Q(λ) = (λ − R(λi ))mi
i=1
λ
P 1 k
où R(λ) = (e − 1)/λ = Si on tronque R, alors Q ◦ R est un polynôme divisible
k≥0 (k+1)! λ .
par P car
Yr Yr Yr
Q(R(λ)) = (R(λ) − R(λi ))mi = (λ − λi )mi (Ri (λ))mi = P (λ) (Ri (λ))mi
i=1 i=1 i=1
où Ri (λ) est un polynôme tel que R(λ) − R(λi ) = (λ − λi )Ri (λ). Cette dernière égalité est encore
λ
vraie pour R non tronquée donc si λi est valeur propre de multiplicté mi de adΩ , alors e λi −1 i
est
valeur propre de d expΩ avec la même multiplicité mi . L’espace Mn (C) étant de dimension finie,
on obtient bien toutes les valeurs propres de d expΩ comptées avec leur multiplicité. L’opérateur
λ
d expΩ est donc inversible si pour tout λi de Sp(adΩ ), e λi −1
i
6= 0. Une condition nécessaire et
suffisante pour que d expΩ soit inversible est donc :
Sp(adΩ ) ∩ 2iπZ∗ = ∅
6. FORMULE DE BAKER-CAMPBELL-HAUSDORFF (BCH) 55
La formule de l’opérateur inverse fait intervenir les nombres de Bernoulli définis par la donnée de
B0 = 1 et la relation récurrente suivante, pour n ≥ 1 :
n
X
n+1
Bk = 0.
k=0
k
On a par exemple
2B1 + B0 = 0
3B2 + 3B1 + B0 = 0
4B3 + 6B2 + 4B1 + B0 = 0
5B4 + 10B3 + 10B2 + 5B1 + B0 = 0
donc B1 = − 21 , B2 = 61 , B3 = 0 et B4 = −130
. D’une manière générale, Bk , k ≥ 1, est défini en
fonction des Bi pour i ≤ k − 1. On considère alors la série entière
∞
X Bk
S(z) = zk
k=0
k!
pour z ∈ C. On peut montrer que le rayon de convergence de S(z) est 2π et que l’on a la relation
z
∀z ∈ C∗ , S(z) =
ez −1
En effet, en multipliant les deux séries suivantes (par produit de Cauchy) on obtient la relation
! ! n−1
!
X Bn X 1 X Bi X X Bk
zn z n−1 = z i+j−1 = z n−1 = B0 = 1.
n≥0
n! n≥1
n! i≥0,j≥1
i!j! n≥1 k=0
k!(n − k)!
Donc S(z)(ez − 1) = z. En outre, pour 0 < |z| < 2π, on a en particulier z ∈ / 2iπZ∗ , donc
ez − 1 6= 0. La série entière de la fonction holomorphe z/(ez − 1) est donc bien définie pour
|z| < 2π, car en z = 0 on a par continuité z/(ez − 1) = 1. Le rayon de convergence de S(z) est
donc bien 2π.
P Bk k
Notons que la fonction z 7→ S(z) − 1 − B1 z étant paire, la série k≥2 k! z ne comporte que
des termes pairs, de sorte que les B2k+1 sont nuls pour k ≥ 1.
Théorème 6.7 On suppose que kΩk < π. Alors l’opérateur d expΩ est inversible et son inverse
est donné par la formule suivante :
X Bk
d exp−1
Ω = adkΩ . (6.77)
k≥0
k!
56
Preuve. Pour kΩk < π, l’opérateur adΩ est borné strictement par 2π, et donc d expΩ est inversible
(d’après la question 3). En outre, la série
X Bk
adkΩ
k≥0
k!
Preuve. On a
|t|2 2 2 3 3 3 1 2 2
kX(t)k ≤ |t|(kAk+kBk)+ (kAk +kBk +2kAkkBk)+C|t| kAk + kBk + (kAk kBk + kAkkBk )
2 2
donc pour |t| < ǫ, ǫ suffisamment petit, kX(t)k < 1 et la série C(t) (du logarithme) converge. On
a bien sûr que exp(C(t)) = exp(log(I + X(t))) = I + X(t) = I + (exp(tA) exp(tB) − I) =
exp(tA) exp(tB).
Premiers termes. Il vient alors
t2
C(t) = t(A + B) + A2 + 2AB + B 2 − (A + B)2 + O(t3 ),
2
de sorte que C(t) = tC1 + t2 C2 + O(t3 ) avec C1 = (A + B) et C2 = A2 + 2AB + B 2 − (A + B)2 .
De la même manière, pour |s| ≤ ǫ et |t| ≤ ǫ, il existe une fonction (s, t) 7→ Z(s, t) à valeurs
dans Mn (C), différentiable et telle que
Il suffit en effet de poser X(s, t) = exp(sA) exp(tB) − I = O(max(|s|, |t|) et Z(s, t) = log(I +
X(s, t)), puis de remarquer que Z(s, t) est différentiable car log est une série entière et X(s, t) est
différentiable par rapport à s et t.
∂Z 1 X Bk
(s, t) = d exp−1
Z(s,t) (A) = A − adZ(s,t) (A) + adkZ(s,t) (A),
∂s 2 k≥2
k!
et
∂Z 1 X Bk
(s, t) = d exp−1
−Z(s,t) (B) = B + adZ(s,t) (B) + adkZ(s,t) (B).
∂t 2 k≥2
k!
Preuve. En dérivant par rapport à s l’égalité exp(sA) exp(tB) = exp(Z(s, t)), il vient
d exp Z(s, t) ∂Z
A exp(sA) exp(tB) = (s, t)
dZ(s, t) ∂s
∂Z
= d expZ(s,t) (s, t) exp(Z(s, t))
∂s
58
où l’on a utilisé la formule (6.75). Ainsi, en multipliant à gauche les deux membres de l’égalité par
le terme exp(−tB) exp(−sA) = (exp(Z(s, t)))−1 , on obtient
∂Z
A = d expZ(s,t) (s, t)
∂s
puis en appliquant l’opérateur d expZ(s,t)
∂Z
(s, t) = d exp−1
Z(s,t) (A)
∂s
c’est-à-dire, en tenant compte du fait que B0 = 1 et B1 = − 21 :
∂Z 1 X Bk
(s, t) = A − adZ(s,t) (A) + adkZ(s,t) (A).
∂s 2 k≥2
k!
Pour la dérivée par rapport à t, l’astuce consiste à écrire d’écrire tout d’abord l’égalité exp(sA) exp(tB) =
exp(Z(s, t)) sous la forme
exp(−tB) exp(−sA) = exp(−Z(s, t))
par passage aux inverses. En dérivant par rapport à t comme précédemment par rapport à s, puis
en remarquant que
adk−Z(s,t) = (−1)k adZ(s,t)
et en tenant compte du fait que les B2k+1 sont nuls pour k ≥ 1, on obtient la formule de l’énoncé.
La fonction recherchée C(t) peut maintenant s’obtenir comme solution d’un équation différentielle
en remarquant qu’on a bien sûr C(t) = Z(t, t). Il vient ainsi :
∂Z(s, t) ∂Z(s, t)
Ċ(t) = +
∂s s=t ∂t s=t
1 X Bk 1 X Bk
= A − adZ(t,t) (A) + adkZ(t,t) (A) + B + adZ(t,t) (B) + adkZ(t,t) (B)
2 k≥2
k! 2 k≥2
k!
et donc
P Bk k
Ċ(t) = A + B − 21 adC(t) (A − B) + k≥2 k! adC(t) (A + B)
C(0) = 0
Premiers termes. En s’appuyant sur l’identité de Jacobi :
∀(D1 , D2 , D3 ) ∈ Mn (C), adD1 (adD2 (D3 )) + adD3 (adD1 (D2 )) + adD2 (adD3 (D1 )) = 0,
il est facile de vérifier que C(t) = tC1 + t2 C2 + t3 C3 + O(t4 ) avec
C1 = A + B,
1
C2 = adA (B),
2
1
C3 = ad2A (B) + ad2B (A)
12
6. FORMULE DE BAKER-CAMPBELL-HAUSDORFF (BCH) 59
pour s un entier non nul et (a1 , . . . , as ) et (b1 , . . . , bs ) des s-uplets de réels. Ceci revient à approcher
la solution du système différentiel linéaire
ẏ = Ay et ẏ = By.
Définition 6.10 On dit que ψs (t) est d’ordre p, si p est le plus grand entier non nul tel que ψs (t) =
exp(tA) exp(tB) + O(tp+1 ).
Théorème 6.11 Le schéma de Lie est d’ordre 1 et le schéma de Strang est d’ordre 2, i.e.
et
1 1
exp( tA) exp(tB) exp( tA) = exp(t(A + B)) + O(t3 ).
2 2
60
t
Ej (t) = αj A + βj B + γj adA (B) + O(t2 )
2
Le même calcul que pour le schéma de Lie avec A et B remplacées respectivement par a1 A et b1 B
donne
exp(a1 A) exp(b1 B) = exp(tC̃(t)) = exp(tẼ1 (t))
avec C̃(t) = a1 A + b1 B + 2t ada1 A (b1 B) + O(t2 ), soit α1 = a1 , β1 = b1 et γ1 = a1 b1 .
(i) αj = αj−1 + aj
(ii) βj = βj−1 + bj
(iii) γj = γj−1 + aj bj + αj−1 bj − βj−1 aj
Preuve. Supposons la récurrence établie pour i = 1, . . . , j−1. Alors, de ψj = ψj−1 exp(taj A) exp(tbj B),
on tire
ψj = exp(tEj−1 (t)) exp(tC̃(t))
où
t t
C̃(t) = aj A + bj B + aj bj adA (B) + O(t2 ), Ej−1 (t) = αj−1A + βj−1B + γj−1 adA (B) + O(t2 ),
2 2
6. FORMULE DE BAKER-CAMPBELL-HAUSDORFF (BCH) 61
et par application de la BCH, ψj = exp(tEj (t)) avec (modulo les termes en O(t2 ))
t
Ej (t) ≡ Ej−1 (t) + C̃(t) + adEj−1 (t) (C̃(t)),
2
t t t
≡ (αj−1 + aj )A + (βj−1 + bj )B + γj−1adA (B) + aj bj adA (B) + adαj−1 A+βj−1 B (aj A + bj B),
2 2 2
t
≡ (αj−1 + aj )A + (βj−1 + bj )B + (γj−1 + aj bj + bj αj−1 − aj βj−1) adA (B).
2
D’où les relations de récurrence.
Ps
Preuve. Pour j = 2,P . . . , s, on a αj = αj−1 + aj . Comme α1 = a1 , on obtient
Ps α s = Pj=1 aj . De
s s
même, il vient βs = j=1 bj . Enfin, la troisième relation (iii) donne γs = i=1 ai bi + i=1 (αi bi −
βi ai ). Maintenant, on sait que
ψs (t) = exp(t(A + B)) + O(t2 )
Ps
si et seulement si Es (t) = A + B + O(t). Donc ψs (t) est d’ordre au moins 1 ssi
P j=1 aj =
s
j=1 bj = 1. En outre,
ψs (t) = exp(t(A + B)) + O(t3 )
2
Ps
si
Ps et seulement si
Ps Es (t) = A
Ps + B + O(t ). Donc, ψs (t) est d’ordre au moins 2 ssi j=1 aj =
j=1Pbj = 1 et i=1 ai bi + i=1 (αi bi − βi ai ) = 0. En tenant compte des deux premières relations
on a si,j=1 ai bj = 1 de sorte que la troisième
s
X s
X
0= ai bi + (βi ai − αi bi )
i=1 i=1
peut se réécrire !
s j
X X 1
ai bj = .
j=1 i=1
2
Pour atteindre des ordres supérieurs à 3, il convient de calculer les termes d’ordres supérieurs
au moyen de la BCH. Bien que le principe de la procédure soit simple, les calculs deviennent
rapidement compliqués. On reprend donc succinctement la technique de composition décrite plus
avant dans ce cours et on suppose que pour une entier s non nul donné, l’approximation ψs (t) est
d’ordre p, et qu’on peut donc écrire :
ψs (t) = exp(tC̃(t))
62
avec C̃(t) = A + B + tp Cp+1 + O(tp+1 ), où Cp+1 est une matrice non nulle de Mn (C) formée de
commutateurs de A et B, que l’on n’explicitera pas. On considère alors l’approximation suivante :
φ(t) = exp(tD̃(t))
avec
D̃(t) = A + B + tp (µp+1
1 + µp+1
2 + µp+1
3 )Cp+1 + O(t
p+1
).
Preuve. En appliquant la BCH avec les deux matrices µ2 C̃(µ2 t) et µ3 C̃(µ3 t), on a
ψs (µ2 t)ψs (µ3 t) = exp t(µ2 C̃(µ2 t) + µ3 C̃(µ3 t)) + t2 R(t)
où R(t) est une somme de commutateurs comprenant k ≥ 2 matrices parmi A + B, tp Cp+1 et
O(tp+1 ). Si les k matrices sont toutes prises égales à A + B, alors le commutateur est nul. Si au
moins l’une d’entre elles est différente, alors le commutateur est un O(tp ). Donc t2 R(t) = O(tp+2 ).
Maintenant, on a :
µ2 C̃(µ2 t) + µ3 C̃(µ3 t) = (µ2 + µ3 )(A + B) + tp (µp+1
2 + µp+1
3 )Cp+1 + O(t
p+1
)
de sorte que
ψs (µ2 t)ψs (µ3 t) = exp t(µ2 + µ3 )(A + B) + tp+1 (µp+1
2 + µp+1
3 )Cp+1 + O(t
p+2
) .
En répétant l’opération avec ψs (µ1 t) et en tenant compte de la relation µ1 + µ2 + µ3 = 1, on obtient
la formule souhaitée.
µ1 + µ2 + µ3 = 1, µp+1
1 + µp+1
2 + µp+1
3 = 0,
Preuve. Pour avoir D̃(t) = A + B + O(tp+1 ), il suffit de prendre (µ1 , µ2 , µ3 ) comme indiqué.
Si p = 2q + 1 est impair, alors l’équation µ2q+2
1 + µ2q+2
2 + µ2q+2
3 = 0 n’a pas de solution réelle, en
dehors de (0, 0, 0), dont la somme de vaut pas 1. Si p = 2q est pair et µ1 = µ3 , alors l’équation
2µ2q+1
1 + (1 − 2µ1 )2q+1 = 0,
1
admet la solution réelle µ1 = 2−21/(2q+1) , µ3 = µ1 , µ2 = 1 − 2µ1. Pour cette solution, l’approxi-
mation φ(t) est d’ordre au moins p + 1. Si on suppose en outre que ψs (t) satisfait la relation de
symétrie suivante pour tout |t| suffisamment petit :
alors p est nécessairement pair. En effet, en reportant ψs (t) = exp(t(A+B)+tp+1 Cp+1 +O(tp+2 ))
dans (ψs (t))−1 = ψs (−t), on obtient (−1)p+1 = −1. Pour la solution (µ1 , µ2, µ1 ) du théorème
précédent, on a alors également :
(φ(t))−1 = φ(−t),
car
ψs (µ1 t)ψs (µ2 t)ψs (µ1 t)ψs (−µ1 t)ψs (−µ2 t)ψs (−µ1 t) = I
L’ordre de φ(t) est donc pair, donc au moins p+2. En partant de ψ3 (t) = exp( 12 tA) exp(tB) exp( 12 A),
qui est d’ordre 2 et on peut ainsi construire des approximations d’ordres 4, 6, 8...
64
Deuxième partie
65
7. INTRODUCTION 67
Ce cours est une introduction à la théorie de la moyennisation des équations différentielles or-
dinaires dans Rn . Bien qu’il soit possible dans certaines situations d’appliquer (au moins en partie)
les techniques développées ici à des équations ni périodiques, ni quasi-périodiques, on se limitera
à donner les démonstrations dans les cas périodique et quasi-périodique, alors que les cas plus
généraux feront seulement l’objet de remarques. Les résultats présentés pourraient pour la plupart
s’envisager dans le cadre d’une équation différentielle posée dans un espace de Banach (ou de
Hilbert lorsque les aspects géométriques interviennent) sans difficulté conceptuelle majeure, mais
là-encore nous nous limiterons au seul cas de Rn afin de simplifier les développements.
7 Introduction
Considérons une équation différentielle ordinaire de la forme
où D désigne un ouvert connexe d’adhérence compacte. La fonction (y, t, ε) 7→ ftε (y) est supposée
définie sur
Ω = D × R×] − ε0 , ε0 [,
où ε0 est un réel strictement positif. En règle générale, (y, t, ε) 7→ ftε (y) est supposée régulière sur
Ω, mais dans certains cas que l’on précisera, cette hypothèse peut être relaxée. On rappelle qu’en
vertu du théorème de Cauchy-Lipschitz, on a
Théorème 7.1 On considère l’équation différentielle (7.78) et on suppose que la fonction (y, t, ε) 7→
ftε (y) est définie, continue sur Ω, bornée par une constante Mf > 0 sur D × R×] − ε0 , ε0[ et uni-
formément lipschitzienne par rapport à y ∈ D, c’est-à-dire qu’il existe une constante Lf > 0 telle
que
Soit ỹ0 ∈ D. Alors, il existe r > 0 tel que pour tout y0 ∈ Br (ỹ0 ), pour tout |ε| < ε0 , il existe une
unique solution y ε (t) ∈ D de (7.78) de classe C 1 en temps sur un intervalle de la forme [0, Tε [. En
T T
outre il existe T > 0 tel que, pour tout 0 < |ε| < ε0 , Tε > |ε| avec y ε (t) ∈ D pour tout 0 ≤ t ≤ |ε| .
68
On utilisera fréquemment le lemme de Gronwall que l’on rappelle ci-après à toutes fins utiles :
Lemme 7.2 Soit I un intervalle de la forme [a, +∞[, [a, b] ou [a, b[ avec a < b. On considère
trois fonctions α, β et u à valeurs réelles définies et continues sur I. Si β est positive et u satisfait
l’inégalité intégrale Z t
u(t) ≤ α(t) + β(s)u(s) ds, ∀t ∈ I,
a
alors Z Z t
t
u(t) ≤ α(t) + α(s)β(s) exp β(r) dr ds, t ∈ I.
a s
En outre, si α est croissante, alors
Z t
u(t) ≤ α(t) exp β(s) ds , t ∈ I.
a
8. EXAMPLES 69
Equation logistique
2
1.8
1.6
y et z
1.4
1.2
0.8
0 10 20 30 40 50 60 70 80 90 100
Temps
8 Examples
On voit donc ici que ftε (y) = y(1 − y) + sin(t) est périodique de période P = 2π et il est facile de
vérifier que (y, t, ε) 7→ ftε (y) est bornée sur D × R × R où D = BR (y0 ) = {y ∈ R, |y − y0 | < R}
par Mf = (|y0|+R)(|y0|+R+1) et lipschitzienne (ou Lipschitz pour plus de brièveté) de constante
de Lipschitz Lf = 1 + 2(|y0| + R). Il existe donc une solution unique y ε (t) ∈ D de l’équation
R
logistique sur l’intervalle [0, |ε|M f
[ pour tout ε non nul (et bien sûr sur [0, +∞[ pour ε = 0). On
peut ici aisément obtenir hf i(z) = z(1 − z), de sorte que l’équation moyennée (à l’ordre 1) s’écrit
simplement
ż ε = εz ε (1 − z ε ), z ε = y0 ,
R
dont la solution existe (et reste dans D) sur [0, |ε|M f
[. En outre, z ε se calcule explicitement et l’on
a:
y0 eεt
z ε (t) = .
1 + y0 (eεt − 1)
70
2 2
1.5 1.5
1 1
0.5 0.5
0 0
p
p
−0.5 −0.5
−1 −1
−1.5 −1.5
−2 −2
−2.5 −2.5
−2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5
q q
F IGURE 11 – Solutions de l’équation de Van Der Pol originale (rouge) et de sa version moyennée
(bleu)
Les équations du mouvement se réduisent alors à une équation sur l’angle θ du second ordre :
d2 ε g 1 vmax t
2
θ = + cos + ϕ0 sin θε
dt l ε l ε
Afin de simplifier l’étude et les notations, on suppose dans la suite que ϕ0 = 0 et on pose α = g/l
et β = vmax /l. Par changement d’échelle de temps τ = εt , on obtient l’équation
d2 θ ε
2
= θ̈ε = ε(εα + β cos(τ )) sin(θε )
dτ
et en posant q ε = θε et pε = 1ε θ̇ε , on est ramené à un système d’ordre 1 de la forme :
q̇ ε = εpε
ṗε = β cos(τ ) sin(q ε ) + εα sin(q ε )
Finalement, en posant y1ε = q ε et y2ε = −β sin(τ ) sin(q ε ) + pε , il vient
ẏ1ε = εpε = ε(y2ε + β sin(τ ) sin(y1ε))
ẏ2ε = −β cos(τ ) sin(y1ε ) − β sin(τ ) cos(y1ε )ẏ1ε + ṗε
= εα sin(y1ε ) − εβ sin(τ ) cos(y1ε )(y2ε + β sin(τ ) sin(y1ε))
Il s’agit donc d’un système de la forme (7.78) périodique en temps de période P = 2π, avec un
champ de composante
ε ε β2
ft,1 (y) = y2 + β sin(τ ) sin(y1 ), ft,2 (y) = α sin(y1) − sin2 (τ ) sin(2y1) − β sin(τ ) cos(y1 )y2 ,
2
dont la moyenne s’écrit :
β2
hf i1 (z) = z2 , hf i2 (z) = α sin(z1 ) − sin(2z1 ).
4
L’équation moyenne est donc de la forme
ε β2
Z̈ = ε α − cos(Z ) sin(Z ε ).
ε
2
2
Cette équation montre qu’une force d’expression − β2 cos(Z ε ) s’oppose à la gravité. Cette force
permet de stabiliser l’équilibre Z ε = 0 pour β suffisamment grand (voir figure 13).
0.1
0.05
0
θ
−0.05
−0.1
−0.15
−0.2
0 10 20 30 40 50 60 70 80 90 100
Temps
Preuve. [Première preuve] Désignons par E ε (t) = y ε (t) − z ε (t) l’erreur d’approximation. Il vient
Z t
ε
E (t) = ε (fsε (y ε (s)) − hf i(z ε (s)) ds
Z0 t Z t
ε ε ε ε
= ε (fs (y (s)) − fs (z (s))) ds + ε (fsε (z ε (s)) − hf i(z ε (s)) ds
0 0
:= E1ε (t) + E2ε (t).
74
Pour le second terme, on découpe l’intervalle [0, t] en [0, nP ] ∪ [nP, t] où n est la partie entière de
t/P . L’erreur E2ε (t) se décompose ainsi en
n−1 Z (i+1)P
X Z t
ε ε ε ε
E2 (t) = ε (fs (z (s)) − hf i(z (s)) ds + ε (fsε (z ε (s)) − hf i(z ε (s)) ds.
i=0 iP nP
RP
En notant ∆εs (z) = fsε (z) − hf i(z) et en remarquant que 0 ∆εs (z)ds = 0, le terme d’indice i de
la somme peut alors être réécrit de la manière suivante :
Z (i+1)P Z P
ε ε
∆s (z (s))ds = ∆εs (z ε (s + iP ))ds
iP 0
Z P
= (∆εs (z ε (s + iP )) − ∆εs (z ε (iP ))) ds
0
et donc estimé par
Z P Z P
∆εs (z ε (s + iP ))ds ≤ k∆εs (z ε (s + iP )) − ∆εs (z ε (iP ))kds
0 0
Z P
≤ 2Lf kz ε (s + iP ) − z ε (iP )kds.
0
ε ε
Or kz (s + iP ) − z (iP )k ≤ s|ε|Mf , de sorte que
kE2ε (t)k ≤ nP 2 ε2 Mf Lf + (t − nP )2 ε2 Mf Lf ≤ ε2 tMf Lf P
Finalement, on a l’inégalité
Z t
ε
kE (t)k ≤ |ε|Lf kE ε (s)kds + tε2 Mf Lf P
0
2
et le lemme de Gronwall (avec α(t) = tε P Mf Lf fonction croissante et β(t) = |ε|Lf > 0) donne
alors
kE ε (t)k ≤ ε2 tP Mf Lf e|ε|Lf t ≤ |ε|P Mf Lf T̃ eLf T̃ .
où la fonction (x, t, ε) 7→ uεt (x) est supposée de classe C 1 par rapport à x et continue par rapport
à (x, t, ε). On impose en outre que ce changement de variable soit périodique en t de période P .
L’idée est de transformer l’équation (7.78) en une équation différentielle sur x de la forme
ẋε = εhf i(xε ) + ε2 Rtε (xε ), xε (0) = (U0ε )−1 (y0 ) (9.84)
9. MOYENNISATION AU PREMIER ORDRE EN ε 75
où la fonction (x, t, ε) 7→ Rtε (x) est un reste périodique en t, de période P , que l’on va chercher à
borner, ce qui permettra d’estimer z ε − xε , puis y ε − z ε . En remplaçant y ε via (9.83), il vient
de sorte que si (I + ε∂x uεt (xε )) est supposé inversible le long de la solution xε , alors
ẋε = ε(I + ε∂x uεt (xε ))−1 (ftε (xε + εuεt (xε )) − ∂t uεt (xε ))
On voit donc que pour transformer (7.78) en (9.84) via (9.83), uεt doit être solution de l’équation
où κ(x) est une fonction quelconque régulière. Notons que la fonction uεt ainsi définie est bien
périodique de période P en temps car fsε (x) − hf i(x) est de moyenne nulle.
ce qui implique x2 = x1 . La fonction x 7→ Utε (x) est donc injective de D dans Utε (D), donc
bijective. En outre, en appliquant le théorème des fonctions implicites au voisinage d’un point
(x, t, ε) de D × R×] − ε1 , ε1 [, on montre que Vtε = (Utε )−1 est régulière et de la forme :
Vtε (y) = y + εvtε (y).
Etape 2 : Soit donc y ε = xε + εuεt (xε ) : les calculs effectués plus haut donnent
(I + ε∂x uεt (xε )) · ẋε = ε (ftε (xε + εuεt (xε )) − ∂t uεt (xε )))
= εhf i(xε ) + ε ftε (xε + εuεt (xε )) − ftε (xε )
c’est-à-dire
ẋε = εhf i(xε ) + ε2 Rtε (xε ) (9.89)
où
Z 1
Rtε (x) = (I + ε∂x uεt (xε ))−1 ∂x ftε (x + ε(1 − s)uεt (x)) · uεt (x)ds − ∂x uεt (xε ) ε
· hf i(x )
0
Notons que l’inverse de l’application (x, t, ε) 7→ I + ε∂x uεt (x) ∈ L(Rn ) est bien définie car
kε∂x uεt (x)kL(Rn ) ≤ |ε|Lu < 1
dès lors que |ε| < ε1 . En outre, Rtε est clairement périodique de période P en temps, puisque uεt
l’est, et est régulière de classe C 0 en (x, t, ε), bornée sur un sous-domaine de D × R×] − ε1 , ε1[
contenant xε (t) pour t < T̃ε par une constante MR > 0.
Remarque 9.2 L’approximation Utε (z ε (t)) de y ε (t) est dite approximation “améliorée” de y ε(t).
En un sens qui apparaitrait plus clairement aux ordres plus élevés en ε, c’est une meilleure approxi-
mation que z ε (t). Cependant à l’ordre un en ε qui est la cas considéré ici, les deux approximations
sont en fait de même ordre puisque Utε (z) = z + O(ε).
Remarque 9.3 Plusieurs choix de la fonction κ sont possibles. Le choix consistant à rendre uεt
de moyenne nulle est bien souvent considéré comme le choix standard en particulier lorsque les
fonctions périodiques sont écrites sous la forme de leurs développements de Fourier. Néanmoins,
le choix κ ≡ 0 est le seul choix qui permette d’assurer la conservation des structures géométriques
de l’équation de départ. Il possède l’avantage que les conditions initiales des équations en y, x et
z sont identiques et ne dépendent pas de ε, puisqu’on a alors U0ε = id. En outre, pour les temps
dits stroboscopiques t = mP , on a également Utε = id de sorte que les solutions y ε (t) et Utε (xε (t))
coı̈ncident. Pour cette raison, la moyennisation correspondant à κ ≡ 0 est appelée moyennisation
stroboscopique. C’est celle que nous présenterons dans ce cours.
Lorsque cette limite est bien définie et uniforme en y sur tout compact de D, alors il convient de
calculer la fonction suivante, dite fonction d’ordre :
Z t
δ(ε) = sup sup ε (fsε (y) − hf i(y)) ds .
y∈D t∈[0,T /ε) 0
Théorème 9.4 On suppose de que f est de classe C 1 ou au moins Lipschitz. Soit z ε (t) la solution
de l’équation
ż ε = εhf i(z), z ε (0) = x0 .
L’étude des conditions dans lesquelles δ est bien définie et son calcul sortent du cadre de ce cours
introductif. Nous renvoyons le lecteur intéressé à l’ouvrage [19].
78
Il est très facile de vérifier que ces propriétés sont conservées par hf i. On a en effet :
— si (y, t, ε) 7→ ftε (y) est hamiltonien, alors
Z Z 1 Z P
1 P ε 1 P −1 −1
hf i(z) = fs (z)ds = ε
J ∇z Hs (z)ds = J ∇z Hsε (z)ds = J −1 ∇z hHi(z)
P 0 P 0 P 0
R P
avec hHi(z) = P1 0 Hsε (z)ds.
— si (y, t, ε) 7→ Itε (y) est un invariant du champ (y, t, ε) 7→ ftε (y) alors I0ε est un invariant à
ε-près de hf i, c’est-à-dire que l’on a
∀z ∈ Rn , ∂z I0ε (z) · hf i(z) = O(ε)
Il vient en effet, pour tout z ∈ Rn
Z P Z P
1 ε
ε∂z I0ε (z) · hf i(z) = ε∂z I0ε (z) · fsε (z)ds
∂z I0ε (z) · fsε (z)ds
=
0 P 0 P
Z P Z
ε ε ε ε 1 P
= (∂z I0 (z) − ∂z Is (z)) · fs (z)ds + ε∂z Isε (z) · fsε (z)ds
P 0 P 0
Z Z Z
ε P s ε ε 1 P
= ∂µ ∂z Iµ · fs (z)dµds + − (∂s Isε (z)) ds
P 0 0 P 0
Z Z
ε2 P s 2 ε ε (I ε (z) − IPε (z))
= ∂z Iµ fs (z), fµε (z) dµds + 0 = O(ε2 )
P 0 0 P
10. MÉTHODE DE MOYENNISATION À UN ORDRE QUELCONQUE : CAS PÉRIODIQUE79
— si (y, t, ε) 7→ ftε (y) est à divergence nulle, alors il vient là-encore très aisément
Z
1 P
div(hf i(z)) = div(fsε (z))ds = 0.
P 0
Dans un premier temps, nous allons déterminer de manière formelle les équations satisfaites par
Φεθ et Ψεt . En différentiant les deux membres de l’équation (10.92) par rapport à t et en utilisant
l’équation (10.91) , on obtient immédiatement
∂Φεt ∂Φε
(Ψεt (y0 )) + ε t (Ψεt (y0 )) F ε (Ψεt (y0 )) = ε ft (Φεt ◦ Ψεt (y0 )) . (10.93)
∂t ∂y
En considérant (10.93) pour y0 = Ψε−t (y) et en remplaçant t par θ ∈ T (rappelons que ft et Φεt
sont périodiques en temps) il vient
∂Φεθ ∂Φε
(y) + ε θ (y) F ε (y) = ε fθ (Φεθ (y)) . (10.94)
∂θ ∂y
En prenant la moyenne des deux membres de l’équation (10.94), on obtient alors
∂hΦε i
(y) F ε (y) = hf ◦ Φε i (y) ,
∂y
où l’on a tenu compte du fait que la moyenne de la dérivée d’une fonction périodique est nulle, ici
∂hΦε i
h∂θ Φεθ i = 0. En supposant que l’opérateur ∆y 7→ (y) ∆y est inversible pour tout y, on peut
∂y
alors écrire
−1
ε ∂hΦε i
F (y) := (y) hf ◦ Φε i (y). (10.95)
∂y
En d’autres termes, nous avons obtenu ici l’expression du champ moyen F ε et donc du flot Ψεt
associé. En insérant cette expression dans l’équation (10.94), il s’ensuit que
−1
∂Φεθ ∂Φεθ ∂hΦε i
(y) + ε (y) (y) hf ◦ Φε i (y) = ε fθ ◦ Φεθ (y), (10.96)
∂θ ∂y ∂y
Il s’agit là dune équation fermée sur Φεθ , bien que non-linéaire et non-locale. En résumé, nous
avons obtenu de manière formelle les équations suivantes :
La première est le pendant à un ordre quelconque de l’équation moyennée (7.81) à l’ordre 1. C’est
elle qui régit la dynamique du système sur les temps longs (en O(1/ε)). La seconde équation four-
nit l’expression du champ moyen en fonction du changement de variables. Résoudre la troisième
équation permet de déterminer le changement de variables lui-même à partir de la donnée f du
problème. C’est cette équation qui concentre la quasi-totalité des difficultés du problème. On sait
en particulier qu’elle ne possède pas de solution générique : il existe en effet des contre-exemples
dans la littérature (voir [17, 8]). En fait, il est impossible d’obtenir mieux qu’une solution ap-
prochée (c’est-à-dire ici à un terme d’erreur près de taille O(exp(−c/ε))), et c’est donc l’objectif
visé ici.
où (·, ·) désigne le produit scalaire sur Rn . Notons que pour tout (y, z) ∈ Rn × Rn , on a (y, z)Cn =
(y, z) et kykCn = kyk, de sorte que nous noterons identiquement (·, ·) et k · k les produit scalaire
et norme sur Rn et Cn .
Si K désigne un ouvert borné de Rn , on considère pour ρ > 0 l’extension ouverte de K dans
n
C définie par
Kρ = {y + ∆y : (y, ∆y) ∈ K × Cn , k∆yk < ρ},
y + ∆y
k∆yk < ρ
y
×
K ⊂ Rn
ρ
Kρ ⊂ C n
Définition 10.1 Soit une fonction (y, θ) ∈ Kρ × T 7→ gθ (y) ∈ Cn . L’application gθ est dite
analytique sur Kρ ⊂ Cn si elle est continûment differentiable sur Kρ , c’est-à-dire si il existe une
application linéaire continue
Kρ × T → L(Cn )
(y, θ) 7→ (∂y gθ )(y)
où L(Cn ) désigne l’ensemble des applications linéaires de Cn dans Cn , telle que
Si gθ est analytique sur Kρ au sens de la définition 10.1, la formule de Cauchy permet d’écrire pour
tout 0 < δ < ρ, tous y et ∆y tels que y ∈ Kρ−δ et ∆y ∈ Cn avec k∆yk = 1, et tout µ ∈ C avec
|µ| < δ Z
1 gθ (y + ξ∆y)
gθ (y + µ∆y) = dξ.
2iπ |ξ|=δ ξ−µ
Cette formule permet de recourir aux estimations dites “de Cauchy” au voisinage d’un point y ∈
Kρ . Par exemple, étant donné 0 < δ < ρ, on peut estimer la norme de sa dérivée première ∂y gθ par
1
k∂y gkρ−δ ≤ kgkρ, (10.98)
δ
inégalité qui découle de la relation
Z
d 1 gθ (y + ξ∆y)
(∂y gθ )(y)∆y = gθ (y + µ∆y) = dξ.
dµ µ=0 2iπ |ξ|=δ ξ2
Au passage, on retrouve le fait connu que ∂y gθ est également analytique au sens de la définition
10.1 comme fonction de Kρ dans L(Cn ). On peut désormais énoncer précisément les hypothèses
de régularité faites sur f :
Hypothèse 10.2 La fonction (y, θ) 7→ fθ (y) est C 0 et périodique en θ. En outre, (y, θ) 7→ fθ (y)
est réelle-analytique en y, c’est-à-dire qu’il existe
R > 0, CK > 0,
tels que, pour tout θ ∈ T, y 7→ fθ (y) est analytique sur K2R au sens de la définition 10.1, et
(y, θ) 7→ fθ (y) est bornée par CK sur K2R × T
Remarque 10.3 Dès lors que ϕθ est périodique, Λ(ϕ)θ est de moyenne nulle, de sorte que Γε (ϕ)θ
est à son tour périodique.
comme une équation de point fixe, que l’on cherche à résoudre par itération (Γε )k , k = 0, 1, . . . , n.
De manière équivalente, sous sa forme différentielle (10.96), elle peut s’interprèter comme
Remarque 10.4 Par construction (cela sera justifié de manière rigoureuse au moyen du théorème
10.12), pour tous n, k ≥ 0,
1 dk F [n+k]
(y) = Fk+1 (y).
k! dεk ε=0
84
Avant de conclure ce paragraphe avec l’énoncé du théorème principal, on rassemble ci-dessous les
notations utilisées :
1 ε0 P
1. ε0 := 8CRK P , ε1 := ε0 /2, ε2 = ε0 min 28 , T ,
2. C0 = 16CK , C1 = 2CK , η = 2/ε0,
R
3. rn := n+1
, Rk = 2R − krn pour k = 1, . . . , n − 1, de sorte que R0 = 2R et Rn+1 = R.
Rk = 2R − krn
K R×
K2R 2R
[k]
Lemme 10.5 Etant donné un entier n ∈ N, les applications Φθ et Fk+1 avec 0 ≤ k ≤ n + 1 sont
définies pour tout ε ∈ C tel que |ε| ≤ ε0 /(n + 1), sont de classe C 1 en θ ∈ T, analytiques en
y ∈ KRk , et analytiques en ε pour |ε| < ε0 /(n + 1). En outre, l’estimation suivante est valable
pour tout ε ∈ C tel que (n + 1)|ε| < ε0 :
rn
kΦ[k] − idkRk ≤ . (10.106)
2
Pour 0 < ε < ε2 /2, on définit l’entier nε ∈ N∗ par l’expression (nε + 1) = E(ε0 /(4|ε|)) et on note
n
X
Feε (y) = Fe [nε] (y), avec Fe [n](y) = εk Fk+1 (y), (10.107)
k=0
où les champs de vecteur Fk+1 (y) sont définis par l’équation (10.105).
Théorème 10.6 Soient ε tel que 0 < ε < min(ε∗ , ε2/2) et (nε + 1) = E(ε0 /(4|ε|)). Alors la
fonction Φe ε = Φ[nε ] définie dans le lemme 10.5 et la fonction Fe ε définie par l’équation (10.107)
θ θ
sont telles que :
(i) La solution y ε (t) de (10.90) vérifie
et où (x, θ) 7→ Rεθ (x) est périodique et régulière en θ, analytique sur KR , et admet la borne
supérieure suivante
ε ε0 log(2)
kR kR ≤ 32 CK exp − . (10.110)
4 |ε|
e ε le flot de l’équation différentielle autonome
(ii) Si l’on considère Ψ t
ż ε = εFeε (z ε ), (10.111)
alors la solution y ε(t) du problème hautement oscillant d’origine (10.90) est exponentiellement
proche de la forme composée Φ eε ◦ Ψ
e ε (y0 ). Plus précisément, on a
t t
ε e ε e ε | log(2)|ε0
∀t ∈ [0, T /|ε|], y (t) − Φt ◦ Ψt (y0 ) ≤ 48 R exp − . (10.112)
8 |ε|
kΛ(ϕ)kρ−δ ≤ 4 CK ,
et en intégrant en θ,
kΓε (ϕ) − idkρ−δ ≤ 4 CK P |ε|.
L’analyticité des fonctions Λ(ϕ)θ et Γε (ϕ)θ sur Kρ−δ est assurée par les théorèmes de composition
habituels.
Le lemme 10.7 montre que, partant d’une fonction (y, θ) ∈ K2R ×T 7→ ϕθ (y) ∈ Cn et pour peu
que ε soit suffisamment petit, il est possible de construire la suite (Γε )k (ϕ)θ , mais au prix d’une
contraction graduelle des domaines de définition (et d’analyticité). Afin d’estimer la convergence
de cette suite, on montre maintenant la propriété de contraction suivante :
Lemme 10.8 Soient δ et ρ tels que 0 < δ < ρ ≤ 2R. On considère deux fonctions (y, θ) ∈
Kρ × T 7→ ϕθ (y) et (y, θ) ∈ Kρ × T 7→ ϕ̂θ (y), analytiques sur Kρ et telles que
δ δ
kϕ − idkρ ≤ et kϕ̂ − idkρ ≤ .
2 2
Alors les inégalités suivantes sont satisfaites pour |ε| ≥ 0 :
C0 C0 P |ε|
kΛ(ϕ) − Λ(ϕ̂)kρ−δ ≤ kϕ − ϕ̂kρ et kΓε (ϕ) − Γε (ϕ̂)kρ−δ ≤ kϕ − ϕ̂kρ .
δ δ
Preuve du lemme 10.8. Afin d’alléger la présentation, on note Aθ (y) = ∂y ϕθ (y) et Âθ (y) =
∂y ϕ̂θ (y). On a tout d’abord
En procédant comme dans le lemme 10.7, on obtient kAkρ−δ ≤ 23 , ainsi que khAi−1 kρ−δ ≤ 2,
et des estimations similaires pour Â. L’identité hAi−1 − hÂi−1 = hAi−1 h − AihÂi−1 , et une
estimation de Cauchy permettent d’obtenir
4
khAi−1 − hÂi−1 kρ−δ ≤ khAi−1 kρ−δ khÂi−1 kρ−δ kA − Âkρ−δ ≤ kϕ − ϕ̂kρ .
δ
10. MÉTHODE DE MOYENNISATION À UN ORDRE QUELCONQUE : CAS PÉRIODIQUE87
Finalement, étant donné que kϕ−idkρ ≤ δ/2 et kϕ̂−idkρ ≤ δ/2, il est clair que, dès que y ∈ Kρ−δ ,
ϕθ (y) ∈ Kρ−δ/2 et ϕ̂θ (y) ∈ Kρ−δ/2 , de sorte que
khf ◦ ϕ̂ikρ−δ ≤ CK
et
2CK
kf ◦ ϕ − f ◦ ϕ̂kρ−δ ≤ k∂y f kρ−δ/2 kϕ − ϕ̂kρ−δ ≤ kϕ − ϕ̂kρ−δ .
δ
En regroupant tous les termes, il vient finalement
2CK 6CK 6CK 2CK
kΛ(ϕ) − Λ(ϕ̂)kρ−δ ≤ + + + kϕ − ϕ̂kρ
δ δ δ δ
4 CK P ε 0 R rn
kΦ[k+1] − idkRk+1 ≤ 4CK P |ε| < = = .
n+1 2(n + 1) 2
Comme composition de fonctions analytiques, elle est de plus analytique en ε pour |ε| < ε0 /(n +
1). La régularité en θ est également claire, ce qui complète la preuve par induction pour Φ[k] .
Par définition de (10.104), F [k] est alors aussi analytique sur KRk de sorte que Fk+1 , définie
par(10.105), également.
[k]
et la suite associée des défauts (δθ )k=0,...,n , définie par
[k]
[k] ∂Φθ
εδθ (y) := (y) − εΛ(Φ[k] )θ (y) (10.113)
∂θ
[k] [k] −1
∂Φθ ∂Φθ ∂hΦ[k] i [k]
= (y) + ε (y) (y) hf ◦ Φ[k] i(y) − εfθ ◦ Φθ (y).
∂θ ∂y ∂y
Alors, les assertions suivantes sont vérifiées :
[n] [n]
(i) Les fonctions Φθ et δθ sont C 1 en θ, analytiques respectivement sur KR+rn et KR , et
analytiques en ε ∈ C pour |ε| < ε0 /(n + 1).
[n] [n]
(ii) Les fonctions Φθ et δθ vérifient les estimations suivantes pour tout |ε| < ε0 /(n + 1) :
rn
kΦ[n] − idkR+rn ≤ , (10.114)
2
kδ [n] kR ≤ C1 (η(n + 1)|ε|)n . (10.115)
[n]
(iii) Pour tout θ ∈ T, l’application y 7→ Φθ (y) possède un inverse défini sur KR à valeurs
dans KR+rn . Cet inverse (Φ[n] )−1 est analytique. En outre, on a k(Φ[n] )−1 kR−rn ≤ R
[n]
Remarque 10.10 En d’autres termes, le n-ième itéré Φθ satisfait l’équation (10.96) à un reste
près de taille C1 ε(η(n + 1)ε)n .
Preuve du théorème 10.9. Le point (i) et l’estimation (10.114) découlent aisément du lemme 10.5,
possiblement à une singularité près en ε = 0, qui peut être éliminée en raison de l’estimation
[n]
(10.115), que nous allons maintenant prouver. D’une part, en différentiant Φθ = Γε (Φ[n−1] )θ , il
vient pour n ≥ 1 :
[n]
εδθ = ε Λ(Φ[n−1] )θ − Λ(Φ[n] )θ .
En conséquence, la propriété de contraction établie dans le lemme 10.8 conduit, pour δ = rn et
ρ = R + δ = Rn , à l’inégalité :
C0 |ε| [n] 1 C0 P |ε| ε [n−1]
|ε|kδ [n]kR ≤ kΦ − Φ[n−1] kRn = Γ Φ − Γε Φ[n−2] Rn
rn P rn
2
1 C0 P |ε|
≤ Γε Φ[n−2] − Γε Φ[n−3] Rn−1 ≤ . . .
P rn
n
1 C0 P |ε|
≤ kΦ[1] − Φ[0] kR1 .
P rn
[1]
D’autre part, par définition de Φθ , on a kΦ[1] − Φ[0] kR1 ≤ 2|ε| P CK , ce qui prouve (10.115).
[n]
En ce qui concerne le point (iii), il est clair que si l’on prend y1 , y2 ∈ KR tels que Φθ (y1 ) =
[n]
Φθ (y2 ), alors on a :
1
ky1 − y2 k ≤ k∂u Φ[n] − IdkR ky1 − y2 k ≤ kΦ[n] − idkR+rn ky1 − y2 k
rn
1
≤ ky1 − y2 k
2
10. MÉTHODE DE MOYENNISATION À UN ORDRE QUELCONQUE : CAS PÉRIODIQUE89
de sorte que y1 = y2 , ce qui prouve l’injectivité. Afin de démontrer l’existence d’un antécédent,
étant donnés (y, ỹ) ∈ K × Cn avec ρ := kỹkCn < R, on considère la suite vk définie par
[n]
v0 = y + ỹ ∈ KR , vk+1 = vk − Φθ (vk ) + y + ỹ.
Tout d’abord, pour vk ∈ K(R+ρ+rn )/2 on a
[n] R + ρ + rn
kvk+1 − yk ≤ ρ + kΦθ − idkR+rn < ,
2
de sorte qu’on peut montrer par récurrence que l’ensemble des éléments de la suite (vk )k∈N appar-
tiennent à K(R+ρ+rn )/2 . On peut alors écrire
[n] [n]
kvk+1 − vk k = k(Φθ − id)(vk ) − (Φθ − id)(vk−1 )k
1
≤ k∂u Φ[n] − Idk(R+ρ+rn )/2 kvk − vk−1 k ≤ kvk − vk−1 k,
1 + (R − ρ)/rn
et donc la suite vk converge vers une limite v ∈ K(R+ρ+rn )/2 ⊂ KR+rn . Si ρ < R − rn , on a de
[n]
plus kvk < R. L’analyticité de (Φθ )−1 est une consequence directe du théorème des fonctions
implicites.
[k]
Remarque 10.11 A partir de la définition de la suite(Φθ )k=0,··· ,n+1, il est possible de montrer par
[n] [n]
récurrence que si y appartient à KR ∩ Rn et ε ∈ R, alors Φθ (y), F [n] (y) et δθ (y) appartiennent
à Rn également.
[n]
Le théorème suivant établit que les n premiers termes du développement de Φθ en puissances
de ε, sont independants du mode de construction utilisé (une propriété similaire prévaut claire-
ment pour F [n] aussi). Cette indépendance valide le caractère “canonique” des équations (10.96)
et (10.95).
Théorème 10.12 [Unicité des quasi-solutions de (10.96)] Pour n ∈ N fixé, on considère la
fonction (y, θ) 7→ Φ b θ (y), qui est C 1 en θ ∈ T, analytique sur KR+rn , analytique en ε pour
|ε| < ε0 /(n + 1), et satisfait
b 0 = id,
Φ kΦb − idkR+rn ≤ rn . (10.116)
2
On suppose que le défaut associé à Φ b θ , défini par
∂Φbθ
b
εδθ := − εΛ Φ b
∂θ θ
bR ≤ C
vérifie, pour tout |ε| < ε0 /(n + 1), l’inégalité kδk b |ε|n pour une constante C
b > 0
indépendante de ε. Alors on a nécessairement pour tout |ε| < ε0 /(4(n + 1))
b − Φ[n] krn ≤ C3 (n) |ε|n+1,
kΦ
[n] b + 20CK (η (n + 1))n .
où Φθ est la fonction définie dans le théorème 10.9 et C3 (n) = P 2C
90
b − idkR+rn −krn ≤ rn ,
k (Γε )k (Φ) k (Γε )k (Φ[n] ) − idkR+rn −krn ≤
rn
, (10.117)
2 2
rn
où on a utilisé le fait que 4 CK P |ε| ≤ 2
. Ces inégalités permettent d’appliquer le lemme 10.8 et
on obtient, pour tout k ≤ n,
k
ε k+1 b − (Γε )k (Φ)k
b R−krn ≤ C0 P |ε| b R ≤ 1 kΓε (Φ)
b − Φk b − Φk
b R
k (Γ ) (Φ) kΓε (Φ)
rn 2k
ainsi qu’une inégalité similaire pour Φ[n] , où l’on a eu recours à l’inégalité C0 P |ε|/rn ≤ 1/2. Par
sommation, il vient alors
n
X
k (Γ ) ε n+1 b − Φk
(Φ) b rn ≤ b − (Γε )k (Φ)k
k (Γε )k+1 (Φ) b − Φk
b rn ≤ 2 kΓε (Φ) b R, (10.118)
k=0
et de même pour Φ[n] . D’autre part, en utilisant la borne (10.117) et en appliquant une fois encore
n + 1 fois le lemme de contraction 10.8, on aboutit à
n+1 n+1
ε n+1 b − (Γε )n+1 (Φ[n] )krn ≤ C0 P |ε| b − Φ[n] kR+rn ≤ C0 P |ε|
k (Γ ) (Φ) kΦ rn ,
rn rn
(10.119)
où la derrière inégalité repose sur (10.116) (de même pour Φ[n] ). Finalement, de l’inégalité (10.118)
portant sur Φ b (et sur Φ[n] ), et de (10.119), on déduit
b − Φ[n] krn
kΦ
≤ kΦb − (Γε )n+1 (Φ)k b − (Γε )n+1 (Φ[n] )krn
b rn + kΦ[n] − (Γε )n+1 (Φ[n] )krn + k (Γε )n+1 (Φ)
n
≤ 2kΦ b R + 2kΦ[n] − Γε (Φ[n] )kR + C0 P C0 P (n + 1)
b − Γε (Φ)k |ε|n+1 (10.120)
R
n
b n+1 n C0 P
≤ 2C P |ε| + 2C1 P (η(n + 1)|ε|) ε + C0 P (n + 1) |ε|n+1,
R
ε0
Preuve. Partie (i). Puisque (nε + 1)|ε| ≤ 4
< ε0 , le théorème 10.9 s’applique avec n = nε et on a
e εθ (y) = ε fθ ◦ Φ
∂θ Φ e εθ (y) F [nε] (y) + ε δ [nε] (y),
e εθ (y) − ε ∂y Φ
θ
dxε (t)
= ε F̃ ε (xε (t)) + Rεt (xε (t)) , xε (0) = y0 ,
dt
−1
[n ]
avec Rεt (x) = F [nε] (x) − Feε (x) − ∂x Φ e εt (x) e εt (xε (t)) vérifie
δt ε (x), la fonction y ε (t) := Φ
clairement y ε (0) = y0 ainsi que l’équation
dy ε e ε (xε (t)) · F [nε ] (xε (t)) + ε δt[nε] (xε (t))
e ε )(xε (t)) − ε ∂x Φ
(t) = ε (ft ◦ Φ t t
dt
−1
e (x (t)) · F ε (x (t)) − ∂x Φ
ε ε [n ] ε e (x (t))
ε ε [n ε ] ε
+ ε ∂x Φ t t δt (x (t))
comme souhaité. Donc, y ε (t) coı̈ncide pour tout temps t ∈ [0, T /|ε|] avec la solution de (10.90).
[n ]
Par ailleurs, le théorème 10.9 et le choix de nε assurent que le défaut δθ ε vérifie
nε
[nε ] nε 1
kδ kR ≤ C1 (η(nε + 1)ε) ≤ C1 .
2
En outre, l’analyticité de F [nε ] par rapport à ε et la formule de Cauchy permettent d’écrire, pour
tout y ∈ KR et δ := 2(nεε0+1) , que
X εk dk F [nε]
kF [nε] (y) − Fe ε (y)k = (y)
k≥n +1
k! dεk ε=0
ε
X |ε|k k! nε
[nε ] (|ε|/δ)nε+1 [nε ] 1
≤ k
sup kF (y)k ≤ sup kF kR ≤ C1 ,
k! δ |ε|<δ 1 − (|ε|/δ) |ε|<δ 2
k≥nε +1
−1
où l’on a utilisé l’inégalité ∂y hΦ[nε ] i ≤ 2 pour borner 6 F [nε ] sur KR par C1 et |ε|/δ ≤ 21 .
R
6. En utilisant le lemme 10.7-(i) et l’inégalité kΦ[nε ] − idkR+rnε ≤ rnε /2 (avec rnε = R/(nε + 1)), comme
établie dans le théorème 10.9.
92
dz ε
= ε Fe ε (z ε ).
dt
e ε (y0) soit défini pour tout 0 ≤ t ≤ T1 /ε, car Feε est analytique sur KR ,
Il existe T1 > 0 tel que Ψ t
donc Lipschitz et continue sur tout Kρ pour ρ < R. On a d’une part
dy ε (t)
= εft (y ε (t))
dt
et d’autre part, pour ỹ ε (t) = Φ eε ◦ Ψ
e ε (y0 )
t t
dỹ ε (t)
= εft (ỹ ε (t)) − ε ∂y Φ e ε ◦ (Φ e ε )−1 (ỹ ε (t))
e ε )−1 (ỹ ε (t)) · Rε ◦ (Φ
t t t t
dt
e εt )−1 (ỹ ε(t)) restent dans KR . Si L = CK désigne une constante de Lipschitz de f
tant que ỹ ε et (Φ R
sur KR , une variante du lemme de Gronwall conduit à
3 eεLt − 1 ε0 | log(2)|
ky ε (t) − ỹ ε (t)kCn ≤ kRε kR ≤ 48R eεLt− 4 |ε| := M(t, ε)
2 L
e ε − idkR ≤ 1/2.
où l’on a utilisé l’inégalité k∂u Φ t
On rappelle que par hypothèse y ε (t) existe et appartient à K pour 0 < |ε| < ε∗ et 0 ≤ t ≤ T /ε.
En particulier, en choisissant |ε| < ε2 /2 avec ε2 tel que M(T /|ε|, ε) ≤ R − rnε , on est assuré que
e ε )−1 (ỹ ε(t)) restent dans KR pour t ≤ T /|ε| (donc T1 ≥ T ). On peut affirmer que pour ce
ỹ ε et (Φ t
choix ε2 P ε
0
0 < ε < min ε∗ , 0 , , (10.121)
2T 56
| log(2)| εLt− | log(2)|
on obtient l’inégalité (10.112). En effet, pour |ε| ≤ 4ηLT
= ε20 | log(2)| PT , on a e 2|ε|η ≤
− | log(2)|
e 4|ε|η sur [0, T /|ε|], de sorte que pour tout 0 ≤ t ≤ T /|ε|
| log(2)|
M(t, ε) ≤ 48Re− 4εη .
La quantité M(T /ε, ε) est alors bornée par R/2 ≤ R − rnε (notons que par définition de nε et ε2 ,
on a nε ≥ 1) si en outre
ε0 − log(2)
ε< .
8 log (96)
10. MÉTHODE DE MOYENNISATION À UN ORDRE QUELCONQUE : CAS PÉRIODIQUE93
L’agrégation des bornes sur ε (10.121) fournit une condition telle que
| log(2)|ε0
∀t ∈ [0, T /|ε|], ky ε(t) − ỹ ε (t)kCn ≤ 48 R e− 8ε .
fθ (y) ≡ Aθ y,
où Aθ est un élément de L(Rn ). Dans cette situation spécifique, il s’avère que le changement
de variable solution de (10.96) peut être construit exactement, en vertu du fait que la procédure
itérative mise en place précédemment est convergente. L’hypothèse 10.2 est ici formulée de la
manière suivante :
d ε
y (t) = εAt y ε (t), y ε (t) ∈ Rn , (10.122)
dt
y ε (0) = y0 , y 0 ∈ Rn .
Théorème 10.14 [Moyennisation √ exacte dans leR cas linéaire] On considère y ε (t) la solution de
(10.122) et on note εl = α1 ( 23 − 2), avec α = T kAθ kdθ. Alors, pour tout 0 ≤ ε < εl , il existe
une application θ 7→ Φεθ ∈ L(Rn ) avec les propriétés suivantes :
(i) La fonction Φεθ est P -periodique et C 1 en θ, et vérifie Φθ=0 = id.
(ii) Pour tout θ ∈ T, la matrice Φεθ est inversible.
(iii) Pour toute valeur initiale y0 ∈ Rn , la solution de (10.122) peut s’écrire sous la forme
ε
∀t ∈ R+ y ε (t) = Φεt eεt F y0 ,
Preuve du théorème 10.14. Dans le cadre linéaire, l’équation (10.96) revêt la forme suivante
dΦεθ
+ εΦεθ hΦε i−1 hAΦε i = εAθ Φεθ , ou de manière équivalente Φεθ = Γε (Φε )θ ,
dθ
où l’application Γε agit sur l’ensemble des functions de C(T, L(Rn )) inversibles pour tout θ, et est
définie, pour toute fonction Φ de la sorte, par
Z θ
ε
Γ (Φ)θ = Id + ε Aξ Φξ − Φξ hΦi−1 hAΦi dξ. (10.124)
0
[0] [n+1]
On introduit alors la suite de fonctions Φθ = Id, Φθ = Γε (Φ[n] )θ pour tout n ∈ N. On peut
affirmer que cette suite est bien définie pour tout n ∈ N et satisfait l’inégalité suivante
r
1 1
kΦ[n] − Idk < d∗ := − εα − ε2 α2 − 3εα + . (10.125)
2 4
Notons que le terme ε2 α2 − 3εα + 14 est positif en vertu de l’hypothèse ε < εl . Nous allons montrer
(10.125) par récurrence (l’initialisation est claire). Supposons donc que (10.125) soit satisfaite pour
un certain n, et posons dn := kΦ[n] − Idk. Alors, on déduit de 0 < ε < εl et de l’hypothèse de
récurrence
1 √
dn < d∗ < + εα ≤ 2 − 2 < 1.
2
[n] [n+1]
Il est alors clair (Neumann) que Φθ est inversible pour tout θ. En particulier, Φθ = Γε (Φ[n] )θ
est bien définie et il vient alors en utilisant (10.124) :
Z
ε [n] [n] kΦ[n] k
kΓ (Φ ) − Idk ≤ ε kΦ k 1+ kAθ kdθ,
1 − kΦ[n] − Idk T
Il est facile de vérifier grâce à l’hypothèse 0 < ε < εl , que si 0 ≤ d ≤ d∗ , alors 0 ≤ f (d) ≤
f (d∗ ) = d∗ où f (d) = 2εα 1+d 1−d
. Il résulte que dn+1 ≤ f (dn ) < d∗ et la majoration (10.125)
est ainsi prouvée. Afin de compléter la preuve du théorème, il reste à vérifier que la constante de
Lipschitz de Γε est strictement inférieure à 1 sur le domaine {Φ s.t. kΦ − Idk ≤ d∗ }. Pour ce faire,
on considère Φ ∈ C(T, L(Rn )) et Φ b ∈ C(T, L(Rn )) telles que kΦ − Idk < d∗ et kΦ b − Idk < d∗ :
pour tout θ ∈ T les applications Φθ et Φ b θ sont inversibles et leurs inverses peuvent être majorés
10. MÉTHODE DE MOYENNISATION À UN ORDRE QUELCONQUE : CAS PÉRIODIQUE95
Théorème 10.17 Etant donné l’équation (10.90), il existe un champ de vecteur (formel)
ε F ε (u) = ε F1 (u) + ε2 F2 (u) + ε3 F3 (u) + · · · (10.127)
tel que, pour toute valeur initiale y0 , la solution y (t) de (10.90) et la solution (formelle) z ε (t) du
ε
Exemple 10.18 On reprend ici l’exemple de l’équation de Van Der Pol considéré dans la première
section de ce cours, qu’on rappelle ci-après :
u̇ε = Jy ε + εg(uε)
avec
0 1 0
J= et g(u) = ,
−1 0 (1 − u21 )u2
système que l’on peut donc réécrire sous la forme
ε ε −tJ tJ ε tJ cos(t) sin(t)
ẏ = εft (y ) := εe g(e y ) avec e = ,
− sin(t) cos(t)
c’est-à-dire encore
ẏ ε = εξt (y ε )Vt
avec
− sin(t) 2
Vt = et ξt (y) = 1 − (cos(t)y1 + sin(t)y2 ) (− sin(t)y1 + cos(t)y2 ).
cos(t)
D’après ce qui précède, on peut calculer les deux premiers termes du développement de F ε de la
manière suivante :
Z 2π " #
1 −1/8 z1 (z2 2 + z1 2 − 4) (kzk22 − 4)
F1 (z) = ξτ (z)Vτ dτ = = − z
2π 0 −1/8 z2 (z2 2 + z1 2 − 4) 8
Z Z
−1 2π t
F2 (z) = (ξt (∇z ξτ , Vt ) Vτ − ξτ (∇z ξt , Vτ ) Vt ) dτ
4π 0 0
" 1
#
− 256 z2 (32 − 24 z2 2 + 5 z2 4 − 88 z1 2 + 21 z1 4 + 10 z1 2 z2 2 )
= 1
256 1
z (21 z1 4 + 32 − 88 z1 2 + 40 z2 2 + 10 z1 2 z2 2 + 5 z2 4 )
(kz ε k22 − 4) ε
ż ε = −ε z
8
et il est facile de voir par dérivation que N = kz ε k2 satisfait l’équation
N(N − 4)
Ṅ = −ε
4
dont un point fixe stable est N = 4. C’est l’approximation au premier ordre (cercle de rayon 2
et de centre (0, 0)) du cycle limite de l’équation de Van Der Pol. On peut obtenir une meilleure
98
Proposition 10.19 (Cas particulier de Magnus) Etant donné un champ de vecteur (y, t) 7→ ft (y),
il existe une unique paire de séries formelles
X [i]
X [i]
ε
Vs,t (y) = εi−1 Vs,t (y), Wtε (y) = εi−1 Wt (y), (10.132)
i≥1 i≥1
telles que
ε
V0,t (y) = 0, W0ε (y) = 0 et ε
V1,t (y) = ft (y), (10.133)
ε (εs)2 ε (εs)3 ε
ε Vs,· = (εs) Rε + R ⊳ Rε + R ⊳ (Rε ⊳ Rε ) + · · · (10.136)
2 3!
Il est immédiat de vérifier que pour y0 quelconque, z(1, t) = y ε (t) est l’unique solution de (10.90),
et que Y (t) = z(t/P, P ) est la solution du problème de valeur initiale
d ε
Y (t) = WPε (Y (t)), Y (0) = y0 ,
dt P
10. MÉTHODE DE MOYENNISATION À UN ORDRE QUELCONQUE : CAS PÉRIODIQUE99
1
et donc que F ε (Y ) := P
WPε (Y ) en vertu du fait que y ε (P ) = z(1, P ) = Y (P ).
Preuve de la proposition 10.19 : Supposons dans un premier temps que les séries (10.132)
vérifient les équations (10.133) et (10.134). En dérivant respectivement la première et la seconde
équation de (10.134) par rapport à respectivement s et t, on obtient
∂s Vs,t(y) − ∂t Wt (y) = ε ∂y Wt (y) Vs,t(y) − ε ∂y Vs,t (y) Wtε (y)
ε ε ε ε ε
(10.137)
évaluée au point (s, t, y) = (s, t, z(s, t)), où z(s, t) désigne l’unique solution de (10.134) corres-
pondant à la valeur initiale y0 . Or y0 est quelconque, donc l’équation (10.137) doit être satisfaite
pour des valeurs arbitraires de (s, t, y). En considérant maintenant
Rtε (y) = ∂t Wtε (y), (10.138)
la condition (10.137) se réécrit (en prenant en compte le fait que W0ε (y) ≡ 0) comme
hZ t i
ε ε ε ε ε
∂s Vs,t = Rt + ε [Wt , Vs,t] = Rt + ε ∂τ Wτε dτ, Vs,t
ε
0
hZ t i
= Rtε + ε Rτε dτ, Vs,t
ε
= Rtε + (R·ε ⊳ Vs,·
ε
)t , (10.139)
0
ce qui conduit à (10.136). Il résulte que les équations (10.132), (10.133) et (10.134) impliquent
(10.136) avec (10.135).
Réciproquement, supposons qu’un couple de fonctions (V ε , W ε ) de la forme (10.132) satis-
fasse (10.133), (10.135) et (10.136). Alors les functions R[i] sont déterminées de manière unique
par l’équation (10.136) pour s = 1, de telle sorte que les W [i] sont également déterminées de
manière unique par (10.135). L’équation (10.137) est alors satisfaite. Il reste Pdonc à montrer que le
couple de fonctions (V ε , W ε ) vérifie (10.134). Soit en effet z(s, t) = y0 + i≥1 εi zi (s, t) l’unique
ε
solution de la première équation de (10.134) : puisque ∂t z(0, t) = ε V0,t (y) = 0, il résulte que
z(0, t) = y0 , et il suffit de prouver par récurrence sur i que, pour tout i ≥ 0,
∂s z(s, t) = ε Wtε (z(s, t)) + O(εi+1 ).
L’initialisation est claire. Supposons donc la relation vérifiée pour un certain i ≥ 0. En dérivant
par rapport à s la première équation de (10.134) et en utilisant l’hypothèse de récurrence, il vient
2 ε ε
∂t,s z(s, t) = ε(∂s Vs,t )(z(s, t)) + ε ∂y Vs,t (z(s, t)) ∂s z(s, t)
ε
= ε(∂s Vs,t )(z(s, t)) + ε2 ∂y Vs,t
ε
(z(s, t)) Wtε (z(s, t)) + O(εi+2 ),
ε 2 ε ε
= ε(∂t Wt )(z(s, t)) + ε ∂y Wt (z(s, t)) Vs,t (z(s, t)) + O(εi+2 ),
et finalement, puisque ∂s z(s, 0) ≡ 0,
Z t
∂s z(s, t) = ε (∂τ Wτε )(z(s, τ )) + ε∂y Wτε (z(s, τ )) Vs,τ
ε
(z(s, τ )) dτ + O(εi+2 )
0
= εWtε (z(s, t)) − εW0ε (z(s, 0)) + O(εi+2 )
= εWtε (z(s, t)) + O(εi+2 ).
100
Qt (y ε (t)) ≡ Q0 (y0 ).
La question que l’on se pose ici est de savoir si le champ de vector moyenné F [n] possède lui aussi
un invariant. Il s’avère que la réponse à cette question est positive et que la preuve de ce résultat
est fort simple. Le point clé est que Qθ ◦ Φθ est (presque) indépendant de θ, tandis que Q0 est
presque invariant pour le système moyenné (invariant à des termes d’erreur petits en ε). Avant de
poursuivre, il convient de définir précisément la notion d’invariant utilisée ici. En différentiant la
relation Qt (y ε (t)) ≡ Q0 (y0 ), il vient
pour (y, θ) = (y ε (t), t). Dans la suite ce paragraphe, on suppose que cette relation est en fait
satisfaite pour tout θ ∈ T et tout y ∈ K.
[n]
En particulier, on a (d/dt)Q0(Ψt (y0 )) = O(εn+1), pour t ∈ [0, T /ε].
Remarque 10.21 Si l’invariant Qθ de dépend pas de ε, alors, puisque F [n+1] − Fe[n] = O(εn+1)
et en remarquant que Fe[n] est un polynôme de degré n en ε, il est possible de conclure de (10.141)
que
(∂y Q0 )(y) Fe[n] (y) = 0,
c’est-à-dire que Q0 est exactement préservé par le système autonome (10.111).
[n]
conduit, en pre-multipliant par (∂y Qθ ) ◦ Φθ et en utilisant le fait que Q est un invariant de εfθ , à
[n] [n] [n]
(∂y Qθ ) ◦ Φθ (y) · ∂θ Φθ (y) + ε∂y (Qθ ◦ Φθ )(y) F [n](y)
[n] [n]
= ε ∂y Qθ ◦ Φθ (y) fθ ◦ Φθ (y) + O(εn+1 )
[n]
= − (∂θ Qθ ) ◦ Φθ (y) + O(εn+1 ),
pour y ∈ K. Notons que le terme d’erreur O(εn+1) est à entendre au sens des fonctions analytiques,
c’est-à-dire par exemple pour toute norme k.kρ′ avec 0 < ρ′ < ρ. En conséquence, il vient
[n] [n]
∂θ Qθ ◦ Φθ (y) + ε∂y Qθ ◦ Φθ (y) F [n] (y) = O(εn+1 ). (10.142)
Le résultat énoncé dans ce théorème découle alors d’un argument de récurrence. Supposons que
pour un certain k < n, et pour tout θ ∈ T, on ait
[n]
Qθ ◦ Φθ (y) = Q0 (y) + O(εk+1 ).
Notons que cette propriété est claire pour k = 0 car Φ[n] = Id + O(ε) et Qθ (y) = Q0 (y) + O(ε)
par intégration de (10.140). Il vient alors
[n]
Qθ ◦ Φθ = Q0 + O(εk+1) = hQ ◦ Φ[n] i + O(εk+1),
et par intégration en θ, à
[n] [n]
Qθ ◦ Φθ = Q0 ◦ Φ0 + O(εk+2 ) = Q0 + O(εk+2),
Définition 10.22 Le champ de vecteur (y, θ) 7→ fθ (y) dans l’hypothèse 10.2 est dit hamiltonien
s’il existe une fonction (y, θ) 7→ Hθ (y) analytique au sens de la définition 10.1, telle que
Nous allons établir dans cette partie que le champ de vecteur F [n] défini dans le théorème 10.6 est
lui-aussi hamiltonien, en tout cas à un terme d’erreur εn+1 près. Avant cela, le lemme suivant est
requis :
[n]
Lemme 10.24 Sous les hypothèses du théorème 10.6, on suppose que Φθ est symplectique à
εk+1-près, avec 0 ≤ k ≤ n, c’est-à-dire que
[n] [n]
(Sk ) ∀y ∈ K, (∂y Φθ (y))T J(∂y Φθ (y)) = J + O(εk+1).
avec
[n]
Preuve : L’équation (10.96) avec les notations Φθ ≡ Φθ et F ≡ F [n] peut être mise sous la forme
et d’autre part
[n]
Théorème 10.25 Sous les hypothèses du théorème 10.6, Φθ est symplectique à εn+1-près et F [n]
est hamiltonien à εn+1 -près, c’est-à-dire qu’on a :
[n] [n]
∀y ∈ K, (∂y Φθ (y))T J(∂y Φθ (y)) = J + O(εn+1 ),
F [n] (y) = J −1 ∇y H [n] (y) + O(εn+1 ),
1 D [n] E
H [n] (y) = hH ◦ Φ[n] (y)i − (Φ (y))T J(∂θ Φ[n] (y)) .
2ε
-
104
[n]
Preuve : Par construction, Φθ (y) = y+O(ε), de sorte que (S0 ) et (H0 ) sont satisfaites. Supposons
maintenant que (Sk ) soit vérifiée pour un certain 0 ≤ k ≤ n − 1. Alors, en vertu du lemme 10.24,
(Hk ) est également vérifiée. L’équation (10.96) peut alors être réécrite, pour θ petit, comme
d [n]
−1
[n]
Φθ ◦ Ψθ = εJ (∇y Hθ ) ◦ Φθ ◦ Ψθ + O(εk+2 )
dθ
et donc T
[n] [n]
∂y Φθ ◦ Ψθ J ∂y Φθ ◦ Ψθ = J + O(εk+2).
d
La relation Ψ
dτ τ
= εJ −1 (∇y H [n] ) ◦ Ψτ + O(εk+2) permet finalement d’écrire
T
∂y Ψθ J ∂y Ψθ = J + O(εk+2),
et donc d’affirmer que (Sk+1 ) est vérifiée. Une récurrence permet de conclure.
ẏ ε = εftω (y ε ), (11.147)
y ε (0) = y0 ∈ Rn , (11.148)
où ε est un petit paramètre, la fonction (y, θ) ∈ Rn × Td 7→ fθ (y) est régulière (on précisera plus
avant le degré de régularité nécessaire), 2π-périodique par rapport à chaque composante 7 θj de
θ ∈ Rd , et où ω ∈ Rd est un vecteur constant de fréquences angulaires. On fait en outre l’hypothèse
fondamentale à travers toute cette partie, que ω est non-résonant, c’est-à-dire que, pour tout multi-
indice k ∈ Zd , k · ω 6= 0. Notons que lorsque le problème est résonant, il peut toujours se réécrire
sous la forme d’un problème non-résonant en diminuant le nombre de fréquences. De même que
dans le cas périodique, il s’agit ici d’étudier la solution de (11.147)–(11.148) sur un intervalle de
la forme 0 ≤ t ≤ L/ǫ de sorte que la solution y ε varie d’une quantité O(1). Le cas d = 1 est le cas
périodique que nous avons déjà étudié. Pour d > 1, le problème est beaucoup plus difficile à traiter,
en raison du problème dit des “petits dénominateurs” (voir la section concernant les estimations
d’erreur).
i∂t ψ ε (t, x) = −∆ψ ε (t, x) + ε|ψ ε (t, x)|2 ψ ε (t, x), x ∈ [0, 2π] × [0, 2πa], t ∈ [0, T /ε[,
7. Les exposants correspondent ici aux composantes des vecteurs considérés.
11. MÉTHODE DE MOYENNISATION À UN ORDRE QUELCONQUE : CAS QUASI-PÉRIODIQUE105
I = {k ∈ Z2 , −N ≤ k1 , k2 ≤ N}
2
et en effectuant le changement de variables ψ̂lε (t) = e−i|l| t ξlε l(t), on obtient un système de la forme
(11.147) :
X 2 2 2 2
∀ l ∈ I, ξ˙lε (t) = iε ei(|k1 | −|k2 | +|k3 | −|l| )t ξkε 1 (t) ξkε 2 (t) ξkε 3 (t)
k1 , k2 , k3 ∈ I
k1 − k2 + k3 = l
Dans un premier temps, et afin de travailler de manière formelle, nous supposons que les fonctions
fk sont de classe C ∞ et nulles sauf pour un nombre fini d’indices k ∈ Zd . Nous verrons par la suite
qu’il est possible de travailler sous des hypothèses beaucoup moins restrictives. A chaque fonction
fk , k ∈ Zd , on peut associer l’opérateur linéaire suivant :
∂·
Ek := fk
∂y
(les exposants correspondent ici aux composantes des vecteurs considérés). L’opérateur Ek agit
sur l’espace C ∞ (Rn ; Rn ) de la manière suivante : pour g ∈ C ∞ (Rn ; Rn ), Ek [g] est défini comme
la fonction de C ∞ (Rn ; Rn ) telle que
n
!
X ∂
∀y ∈ Rn , Ek [g](y) = fkj (y) j g i (y) = g ′ (y)fk (y).
j=1
∂y
i=1...n
106
On considère maintenant l’opérateur Φt agissant sur les fonctions g ∈ C ∞ (Rd ; Rd ) tel que
pour toute solution y ε de (11.147) et toute fonction g. Clairement Φ0 = I, où I désigne l’opérateur
identité, c’est-à-dire l’opérateur tel que I[g] = g pour tout g ∈ C ∞ (Rd ; Rd ). On a alors :
d X
g(y ε(t)) = g ′ (y(t))ftω (y ε (t)) = ε eik·ωt Ek [g](y ε (t)),
dt d k∈Z
ou de manière équivalente
d X
Φt [g](y(0)) = ε eik·ωt Φt Ek [g](y(0)),
dt d k∈Z
car Ek [g](y ε (t)) = (g ′ fk )(y ε (t)) = Φt [g ′fk ](y(0)) = Φt Ek [g](y(0)). Cette dernière équation
montre que l’opérateur Φt est solution de l’équation différentielle linéaire
d X
Φt = ε eik·ωt Φt Ek , Φ0 = I,
dt d k∈Z
que nous pouvons résoudre par itération de Picard. Partant de Φ[0] ≡ I, on construit la suite
Z t X
[l+1]
Φt =I +ε eikl+1 ·ωtl+1 Φtl+1 Ekl+1 dtl+1 , l = 0, . . . , +∞.
0
kl+1 ∈Zd
[1]
X Z t
ik1 ·ωt1
Φt = I + ε e dt1 Ek1 ,
k1 0
[2]
XZ t
[1]
Φt = I +ε eik2 ·ωt2 Φt2 Ek2 dt2
k2 0
!
XZ t X Z t2
ik2 ·ωt2 ik1 ·ωt1
= I +ε e I +ε e dt1 Ek1 Ek2 dt2
k2 0 k1 0
X Z t X Z t Z t2
ik2 ·ωt2 2 ik2 ·ωt2 ik1 ·ωt1
= I +ε e dt2 Ek2 + ε e e dt1 dt2 Ek1 Ek2 dt2 ,
k2 0 k1 ,k2 0 0
Par définition de l’opérateur Φt , pour toute fonction g ∈ C ∞ (Rd ; Rd ) et toute solution y ε (t) de
(11.147), on a la relation :
∞
X X
ε ε
g(y (t)) = g(y (0)) + εj αk1 ···kj (t) Ek1 · · · Ekj [g](y ε(0)),
j=1 k1 ,...,kj ∈Zd
relation qui, dans le contexte de la théorie du contrôle non-linéaire, est connue sous le nom de
développement de Chen-Fliess. En prenant pour g la fonction y 7→ y, la solution y ε (t) de (11.147)–
(11.148) peut alors être représentée de manière formelle par la formule
∞
X X
ε
y (t) = y0 + εj αk1 ···kj (t) fk1 ···kj (y0 ), (11.151)
j=1 k1 ,...,kj ∈Zd
avec
fk1 ···kj (y) := fk′ 2 ···kj (y)fk1 (y). (11.152)
Il est en effet facile de vérifier par récurrence que
Ek1 · · · Ekj [id] = fk1 ···kj .
Le coefficient X
αk1 ···kj (t) fk1 ···kj (y0 ) (11.153)
k1 ,...,kj ∈Zd
de εi dans (11.151) est constitué de deux ingrédients : les fonctions mots 8 fk1 ···kj (y) définies par
l’équation (11.152) et les coefficients scalaires αk1 ···kj dans (11.150). Il est important de noter que
ces coefficients dépendent du vecteur de fréquences ω mais sont complètement indépendants de la
fonction (y, θ) 7→ fθ (y) dans (11.149). D’un autre coté, les fonctions mots dépendent uniquement
des coefficients de Fourier fk ; ce sont des combinaisons de leurs dérivées partielles. On a par
exemple :
flk (y) = fk′ (y) fl(y), (11.154)
′
fmlk (y) = flk (y) fm(y)
= fk (y) (fl (y), fm (y)) + fk′ (y) fl′(y) fm (y).
′′
(11.155)
8. La terminologie “fonction mots” sera motivée ci-après.
108
Ici, fk′ (y) désigne la dérivée première de fk au sens de Fréchet, évaluée au point y, et fk′ (y) fl(y)
est le produit d’une matrice de Rn × Rn par un vecteur de Rn . De même, fk′′ (y) désigne la dérivée
seconde (toujours au sens de Fréchet) de fk au point y et fk′′ (y) (fl (y), fm (y)) représente son
action sur les deux vecteurs de Rn que sont fl (y) et fm (y) ; ainsi, la r-ième composante de
fk′′ (y) (fl (y), fm (y)) s’écrit
X n
∂2
f (y) fli (y)fm
r
i ∂y j k
j
(y).
i,j=1
∂y
Les fonctions fk′ fl , fk′′ (fl , fm ), . . ., apparaissant dans (11.154)–(11.155) sont appelées éléments
différentiels relatifs à la fonction f dans (11.147). L’élément différentiel fk′ fl est d’ordre 2 (il
contient deux facteurs) et les éléments différentiels fk′′ (fl , fm ) et fk′ fl′ fm sont tous les deux
d’ordre 3. D’une manière générale, chaque fonction mot fk1 ···kl est une combinaison linéaire d’éléments
différentiels d’ordre l. Notons qu’un élément différentiel peut entrer dans la définition de plusieurs
fonctions mots. Par exemple, fk′′ (fl , fm ), l 6= m est un terme à la fois de fmlk et de flmk .
Les notations dans (11.153) et dans les expressions à suivre peuvent être simplifiées en considérant
la notion de mots k1 · · · kl , constitués des lettres kr , r = 1, . . . , l, prises parmi l’alphabet Zd . Si les
Wl , l = 1, 2, . . . , désignent les ensembles de mots à l lettres, alors la série (11.153) peut s’écrire :
X
αw (t) fw (y0). (11.156)
w∈Wl
où δ est une application qui, à chaque mot w ∈ W, associe un nombre complexe δw . Les séries de
la forme (11.157) sont appelées séries indexées par des mots (relatives à l’équation (11.147)). Le
cas particulier où δ est l’application telle que δ∅ = 1 et δw = 0 pour tout w 6= ∅ correspond à la
série identité W (δ, y) ≡ y.
est encore une série du même type S(δ ⋆ δ̂) = S(δ)S(δ̂) où leproduit des coefficients β = δ ⋆ δ̂ est
défini par
X
∀w ∈ W, (δ ⋆ δ̂)w = δu δ̂v = δ∅ δ̂w + · · · + δw δ̂∅ (11.158)
(u, v) ∈ W 2
uv = w
L’élément neutre de cette multiplication est la série S(1) = I associée aux coefficients (1 est
l’élément neutre de ⋆)
1 si w=∅
1w =
0 si w ∈ W\{∅}
Etant donnée une application δ : W → C satisfaisant δ∅ = 1, on peut également définir l’inverse
δ −1 de δ comme étant l’unique application de W dans C telle que S(δ −1 )S(δ) = S(δ)S(δ −1) = I.
Il suffit pour cela de contruire δ −1 par récurrence sur le nombre de lettres à partir de la relation
δ −1 ⋆ δ = δ ⋆ δ −1 = 1. On a en effet, par définition (11.158), que δ −1 ⋆ δ = 1 si et seulement si
δ∅−1 δ∅ = 1,
X
∀l ≥ 1, ∀w ∈ Wl , δu−1 δv = 0.
(u, v) ∈ W 2
uv = w
La première relation impose δ∅−1 = 1. Supposons que la seconde soit vérifiée pour l ≥ 1 donné et
considérons un mot quelconque w ∈ Wl+1 . Alors nécessairement
X
δw−1 = − δu−1 δv
(u, v) ∈ W 2 , |u| ≤ l,
uv = w
où les mots u apparaissant dans le membre de droite comporte au plus l lettres. Ce qui définit δw−1
de manière unique. Par symétrie de la relation (11.158), l’inverse à droite est également inverse à
gauche.
Il est intéressant de noter que la série Φt = S(α(t)) est la solution de l’équation différentielle
d
S(α(t)) = S(α(t))S(β(t))
dt
où
X
itk·ω eitk·ω si w = k ∈ W1 = Zd ,
S(β(t)) = ε e Ew avec βw (t) =
0 sinon.
k∈W1
avec la condition initiale S(α(0)) = I. Traduite en termes de relations sur les coefficients, cette
équation devient
d
α(t) = α(t) ⋆ β(t), α(0) = 1. (11.159)
dt
110
d X
αk1 ...kl (t) = αu (t)βv (t) = αk1 ...kl−1 (t) βkl (t) = αk1 ...kl−1 (t) eitkl ·ω
dt uv=k ...k
1 l
ε3 ′′
εfk (y + fk (y)) = εfk (y) + ε2 fk′ (y)fk (y) + f (y) (fk (y), fk (y)) + O(ε4 )
2 k
il n’y a pas d’identification possible avec les premiers termes d’une série indexée par des mots,
puisque comme nous l’avons détaillé ci-avant, on a
En revanche, si
y(t) = W (δ(t), y(0)) = S(δ(t))[id](y(0))
est la solution d’une équation différentielle associée à un champ de vecteur
c’est-à-dire si X
ẏ(t) = ε|w|β̃w (t)Ew [id](y(t))
w∈W
11. MÉTHODE DE MOYENNISATION À UN ORDRE QUELCONQUE : CAS QUASI-PÉRIODIQUE111
alors S(α(t)) est l’opérateur solution de l’équation différentielle (le processus est le même que
pour Φt )
d
S(δ(t)) = S(δ(t))S(β̃(t)), S(δ(0)) = I,
dt
ce qui équivaut, par définition du produit S(δ(t))S(β̃(t)) de séries d’opérateurs, à l’équation
différentielle
δ̇w (t) = δ(t) ⋆ β̃(t), δ(0) = 1.
Ainsi, pour toute fonction g ∈ C ∞ (Rn ; Rn ), on a
de sorte que si la fonction g est elle-même une série indexée par des mots de la forme g(y) =
W (ξ, y), alors la composition (11.160) peut s’écrire
W (ξ, W (δ(t), y(0))) = g W (δ(t), y(0)) = S(δ(t))[g](y(0))
h i i
= S(δ(t)) S(ξ)[id] (y(0)) = S(δ(t))S(ξ) [id] (y(0))
= S(δ(t) ⋆ ξ)[id](y(0)) = W (δ(t) ⋆ ξ, y(0)).
Cette relation étant valable pour tout y(0), on a donc finalement le résultat suivant :
Lemme 11.2 Soient δ et ξ deux applications de W dans C et W (δ, y), W (ξ, y) les séries indexées
par des mots associées. Si il existe une application β̃(t) de W dans C, continue en t et vérifiant
β̃∅ ≡ 0, telle que δ coı̈ncide avec la valeur en un instant t de la solution de l’équation différentielle
˙
α̃(t) = α̃(t) ⋆ β̃(t), α̃(0) = 1,
alors la composition W (ξ, W (δ, y)) est encore une série indexée par des mots et on a
On peut ainsi se convaincre aisément que les coefficients α sont des polynômes en les 2d + 1
variables t, eiω1 t , . . ., eiωd t , e−iω1 t , . . ., e−iωd t . On montre alors par récurrence que pour tout mot w,
il existe une unique fonction γw (t, θ), t ∈ R, θ ∈ Td , polynômiale en les variables t, eiθ1 , . . ., eiθd ,
e−iθ1 , . . ., e−iθd telle que
αw (t) = γw (t, tω). (11.161)
On a en particulier γ∅ (t, θ) ≡ 1.
Définition 11.3 On dit qu’une fonction µ de R × Td à valeurs dans C est polynomiale s’il existe
un polynôme P ∈ C[X0 , X1 , . . . , X2d ] à 2d + 1 variables X0 , . . . , X2d , tel que
1 d 1 d
∀(t, θ) ∈ R × Td , µ(t, θ) = P (t, eiθ , . . . , eiθ , e−iθ , . . . , e−iθ )
Lemme 11.5 Etant donnée une fonction polynomiale µ : R×Td → C, il existe une unique solution
polynomiale de l’équation
où, par hypothèse (µ polynomiale), I est un sous-ensemble fini de Zd et les µ̂k (t) et ν̂k (t) sont des
polynômes en t. L’équation (11.162) s’écrit alors
d µ̂k (t) si k ∈ I
ν̂k (t) + i(k · ω)ν̂k (t) = (11.163)
dt 0 sinon
Proposition 11.6 Pour tout w ∈ W et tout k ∈ Zd , γwk est l’unique fonction polynomiale telle
avec γ∅ ≡ 1.
11. MÉTHODE DE MOYENNISATION À UN ORDRE QUELCONQUE : CAS QUASI-PÉRIODIQUE113
Preuve : L’unicité d’une solution polynomiale est une conséquence immédiate du lemme précédent
et d’une récurrence sur le nombre de lettres. Considérons maintenant la fonction γ définie par
(11.161) : pour θ = tω, on a γwk (t, θ) = αwk (t) et γw (t, θ) = αw (t), de sorte que
(∂t + ω · ∇θ )γwk (t, tω) = α̇wk (t)
et d’après (11.159)
α̇wk (t) = αw (t)eitk·ω = γw (t, tω)eitk·.ω
Par ailleurs, γwk (0, 0) = αwk (0) = 0. La fonction γ définie par (11.150) satisfait donc (11.164)
pour tout (t, tω) ∈ R × Td : elle coı̈ncide donc avec l’unique solution de (11.164) pour tout
(t, tω) ∈ R × Td , et donc, en vertu du lemme 11.4, pour tout (t, θ) ∈ R × Td .
L’équation de transport (11.164) joue ici le rôle de l’équation dite homologique dans l’ap-
proche usuelle de la moyennisation (c’est-à-dire l’approche basée sur une construction récursive
de changements de variables). C’est sur cette équation que repose les résultats suivants.
Proposition 11.7 Pour tous t, t′ ∈ R et tout y ∈ Rn
W γ(t′ , 0), W (γ(t, 0), y) = W (γ(t + t′ , 0), y), (11.165)
et pour tout t ∈ R , tout θ ∈ Td et tout y ∈ Rn
W γ(0, θ), W (γ(t, 0), y) = W (γ(t, θ), y). (11.166)
Preuve : Soient t′ ∈ R fixé et γ̃(t, θ) l’application de R × Td dans C définie par
γ̃(t, θ) = γ −1 (t′ , 0) ⋆ γ(t + t′ , θ).
Pour w ∈ W et k ∈ Zd , on a :
X
(∂t + ω · ∇θ )γ̃wk (t, θ) = γu−1 (t′ , 0) ′
(∂t + ω · ∇θ )γv (t + t , θ)
(u, v) ∈ W 2
uv = wk
Les couples de mots tels que uv = wk sont, soit de la forme (u, v ′k) avec uv ′ = w, soit de la
forme (wk, ∅). En tenant compte du fait que γ∅ (t, θ) ≡ 1, on a (∂t + ω · ∇θ )γ∅ (t, θ) = 0 et il vient
donc
X
−1 ′ ′
(∂t + ω · ∇θ )γ̃wk (t, θ) = γu (t , 0) (∂t + ω · ∇θ )γvk (t + t , θ)
(u, v) ∈ W 2
uv = w
X
= γu−1 (t′ , 0) γv (t + t′ , θ)eik·θ
(u, v) ∈ W 2
uv = w
X
= eik·θ γu−1 (t′ , 0) γv (t + t′ , θ)
(u, v) ∈ W 2
uv = w
où l’on a utilisé le fait que γ(t+t′ , θ) satisfait l’équation de transport (11.164) évaluée en (t+t′ , θ).
On a en outre γ̃(0, 0) = γ −1 (t′ , 0) ⋆ γ(t′ , 0)) = 1. Par unicité de la solution de l’équation de
transport (11.164), on peut conclure que γ̃(t, θ) = γ(t, θ), soit encore
y ∈ Rn 7→ W (ᾱ(t), y) ∈ Rn .
avec X
Fl (z) := β̄w fw (z) l = 1, 2, . . . (11.171)
w∈Wl
En particulier, il est facile de voir que ᾱk (t) = 0 pour k 6= 0 et que ᾱ0 = t, de telle sorte que, par
définition de β̄, β̄k (t) = 0 pour k 6= 0 et β̄0 = 1. Il vient donc
Z
1
F1 (z) = f0 (z) = fθ (z)dθ, (11.172)
(2π)d Td
où X
[l]
uθ (z) = κw (θ) fw (z), l = 1, 2, . . . (11.175)
w∈Wl
ε
et où z (t) désigne la solution de l’équation autonome (moyennée) (11.170) avec condition initiale
z ε (0) = y0 .
116
i
γk (t, θ) = (1 − eik·θ ),
k·ω
tl
γ0l (t, θ) = ,
l!
i
γ0l k (t, θ) = γ0l−1 k (t, θ) − γ0l (t, θ)eik·θ ,
k·ω
i
γkl1 ···ls (t, θ) = γl1 ···ls (t, θ) − γ(k+l1 )l2 ···ls (t, θ) ,
k·ω
i
γ0l kl1 ···ls (t, θ) = γ0l−1 kl1 ···ls (t, θ) − γ0l (k+l1 )l2 ···ls (t, θ) .
k·ω
β̄k = 0,
β̄0 = 1,
β̄0l+1 = 0,
i
β̄0l k = (β̄ l−1 − β̄0l ),
k·ω 0 k
i
β̄kl1 ···ls = (β̄l ···l − β̄(k+l1 )l2 ···ls ),
k·ω 1 s
i
β̄0l kl1 ···ls = (β̄ l−1 − β̄0l (k+l1 )l2 ···ls ).
k · ω 0 kl1 ···ls
La valeur des coefficients κw dans (11.173), nécessaires pour exprimer le changement de va-
riables (11.174)–(11.175), peuvent également être calculés à partir de la proposition 11.9.
Remarque 11.11 Dans l’approche de la moyennisation retenue ici, le champ moyenné et le chan-
gement de variables peuvent être déterminés (et calculés) de manière indépendante.
11. MÉTHODE DE MOYENNISATION À UN ORDRE QUELCONQUE : CAS QUASI-PÉRIODIQUE117
Remarque 11.13 L’hypothèse faite sur f , à savoir que son développement en série de Fourier est
fini, a été utilisée afin d’assurer que les séries (11.153), (11.171), (11.175) convergent, de sorte que
(11.151), (11.170), (11.174) puissent être considérées comme des séries formelles en puissances
de ε. En tronquant ces dernières séries, on obtiendrait bien sûr des développements de Taylor poly-
nomiaux, et les hypothèses de régularité de f pourraient être affaiblies. Il est cependant aussi pos-
sible d’effectuer la même construction sans cette hypothèse de finitude du développement de Fou-
rier, tout en assurant la convergence des séries (11.153), (11.171), (11.175). C’est le choix retenu
dans le paragraphe qui détaille les estimations d’erreur (f est alors supposée réelle-analytique en
y mais aussi en θ).
avec F1 = f0 , et
X i
F2 = ([f−k , fk ] + [fk − f−k , f0 ])
k>0
k · ω
X 1 1
F3+ = [f 0 , [f ,
0 kf ]] + [fk , [f , f
k −k ]] + [f k , [f ,
0 kf ]] + 2[f−k , [f ,
k 0f ]]
k>0
(k · ω)2 2
X 1
+ [f−l , [fk , fl−k ]] − [fl , [fk , f−k−l ]] ,
k>l>0
(k · ω)(l · ω)
où F3− est obtenu en remplaçant dans l’expression de F3+ , (k, l) par (−k, −l). Ici, la notation >
désigne une relation d’ordre total sur Zd compatible avec l’addition.
Théorème 11.15 (i) Le système autonome (11.170)–(11.171) obtenu par moyennisation quasi-
stroboscopique peut s’écrire sous la forme géométrique suivante
∞
d ε X X εl
z = β̄k1 ···kl [[· · · [[fk1 , fk2 ], fk3 ] · · · ], fkl ](z ε ). (11.177)
dt l=1
l
k1 ,...,kl ∈Z d
(ii) Supposons que le champ de vecteur f dans l’équation hautement oscillante d’origine soit
à divergence nulle , alors le champ de vecteur moyenné (11.177) est également à divergence nulle.
(iii) Supposons que le champ de vecteur f dans l’équation hautement oscillante d’origine soit
hamiltonien dans une espace de dimension paire (n = 2D), associé à l’hamiltonien Hθ (y)
X
Hθ (y) = eik·θ Hk (y)
k∈Zd
(de sorte que pour tout k ∈ Zd , fk (y) soit le champ de vecteur hamiltonien associé à Hk ), alors le
champ de vecteur moyenné (11.177) est également hamiltonien et est associé à l’hamiltonien 9
∞
X X εr
H̄ = β̄k ···k {{· · · {{Hk1 , Hk2 }, Hk3 } · · · }, Hkr }. (11.178)
r=1 k1 ,...,kr ∈Z d
r 1 r
Preuve : Le point (i) repose, comme mentionné en préambule de ce paragraphe, sur le théorème
de Dynkin-Specht-Wever. Le point (ii) découle du calcul suivant : pour g et h deux champs de
9. Le crochet de Poisson de deux hamiltoniens H et K est défini par la formule
vecteur, on a :
∂ X ∂hi j ∂g i j
Xn n
div[g, h] = g − jh
i=1
∂y i j=1 ∂y j ∂y
n
X ∂ 2 hi j ∂ 2 g i j ∂hi ∂g j ∂g i ∂hj
= g − h + −
i,j=1
∂y j ∂y i ∂y j ∂y i ∂y j ∂y i ∂y j ∂y i
n
X ∂(divh) ∂(divg) j
= gj − h
j=1
∂y j ∂y j
de sorte que si divh = divg = 0, alors div[g, h] = 0. Le point (iii) s’obtient grâce à un calcul
similaire qui montre que le crochet de Lie de deux champs hamiltoniens est encore hamiltonien.
En effet, si g = J −1 ∇y K et h = J −1 ∇y L, alors
Kρ = {y + z ∈ Cn : y ∈ K, kzk ≤ ρ}.
Remarque 11.17 Sous l’hypothèse A, la série de Fourier (11.149) converge absolument et uni-
formément pour y ∈ KR , θ ∈ Tdµ , où l’ensemble Tdµ est défini par
Ainsi, f peut être étendue à une fonction analytique sur KR × Tdµ . Remarquons que
∀θ ∈ Td , kf (·, θ)kR ≤ M.
120
La décroissance exponentielle des coefficients de Fourier , telle que supposée dans (11.179), per-
met de contrer l’influence des petits dénominateurs. Dans le cas d = 1, il est possible de prendre
µ = 0 dans l’hypothèse A : f n’a alors pas besoin d’être supposée analytique en θ. Cependant,
afin de traiter ces petits dénominateurs, une hypothèse supplémentaire est nécessaire : on suppose
donc que le vecteur ω ∈ Rd , d > 1, satisfait une condition de non-résonance forte :
y = z + εǓtε,N (z)
où
[1] [2] [N ]
Ǔθε,N (z) := uθ (z) + εuθ (z) + · · · + ǫN −1 uθ (z)
(les fonctions ul sont définies comme dans (11.174)). Si ε ∈ C est tel que
R
|ε| ≤ ε0 , ε0 = ε0 (N) := , (11.181)
4LN ν+1
alors
1. Pour tout θ ∈ T, l’application z ∈ KR/2 7→ z + εǓθε,N (z) est analytique à valeurs KR .
2. Pour tous θ ∈ T, k(ω · ∂θ )Ǔ ε,N − θ(·)kR/2 ≤ 3M/2.
ε,N )
3. Pour tout θ ∈ T and z ∈ KR/2 , la matrice jacobienne I + ε∂z žθ (z) est inversible et
Le théorème savant établit les estimations d’erreur recherchées : notins que le terme d’erreur en
exp(−c̃|ε|−1/(ν+1) ) dans (11.182) s’accroı̂t avec ν, c’est-à-dire que la qualité de l’approximation
par moyennisation se détériore à mesure que l’on se rapproche d’une valeur résonante de ω. Dans
le cas périodique, la borne sur Rε,N dans (11.182) était simplement de la forme C exp(−c̃/|ε|).
11. MÉTHODE DE MOYENNISATION À UN ORDRE QUELCONQUE : CAS QUASI-PÉRIODIQUE121
5|ε/ε0|N
kRε,N kR/2 ≤ M.
1 − |ε/ε0|
ν
En particular, si |ε| ≤ R/(4eL) (L est donné par la relation L = 2M ν
cµν eν
), et si N est choisi comme
1/(ν+1)
la partie entière de R/(4eL|ε|) ≥ 1, alors l’estimation exponentielle suivante prévaut :
1/(ν+1)
ε,N 5e2 c̃ R
kR kR/2 ≤ M exp − 1/(ν+1) , c̃ = . (11.182)
e−1 |ε| 4eL
122
Bibliographie
[1] V. I. Arnold, Geometrical Methods in the Theory of Ordinary Differential Equations, 2nd ed.,
Springer, New York, 1988.
[2] V. I. Arnold, Mathematical Methods of Classical Mechanics, 2nd ed., Springer, New York,
1989.
[3] M.P. Calvo, Ph. Chartier, A. Murua, J.M. Sanz-Serna, A stroboscopic method for highly
oscillatory problems. In : B. Engquist, O. Runborg, R. Tsai, (eds.) Numerical Analysis of
Multiscale Computations, Lect. Notes Comput. Sci. Eng. 82, pp. 73–87. Springer, Berlin,
2011.
[4] M.P. Calvo, Ph. Chartier, A. Murua, J.M. Sanz-Serna, Numerical stroboscopic averaging
for ODEs and DAEs, Applied Numerical Mathematics, Appl. Numer. Math. 61, 1077–1095
(2011).
[5] F. Castella, P. Chartier, F. Méhats and A. Murua, Stroboscopic averaging for the nonlinear
Schrödinger equation, hal-00732850, http ://hal.archives-ouvertes.fr/hal-00732850, FOCM
(2014).
[6] P. Chartier, A. Murua, J.M. Sanz-Serna, Higher-Order averaging, formal series and numerical
integration I : B-series, Found. Comput. Math. 10, 695–727 (2010).
[7] P. Chartier, A. Murua, J.M. Sanz-Serna, Higher-Order averaging, formal series and numerical
integration II : the quasi-periodic case, Found. Comput. Math. 12, 471-508 (2012).
[8] P. Chartier, A. Murua, J.M. Sanz-Serna, A formal series approach to averaging : exponentially
small error estimates, DCDS A 32, 3009-3027 (2012).
[9] G. Contopoulos, C. Efthymiopoulos and A. Giorgilli, Non-convergence of formal integrals
of motion, J. Phys. A 36, 8369–8660 (2003).
[10] A. Durga Devi, R. Gladwin Pradeep, V. K. Chandrasekar, M. Laksmanan, Method of genera-
ting N-dimensional nonsigular Hamiltonian systems, J. Nonlinear Math. Phys. to appear.
[11] C. Efthymiopoulos, A. Giorgilli and G. Contopoulos, Non-convergence of formal integrals :
II. Improved estimates for the optimal order of truncation, J. Phys. A 37, 10831–10858
(2004).
[12] E. Hairer, Ch. Lubich and G. Wanner, Geometric Numerical Integration, 2nd ed., Springer,
Berlin, 2006.
123
124 BIBLIOGRAPHIE
[13] E. Hairer, S. P. Nørsett and G. Wanner, Solving Ordinary Differential Equations I, Nonstiff
Problems, 2nd ed., Springer, Berlin, 1993.
[14] P. Lochak, C. Meunier, Multiphase Averaging for Classical Systems, with Applications to
Adiabatic Theorems, Springer, New York, 1988.
[15] R.I. McLachlan, On the numerical integration of ordinary differential equations by symmetric
composition methods, SIAM J. Sci. Comput. 16 (1995) 151–168.
[16] A. Murua, The Hopf algebra of rooted trees, free Lie algebras and Lie series, Found. Comput.
Math. 6, 387–426 (2006).
[17] A.I. Neishtadt, The separation of motions in systems with rapidly rotating phase, J. Appl.
Math. Mech. 48, 133–139 (1984).
[18] L.M. Perko, Higher order averaging and related methods for perturbed periodic and quasi-
periodic systems, SIAM J. Appl. Math. 17, 698–724 (1968).
[19] J. A. Sanders, F. Verhulst and J. Murdock, Averaging Methods in Nonlinear Dynamical
Systems (2nd. ed.), Springer, New York, 2007.
[20] J. M. Sanz-Serna and M. P. Calvo, Numerical Hamiltonian Problems, Chapman and Hall,
London, 1994.
[21] C. Simó, Averaging under fast quasi-periodic forcing. In : I. Seimenis (ed.) Proceedings of
the NATO-ARW Integrable and Chaotic Behaviour in Hamiltonian Systems, Torun, Poland,
1993, pp. 13–34. Plenum, New York, 1994.
[22] M. Suzuki, General theory of higher-order decomposition of exponential operators and sym-
plectic integrators, Phys. Lett. A 165 (1992) 387–395.
[23] H. Yoshida, Construction of higher order symplectic integrators, Phys. Lett. A 150 (1990)
262–268.