Maple - le pendule simple

Pourquoi revenir sur le pendule simple

Nous avons déjà abordé la physique du pendule simple à l'aide de SciLab, alors pourquoi revenir sur le sujet ? Parce que j'avais envie d'y introduire la notion d'analyse spectrale, qui est maintenant abordée en TS, mais aussi parce que j'ai eu aussi envie de l'aborder avec Maple plutôt que SciLab. Ce n'est pas un changement d'allégence (je n'ai d'allègence à aucun langage...) mais plutôt la prise en compte des besoins des élèves de prépa, qui étudieront Maple au cours de leurs études. Et puis, à la pratique, je constate que Maple est bien puissant! Pas adapté à tout mais sympa quand même...

Avouons aussi que le pendule simple, et plus généralement l'oscillateur harmonique est un sujet aussi étendu qu'intéressant et que j'aurais sans doute l'occasion d'y revenir, avec Maple, SciLab ou Python !

La modélisation du pendule simple

Je rappelle brièvement l'équation diffréentiell de second ordre non linéaire que le PFD permet d'établir et qui décrit le mouvement du pendule simple:

d²θ/dt² + ω₀²sinθ = 0

avec ω₀ = sqrt(g/l).
En TS, on étudie l'approximation harmonique de ce mouvement, en posant sinθ ≈ θ. Dans ce cas, la solution est simple: on obtient une sinusoïde de pulsation ω₀ et de période T₀ = 2*π/ ω₀. Mais vous êtes tous prévenus des limites de cette approximation: au délà d'une certaine valeur d'angle initial, le mouvement réel différe totalement du mouvement calculé. Le modèle ne rend plus compte de la réalité. Je vous propose donc d'aborder ici quelques aspects du mouvement réel du pendule simple, sachant que je néglige tout frottement et donc que le qualificatif de "réel" est un peu exagéré...

Nous allons commencer par adimensionner cette équation pour la rendre plus "lisible". Je poserai donc τ = t/T₀ et y(τ) = θ(T₀τ), ce qui nous donne l'équation:

d²y/dτ² + 4*π²sin(y) = 0

L'unité de temps est donc devenue T₀. Nous travaillerons dans la suite avec cette équation, en gardant toujours présent à l'esprit l'unité de temps.

Le programme Maple

Vous pouvez télécharger le code du programme FFTPendule.mpl. Je vous recommande d'ailleurs de le faire tout de suite, d'ouvrir Maple et de charger le code (read("FFTpendule.mpl")) pour pouvoir faire les manips qui vont suivre.
Le code présente plusieurs particularités:

la résolution par dsolve d'une EDO avec des paramètres, qui me permet de modifier facilement les conditions initiales de la trajectoire du pendule.
l'utilisation de la fonction DEplot (package DEtools) pour tracer simplement le portrait de phase d'un système différentiel.
l'utilisation de la fonction FFT pour calculer la transformée de Fourier discrète du "signal" du pendule.

Nous verrons dans la suite comment fonctionnent ces différents composants.

Trajectoire et approximation harmonique

A l'aide du programme, nous allons tracer, pour plusieurs valeurs de l'angle initial, la trajectoire du pendule sur un domaine temporel de 0 à 5 unités (5 T₀ donc) . Pour pouvoir visualiser l'écart entre l'approximation harmonique et la trajectoire réelle, j'ai superposé la trajectoire calculée en bleu, et l'approximation harmonique en rouge. Pour modifier la valeur de l'angle initial, vous devez modifier la ligne:

s(parameters = [1.0,0.0]):

Dans cet exemple, l'angle initial vaut 1. rd et la vitesse initiale 0. rd.s^-1. Lançons le programme avec cette valeur d'angle initial. Nous obtenons la courbe suivante:

Pour un angle d'un radian, nous observons la dérive très nette de l'approximation harmonique (en rouge) par rapport à la trajectoire calculée avec l'équation horaire complète. Vous noterez également que la dérive est cumulative,ce qui est due aux méthodes de résolution des EDO. Pour information, Maple utilise ici une méthode RK4 un peu améliorée.
Voyons maintenant ce que cela donne avec un angle plus grand, pour lequel l'approximation harmonique n'est plus acceptable. Pour cela, nous allons modifier la valeur de l'angle dans la ligne, en indiquant 3.1 rd pour la valeur de l'angle initial (proche de π) :

s(parameters = [3.1,0.0]):

Nous obtenons la trajectoire suivante:

Elle n'a plus rien de sinusoïdale ! D'après le tracé, sa période est d'environ 3,3 unités de temps. Vous comprenez mieux pourquoi il ne vaut mieux pas employer l'approximation harmonique (dite aussi "des petits angles"...) pour des angles qui dépassent le radian! Et même pour un radian, on a vu que ce n'était pas terrible.

Portrait de phase

Le portait de phase, incontournable dans l'étude des systèmes dynamiques est la courbe paramétrique x= y(t) et y = dy(t)/dt, les grandeurs y(t) et dy(t)/dt étant fournies par la résolution de l'équation différentielle de trajectoire.
Maple dispose d'une fonction puissante pour traiter le problème, qui nous évite bien des calculs. Il s'agit de DEplot, que j'utilise comme suit:

DEplot(sys,[y(t),v(t)],t=0..2.,

[[y(0) = 0, v(0)= 1.], #description des différentes conditions initiales

[y(0) = 0, v(0)= 1.5],

[y(0) = 0, v(0)= 1.99],

[y(0) = 0, v(0)= 2.01],

[y(0) = 0, v(0)= 2.],

[y(0) = 0, v(0)= 2.5]],

stepsize = 0.01, linecolor = blue,thickness=1,arrows=none,

scaling=constrained,title = "Portrait de phase");

La variable sys contient la définition du système différentiel. t définit le domaine temporel de tracé. Les lignes suivantes définissent les différentes valeurs de couples (y₀,v₀) utilisées pour tracer le portait de phase. Vous pouvez bien sur les modifier! Voyons ce que cela donne:

Nous observons plusieurs choses intéressantes:

pour v₀ égale à 1, la courbe est pratiquement un cercle. Il s'agit d'un mouvement oscillatoire proche de l'approximation harmonique.
pour v₀ égale à 1.5, le cercle est déformé, nous nous éloignons de l'approximation harmonique.
pour v₀ autour de 2.0 mais sans prendre la valeur 2.0, les courbes prennent des formes très différentes. Pour une vitesse de 1.99, le mouvement est toujours oscillatoire car la courbe se referme. Par contre, pour 2.01, la courbe est ouverte, le pendule entre en rotation...
la valeur v₀ égale à 2 est la valeur critique de vitesse du système. La courbe semble s'arrêter au point (π,0). Le pendule arrive en position haute de sa trajectoire avec une vitesse nulle, position à l'équilibre hautement instable! Les variations de trajectoire sur le portrait de phase autour de v₀ = 2. rd.s^-1 témoignent de la sensibilité du système aux conditions initiales.
pour v₀ égale à 2.5, la courbe est ouverte. Le pendule est en rotation.

FFT et analyse spectrale

Il est intéressant de comparer le spectre des oscillations du pendule alors que sa trajectoire est approximativement harmonique et lorsqu'elle en est très éloignée. Là aussi, Maple possède des outils assez simples à mettre en oeuvre. Dans le programme, j'ai utilisé la fonction FFT de Maple, sur un vecteur de 256 points (2^8). Après avoir initialisé les paramètres et les vecteurs réels et complexes, le calcul du spectre tient en deux lignes:

evalf(FFT(p,uR,uI)):

sp := seq([[i/Tacq,0],[i/Tacq,sqrt(uR[i]^2+uI[i]^2)]],i=1..N/2):

La première ligne permet d'évaluer la transformée de Fourier discrète du signal contenu dans les vecteurs uR et uI sur p points. La seconde ligne permet de construire la séquence de points nécessaires à construire les bâtons du spectre. Notez le calcul de la magnitude de chaque bâton. Ah oui, vous aurez remarqué que dans la construction de la séquence de calcul du spectre, je n'ai pris que la première moitié des points. c'est pour éviter certains aspects déplaisants de l'algorithme de FFT, en particulier le repliement du spectre... Pour plus de détails voir votre cours de FFT.
Voici ce que donne le spectre des oscillations pour un angle initial de 1. rd, c'est à dire dans le cadre de l'approximation harmonique:

Sans surprise, nous obtenons un pic unique centré sur f = 1/T₀, ce qui vaut 1 dans notre système unité. Quand je dis unique, je néglige les petits pics parasites qui sont dus à l'algorithme de FFT.

Voyons maintenant ce que donne le spectre des grandes oscillations, pour un angle initial de 3.1 rd:

Lors de l'analyse de la trajectoire, nous avions estimé que la période était de 3.3 unité de temps, ce qui nous donne une fréquence fondamentale de 0.3 dans notre système d'unités. C'est bien pour cette valeur que nous voyons le pic principal.
On observe aussi un pic pour l'harmonique 3 et pour l'harmonique 5. L'amplitude des pics décroît fortement avec le rang de l'harmonique: 1/5 pour l'harmonique 3 et 1/10 pour l'harmonique 5.
Les harmoniques paires sont absentes. A noter les mêmes pics parasites dus à l'algorithme de FFT.

A vous de jouer maintenant.Le programme vous permet d'explorer le portait de phase et le spectre autour du point critique ou en tout autre point. Vous pouvez aussi l'améliorer pour étudier la variation de période en fonction de l'angle initial. Bon travail !