Modélisation du mouvement d’une

plateforme en mer
Sujet du jeudi 2 mai 2019, Concours Commun INP, filière PC, modélisation

On s’intéresse à la résolution d’équations du mouvement dans une approche classique de la mécanique afin
d’étudier le mouvement d’une plateforme en mer. Le modèle envisagé est un système à un degré de liberté
considéré comme oscillateur harmonique : une masse est reliée à un ressort, avec ou sans amortissement, et
peut être soumise à une excitation externe.
La résolution est tout d’abord abordée de façon analytique puis de façon numérique, avant enfin de comparer
les résultats obtenus.
Les résolutions analytiques et numériques sont largement indépendantes.
Dans la suite de l’énoncé, toutes les grandeurs vectorielles sont indiquées en gras. Un aide-mémoire sur
numpy est donné en Annexe.
On considère le mouvement d’une plateforme en mer soumise à un courant marin. Sa partie supérieure de
masse 𝑚 = 110 tonnes est considérée comme rigide et le mouvement principal de la plateforme a lieu suivant 𝑥
(figure 1a).
Afin d’étudier le mouvement de cette plateforme, on la représente par une masse 𝑚, liée à un ressort de
constante de raideur 𝑘 et à un amortisseur de constante d’amortissement 𝛾, pouvant subir une excitation externe
# 
de force 𝐹𝑒𝑥𝑐 , et se déplaçant sur un support (figure 1b). Le ressort représente la rigidité de l’ensemble du support
de la plateforme. L’amortisseur permet de prendre en compte l’effet de l’eau environnante et la force d’excitation
externe celui des vagues qui frappent périodiquement la plateforme. La masse est supposée se déplacer selon
une seule direction parallèle à l’axe 𝑂𝑥 en fonction du temps 𝑡.
Les projections sur l’axe 𝑂𝑥 de la position, de la vitesse et de l’accélération de la masse en fonction du temps
#
sont notées respectivement 𝑥(𝑡), 𝑥(𝑡) ̇ et 𝑥(𝑡).
̈ La force totale 𝐹𝑡𝑜𝑡 agissant sur la masse correspond à la réaction
# # # #
normale 𝑅𝑁 de la base horizontale, à la force de frottement 𝐹𝑑 , à la force de rappel 𝐹𝑘 du ressort, au poids 𝑃
# 
de la masse et à la force 𝐹𝑒𝑥𝑐 d’excitation externe. La position d’équilibre de la masse sera choisie à 𝑥 = 0. En
l’absence d’action de l’amortisseur, la masse se déplace sur la base horizontale sans frottements.

0.1 Résolution analytique et détermination des paramètres pour la modélisation

# #
Q1 — En effectuant une projection sur l’axe 𝑂𝑥, montrer que 𝑃 et 𝑅𝑁 n’interviennent pas dans le bilan des
forces.

Solution : On applique le PFD dans le référentiel terrestre, galiléen. Le bilan des forces est donné par
l’énoncé, on obtient :
# #  # # # # 
𝑚 #
𝑎 = ∑ 𝐹 = 𝑅𝑁 + 𝐹𝑑 + 𝐹𝑘 + 𝑃 + 𝐹𝑒𝑥𝑐
#
Le mouvement est suivant 𝑢#𝑥 uniquement, la projection sur 𝑂𝑧 la verticale indique que le poids et 𝑅𝑁

masse m

support k
Fexc
vague m
γ

x
0
base (b) système masse (𝑚), ressort (𝑘), amortisseur (𝛾) et ex-
#
x citation externe (𝐹𝑒𝑥𝑐 )
(a) Plateforme en mer soumise aux vagues marines

Figure 1 : Schématisations du problème modélisé.

 s’annulent, et suivant l’horizontale on obtient l’équation du mouvement :

𝑚𝑥 ̈ = 0 + 𝐹𝑑 + 𝐹𝑘 + 0 + 𝐹𝑒𝑥𝑐

Le poids et la réaction normale n’interviennent pas dans cette équation du mouvement, mais font cependant
bel et bien partie du bilan des forces…

0.1.1 Ressort sans amortissement et sans excitation

Q2 —Démontrer que l’équation du mouvement de la masse correspond à l’équation différentielle du second
ordre suivante :
𝑚𝑥 ̈ + 𝑘𝑥 = 0 (1)

Solution : 𝐹𝑑 = 0 = 𝐹𝑒𝑥𝑐 traduit les conditions de travail de cette partie, et la force élastique est une force
de rappel vers la position d’équilibre. Dans la limite d’élasticité de la structure, 𝐹𝑘 = −𝑘𝑥 puisque la position
d’équilibre est choisie comme origine de l’axe horizontal. On a donc bien finalement 𝑚𝑥 ̈ + 𝑘𝑥 = 0

Q3 — La solution de cette équation prend la forme générale suivante

𝑥(𝑡) = 𝐴0 sin(𝜔0 𝑡) + 𝐵0 cos(𝜔0 𝑡) (2)

avec 𝐴0 et 𝐵0 deux coefficients réels. Exprimer 𝜔0 en fonction des grandeurs caractéristiques du système
et donner sa signification physique. De plus, en remarquant qu’à 𝑡 = 0 : 𝑥(𝑡) = 𝑥0 et 𝑥(𝑡)
̇ = 𝑥0̇ , déterminer les
expressions de 𝐴0 et de 𝐵0 en fonction de 𝑥0 , 𝑥0̇ et de 𝜔0 .

Solution : L’équation différentielle obtenue est celle d’un oscillateur harmonique de pulsation propre
𝑘
𝜔0 = , cette pulsation est celle des oscillations libres de ce système. La condition initiale sur la position
√𝑚
𝑥0̇ 𝑥0̇
donne 𝑥(𝑡 = 0) = 𝑥0 = 𝐵0 et celle sur la vitesse 𝜔0 𝐴0 = 𝑥0̇ , soit 𝐴0 = et donc 𝑥(𝑡) = 𝑥0 cos 𝜔0 𝑡 + sin 𝜔0 𝑡
𝜔0 𝜔0

Q4 —On cherche à reformuler l’équation précédente sous une forme plus compacte du type :

𝑥(𝑡) = 𝑅0 cos(𝜔0 𝑡 − 𝜑0 ) (3)

Donner les expressions de 𝑅0 et de 𝜑0 en fonction de 𝑥0 , 𝑥0̇ et de 𝜔0 .

Solution :
𝑥(𝑡) = 𝑅0 cos(𝜔0 𝑡 − 𝜑0 ) = 𝑅0 (cos(𝜔0 𝑡) cos(𝜑0 ) + sin(𝜔0 𝑡) sin(𝜑0 )
Par identification avec 𝑥(𝑡) = 𝐴0 sin(𝜔0 𝑡) + 𝐵0 cos(𝜔0 𝑡), il vient :

𝐴0 = 𝑅0 sin(𝜑0 )
{𝐵0 = 𝑅0 cos(𝜑0 )

Ce qui permet d’obtenir 𝑅0 et 𝜑0 en réinjectant les résultats de la question précédente :

⎧⎪
⎪⎪𝑅 = 𝐴2 + 𝐵2 = 𝑥 2 + 𝑥0̇ 2
⎪ 0 √ 0 0 0 (𝜔 )
⎨⎪ √ 0
⎪⎪tan(𝜑 ) = 𝐴0 = 𝑥0̇
⎪⎩ 0 𝐵0 𝜔0 𝑥0

Page 2
 Q5 —Représenter qualitativement 𝑥(𝑡) en fonction de 𝑡 et indiquer sur le tracé 𝑅0 , 𝑥0 et 2𝜋/𝜔0 .

Solution : 𝑅0 est l’amplitude de la sinusoïde, 𝑥0 sa valeur à 𝑡 = 0 et 2𝜋/𝜔0 la période :

2𝜋
𝑥 Δ𝑡 =
𝜔0
𝑅0
𝑥0
𝑡

Q6 —En utilisant les expressions des énergies cinétique 𝐾 (𝑡) et potentielle 𝑈 (𝑡) du système, montrer que
l’énergie totale 𝐸(𝑡) du système est alors :

𝑘𝑅02
𝐸(𝑡) = (4)
2
Justifier le résultat obtenu.

1 1 2 1
Solution : Par définition de l’énergie cinétique, 𝐾 (𝑡) = 𝑚𝑣 2 = 𝑚 (−𝑅0 𝜔0 sin(𝜔0 𝑡 − 𝜑0 )) = 𝑚𝑅02 𝜔02 sin2 (𝜔0 𝑡−
2 2 2
𝜑0 ). Pour un système élastique de constante de raideur 𝑘, l’énergie potentielle élastique s’exprime sous la
1 1
forme 𝑈 = 𝑘𝑥 2 , soit 𝑈 (𝑡) = 𝑘 2 𝑅02 cos2 (𝜔0 𝑡 − 𝜑0 ). Avec 𝜔02 = 𝑘/𝑚, on en déduit 𝐸(𝑡) = 𝐾 (𝑡) + 𝑈 (𝑡) =
2 2
1 1 2 2 1 2
𝑘𝑅02 sin2 (𝜔0 𝑡 − 𝜑0 ) + 𝑘 𝑅0 cos2 (𝜔0 𝑡
𝑘𝑅 .− 𝜑0 ) =
2 2 2 0
𝐸 est une constante car il n’y a que des forces conservative. Lorsque l’énergie cinétique est nulle, le
déplacement esst maximal, vaut ±𝑅0 et toute l’énergie est sous forme potentielle élastique, d’où cette ex-
pression de cette constante.

Q7 —Représenter qualitativement 𝐸(𝑡), 𝐾 (𝑡) et 𝑈 (𝑡) en fonction de 𝑡.

Solution :
𝐸, 𝑈 , 𝐾 2𝜋
Δ𝑡 =
𝜔0
1
𝐸 = 𝑘𝑅02
2 𝑈 (𝑡)
1
𝑘𝑥 𝐾 (𝑡)
2 0 𝑡

0.1.2 Ressort avec amortissement et sans excitation

Q8 —La force de frottement que l’amortisseur exerce sur la masse est considérée comme linéaire, c’est-à-
#
dire proportionnelle au vecteur vitesse #
𝑣 de celle-ci : 𝐹𝑑 = −𝛾 #
𝑣 , avec 𝛾 la constante d’amortissement, positive.
En considérant une projection sur l’axe 𝑂𝑥, démontrer que la position de la masse en fonction du temps suit
l’équation du mouvement ci-après

𝑥 ̈ + 2𝜁 𝜔0 𝑥 ̇ + 𝜔02 𝑥 = 0 (5)

avec 𝜔0 défini en question 3 et 𝜁 à exprimer en fonction de 𝛾, 𝑘 et 𝑚.

Page 3
 𝛾 𝑘
Solution : On repart de la question 1 avec 𝐹𝑑 = −𝛾𝑥 ̇ et toujours 𝐹𝑒𝑥𝑐 = 0, on obtient 𝑥 ̈ + + 𝑥=0.
𝑚 𝑚
Par identification avec la forme de l’équation ??, on obtient bien 𝜔02 = 𝑘/𝑚 comme précédemment, et
𝛾 𝛾
2𝜁 𝜔0 = , ce qui donne 𝜁 =
𝑚 2√𝑚𝑘

Q9 —Dans le cas où 𝜁 < 1, 𝑥(𝑡) prend la forme suivante :

𝑥(𝑡) = 𝑒 −𝜁 𝜔0 𝑡 (𝐴𝑑 cos(𝜔𝑑 𝑡) + 𝐵𝑑 sin(𝜔𝑑 𝑡)) (6)

Déterminer les deux coefficients réels 𝐴𝑑 et 𝐵𝑑 en fonction de 𝑥0 , 𝑥0̇ , 𝜁, 𝜔0 et 𝜔𝑑 = 𝜔0 · √1 − 𝜁 2 . On utilisera

pour cela les mêmes conditions initiales que celles utilisées en question 3.

Solution : À 𝑡 = 0, 𝑥(𝑡 = 0) = 𝑥0 = 𝐴𝑑 𝑒 0 , donc 𝐴𝑑 = 𝑥0 . On calcule ensuite 𝑥 ̇ pour pouvoir utiliser la

seconde condition initiale :

̇ = 𝑒 −𝜁 𝜔0 𝑡 ((−𝜁 𝜔0 𝑥0 + 𝐵𝑑 𝜔𝑑 ) cos(𝜔𝑑 𝑡) + (−𝜁 𝜔0 𝐵𝑑 − 𝑥0 𝜔𝑑 ) sin(𝜔𝑑 𝑡))

𝑥(𝑡)

𝑥0̇ + 𝜁 𝜔0 𝑥0
d’où 𝑥(𝑡
̇ = 0) = −𝜁 𝜔0 𝑥0 + 𝐵𝑑 𝜔𝑑 = 𝑥0̇ et donc 𝐵𝑑 =
𝜔𝑑

Q10 —Montrer alors que l’on peut obtenir une forme du type

𝑥(𝑡) = 𝑅𝑑 𝑒 −𝜁 𝜔0 𝑡 cos(𝜔𝑑 𝑡 − 𝜙𝑑 ) (7)

avec 𝑅𝑑 et 𝜙𝑑 à préciser.

Solution : On développe cos(𝜔𝑑 𝑡 − 𝜙𝑑 ) = cos(𝜔𝑑 𝑡) cos(𝜙𝑑 ) + sin(𝜔𝑑 𝑡) sin(𝜙𝑑 ), ce qui permet d’identifier les
𝐴𝑑 = 𝑅𝑑 cos(𝜙𝑑 ) 𝑅𝑑 = √𝐴2𝑑 + 𝐵𝑑2
termes entre les deux formes ce qui donne On pourrait réinjecter
{𝐵𝑑 = 𝑅𝑑 sin(𝜙𝑑 ) {tan(𝜙𝑑 ) = 𝐵𝑑 /𝐴𝑑
les valeurs de 𝐴𝑑 et 𝐵𝑑 en fonction des conditions initiales, mais ce n’est pas explicitement demandé et sans
intérêt majeur.

Q11 —Représenter qualitativement 𝑥(𝑡) en fonction de 𝑡 et indiquer sur le tracé 𝑅𝑑 𝑒 −𝜁 𝜔0 𝑡 , 𝑥0 et 2𝜋/𝜔𝑑 .

Solution :
𝑥 2𝜋
Δ𝑡 =
𝜔𝑑

𝑅𝑑

𝑥0 𝑅𝑑 𝑒 −𝜁 𝜔0 𝑡

𝑡
1/𝜁 𝜔0
𝑅𝑑 𝑒 −𝜁 𝜔0 𝑡 cos(𝜔𝑑 𝑡 − 𝜙𝑑 )

−𝑅𝑑

Page 4
 Q12 —Donner l’expression de 𝐸(𝑡) et commenter les cas où 𝜁 = 0 et 𝜁 = 1.

Solution : Pour 𝜁 = 0, l’énergie 𝐸 est constante : c’est le cas étudié précédemment. Pour 𝜁 = 1, 𝜔𝑑 = 0 les
expressions obtenues précédemment ne sont plus valables car cela correspond au régime critique et donc
plus au régime pseudo-périodique sur lequel les expressions fournies par l’énoncé sont valables. On a alors
un retour à l’équilibre particulièrement rapide. Pour 1 > 𝜁 > 0, on peut essayer de se laisser guider par le
calcul précédent :

1 1 1 𝑥̇2 1
𝐸(𝑡) = 𝐾 (𝑡) + 𝑈 (𝑡) = 𝑚𝑥 ̇ 2 + 𝑘𝑥 2 = 𝑘 + 𝑘𝑥 2
2 2 2 𝜔02 2
2
1 d (𝑒 −𝜁 𝜔0 𝑡 cos(𝜔𝑑 𝑡−𝜙𝑑 )) 1
= 𝑘𝑅02 + 𝑘 2 𝑅02 𝑒 −2𝜁 𝜔0 𝑡 cos2 (𝜔𝑑 𝑡 − 𝜙𝑑 )
2 ( 𝜔0 d 𝑡 ) 2
1 𝜔 2
= 𝑘𝑅02 𝑒 −2𝜁 𝜔0 𝑡 (−𝜁 cos(𝜔𝑑 𝑡 − 𝜙𝑑 ) − 𝑑 sin(𝜔𝑑 𝑡 − 𝜙𝑑 )) + 𝑒 −2𝜁 𝜔0 𝑡 cos2 (𝜔𝑑 𝑡 − 𝜙𝑑 )
2 [ 𝜔0 ]
1 2 −2𝜁 𝜔 𝑡 2 2 2 2
= 𝑘𝑅0 𝑒 0
(𝜁 cos (𝜔𝑑 𝑡 − 𝜙𝑑 ) + (1 − 𝜁 ) sin (𝜔𝑑 𝑡 − 𝜙𝑑 ) + …
2
… +2𝜁√1 − 𝜁 2 cos(𝜔𝑑 𝑡 − 𝜙𝑑 ) sin(𝜔𝑑 𝑡 − 𝜙𝑑 ) + cos2 (𝜔𝑑 𝑡 − 𝜙𝑑 ))
1
= 𝑘𝑅02 𝑒 −2𝜁 𝜔0 𝑡 (1 + 𝜁 2 (cos2 (𝜔𝑑 𝑡 − 𝜙𝑑 ) − sin2 (𝜔𝑑 𝑡 − 𝜙𝑑 )) + 2𝜁√1 − 𝜁 2 cos(𝜔𝑑 𝑡 − 𝜙𝑑 ) sin(𝜔𝑑 𝑡 − 𝜙𝑑 ))
2
1
= 𝑘𝑅02 𝑒 −2𝜁 𝜔0 𝑡 (1 + 𝜁 2 cos(2𝜔𝑑 𝑡 − 2𝜙𝑑 ) + 𝜁√1 − 𝜁 2 sin(2𝜔𝑑 𝑡 − 2𝜙𝑑 ))
2

Pas certain de l’intérêt de ce calcul. On retrouve bien la valeur constante pour 𝜁 = 0. Pour 𝜁 → 1 dans
cette expression, le terme en sinus disparaît, et le cosinus est constant le temps de la décroissance de l’ex-
ponentielle, puisque 𝜔𝑑 = 0. On a alors principalement une décroissance exponentielle de l’énergie selon
𝐸(𝑡) = 𝐸0 𝑒 −2𝜔0 𝑡 .

Q13 —Montrer de façon simple que 𝐸 est une fonction décroissante de 𝑡. À quoi cela est-il dû ?

Solution : Le théorème de l’énergie mécanique écrit sur les puissances nous indique que la dérivée tem-
porelle de 𝐸 est égale à la puissance des forces non conservatives, c’est-à-dire ici la puissance de la force de
frottement :
d𝐸
= −𝛾 #
𝑣 ⋅ #
𝑣 = −𝛾𝑥 ̇ 2
d𝑡
𝛾 étant une constante positive, cette expresssion est nécessairement toujours négative et l’énergie totale
𝐸 est bien une fonction décroissante du temps. C’est dû à la caractéristique principale de toute force de
frottement, qui s’oppose toujours a déplacement (ce qui ce traduit dans les équations en imposant 𝛾 > 0).

Q14 — On envisage deux temps successifs 𝑡1 et 𝑡2 pour lesquels les déplacements sont 𝑥1 et 𝑥2 , tels que
𝑠2 > 𝑡1 et 𝑡2 − 𝑡1 = 𝜏𝑑 , avec 𝜏𝑑 : période des oscillations amorties. En utilisant l’équation 7 et en considérant que
𝜁 ≪ 1, montrer que :
ln(𝑥1 /𝑥2 ) ≈ 2𝜋𝜁 . (8)

Solution : Les deux instants étant séparés d’une période de la fonction sinusoïdale, les cosinus ont la
même valeur et se simplifient en faisant le rapport, tout comme l’amplitude 𝑅𝑑 . Dès lors il ne reste dans ce
rapport que les termes exponentiels :

𝑒 −𝜁 𝜔0 𝑡1 2𝜋 2𝜋𝜁
ln(𝑥1 /𝑥2 ) = ln = 𝜁 𝜔0 (𝑡2 − 𝑡1 ) = 𝜁 𝜔0 =
( 𝑒 −𝜁 𝜔0 𝑡2 ) 𝜔𝑑 √1 − 𝜁 2

Dans la limite 𝜁 ≪ 1, √1 − 𝜁 2 ≈ 1 et on obtient bien ln(𝑥1 /𝑥2 ) ≈ 2𝜋𝜁

Q15 — Le relevé du déplacement horizontal de la plateforme en fonction du temps est représenté en figure 2.
En utilisantl’ les deux points qui sont indiqués sur la figure, déterminer 𝑘, 𝜁 et 𝛾. Comment ce tracé serait
modifié en fonction de la valeur de 𝜁 ?

Page 5
 0,03
(8,008008 ; 0,010661)
0,02
(4,004004 ; 0,014602)

0,01

x(m)
0,00

− 0,01
t1 t2
− 0,02
0 2 4 6 8 10 12 14
t( s )

Figure 2 : Relevé du déplacement horizontal x (en m) de la plateforme de masse m = 110 tonnes en fonction du
temps t (en s). Les deux temps t1 et t2 mentionnés en question 14 sont indiqués.

Solution : La figure nous indique 2 points séparés d’une période de la pseudo-oscillation, on utilise pour ces
1
points le résultat précédent. Avec 𝑥1 = 1,460 2 cm et 𝑥2 = 1,066 1 cm, 𝜁 = ln(𝑥1 /𝑥2 ) = 5,01 × 10−2 = 𝜁 . La
2𝜋
condition 𝜁 ≪ 1 est raisonnablement vérifiée. En conservant la même approximation, on a 𝜔𝑑 ≈ 𝜔0 , qui peut
2
2𝜋 2𝜋
être évalué via les valeurs de 𝑡1 et 𝑡2 : 𝜔0 ≈ . Comme 𝑘 = 𝑚𝜔02 , 𝑘 = 𝑚 = 2,71 × 105 kg ⋅ s−2 .
𝑡2 −𝑡1 ( 𝑡2 − 𝑡 1 )

Enfin, 𝛾 = 2𝜁√𝑚𝑘 = 1,73 × 104 kg ⋅ s−1 .

Lorsque le coefficient d’amortissement 𝜁 augmente, la décroissance est plus rapide, et lorsqu’on quitte
la condition 𝜁 ≪ 1, la période des pseudo-oscillation augmente (𝜔𝑑 décroit avec la croissance de 𝜁)

0.1.3 Ressort avec amortissement et avec excitation

On envisage enfin le cas où le système est soumis à la fois aux esffets d’amortissement et d’excitation.
On se limite ici à la réponse à une excitation harmonique sinusoïdale de fréquence 𝜔 produite par une force
extérieure au système
# 
𝐹𝑒𝑥𝑐 (𝑡) = 𝐹0 cos(𝜔𝑡)𝑒#𝑥 (9)
avec 𝑒#𝑥 vecteur unitaire sur l’axe 𝑂𝑥 et on se place dans le cas traité précédemment pour l’étude de l’amor-
tisseur, c’est-à-dire 𝜁 < 1 (partie 0.1.2).
On admet de plus dans ce qui suit que la réponse du système dans le cas où amortisseur et excitation sont
pris en compte peut s’écrire comme somme de la solution donnée par l’équation 6 et de la contribution due à
l’excitation :
𝑥𝑒𝑥𝑐 (𝑡) = 𝑋 cos(𝜔𝑡 − 𝜙). (10)
Q16 —Montrer que l’équation différentielle caractérisant le système devient alors :
𝐹0
𝑥 ̈ + 2𝜁 𝜔0 𝑥 ̇ + 𝜔02 𝑥 = cos(𝜔𝑡). (11)
𝑚

Solution : On reprend le PFD établit à la question 1, et on obtient pour la projection de ce PFD sur l’axe
𝑂𝑥 𝑚𝑥 ̈ = −𝛾𝑥 ̇ − 𝑘𝑥 + 𝐹0 cos(𝜔𝑡), ce qui se réécrit :

𝛾 𝑘 𝐹0
𝑥̈ + 𝑥̇ + 𝑥 = cos(𝜔𝑡)
𝑚 𝑚 𝑚

𝑘 𝛾
On a bien l’expression voulue en gardant comme précédemment 𝜔02 = et 2𝜁 𝜔0 =
𝑚 𝑚

Page 6
 Q17 —En utilisant l’équation 10 et en privilégiant une représentation complexe, vérifier que :
⎧⎪𝑋 = 𝐹0 ⋅ 1
⎪⎪
⎪ 𝑚 𝜔 2
−𝜔 2 2 + 2𝜁 𝜔 𝜔 2
⎨⎪ √( 0 ) ( 0 ) (12)
⎪⎪tan 𝜙 = 2𝜁 𝜔0 𝜔
⎪⎩ 𝜔02 −𝜔 2

Solution : L’équation différentielle étant linéaire, traiter la solution particulière et la solution de l’équation
homogène séparément pour ensuite les additionner et obtenir l’ensemble des solutions ne surprend guère.
On introduit les représentations complexes des grandeurs évoluant sinusoïdalement, la grandeur com-
plexe est soulignée :
(𝑡) = 0 cos(𝜔𝑡 + 𝜙) ↔ (𝑡) = 0 𝑒 𝑗𝜔𝑡+𝜙 = 0 𝑒 𝑗𝜔𝑡
Pour les représentations complexes des grandeurs, l’équation 11 devient, les exponentielles se factori-
sant puis pouvant être simplifiées (une dérivée temporelle devient une multiplication par 𝑗𝜔) :
𝐹0
𝑋 (−𝜔 2 + 2𝜁 𝜔0 𝑗𝜔 + 𝜔02 ) =
𝑚
D’où :
𝐹0
𝑚
𝑋=
−𝜔 2 + 2𝜁 𝜔0 𝑗𝜔 + 𝜔02
𝑋 et 𝜙 sont respectivement le module et l’argument de 𝑋, on en déduit bien les résultats annoncés :
⎧⎪ 𝐹0
⎪⎪𝑋 = |𝑋| = 𝑚
⎪⎪ 2 2 2
2
⎨⎪ √(𝜔02−𝜔 2) +(2𝜁 𝜔0 𝜔)
⎪⎪ 𝐼 𝑚(𝜔0 −𝜔 +𝑗⋅2𝜁 𝜔0 𝜔) 2𝜁 𝜔0 𝜔
⎪⎪tan 𝜙 = − =
⎩ 𝑅𝑒(𝜔02 −𝜔 2 +𝑗⋅2𝜁 𝜔0 𝜔) 𝜔02 −𝜔 2

𝑋
Q18 —Exprimer la grandeur 𝑀 = en fonction de 𝑟 = 𝜔/𝜔0 et expliciter le sens physique de 𝑀.
𝐹0 /𝑘

𝑘
𝑚 1 1
Solution : 𝑀 = = =
𝜔 2
−𝜔 2 2 + 2𝜁 𝜔 𝜔 2 2 2
2
2
2
1−𝑟 2 +4𝜁 2 𝑟 2
√ ( 0 ) ( 0 ) (𝜔0 −𝜔 ) (2𝜁 𝜔0 𝜔) √( )
4 + 4
𝜔0 𝜔0
√
𝑀 est un facteur sans dimension qui caractérise la réponse fréquentielle du système, ce qui ramène
au cas classique d’une fonction de transfert lorsque l’excitation et la réponse sont mesurées sur la même
grandeur physique.

Q19 —Trouver la condition sur 𝑟 puis sur 𝜔 pour laquelle 𝑀 est maximale.

2
Solution : 𝑀(𝑟) atteindra son maximum lorsque (1 − 𝑟 2 ) + 4𝜁 2 𝑟 2 = 𝑟 4 + (4𝜁 2 − 2)𝑟 2 + 1 sera minimal.
Cette expression est un polynôme du second degré en 𝑅 = 𝑟 2 , le minimum est atteint en « −𝑏/2𝑎 », soit
1
pour 𝑅 = 1 − 2𝜁 2 , soit 𝑟 = √1 − 2𝜁 2 si 𝜁 < . Dans la limite 𝜁 ≪ 1 du problème, cette condition est bien
√2
vérifiée. On en déduit la pulsation de résonance pour laquelle 𝑀(𝜔) est maximale, 𝜔𝑟 = 𝜔0 ⋅ 𝑟 = 𝜔0 √1 − 2𝜁 2

Q20 —Si l’on considère une période moyenne des vagues en mer de 8 s, que peut-on conclure sur le mouve-
ment de la plateforme ?

Solution : La période moyenne des vagues en mer est de 8 s, la valeur du coefficient d’amortissement et de
la période propre obtenus en exploitant la figure 2 nous permettent de dire que la fréquence de résonance

Page 7
 sera proche de la fréquence propre, puisque 𝜔𝑟 = 𝜔0 √1 − 2𝜁 2 ≈ 𝜔0 pour 𝜁 ≪ 1. La plateforme peut entrer
en résonance pour une excitation de période proche de 4 s, pour une excitation avec une période de 8 s
on est suffisamment loin de la résonance pour ne pas avoir le développement de mouvements de grande
amplitude.

0.2 Modélisation : codage

On souhaite maintenant obtenir 𝑥(𝑡) et 𝐸(𝑡) de façon numérique et comparer les résultats obtenus à ceux
fournis par les solutions analytiques précédentes pour 𝑥(𝑡). On rappelle que 𝑥(𝑡) et 𝐸(𝑡) représentent respecti-
vement la position de la masse et l’énergie mécanique totale en fonction du temps.
Pour cela, le temps est discrétisé en 𝑁 points 𝑡 = 0, ∆𝑡, 2∆𝑡, ..., (𝑁 − 1)∆𝑡 avec un pas de temps constant ∆𝑡.
Les 𝑁 − 1 pas sont effectués pendant la simulation de durée totale 𝑡𝑚𝑎𝑥 . On note respectivement 𝑥𝑛 , 𝑣𝑛 , 𝑎𝑛 , 𝐸𝑛
# 
et 𝐹𝑛 les valeurs de 𝑥(𝑡), 𝑥(𝑡),
̇ 𝑥(𝑡),
̈ 𝐸(𝑡) et 𝐹𝑒𝑥𝑐 (𝑡) à 𝑡 = 𝑛∆𝑡.
À chaque pas, les équations du mouvement reliant 𝑥𝑛+1 et 𝑣𝑛+1 à 𝑥𝑛 et 𝑣𝑛 sont utilisées afin d’obtenir
les valeurs de 𝑥, 𝑣 et 𝐸. Les conditions initiales 𝑥0 et 𝑣0 sont connues et permettent de démarrer le processus
d’intégration numérique. Deux algorithmes distincts (Euler et Leapfrog) vont être utilisés dans la suite.
Pour l’écriture du code, on se place dans le cas le plus général, c’est-à-dire avec amortissement et excitation
harmonique externe. Les variables et tableaux choisis sont notamment listés dans la table 1

N : nombre 𝑁 de points sur l’axe des temps utilisés pendant toute la simulation
t[] : tableau des temps 𝑡 (s), de dimension 𝑁
x[] : tableau des positions 𝑥 (m), de dimension 𝑁
v[] : tableau des vitesses 𝑣 (m ⋅ s−1 ), de dimension 𝑁
E[] : tableau des énergies totales (J), de dimension 𝑁
F[] : tableau des forces d’excitation (N), de dimension 𝑁
dt : pas de temps Δ𝑡 (s)
tmax : temps total de la simulation 𝑡𝑚𝑎𝑥 (s)
k : constante de raideur 𝑘 du ressort (N ⋅ m−1 )
m : masse 𝑚 du système (kg)
om0 : 𝜔0 (s−1 )
zeta : 𝜁 (sans unité)

Table 1 : Principaux tableaux et principales variables utilisés pour la résolution numérique

On rappelle qu’un aide-mémoire sur numpy est fourni en Annexe, page 7.

Q21 —Écrire les lignes de code permettant de définir l’entier 𝑁. On suppose 𝑡𝑚𝑎𝑥 et ∆𝑡 connus et fixés par
l’utilisateur en début de code.

Solution : Le nombre de points est le nombre d’intervalles plus un, on divise le temps total par la durée
d’un intervalle pour avoir le nombre d’intervalles. Le code étant exécuté pour la rédaction de ce document,
les réponses fournies contiennent un peu plus que ce qui est demandé afin d’illustrer les commande sans
obtenir de message d’erreur. Ainsi pour cette question, la ligne 5 est le seul code demandé.

1 >>> from math import *

2 >>> import numpy as np
3 >>> tmax=5
4 >>> dt=1
5 >>> N=int(tmax/dt + 1)
6 >>> N
7 6

Q22 —Écrire alors l’instruction permettant de définir le tableau t qui contient toutes les valeurs de 𝑡 telles
que : 0𝑠 ≤ 𝑡 ≤ 𝑡𝑚𝑎𝑥 . On rappelle que le temps est discrétisé en 𝑁 points 𝑡 = 0, ∆𝑡, 2∆𝑡, ..., (𝑁 − 1)∆𝑡 avec un pas de
temps constant ∆𝑡 et que 𝑁 − 1 pas sont effectués.

Page 8
 (𝑥0 ; 𝑣0 ) (𝑥1 ; 𝑣1 ) (𝑥2 ; 𝑣2 )

𝑡0 𝑡1 𝑡2

Figure 3 : Discrétisation en temps utilisée dans le cas de l’algorithme d’Euler. Les conditions initiales corres-
pondent à (𝑥0 ; 𝑣0 )

Solution :
1 >>> t=np.linspace(0,tmax,N)
2 >>> t
3 array([0., 1., 2., 3., 4., 5.])

Q23 —En effectuant des développements de Taylor de 𝑥 et 𝑣 tronqués à l’ordre 1, on obtient l’algorithme
d’Euler, où 𝑥 et 𝑣 sont évalués au même temps 𝑡 selon le schéma donné figure 3 :

𝑥𝑛+1 ≈ 𝑥𝑛 + 𝑣𝑛 ⋅ Δ𝑡
(13)
{𝑣𝑛+1 ≈ 𝑣𝑛 + 𝑎𝑛 ⋅ Δ𝑡
Donner les expressions de 𝑥𝑛+1 et 𝑣𝑛+1 en fonction de 𝑥𝑛 , 𝑣𝑛 et 𝐹𝑛 .

Solution : On remplace dans la formule de récurrence l’accélération 𝑎𝑛 par son expression obtenue à partir
𝐹
de l’équation différentielle du système, 𝑎𝑛 = 𝑛 − 2𝜁 𝜔0 𝑣𝑛 − 𝜔02 𝑥𝑛 . On en déduit :
𝑚

𝑥𝑛+1 ≈ 𝑥𝑛 + 𝑣𝑛 ⋅ Δ𝑡
𝐹𝑛
{𝑣𝑛+1 ≈ 𝑣𝑛 + ( − 2𝜁 𝜔0 𝑣𝑛 − 𝜔02 𝑥𝑛 ) ⋅ Δ𝑡
𝑚

Q24 —Écrire la boucle en i permettant d’obtenir toutes les valeurs de x[i+1], v[i+1] et E[i+1] où i
correspond à un point sur l’axe des temps. On précisera les variables éventuellement introduites en supposant
qu’elles ont été définies dans le code.

Solution :
1 >>> ## les paramètres expérimentaux doivent être rentrés: zeta, w0, la force
2 >>> zeta=5e-2
3 >>> k=2.7e5
4 >>> om0=2*np.pi/4 #T0=4s
5 >>> om=2*np.pi/8 #Tv=8s pour la période des vagues
6 >>> F0=1e5
7

8 >>> F=np.linspace(-F0,F0,N) #arbitrairement, pour exister et F/m raisonnable

9 >>> m=110e3 #110 tonnes
10

11 >>> # initialisation des tableaux, cf Annexe

12 >>> x=np.zeros(N)
13 >>> v=np.zeros(N)
14 >>> E=np.zeros(N)
15

16 >>> # valeurs initiales

17 >>> x0=0.02

Page 9
 𝑣−1/2 𝑥0 𝑣1/2 𝑥1 𝑣3/2 𝑥2 𝑣5/2

𝑡−1/2 𝑡0 𝑡1/2 𝑡1 𝑡3/2 𝑡2 𝑡5/2

Figure 4 : Discrétisation en temps utilisée dans le cas de l’algorithme de Leapfrog. Les conditions initiales
correspondent à (𝑥0 ; 𝑣−1/2 )

18 >>> v0=0
19 >>> E0=0.5*k*x0**2+0.5*m*v0**2
20

21 >>> ## la question commence là:

22 >>> x[0]=x0
23 >>> v[0]=v0
24

25 >>> for i in range(N-1):

26 ... x[i+1]=x[i]+v[i]*dt
27 ... v[i+1]=v[i]+(F[i]/m-2*zeta*om0*v[i]-x[i]*om0**2)*dt
28 ... E[i+1]=0.5*k*x[i+1]**2+0.5*m*v[i+1]**2
29 ...
30 >>> x
31 array([ 0.02 , 0.02 , -0.93843893, -2.34112919, -1.38979832,
32 5.37042079])
33 >>> v
34 array([ 0. , -0.95843893, -1.40269026, 0.95133088, 6.7602191 ,
35 9.67297081])
36 >>> E
37 array([ 0. , 50577.28515526, 227104.82831241,
38 789696.27126188, 2774288.74267099, 9039741.6627917 ])

Q25 —Dans l’algorithme de Leapfrog, les 𝑥 sont évalués aux temps entiers, c’est à dire à 𝑡 = 0, Δ𝑡, 2Δ𝑡, ..., (𝑁 −
1)Δ𝑡, alors que les 𝑣 sont évalués à 𝑡 = −∆𝑡/2, ∆𝑡/2, 3∆𝑡/2, ..., (𝑁 − 1)∆𝑡 − ∆𝑡/2 selon le schéma donné figure 4.
Ainsi, pour cet algorithme, x[i] représente de façon approchée la position à l’instant 𝑖∆𝑡, et v[i] représente
de façon approchée la vitesse à l’instant (𝑖 − 1/2)∆𝑡.
Pour le système considéré, montrer alors que 𝑥𝑛+1 et 𝑣𝑛+1/2 prennent les formes suivantes :
⎧⎪𝑥
⎪⎪ 𝑛+1 ≈ 𝑥𝑛 + 𝑣𝑛+1/2 ⋅ Δ𝑡
⎨⎪ 𝜔02 Δ𝑡 1−𝜁 𝜔0 Δ𝑡 𝐹𝑛 (14)
⎪⎪𝑣𝑛+1/2 ≈ − 𝑥 + 𝑣 + Δ𝑡
⎩ 1+𝜁 𝜔0 Δ𝑡 𝑛 1+𝜁 𝜔0 Δ𝑡 𝑛−1/2 𝑚(1+𝜁 𝜔0 Δ𝑡)

Solution : Le principe du décalage se comprend en considérant que pour évaluer au mieux la valeur d’une
fonction Δ𝑡 plus tard, il faut connaitre la valeur moyenne de la vitesse sur cet intervalle, et que celle-ci est
plus probablement proche de celle du milieu de l’intervalle qu’à une extrémité. On a alors comme équivalent
de l’équation 13 :
𝑥𝑛+1 ≈ 𝑥𝑛 + 𝑣𝑛+1/2 ⋅ Δ𝑡
{𝑣𝑛+1/2 ≈ 𝑣𝑛−1/2 + 𝑎𝑛 ⋅ Δ𝑡
𝑎𝑛 peut être obtenu par l’équation différentielle, mais celle-ci fait intervenir 𝑣𝑛 dont on ne dispose pas.
On le remplace dans l’équation différentielle par la valeur moyenne entre les valeurs de 𝑣𝑚±1/2 , l’équation
différentielle permet donc d’écrire :
𝐹𝑛 𝑣𝑛−1/2 + 𝑣𝑛+1/2
𝑎𝑛 = − 2𝜁 𝜔0 − 𝜔02 𝑥𝑛
𝑚 2

Page 10
 L’équation de récurrence pour les vitesses devient alors :

𝐹𝑛
(1 + 𝜁 𝜔0 Δ𝑡) 𝑣𝑛+1/2 ≈ (1 − 𝜁 𝜔0 Δ𝑡) 𝑣𝑛−1/2 + ( − 𝜔02 𝑥𝑛 ⋅ Δ𝑡
𝑚 )

On obtient alors bien :

⎧⎪𝑥
⎪⎪ 𝑛+1 ≈ 𝑥𝑛 + 𝑣𝑛+1/2 ⋅ Δ𝑡
⎨⎪ 𝜔02 Δ𝑡 1−𝜁 𝜔0 Δ𝑡 𝐹𝑛
⎪⎪𝑣𝑛+1/2 ≈ − 𝑥 + 𝑣 + Δ𝑡
⎩ 1+𝜁 𝜔0 Δ𝑡 𝑛 1+𝜁 𝜔0 Δ𝑡 𝑛−1/2 𝑚(1+𝜁 𝜔0 Δ𝑡)

Q26 —Quel problème pose l’évaluation de E[i+1] ? Dans la suite, on préfèrera ainsi évaluer E[i].

Solution : Pour évaluer E[i+1], il faut la position et la vitesse à l’instant (𝑖 + 1)Δ𝑡. La vitesse n’étant pas
calculée sur les temps entiers, il faut l’évaluer par sa moyenne entre les valeurs décalées de Δ𝑡/2 avant et
après, ce qui demande d’avoir accès à la valeur qui sera stockée dans v[i+2] et qui ne sera évaluée qu’à
l’itération suivante de la boucle. Calculer E[i] par contre est réalisable dans la boucle qui puisqu’on aura
accès à toutes les valeurs de 𝑥 et de 𝑣 utiles.

Q27 —Comment obtenir E[i] à partir de x[i] et v[i] ?

1 1 1 1 𝑣𝑛−1/2 +𝑣𝑛+1/2 2
Solution : 𝐸𝑛 = 𝑘𝑥𝑛2 + 𝑚𝑣𝑛2 = 𝑘𝑥𝑛2 + 𝑚 ( ) ce qui ce code en :
2 2 2 2 2

E[i]=0.5*k*x[i]**2+0.125*m*(v[i]+v[i+1])**2

1
Q28 — En introduisant deux variables fac1 et fac2 correspondant respectivement à 1 − 𝜁 𝜔0 ∆𝑡 et ,
1+𝜁 𝜔0 ∆𝑡
écrire la boucle en i permettant d’obtenir toutes les valeurs de x[i+1] et v[i+1].

Solution :
1 >>> ### les facteurs de proportionnalité
2 >>> fac1=1-zeta*om0*dt
3 >>> fac2=1/(1+zeta*om0*dt)
4

5 >>> ### initialisation: les tableaux existent déjà, je réinitialise juste les
↪ valeurs initiales
6

7 >>> x[0]=x0
8 >>> v[0]=v0
9

10 >>> ### la boucle demandée:

11

12 >>> for i in range(N-1):

13 ... v[i+1]=-fac2*om0**2*dt*x[i]+fac1*fac2*v[i]+F[i]*fac2*dt/m
14 ... x[i+1]=x[i]+v[i+1]*dt
15 ...
16 >>> x
17 array([ 0.02 , -0.86864492, -0.14638139, 0.63699303, 0.01759004,
18 -0.04610936])
19 >>> v
20 array([ 0. , -0.88864492, 0.72226353, 0.78337442, -0.61940299,
21 -0.0636994 ])

Page 11
 Q29 —Compléter la boucle de la question 28 avec le calcul du terme d’énergie totale.

Solution :
1 >>> for i in range(N-1):
2 ... v[i+1]=-fac2*om0**2*dt*x[i]+fac1*fac2*v[i]+F[i]*fac2*dt/m
3 ... x[i+1]=x[i]+v[i+1]*dt
4 ... E[i]=0.5*k*x[i]**2+0.125*m*(v[i]+v[i+1])**2
5 ...
6 >>> ### Il reste le dernier élément du tableau E à remplir, pour lequel il faut un
↪ calcul de vitesse supplémentaire:
7

8 >>> vsup=-fac2*om0**2*dt*x[N-1]+fac1*fac2*v[N-1]+F[N-1]*fac2*dt/m
9 >>> E[N-1]=0.5*k*x[N-1]**2+0.125*m*(v[N-1]+vsup)**2
10

11 >>> E
12 array([ 10912.23471717, 102244.07814066, 34063.21676746, 55147.30811376,
13 6457.91733222, 9765.20164435])

Q30 —Écrire une fonction integration(F) qui prend en argument le tableau F[] des forces d’excitation
et renvoie les tableaux x[], v[] et E[] complétés.
On introduira pour cela une variable algo supposée définie en global, permettant d’appliquer l’algorithme
d’Euler si algo==0 ou de Leapfrog si algo==1.
On considèrera également le cas 𝑡 = 0.

Solution :
1 >>> def integration(F):
2 ... ''' À partir d'une excitation donnée par un tableau numpy, retourne les
↪ tableaux x, v et E, initialisés à 0.'''
3 ... global algo # 0 pour Euler, 1 pour Leapfrop
4 ... ### Initialisation. Le calcul de E[0] est omis car tout est au repos
5 ... x=np.zeros(len(F))
6 ... v=np.zeros(len(F))
7 ... E=np.zeros(len(F))
8 ... ## Valeurs initiales
9 ... x[0],v[0]=x0,v0
10 ... ### Boucles
11 ... if algo==0:
12 ... # Euler
13 ... for i in range(N-1):
14 ... x[i+1]=x[i]+v[i]*dt
15 ... v[i+1]=v[i]+(F[i]/m-2*zeta*om0*v[i]-x[i]*om0**2)*dt
16 ... E[i+1]=0.5*k*x[i+1]**2+0.5*m*v[i+1]**2
17 ... elif algo==1:
18 ... # Leapfrog
19 ... for i in range(N-1):
20 ... v[i+1]=-fac2*om0**2*dt*x[i]+fac1*fac2*v[i]+F[i]*fac2*dt/m
21 ... x[i+1]=x[i]+v[i+1]*dt
22 ... E[i]=0.5*k*x[i]**2+0.125*m*(v[i]+v[i+1])**2
23 ... # dernier élément du tableau E
24 ... vsup=-fac2*om0**2*dt*x[N-1]+fac1*fac2*v[N-1]+F[N-1]*fac2*dt/m
25 ... E[N-1]=0.5*k*x[N-1]**2+0.125*m*(v[N-1]+vsup)**2
26 ... else:
27 ... print("Erreur: pas d'algo reconnu!")

Page 12
28 ... return x,v,E
29 ...
30 >>> algo=1
31 >>> x,v,E=integration(F)
32 >>> x
33 array([ 0.02 , -0.86864492, -0.14638139, 0.63699303, 0.01759004,
34 -0.04610936])
35 >>> v
36 array([ 0. , -0.88864492, 0.72226353, 0.78337442, -0.61940299,
37 -0.0636994 ])
38 >>> E
39 array([ 10912.23471717, 102244.07814066, 34063.21676746, 55147.30811376,
40 6457.91733222, 9765.20164435])

Q31 — Écrire une fonction force(f,t,w) qui prend en argument 𝐹0 , 𝑡, et 𝜔 de l’équation 9 et retourne la
# 
valeur de 𝐹𝑒𝑥𝑐 pour un temps t donné.

Solution :

1 >>> def force(f,t,w):

2 ... ''' Calcul et retourne la force suivant x à un instant donné t, en
↪ fonction de son amplitude et de sa pulsation.'''
3 ... return f*np.cos(w*t)
4 ...
5 >>> force(2,3,2*np.pi/6)
6 -2.0

Q32 —Écrire une fonction force_exc() qui complète et retourne le tableau F[] des forces d’excitation
en fonction du booléen exc défini globalement valant True si une excitation est appliquée au système, False
sinon.
Dans le cas où exc==True, on appellera la fonction force définie précédemment en question 31. Dans le cas
où exc==False, on prendra alors : F[i]=0, ∀i.

Solution : C’est étrange de ne pas passer les caractéristiques de la force en argument, mais les arguments
sont toujours spécifiés dans ce sujet.Je considère que l’amplitude de la force F0 et sa pulsation om sont
définis par ailleurs avant l’appel et que Python ira les chercher en les considérant comme des variables
globales.

1 >>> def force_exc():

2 ... ''' Retourne le tableau décrivant la force d'excitation'''
3 ... global exc
4 ... F=np.zeros(N)
5 ... if exc:
6 ... for i in range(len(F)):
7 ... F[i]=force(F0,i*dt,om)
8 ... return F
9 ...

Q33 —Donner alors les lignes de code permettant de réaliser la simulation numérique à partir des fonctions
précédentes.

Solution :

Page 13
1 >>> ### on reprend la définition des paramètres du problème:
2 >>> zeta = 5e-2
3 >>> k = 2.7e5
4 >>> om0 = 2*np.pi/4 #T0=4s
5 >>> om = 2*np.pi/8 #Tv=8s pour la période des vagues
6 >>> F0 = 5e3 #N, au pif, mais tel que x soit proche des 2cm en statique d'après la
↪ valeur de k et celle de x0...
7 >>> m = 110e3 #110 tonnes
8 >>> x0 = 0.02 #m
9 >>> v0 = 0 #m.s^-1
10

11 >>> tmax=10 #s d'après partie C

12 >>> dt=0.05 #s d'après le tableau fourni
13

14 >>> N=int(tmax/dt + 1)
15

16 >>> algo=1
17 >>> exc=0
18

19 >>> x,v,E=integration(force_exc())
20

21 >>> x[:10]
22 array([0.02 , 0.01988561, 0.01967415, 0.01938097, 0.01901964,
23 0.01860216, 0.01813908, 0.01763971, 0.01711218, 0.0165636 ])
24 >>> len(E)
25 201

Q34 —On souhaite désormais estimer la qualité des résultats numériques par rapport aux données analy-
tiques de référence. Écrire une fonction ema(d,dref) qui calcule l’erreur maximale absolue entre un jeu de
données numériques dn et analytiques da. On supposera que les tableaux dn et da sont de même dimension 𝑛.

Solution :

1 >>> # pour le test

2 >>> dn=np.array([1.1,1.9,3.2,3.6])
3 >>> da=np.array([1,2,3,4])
4

5 >>> def ema(d,dref):

6 ... return np.max(abs(d-dref)) # par défaut sur un tableau numpy les
↪ opérations se font sur chaque élément.
7 ...
8 >>> ema(dn,da)
9 0.3999999999999999

Q35 —Pour le problème particulier qui nous intéresse, si l’on souhaite appliquer cette fonction aux tableaux
contenant les données numériques et analytiques pour 𝐸, quel est l’indice maximal de ces tableaux à considérer ?
Réécrire alors la fonction ema pour qu’elle soit applicable aux deux algorithmes considérés.

Solution : À lire cette question, il n’était probablement pas attendu de calculer le dernier élément du
tableau des énergies : l’ayant fait, la fonction précédente fonctionne pour les deux algorithmes. Si ce n’était
pas le cas, il faudrait dans le cas d’un calcul avec l’algorithme de Leapfrog ne pas tenir compte du dernier
point. On peut tester la valeur de la variable globale pour essayer de savoir quel algorithme, ou poser la
question en le demandant en paramètre. Pour allier les deux possibilités, on peut ajouter un paramètre
optionnel qu’on initialise à la valeur de la variable globale s’il n’est pas renseigné.

Page 14
 EMA
∆𝑡 (s) Euler Leapfrog
0,050 128,0203160 0,0837314
0,010 15,1373529 0,0033484
0,005 7,1159053 0,0008371
0,001 1,3556230 0,0000335

Table 2 : EMA obtenues pour 𝐸 (en J) par rapport à la valeur analytique, pour différents pas de temps ∆𝑡, dans
le cas où il n’y a pas d’amortissement et d’excitation externe.

1 >>> def ema2(d,dref,algorithme=None):

2 ... global algo
3 ... if algorithme==None:
4 ... algorithme=algo
5 ... if algorithme==1:
6 ... return np.max(abs(d[:-1]-dref[:-1]))
7 ... else:
8 ... return np.max(abs(d-dref))
9 ...
10 >>> print(ema2(da,dn),'en utilisant la variable globable par défaut,',algo)
11 0.20000000000000018 en utilisant la variable globable par défaut, 1
12 >>> algo=0
13 >>> print(ema2(da,dn),'en utilisant la variable globable par défaut,',algo)
14 0.3999999999999999 en utilisant la variable globable par défaut, 0
15 >>> print(ema2(da,dn,1),'en précisant utiliser Leapfrog')
16 0.20000000000000018 en précisant utiliser Leapfrog
17 >>> print(ema2(da,dn,0),'en précisant utiliser Euler')
18 0.3999999999999999 en précisant utiliser Euler

0.3 Modélisation : analyse des résultats d’un cas simple

Toutes les données numérique suivantes ont été obtenues avec : 𝑡𝑚𝑎𝑥 = 10 s, 𝑚 = 110 tonns, 𝑥0 = 0,02 m et
𝑣0 = 0 m ⋅ s−1 et la valeur de 𝑘 obtenue en question 15, dans le cas d’un système sans amortissement et sans
excitation externe. Le tableau 2 présente les erreurs maximales absolues (EMA) calculées avec la fonction ema
des énergies obtenues numériquement par rapport à celles obtenues analytiquement, pour les deux algorithmes
envisagés et divers pas de temps.
Q36 —Justifier l’ordre de grandeur des ∆𝑡 considérés du tableau 2 pour la discrétisation en temps utilisée.

Solution : Les Δ𝑡 doivent être petits par rapport aux durées caractéristiques des évolutions pour pouvoir
les suivre, donc petits par rapport à la plus petite période caractéristique du système. Ici, c’est 𝑇0 = 4 s, on
a bien Δ𝑡 ≪ 𝑇0 pour toutes les valeurs, tout en gardant un nombre de points facilement gérable pour les
calculs à faire (𝑁 ≈ 1 × 104 sur 10 s).

Q37 —En utilisant les données numériques du tableau 2 donner l’ordre approximatif de l’erreur globale sur
𝐸 des deux algorithmes considérés. Justifier votre réponse.

Solution : On cherche un coefficient 𝑛 tel que 𝐸𝑀𝐴 ≈ 𝐴(Δ𝑡)𝑛 , on peut donc regarder s’il y a un nombre 𝑛 =
ln(𝐸𝑀𝐴(Δ𝑡1 )/𝐸𝑀𝐴(Δ𝑡2 )
qui se dégage. En prenant Δ𝑡1 = 0,001 s comme référence, on obtient 𝑛 = 1,16, 1,045, 1,040
𝑙𝑛(Δ𝑡1 /Δ𝑡2 )
et donc 𝑛 = 1 probablement pour l’ordre de l’erreur pour l’algorithme d’Euler et 𝑛 = 2 (remarquablement
précisément !) pour Leapfrog.

Q38 —Dans le cas simple de l’algorithme d’Euler, comment augmente 𝐸 lorsqu’on passe de 𝑡𝑛 à 𝑡𝑛+1 ?

Page 15
 Relier alors ce résultat aux données du tableau 2.

Solution : L’algorithme d’Euler est défini par l’équation 13, évaluons 𝐸𝑛+1 − 𝐸𝑛 à partir de là.
1 2 1 2
𝐸𝑛+1 = 𝑘 𝑥 + 𝑣𝑛 ⋅ Δ𝑡) + 𝑚 (𝑣𝑛 + 𝑎𝑛 ⋅ Δ𝑡)
2 ( 𝑛 2
1 1
𝐸𝑛+1 = 𝑘 𝑥 2 + 2𝑥𝑛 𝑣𝑛 Δ𝑡 + 𝑣𝑛2 (Δ𝑡)2 ) + 𝑚 (𝑣𝑛2 + 2𝑣𝑛 𝑎𝑛 Δ𝑡 + 𝑎𝑛2 (Δ𝑡)2 )
2 ( 𝑛 2
1 1
𝐸𝑛 = 𝑘𝑥𝑛2 + 𝑚𝑣𝑛2
2 2
1 1
𝐸𝑛+1 − 𝐸𝑛 = 𝑘 (2𝑥𝑛 𝑣𝑛 Δ𝑡 + 𝑣𝑛2 (Δ𝑡)2 ) + 𝑚 (2𝑣𝑛 𝑎𝑛 Δ𝑡 + 𝑎𝑛2 (Δ𝑡)2 )
2 2
1 2 1 2 2
𝐸𝑛+1 − 𝐸𝑛 = 𝑘𝑥
( 𝑛 𝑛𝑣 + 𝑚𝑣 𝑎
𝑛 𝑛) Δ𝑡 + ( 2 𝑛 + 2 𝑘𝑎𝑛 ) (Δ𝑡)
𝑘𝑣

Le premier terme 𝑃𝑛 = 𝑘𝑥𝑛 𝑣𝑛 + 𝑚𝑣𝑛 𝑎𝑛 ) correspond à ce qu’on attend pour l’évaluation des puissances dans
l’algorithme d’Euler, la méthode utilisée ici s’écarte de ce calcul en ajoutant un terme toujours positif
1 2 1 2 2
( 2 𝑘𝑣𝑛 + 2 𝑘𝑎𝑛 ) (Δ𝑡) sur chaque pas de temps.
Si on regarde l’effet sur une intégration sur une durée finie 𝑡𝑚𝑎𝑥 divisée en 𝑁 pas, l’effet cumulatif
𝑡 2
donne une erreur en 𝑁 × 𝑐 ( 𝑚𝑎𝑥 ) , on retrouve une erreur cumulée proportionnelle à 1/𝑁, le coefficient de
𝑁
proportionnalité 𝑐 correspondant à la moyenne d’une somme de termes tous positifs n’ayant aucune raison
ni de tendre vers 0 ni de diverger.

Une propriété importante que devrait vérifier un algorithme d’intégration est la réversibilité dans le temps :
en partant des positions et vitesses d’un temps 𝑡 + ∆𝑡 et en appliquant un pas de temps −∆𝑡, un algorithme
réversible en temps devrait redonner les positions et vitesses du temps 𝑡.
En pratique, en partant d’un couple (𝑥𝑛 ; 𝑣𝑛 ), on applique donc tout d’abord un pas ∆𝑡 pour déterminer
(𝑥𝑛+1 ; 𝑣𝑛+1 ), puis on applique un pas −∆𝑡 afin d’obtenir (𝑥𝑛̅ ; 𝑣𝑛̅ ). Si (𝑥𝑛̅ ; 𝑣𝑛̅ ) = (𝑥𝑛 ; 𝑣𝑛 ), alors l’algorithme est
dit réversible en temps.
Q39 —Pourquoi s’agit-il d’une propriété importante à vérifier pour le problème considéré ?

Solution : Les équations de la mécanique sont réversibles, c’est fondamentalement bien que les méthodes
numériques d’intégration vérifient cette propriété.
Pour un système dont le mouvement est périodique, on espère se retrouver après intégration sur une
période dans le même état. Sur chaque demi-période successive dans le problème étudié, on peut considérer
qu’on reproduit le même mouvement en inversant le sens du temps. Si la méthode est réversible, les erreurs
accumulées sur une demi-période pourront annuler celles accumulées sur la demi-période suivante.

Q40 —Donner l’expression de 𝑥𝑛̅ en fonction de 𝑥𝑛 pour les deux algorithmes considérés. Peut-on déjà
conclure sur lsilité en stemps de chacun de ces algorithmes ?

Solution : Pour Euler :

+Δ𝑡 −Δ𝑡
𝑥𝑛 −−−−−−−−−→ 𝑥𝑛+1 = 𝑥𝑛 + 𝑣𝑛 Δ𝑡 −−−−−−−−−→ 𝑥𝑛̅ = 𝑥𝑛+1 − 𝑣𝑛+1 Δ𝑡 = 𝑥𝑛 + (𝑣𝑛 − 𝑣𝑛+1 )Δ𝑡

Pour Leapfrog :
+Δ𝑡 𝑣𝑛−1/2 + 𝑣𝑛+1/2 −Δ𝑡 𝑣𝑛+1/2 + 𝑣𝑛−1/2
𝑥𝑛 −−−−−−−−−→ 𝑥𝑛+1 = 𝑥𝑛 + Δ𝑡 −−−−−−−−−→ 𝑥𝑛̅ = 𝑥𝑛+1 − Δ𝑡 = 𝑥𝑛
2 2
L’algorithme d’Euler n’est pas réversible puisque les deux valeurs ne sont pas identiques. Pour l’algorithme
de Leapfrog, la vérification nécessaire pour les positions est validée, reste à vérifier que cet algorithme est
également réversible pour les évaluations des vitesses (il faut que positions et vitesses soient les mêmes
pour se retrouver dans le même état physique, au même point de l’espace des phases.)

Q41 —Donner alors l’expression de 𝑣𝑛̅ en fonction de 𝑣𝑛 lorsque nécessaire et conclure sur la réversibilité en
temps.

Page 16
 Solution : On doit traiter l’algorithme de Leapfrog pour les vitesses :
+Δ𝑡 −Δ𝑡
𝑣𝑛−1/2 −−−−−−−−−→ 𝑣𝑛+1/2 = 𝑣𝑛+1/2 + 𝑎𝑛 Δ𝑡 −−−−−−−−−→ 𝑣𝑛−1/2
̅ = 𝑣𝑛+1/2 − 𝑎𝑛 Δ𝑡 = 𝑣𝑛−1/2

La réversibilité est donc vérifiée pour les vitesses également, il l’est donc pour un état du système : l’algo-
rithme de Leapfrog est réversible.

Q42 —On s’intéresse enfin à une simulation de plus grande durée (𝑡𝑚𝑎𝑥 = 60 s). La figure 5 donne l’évolution
de l’erreur absolue sur 𝑥 entre données numériques et analytiques (|𝑒𝑥 |) pour les deux algorithmes choisis. Que
peut-on mettre en évidence sur cette figure ?

0,0015 Euler
Leapfrog

0,0010
| ex|

0,0005

0,0000
0 10 20 30 40 50
t( s)

Figure 5 : Erreurs absolues calculées sur la position de la masse à ∆𝑡 = 0,001 s, pour les deux algorithmes
considérés.

Solution : Vu la régularité de la figure, on peut penser que le phénomène simuler aboutit à un régime
périodique, très probablement sinusoïdal. Pour l’algorithme de Leapfrog, il semble qu’après chaque pé-
riode, le système se retrouve bien dans son état initial, et qu’en tous cas l’erreur sur une période mesurée
sur les différentes périodes est de moyenne nulle : l’erreur reste contenue. Pour l’algorithme d’Euler, sur
une oscillation l’énergie calculée augmente systématiquement : à chaque passage au niveau de la position
d’équilibre, la vitesse est un peu plus grande qu’attendue. Ces erreurs systématiques s’additionnent période
après période, avec une condition initiale de plus en plus dégradée au départ de la période, la simulation
est de plus en plus mauvaise.

Q43 —Conclure sur les avantages et les inconvénients de chacun de ces algorithmes et sur leur adéquation
pour le traitement numérique de ce problème dans le cas où on envisage une simulation de plusieurs eures.

Solution : L’algorithme d’Euler a pour lui sa simplicité qui le rend très lisible, très simple à implémenter
car la démarche transparait dans chacune des relations de l’implémentation.
L’algorithme d’Euler est donc particulièrement pour une approche didactique de l’intégration numé-
rique.
Les algorithmes plus évolués comme celui de Leapfrog peuvent être beaucoup plus précis, au prix de
quelques calculs préliminaires et d’une ligne de calcul un peu moins lisible pour la mise en place de la
relation de récurrence.
Pour un travail de simulation demandant de la précision et limité par les ressources disponibles, les
algorithmes plus évolués sont à préférer.

Page 17
ANNEXE : Aide-mémoire sur numpy
Les bibliothèques sont importées de la façon suivante :

>>> from math import *

>>> import numpy as np

La création d’un tableau numpy tab à une dimension possédant 𝑛 éléments, tous initialisés à 0, est réalisée
à l’aide de l’instruction :

>>> n=5
>>> tab =np.zeros ( n )

Celle d’un tableau numpy tab à une dimension possédant 𝑛 éléments, uniformément répartis entre deux
valeurs debut et fin, se fait avec :

>>> debut = 0; fin = 10; n = 5

>>> tab=np.linspace(debut,fin,n)
>>> print(tab)
[ 0. 2.5 5. 7.5 10. ]

L’accès à un élément du tableau tab (en lecture ou en écriture) se fait par tab[i], la numérotation des
indices se faisant à partir de 0 :

>>> tab = np.zeros ( 4 )

>>> tab
array([0., 0., 0., 0.])
>>> tab[1]=2; tab[2]=6
>>> print(tab)
[0. 2. 6. 0.]

La sélection de l’ensemble des j premiers éléments du tableau tab est possible avec :

>>> j=3; print(tab[:j])

[0. 2. 6.]

Le maximum des éléments d’un tableau tab s’obtient avec :

>>> np.max ( tab )

6.0

FIN

Page 18