Système hamiltonien à champ magnétique

1. Équations du mouvement

On considère un point matériel de masse m en mécanique classique, soumis à un potentiel V et à une force magnétique (ou force de Coriolis). La dérivée par rapport au temps de la quantité de mouvement est :

\frac{d \vec{p}}{d t} = \vec{Ω} \land \vec{p} - \vec{g r a d} V (1)

Dans une base orthonormée, les composantes des vecteurs sont définis par les matrices colonnes suivantes (position du point et quantité de mouvement) :

\begin{matrix} q = {(q_{1}, q_{2}, q_{3})}^{T} & (2) \\ p = {(p_{1}, p_{2}, p_{3})}^{T} = {(m q_{1}, m q_{2}, m q_{3})}^{T} & (3) \end{matrix}

On définit la matrice antisymétrique :

Ω = (\begin{matrix} 0 & - ω_{3} & ω_{2} \\ ω_{3} & 0 & - ω_{1} \\ - ω_{2} & ω_{1} & 0 \end{matrix}) (4)

Lorsque le champ magnétique est non-uniforme, cette matrice est fonction de q.

Les équations du mouvement s'écrivent :

\begin{matrix} \frac{d q}{d t} = \frac{d K}{d p} = \frac{p}{m} & (5) \\ \frac{d p}{d t} = Ω p - \frac{\partial V}{\partial q} & (6) \end{matrix}

où l'on a introduit la fonction énergie cinétique K(p). Ce système admet un hamiltonien séparable :

H (q, p) = K (p) + V (q) (7)

Contrairement au système considéré pour la méthode de Störmer-verlet, celui-ci n'a pas la forme canonique.

Une écriture plus générale des équations peut être obtenue ([1]) en définissant le vecteur colonne :

z = {(q, p)}^{T} (8)

et la matrice 6x6 suivante :

J (z) = (\begin{matrix} 0 & I_{d} \\ - I_{d} & Ω (z) \end{matrix}) (9)

On obtient ainsi :

\frac{d z}{d t} = J (z) \frac{\partial H}{\partial z} (10)

Dans le cas d'un champ magnétique uniforme (ou d'une force de Coriolis dans un référentiel tournant), la matrice J est constante.

2. Solution sans potentiel dans un champ uniforme

Dans le cas d'un mouvement sans potentiel (V=0) dans un champ uniforme (Ω constante), on obtient un système linéaire à coefficients constants :

\begin{matrix} \frac{d q}{d t} = \frac{p}{m} & (11) \\ \frac{d p}{d t} = Ω p & (12) \end{matrix}

La solution pour la quantité de mouvement s'obtient en utilisant l'exponentielle d'une matrice ([2]) :

p (t) = exp (Ω t) p (0) (13)

Pour une matrice antisymétrique, l'exponentielle s'écrit ([1]) :

exp (Ω t) = I_{d} + \frac{sin (ω t)}{ω} Ω + 2 (\frac{sin (ω t / 2)}{ω})^{2} Ω^{2} (14)

avec :

ω = \sqrt{ω_{1}^{2} + ω_{2}^{2} + ω_{3}^{2}} (15)

Voici par exemple la matrice lorsque le champ magnétique est colinéaire à l'axe 3 :

Omega={{0,-omega,0},{omega,0,0},{0,0,0}};
Id={{1,0,0},{0,1,0},{0,0,1}};
A=Id+Sin[omega*t]/omega*Omega+2*(Sin[omega*t/2]/omega)^2*Omega.Omega

(\begin{matrix} 1 - 2 {sin}^{2} (\frac{ω t}{2}) & - sin (ω t) & 0 \\ sin (ω t) & 1 - 2 {sin}^{2} (\frac{ω t}{2}) & 0 \\ 0 & 0 & 1 \end{matrix})

La position s'obtient par intégration directe :

\begin{matrix} q (t) & = q (0) + \frac{1}{m} \int_{0}^{t} exp (Ω s) p (0) d s & (16) \\ = q (0) + \frac{1}{m} F (t) p (0) & (17) \end{matrix}

où la matrice F(t) est définie par :

F (t) = I_{d} t + \frac{1 - cos (ω t)}{ω^{2}} Ω - \frac{sin (ω t) - ω t}{ω^{3}} Ω^{2} (18)

Si le champ magnétique est colinéaire à l'axe 3, on a :

F=Simplify[Id*t+(1-Cos[omega*t])/omega^2*Omega-(Sin[omega*t]-omega*t)/omega^3*Omega.Omega]

(\begin{matrix} \frac{sin (ω t)}{ω} & \frac{cos (ω t) - 1}{ω} & 0 \\ \frac{1 - cos (ω t)}{ω} & \frac{sin (ω t)}{ω} & 0 \\ 0 & 0 & t \end{matrix})

3. Méthode numérique symplectique

On cherche une solution numérique des équations (1) et (2) avec les conditions initiales à t=0 :

\begin{matrix} q (0) = q_{0} & (19) \\ p (0) = p_{0} & (20) \end{matrix}

Pour obtenir N points sur l'intervalle [0,T], on définit le pas de temps h=T/N et on cherche des valeurs approchées pour les instants :

t_{n} = n h (21)

où l'entier n varie de 0 à N.

Une méthode symplectique est proposée en [1], analogue à la méthode de Störmer-verlet, avec prise en charge du champ magnétique :

\begin{matrix} p_{i} = p_{n} - \frac{1}{2} h \frac{\partial V}{\partial q} (q_{n}) & (22) \\ q_{n + 1} = q_{n} + \frac{1}{m} F (h) p_{i} & (23) \\ p_{n + 1} = exp (Ω h) p_{i} - \frac{1}{2} h \frac{\partial V}{\partial q} (q_{n + 1}) & (24) \end{matrix}

Les quantités de mouvement intermédiaires (vecteur $p_{i}$ ) sont tout d'abord évaluées puis les nouvelles positions sont calculées avec ces quantités de mouvement. Enfin les nouvelles quantités de mouvement sont calculées avec les nouvelles positions.

Lorsque le champ magnétique est uniforme (par exemple une force de Coriolis dans un référentiel en rotation), on aura intérêt à choisir l'axe 3 comme direction du champ et à utiliser les matrices calculées plus haut.

Lorsque le champ magnétique n'est pas uniforme, la méthode est en principe utilisable à condition d'évaluer les matrices F(h) et exp(Ωh) pour les positions q_n.

Les mêmes équations s'appliquent à un système de plusieurs particules soumises au champ magnétique, et dont les forces d'interaction dérivent d'un potentiel. Dans ce cas, les matrices 3x3 F(h) et exp(Ωh) s'appliquent aux colonnes p_i correspondant aux coordonnées des particules. En revanche, le cas d'interactions magnétiques entre particules ne relève pas de ce système.

4. programmation avec Scilab

4.a. Fonctions de calcul

La dérivée de l'énergie potentielle par rapport aux positions (l'opposée de la force) est définie dans une fonction de la forme suivante. L'exemple considéré ici est un champ coulombien attractif.

function [dp]=dV(q)
    dp = q/(q(1)^2+q(2)^2+q(3)^2)^(1.5)        
endfunction

La matrice Ω est aussi définie dans une fonction. Pour l'instant, le champ magnétique est uniforme

function Om=Omega(q)
    Om=[0,-1,0;1,0,0;0,0,0]
endfunction

On fournit aussi une fonction qui calcule l'énergie du système, afin de contrôler la stabilité du calcul :

function e=energie(q,p)
    e=0.5*(p(1)^2+p(2)^2+p(3)^2)-1/(q(1)^2+q(2)^2+q(3)^2)^0.5
endfunction

La fonction suivante effectue un pas élémentaire de la méthode de Stormer-verlet :

function [q,p]=pasVerletMag(dV,Omega,h,qn,pn)
    pi=pn-0.5*h*dV(qn);
    Om=Omega(qn);
    Om2=Om*Om;
    w=(Om(1,2)^2+Om(1,3)^2+Om(2,3)^2)^0.5;
    Id=[1,0,0;0,1,0;0,0,1];
    F=Id*h+(1-cos(w*h))/w^2*Om-(sin(w*h)-w*h)/w^3*Om2;
    q=qn+F*pi;
    E=Id+sin(w*h)/w*Om+2*(sin(w*h*0.5)/w)^2*Om2;
    p=E*pi-0.5*h*dV(q);
endfunction

Cette fonction traite le cas général; il est souhaitable de la simplifier dans le cas d'un champ magnétique uniforme.

Pour stocker les valeurs obtenues sur l'intervalle [0,T], on utilise la convention de la fonction ode de scilab : les valeurs de q sont stockées dans une matrice dont chaque ligne correspond à une composante de q (de même pour p). La fonction suivante effectue le calcul complet, avec une condition initiale (instant 0) définie dans deux matrices colonne.

n₂ est le nombre de points mémorisés, n₁ le nombre de pas entre deux points mémorisés.

function [mq,mp,e,t] = verletMag(q0,p0,T,n1,n2,dV,Omega,energie)
    h=T/(n1*n2);
    mq=zeros(3,n2+1);
    mp=zeros(3,n2+1);
    e=zeros(1,n2+1);
    t=zeros(1,n2+1);
    mq(:,1)=q0;
    mp(:,1)=p0;
    e(:,1)=energie(q0,p0);
    q=q0; 
    p=p0;
    t(1)=0;
    for k2=1:n2,
        for k1=1:n1,
            [q,p]=pasVerletMag(dV,Omega,h,q,p);
        end;
        mq(:,k2+1)=q;
        mp(:,k2+1)=p;
        e(:,k2+1)=energie(q,p);
        t(k2+1)=t(k2)+n1*h;
    end;
endfunction

4.b. Exemple : champ coulombien

Dans cet exemple, une particule chargée est soumise à un champ coulombien attractif et à un champ magnétique uniforme.

function [dp]=dV(q)
    dp =  q/(q(1)^2+q(2)^2+q(3)^2)^(1.5)
endfunction
function Om=Omega(q)
    Om=[0,-1,0;1,0,0;0,0,0]
endfunction
function e=energie(q,p)
    e=0.5*(p(1)^2+p(2)^2+p(3)^2)-1/(q(1)^2+q(2)^2+q(3)^2)^0.5
endfunction
q0=[1;1;1];
p0=[0;0;0];
T=100;
n1=10;
n2=1000;
[mq,mp,e,t]=verletMag(q0,p0,T,n1,n2,dV,Omega,energie);

On trace la projection de la trajectoire dans le plan xy :

figA=scf();
plot2d(mq(1,:),mq(2,:))
xtitle('','x','y')

figA.pdf

Voyons l'énergie en fonction du temps :

figB=scf();
plot2d(t,e)
xtitle('','t','E')

figB.pdf

Références

[1] B. Leimkuhler, S. Reich, Simulating hamiltonian dynamics, (Cambridge university press, 2004)

[2] J.P. Demailly, Analyse numérique et équations différentielles, (P.U.G., 1991)