Retour
Dominique Lefebvre octobre 2009
Retour
Dominique Lefebvre octobre 2009
Simulation du comportement chaotique d'une boussole
La physique de l'expérience |
Expériences
et résultats |
| Le
programme de simulation |
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 q formé entre son axe longitudinal et un axe Ox orienté vers le Nord. Elle est soumise à un couple de freinage égal à K*dq/dt.
Le champ magnétique fixe est de module B0 constant. Le champ magnétique tournant est de module B1, de vitesse de rotation w.
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²q/dt² = -K*dq/dt - M*B0*sin(q) - M*B1*sin(q - wt)
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 t = 2*p*t/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 dq/dt = dq/dt*dt/dt = dq/dt*w, on obtient:
J*w²*d²q/dt² = -K*w*dq/dt - M*B0*sin(q) - M*B1*sin(q - t)
En regroupant les constantes j'obtiens finalement l'équation:
d²q/dt² = -a*dq/dt - b*sin(q) - g*sin(q - t)
avec a = K/J*w b = M*B0/J*w² et g = M*B1/J*w² . 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 (q,dq/dt).
Etudions d'abord rapidement la forme de cette équation:
nous constatons que si a=0 et g=0, nous retrouvons l'équation différentielle d'un oscillateur harmonique de pulsation w. Ce cas ne nous intéresse pas vraiment ici, sauf à visualiser le portait de phase d'un oscillateur harmonique.
si a=0 alors la solution de l'équation est une sinusoïde de la forme q = q0*cos(wt +f), avec dq/dt = -w*q0*sin(wt +f). 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 b et g sur la forme de l'ellipse.
si a>0 la solution de l'équation est de la forme q = q0*exp(-Kt)*cos(wt +f), avec dq/dt = -w*q0*exp(-Kt)*sin(wt +f) - K*q0*exp(-Kt)*cos(wt +f).
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 a, b et g, et aussi pour différentes valeurs des conditions initiales, en particulier la vitesse angulaire initiale.
HOME
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.
//******************************************************************************
// Programme de simulation du comportement d'une boussole dans un champ
// magnétique tournant - Illustration du chaos
//
// Dominique Lefebvre - TangenteX.com Octobre 2009
//******************************************************************************
#include <cstdlib>
#include <iostream>
#include <math.h>
#include <dislin_d.h>
using namespace std;
//******************************************************************************
// Définition des constantes
//******************************************************************************
#define TMAX 50 // durée max de la manip
#define PAS 0.01 // pas de calcul
#define VMAX 5.0 // vitesse angulaire
#define N 3 // nombre d'équations du système différentiel
#define PI 3.1415927
#define DEUXPI 6.2831853
//******************************************************************************
// Définition des variables globales
//******************************************************************************
double alpha, beta, gamma, dt;
//******************************************************************************
// Définition des subroutines de service
//******************************************************************************
void plot_d(double x, double y)
{
rlstrt(x,y);
rlconn(x,y);
return;
}
void InitGraphDislin(void)
{
metafl("XWIN");
disini();
pagera();
pagfll(255);
color("black");
winfnt("arial");
return;
}
//******************************************************************************
// 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;
}
//******************************************************************************
// Routine RK4 - tirée des Numerical Recipes
//******************************************************************************
void RK4(double yin[], double yout[], double h)
{
// La constante N, qui indique le nombre d'équations du système différentiel
// est définie globalement pour tout le programme.
int i;
double dydx[N], dydxt[N], yt[N], k1[N], k2[N], k3[N], k4[N];
Derivee(yin, dydx);
for (i = 0; i < N; i++)
{
k1[i] = h*dydx[i];
yt[i] = yin[i] + 0.5*k1[i];
}
Derivee(yt, dydxt);
for (i = 0; i < N; i++)
{
k2[i] = h*dydxt[i];
yt[i] = yin[i] + 0.5*k2[i];
}
Derivee(yt, dydxt);
for (i = 0; i < N; i++)
{
k3[i] = h*dydxt[i];
yt[i] = yin[i] + k3[i];
}
Derivee(yt, dydxt);
for (i = 0; i < N; i++)
{
k4[i] = h*dydxt[i];
yout[i] = yin[i] + k1[i]/6. + k2[i]/3. + k3[i]/3. + k4[i]/6.;
}
return;
}
//******************************************************************************
// Programme principal
//******************************************************************************
int main(int argc, char *argv[])
{
double theta0, v0, t0, theta, v, t, dt;
double Yin[N], Yout[N];
char buffer[50];
// Déclaration des conditions initiales statiques
t0 = 0.0;
theta0 = 0.0;
// Saisie des paramètres de calcul v0, alpha, beta, gamma
printf("Vitesse angulaire initiale: ");scanf("%lf",&v0);
printf("Valeur de alpha: "); scanf("%lf",&alpha);
printf("Valeur de beta: ");scanf("%lf",&beta);
printf("Valeur de gamma: "); scanf("%lf",&gamma);
// Initialisation des variables
t = t0;
dt = PAS;
theta = theta0;
v = v0;
// Initialisation de la librairie graphique
InitGraphDislin();
// Tracé des axes et des légendes
wintit("Boussole - Espace de phase");
axstyp("CROSS");
name("d(theta)/dt","X"); namjus("RIGHT","X");
name("theta","Y");namjus("RIGHT","Y");
labels("NONE","X");
labels("NONE","Y");
graf(-VMAX,VMAX,-VMAX,1,-PI,PI,-PI,PI/10);
titlin("Diagramme de l'espace de phase (v, theta)",1);
title();
// Affichage des paramètres
sprintf(buffer,"V0: %.2lf",v0);
messag(buffer,2500,1660);
sprintf(buffer,"Alpha: %.2lf",alpha);
messag(buffer,2500,1710);
sprintf(buffer,"Beta : %.2lf",beta);
messag(buffer,2500,1760);
sprintf(buffer,"Gamma: %.2lf",gamma);
messag(buffer,2500,1810);
// tracé des points en rouge
color("red");
// positionnement aux valeurs initiales
rlstrt(v,theta);
// boucle de calcul et d'affichage de la trajectoire par RK4
do
{
// Calcul par RK4
Yin[0] = theta;
Yin[1] = v;
Yin[2] = theta - t;
RK4(Yin, Yout, dt);
theta = Yout[0];
v = Yout[1];
// increment du temps
t += dt;
// contraintes de limitation au domaine [-PI, PI]
if (theta > PI)
theta = theta - DEUXPI;
else if (theta < -PI)
theta = theta + DEUXPI;
// plot du point (v,theta)
plot_d(v,theta);
}
while (t <= TMAX);
// Fin du programme
disfin();
return EXIT_SUCCESS;
}
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.
En cas de question, mon forum
HOME
La découverte de la notion d'attracteur
Dans cette première expérience, nous allons fixer les paramètres a, b et g (respectivement 0.5, 1.5 et 1.0) et faire varier la vitesse initiale dq/dt . Les figures ci-dessous montrent la forme de la trajectoire pour chaque valeur de dq/dt, variant entr 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 dq/dt 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
Nous avons vu plus haut que dans le cas où a = 0 et g = 0, nous avions affaire à un oscillateur harmonique de pulsation w. 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 b en fixant dq/dt à 1.0. Par exemple pour les valeurs 1, 10, 50 et 100:
Vous remarquez que l'ellipse s'aplatit lorsque b augmente. Rappelons que le coefficient sans dimension b est égal à M*B0/J*w² . Les grandeurs M, B0 et J étant des constantes, on note donc que dans le portait de phase, la trajection ellipsoïdale d'un oscillateur harmonique s'aplatit lorsque sa pulsation w diminue.
Analyse du cas a = 0 et des variations relatives de b et g
Conservons une valeur nulle pour a et notre vitesse initiale égale à 1. Nous allons nous intéresser à la variation relation de b et g. Le coefficient sans dimension b est égal à M*B0/J*w² . le coefficient sans dimension g est égal à M*B1/J*w². 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é (b=1, g=1) voilà la trajectoire obtenue:
La trajectoire est chaotique.
Voyons l'évolution de la trajectoire en augmentant le rapport b/g . 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 b/g , 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'a (disons entre 0.01 et 0.5) et des valeurs similaires de b et g (en commençant par b = 1 et g =1) . Conservez la valeur de dq/dt à 1.Vous devriez obtenir une image de ce genre:
Voici, par exemple, la trajectoire obtenue avec a=0.0005, b=2, g=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!
HOME