Boussole chaotique

1. Equation différentielle du mouvement

Une aiguille aimantée de boussole, de moment magnétique M, disposée horizontalement, est libre de tourner sur un pivot vertical. Elle est soumise à l'action de deux champs magnétiques : un champ B₁ constant et un champ B₂ tournant à la pulsation Ω. Si Oxy est un repère centré sur la boussole dont les deux axes sont horizontaux, le champ magnétique total s'écrit :

\vec{B} = (B_{1} + B_{2} cos (Ω t)) {\vec{e}}_{x} + B_{2} sin (Ω t) {\vec{e}}_{y}

Le moment des forces magnétiques par rapport au point O (situé sur l'axe de rotation) est :

\vec{Γ} = \vec{M} \land \vec{B}

Si on note θ l'angle que fait la boussole avec le champ fixe (axe Ox) et J son moment d'inertie par rapport à l'axe de rotation, l'équation différentielle du mouvement est :

J \frac{d^{2} θ}{d t^{2}} + B_{1} M sin θ + B_{2} M sin (θ - Ω t) = 0

Lorsque B2=0 (pas de champ tournant), on retrouve l'équation du pendule conservatif. Dans ce cas, les petites oscillations autour de 0 se font à la pulsation :

ω_{1} = \sqrt{\frac{B_{1} M}{J}}

Lorsque l'amplitude des oscillations augmente, leur pulsation diminue et des harmoniques apparaissent dans θ(t).

Lorsque B1=0 (pas de champ fixe), le changement de variable

θ' = θ - Ω t

permet aussi de se ramener à l'équation du pendule. Dans ce cas, il y a oscillation de la boussole autour du champ magnétique tournant. Si ces oscillations sont de petite amplitude, la pulsation est

ω_{2} = \sqrt{\frac{B_{2} M}{J}}

2. Section de Poincaré

En introduisant la vitesse angulaire de la boussole u et l'angle ϕ entre le champ tournant et l'axe Ox, on obtient le système différentiel du premier ordre :

\begin{matrix} \frac{d θ}{d t} = u \\ \frac{d u}{d t} = - ω_{1}^{2} sin θ - ω_{2}^{2} sin (θ - φ) \\ \frac{d φ}{d t} = Ω \end{matrix}

On a donc un système non linéaire à 3 degrés de libertés (θ,u,ϕ). L'espace des phases a trois dimensions. Pour obtenir une représentation à 2 dimensions, on effectue une section de Poincaré par des plans perpendiculaires à l'axe de ϕ. Autrement dit, on trace les points de coordonnées (θ,u) obtenus lorsque ϕ est multiple entier de 2π. Pour le calcul numérique, on posera :

Ω = 2 π

De cette manière, les points de la section de Poincaré sont obtenus pour des valeurs entières du temps.

On posera également :

\begin{matrix} ω_{1} = 2 π a_{1} \\ ω_{2} = 2 π a_{2} \end{matrix}

Le paramètre de stochasticité est défini par :

s = \frac{2 (ω_{1} + ω_{2})}{Ω} = 2 (a_{1} + a_{2})

3. Fonctions Scilab

3.a. Module de calcul

Module scilab : télécharger.

La fonction ode de scilab utilise le solveur LSODE (Livermore Solver for Ordinary Differentiel Equations). Ce solveur utilise par défaut une méthode prédicteur-correcteur (Adams-Moulton), puis bascule sur la méthode BDF (backward differentiation formula) si le système devient raide. On doit fournir aussi le jacobien du système.

La fonction suivante effectue l'intégration numérique pour une condition initiale donnée. Le temps final et la période d'échantillonnage sont précisés en argument.

function [t,y]=solution(a1,a2,thetaInit,uInit,tmax,te),
    omega1=2*%pi*a1;
    omega2=2*%pi*a2;
    function deriv=systeme(t,y),
        deriv(1)=y(2),
        deriv(2)=-omega1^2*sin(y(1))-omega2^2*sin(y(1)-2*%pi*t),
    endfunction
    function J=jacobien(t,y), 
        J(1,1)=0,
        J(1,2)=1, 
        J(2,1)=-omega1^2*cos(y(1))-omega2^2*cos(y(1)-2*%pi*t), 
        J(2,2)=0,
    endfunction  
    t=[0:te:tmax];
    tolA=1d-7;
    tolR=1d-10;
    y=ode([thetaInit;uInit],0,t,tolR,tolA,systeme,jacobien);
endfunction

La fonction suivante trace la section de Poincaré pour une liste de conditions initiale fournie dans un tableau (une ligne pour chaque condition initiale).

 
function poincare(a1,a2,initial,tmax,rect),
    s=size(initial);
    for k=1:s(1),
        [t,y]=solution(a1,a2,initial(k,1),initial(k,2),tmax,1);
        plot2d(y(1,:),y(2,:),style=0,rect=rect);
    end
endfunction

La fonction suivante trace le spectre en décibel. La période d'échantillonnage sera choisie en fonction de la fréquence maximale désirée.

function spectre(a1,a2,thetaInit,uInit,tmax,te,dbmin,dbmax),
    [t,y]=solution(a1,a2,thetaInit,uInit,tmax,te);
    tfd=fft(y(1,:));
    n=length(y(1,:));
    sp=zeros(1,n);
    f=zeros(1,n);
    for k=1:n,
        sp(k)=20*log10(abs(tfd(k))); 
        f(k)=(k-1)/tmax;
    end; 
    plot2d(f,sp,rect=[0,dbmin,1/te/2,dbmax]); 
endfunction

3.b. Exemple

exec('boussole.sci');
a1=0.5;
a2=0.2;
tmax=500;
initial=[0.1,0;0.5,0;1,0;2,0];

plotA=scf();
poincare(a1,a2,initial,tmax,[-3,-3,3,3])

plotB=scf();
spectre(a1,a2,0.5,0,tmax,0.05,-100,100)

plotB.pdf

4. Package DYNODE pour Python

Le Package DYNODE pour Python permet d'utiliser le solveur SUNDIALS CVODE (méthode à pas multiples), ou un solveur utilisant la méthode de Stormer-Verlet.

On commence par ouvrir un solveur CVODE pour l'équation de la boussole :

import dynode.main as dyn
import numpy
from pylab import *
import math
solver=dyn.CVOde(dyn.OdeBoussole,dyn.OdeAdams,dyn.OdeFunctional)

On fixe les valeurs des constantes :

a1=0.2
a2=0.15
solver.set_cst([(2*math.pi*a1)**2,(2*math.pi*a2)**2,0])

On fixe plusieurs conditions initiales et on effectue le calcul pour chacune d'elles:

reltol=1e-5
abstol=1e-6
y0=[[0.1,0],[0.5,0],[1,0],[2,0]]
data=[]
for k in range(len(y0)):
    solver.init(0.0,y0[k],reltol,[abstol])
    data.append(solver.solve(1,500))

Voici la section de Poincaré :

figure(0)
for k in range(len(y0)):
    plot(data[k][1],data[k][2],marker=".",linestyle=" ",markersize=2)
xlabel('theta')
ylabel('dtheta')

plotC.pdf

Analyse spectrale de l'évolution temporelle de l'angle :

import numpy.fft
solver.init(0.0,[1,0],reltol,[abstol])
tmax=100
data=solver.solve(0.1,tmax)
tfd=numpy.fft.fft(data[1])
n=size(tfd)
f=numpy.zeros(n)
dB=numpy.zeros(n)
for k in range(n):
    f[k] = (k-1)*1.0/tmax
    dB[k] = 20*math.log10(abs(tfd[k]))
figure(1)
plot(f,dB)
xlabel('f')
ylabel('dB')
axis([0,3,-10,60])

plotD.pdf

solver.close()

Boussole chaotique : définition du problème