Simulation du comportement chaotique d'une boussole

La physique de l'expérience

Le principe de l'expérience est assez simple : exposer une boussole à deux champs magnétiques - un champ fixe (le champ magnétique terrestre par exemple) et un champ tournant. Nous examinerons les mouvements de la boussole lorsque elle est exposée à ces deux champs en fonction des caractéristiques de ces champs.
Pratiquement, cette expérience peut être réalisée en plaçant la boussole entre deux bobines de Helmholtz, que l'on alimentera avec un courant alternatif.
L'aiguille de la boussole est placée sur un axe de rotation vertical. Elle présente un moment d'inertie J et un moment magnétique M. On repère son mouvement par l'angle θ formé entre son axe longitudinal et un axe Ox orienté vers le Nord. Elle est soumise à un couple de freinage égal à K*dθ/dt.

Le champ magnétique fixe est de module B0 constant. Le champ magnétique tournant est de module B1, de vitesse de rotation ω.
Etablissons l'équation différentielle de mouvement de l'aiguille de la boussole, ce qui est pratiquement immédiat en appliquant le théorème du moment cinétique:

J*d²θ/dt² = -K*dθ/dt - M*B0*sin(θ) - M*B1*sin(θ - ωt)

Pour faciliter le calcul, je vais procéder à une opération classique de changement de temps, afin que la pulsation du champ tournant soit unitaire et adimensionnée. Il faut donc poser τ = 2*π*t/T, τ étant sans dimension.
Si l'on réécrit l'équation différentielle en tenant compte de ce changement de variable et en remarquant que dθ/dt = dθ/dτ*dτ/dt = dθ/dτω, on obtient:

J*ω²*d²θ/dτ² = -K*ω*dθ/dτ - M*B0*sin(θ) - M*B1*sin(θ - τ)

En regroupant les constantes j'obtiens finalement l'équation:

d²θ/dτ² = - α*dθ/dτ - β*sin(θ) - γ*sin(θ - τ)

avec α = K/Jωβ = M*B0/J*ω² et γ = M*B1/J*ω². Remarquons que ces paramètres sont sans dimension.

Dans la suite, nous nous intéresserons aux trajectoires de cette équation dans l'espace de phase (θ,dθ/dτ).

Etudions d'abord rapidement la forme de cette équation:

nous constatons que si α >= 0 et γ = 0, nous retrouvons l'équation différentielle d'un oscillateur harmonique de pulsation ω. Ce cas ne nous intéresse pas vraiment ici sauf à visualiser le portait de phase d'un oscillateur harmonique.
si α = 0 alors la solution de l'équation est une sinusoïde de la forme θ = θ₀*cos(ωτ + φ), avec dθ/dτ = -ω*q₀*sin(ωτ + φ). La trajectoire dans l'espace de phase est une ellipse (je vous laisse le démontrer). Nous pourrons étudier l'influence de la valeur de β et γ sur la forme de l'ellipse.
si α > 0 la solution de l'équation est de la forme θ = θ₀*exp(-Kτ)*cos(ωτ+φ), avec dθdτ = -ω*θ₀*exp(-Kτ)*sin(ωτ + φ) - K*q₀*exp(-Kτ)*cos(ωτ + φ).

Voilà pour la physique du problème. Nous allons maintenant écrire un programme très simple qui nous permettra de tracer dans l'espace de phase la trajectoire d'une solution de l'équation pour différentes valeurs des paramètres α, β et γ, et aussi pour différentes valeurs des conditions initiales, en particulier la vitesse angulaire initiale.

Le programme de simulation

Le code source du programme est assez simple. Il comporte quelques routines de services graphiques et les routines utilisées pour RK4. Ce schéma numérique est expliqué ici. Dans ce programme, la routine a été adaptée en C. Le corps du programme principal est trivial : saisie des paramètres, configuration graphique, boucle de calcul RK4 et affichage des résultats.
Je vous indique ci-dessous le code de la routine Derivee, qui est la seule qui soit spécifique à ce programme:

//******************************************************************************

// Définition du système différentiel

//******************************************************************************

void Derivee(double Y[], double dY[])

{

dY[0] = Y[1]; // y = d(theta)/dt

dY[1] = -alpha*Y[1]-beta*sin(Y[0])-gamma*sin(Y[2]); // d²(theta)/dt²

dY[2] = -1; // dz/dt = d(theta-t)/dt

return;

}

Le source du programme est disponible ici. Pour compiler et linker le programme avec dislin, notez les points suivants:

Ne pas oublier la déclaration #include <dislin_d.h>, car j'utilise des doubles. En cas d'utilisation de float, inclure dislin.h . Attention, en C++, cette inclusion doit être faite après l'inclusion de <iostream.h>.
Aller dans Project/Project Options/Parameters/Linker et ajouter les commandes -ldismg_d -lgdi32. Attention: si vous utilisez des float, il faut ajouter -ldismg

Ce faisant, la compilation et le link devraient se passer sans problème! Si il y a un problème, revoir votre installation ou/et le type de variables utilisées.

Expériences et résultats

La découverte de la notion d'attracteur

Dans cette première expérience, nous allons fixer les paramètres α, β et γ (respectivement 0.5, 1.5 et 1.0) et faire varier la vitesse initiale dθ/dτ . Les figures ci-dessous montrent la forme de la trajectoire pour chaque valeur de dθ/dτ, variant entre 0.0 et 5.0.

Il est remarquable que malgré la forme initiale de la trajectoire différente, due à la vitesse initiale différente, la trajectoire tende dans tous les cas vers un cycle limite identique. Ce cycle met en évidence une solution stable de l'équation différentielle. Ce cycle limite est un attracteur.
A titre d'expérience, vous pouvez faire varier dθ/dτ entre 1.0 et 2.0 par pas de 0.1 pour visualiser les différentes formes de trajectoires initiales et la "chute" vers l'attracteur.

Le portrait de phase d'un oscillateur harmonique

Analyse pour un couple de freinage et un champ tournant nuls

Nous avons vu plus haut que dans le cas où α= 0 et γ = 0, nous avions affaire à un oscillateur harmonique de pulsation ω. Il s'agit du cas où le couple de freinage est nul (K = 0) et où le champ magnétique tournant est nul.
Examinons son portait de phase en fonction de différentes valeurs de β en fixant dθ/dτ à 1.0. Par exemple pour les valeurs 1, 10, 50 et 100:

Vous remarquez que l'ellipse s'aplatit lorsque β augmente. Rappelons que le coefficient sans dimension β est égal à M*B0/J*ω² . Les grandeurs M, B0 et J étant des constantes, on note donc que dans le portait de phase, la trajectoire ellipsoïdale d'un oscillateur harmonique s'aplatit lorsque sa pulsation ω diminue.

Analyse du cas α = 0 et des variations relatives de β et γ

Conservons une valeur nulle pour α et notre vitesse initiale égale à 1. Nous allons nous intéresser à la variation relation de β et γ. Le coefficient sans dimension β est égal à M*B0/J*ω². Le coefficient sans dimension γ est égal à M*B1/J*ω². Il s'agit donc d'étudier le comportement des trajectoires selon le rapport des intensités du champ fixe par rapport au champ tournant.
Si les champs ont la même intensité (β=1,γ=1) voilà la trajectoire obtenue:

La trajectoire est chaotique.
Voyons l'évolution de la trajectoire en augmentant le rapport β/γ . Les valeurs des paramètres figurent sur les images:

On remarque nettement le passage d'un désordre vers un ordre. A partir d'un rapport donné , autour de 1.7/1.0, la trajectoire tend vers un attracteur ellipsoïdal. Cette tendance est confirmée par les vues suivantes.
Vous pouvez faire varier le rapport β/γ , on obtient des figures très intéressantes. Essayez par exemple autour de 100/50 et 100/70...

Recherche du chaos

Je vous invite maintenant à utiliser le programme "Boussole" pour explorer le chaos. Commencez par explorer le domaine autour des valeurs faibles d'α (disons entre 0.01 et 0.5) et des valeurs similaires de β et γ (en commençant par β = 1 et γ =1) . Conservez la valeur de dθ/dτ à 1.Vous devriez obtenir une image de ce genre:

Voici, par exemple, la trajectoire obtenue avec α=0.0005, β=2, γ=2.5:

Vous pouvez maintenant fixer le jeu des 3 paramètres physiques et faire varier la vitesse initiale. Je vous laisse jouer à votre tour!