Méthode d'Euler symplectique

1. Système différentiel hamiltonien

1.a. Équations de Hamilton

On considère un système dynamique à N_v degrés de libertés. On note q les coordonnées généralisées (un vecteur de N_v variables), et p les variables conjuguées. L'hamiltonien H(q,p) est une fonction des 2N_v variables p et q.

Les équations du mouvement de Hamilton s'écrivent :

\begin{matrix} \frac{d q}{d t} = + \frac{\partial H}{\partial p} & (1) \\ \frac{d p}{d t} = - \frac{\partial H}{\partial q} & (2) \end{matrix}

Pour les systèmes mécaniques conservatifs, l'hamiltonien est séparable en une partie énergie cinétique et une partie énergie potentielle :

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

Les équations de Hamilton s'écrivent alors :

\begin{matrix} \frac{d q}{d t} = + \frac{\partial K}{\partial p} & (3) \\ \frac{d p}{d t} = - \frac{\partial V}{\partial q} & (4) \end{matrix}

1.b. Exemples en mécanique classique

L'hamiltonien de l'oscillateur harmonique à une dimension est :

H = \frac{1}{2} p^{2} + \frac{1}{2} q^{2}

Les équations du mouvement s'écrivent :

\begin{matrix} \frac{d q}{d t} = p & (5) \\ \frac{d p}{d t} = - q & (6) \end{matrix}

Le problème à N corps en interaction gravitationnelle est un autre exemple dont l'hamiltonien est séparable. Par exemple pour 2 corps, si q₁ et q₂ sont les vecteurs position des deux corps :

\begin{matrix} K (p_{1}, p_{2}) = \frac{p_{1}^{2}}{2 m_{1}} + \frac{p_{2}^{2}}{2 m_{2}} & (7) \\ V (q_{1}, q_{2}) = - \frac{G m_{1} m_{2}}{| q_{1} - q_{2} |} & (8) \end{matrix}

2. Méthodes d'Euler symplectique

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

\begin{matrix} q (0) = q_{0} & (9) \\ p (0) = p_{0} & (10) \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 (11)

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

La méthode d'Euler A ([1]) est définie par :

\begin{matrix} q_{n + 1} = q_{n} + h \frac{\partial H}{\partial p} (q_{n + 1}, p_{n}) & (12) \\ p_{n + 1} = p_{n} - h \frac{\partial H}{\partial q} (q_{n + 1}, p_{n}) & (13) \end{matrix}

Cette méthode est d'ordre 1. La première étape est en général implicite. Dans le cas d'un hamiltonien séparable elle devient :

\begin{matrix} q_{n + 1} = q_{n} + h \frac{\partial K}{\partial p} (p_{n}) & (14) \\ p_{n + 1} = p_{n} - h \frac{\partial V}{\partial q} (q_{n + 1}) & (15) \end{matrix}

La méthode est alors explicite. Comparée à la méthode d'Euler simple (explicite), l'évaluation des coordonnées q et des moments conjugés p (les vitesses), se fait en deux étapes successives. Pour cette raison, on parle aussi de méthode d'Euler assymétrique.

La méthode d'Euler B est définie par :

\begin{matrix} q_{n + 1} = q_{n} + h \frac{\partial H}{\partial p} (q_{n}, p_{n + 1}) & (16) \\ p_{n + 1} = p_{n} - h \frac{\partial H}{\partial q} (q_{n}, p_{n + 1}) & (17) \end{matrix}

La seconde étape est en général implicite. Dans le cas d'un hamiltonien séparable elle devient explicite :

\begin{matrix} q_{n + 1} = q_{n} + h \frac{\partial K}{\partial p} (p_{n + 1}) & (18) \\ p_{n + 1} = p_{n} - h \frac{\partial V}{\partial q} (q_{n}) & (19) \end{matrix}

On remarque que l'évaluation des p doit se faire en premier.

Ces méthodes d'Euler sont symplectiques ([1]). Nous allons voir les conséquences de cette propriété sur un exemple.

3. Programmation avec Scilab

3.a. Fonctions de calcul

On se limite au cas d'un hamiltonien séparable. Le système différentiel doit être défini par deux fonctions, une pour la dérivée de l'énergie cinétique, l'autre pour la dérivée de l'énergie potentielle. On considère l'exemple de l'oscillateur harmonique à une dimension :

function [dq]=dK(p),
    dq(1)=p(1)
endfunction
function [dp]=dV(q),
    dp(1)=q(1)
endfunction

La fonction suivante effectue un pas de la méthode Euler-A :

function [q,p]=pasEulerA(dK,dV,nv,h,qn,pn),
     dq=dK(pn);
     for v=1:nv,
        q(v)=qn(v)+h*dq(v)
     end;
     dp=dV(q);
     for v=1:nv,
        p(v)=pn(v)-h*dp
     end;
endfunction

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.

function [mq,mp,t]=eulerA(qi,pi,T,N,dK,dV,nv),
    h = T/N;
    mq = zeros(nv,N+1);
    mp = zeros(nv,N+1);
    t=zeros(1,N+1);
    mq(:,1)=qi;
    mp(:,1)=pi;
    t(1)=0;
    for k=1:N,
        [mq(:,k+1),mp(:,k+1)]=pasEulerA(dK,dV,nv,h,mq(:,k),mp(:,k));
        t(k+1)=t(k)+h;
    end;
endfunction

3.b. Exemple : oscillateur harmonique

On teste le fonctionnement avec l'oscillateur harmonique défini plus haut :

qi=[1];
pi=[0];
nv=1;
N=1000;
T=20;
[mq,mp,t]=eulerA(qi,pi,T,N,dK,dV,nv);

figA=scf();
plot2d(t,mq(1,:));
xtitle('','t','q');

figA.pdf

La solution exacte étant connue, on peut tracer l'erreur en fonction du temps :

T=20;
N=1000;
[mq,mp,t]=eulerA(qi,pi,T,N,dK,dV,nv);

figB=scf();
plot2d(t,abs(mq(1,:)-cos(t)));
xtitle('','t','erreur');

figB.pdf

Pour un même pas de temps, l'erreur est beaucoup plus faible qu'avec la méthode d'Euler.

Pour voir le comportement à long terme, on trace la trajectoire dans le plan de phase :

T=2000;
N=T*50;
[mq,mp,t]=eulerA(qi,pi,T,N,dK,dV,nv);

figC=scf();
plot2d(mq(1,:),mp(1,:),rect=[-1.5,-1.5,1.5,1.5]);
xtitle('','p','q');

figC.pdf

Contrairement aux méthodes d'Euler (explicite et implicite), la trajectoire reste très proche de la solution exacte. On obtient le comportement attendu pour un oscillateur non dissipatif : une trajectoire fermée dans l'espace des phases. Voyons l'énergie mécanique en fonction du temps :

Em=zeros(1,N+1);
for k=1:N+1,
    Em(k)=0.5*mp(1,k)^2+0.5*mq(1,k)^2;
end;

figD=scf();
plot2d(t,Em,rect=[0,0,T,1]);
xtitle('','t','Energie');

figD.pdf

On constate sur cet exemple une propriété très importante des méthodes symplectiques appliquées aux systèmes conservatifs : l'énergie fluctue mais ses variations restent bornées. Les méthodes classiques de type Runge-Kutta ou les méthodes multi-pas présentent généralement des dérives non bornées de l'énergie mécanique.

4. Exemple avec Mathematica

Les méthodes symplectiques sont implémentées dans la fonction NDSolve. On commence par définir les équations et les conditions initiales, par exemple pour l'oscillateur harmonique :

H=0.5*p[t]^2+0.5*q[t]^2;
equations = {q'[t]==D[H,p[t]],p'[t]==-D[H,q[t]]};
init={q[0]==1,p[0]==0};

On choisit le temps final, le pas de temps puis on exécute NDSolve en précisant les variables de positions.

T=2000;
n=T*50;
h=T/n;
sol = NDSolve[{equations,init},{q[t],p[t]},{t,0,T},StartingStepSize->h,MaxSteps->Infinity,
        Method->{"SymplecticPartitionedRungeKutta","DifferenceOrder"->1,"PositionVariables"->{q[t]}}];

ParametricPlot[{q[t],p[t]}/.sol,{t,0,T},PlotRange->{{-1.5,1.5},{-1.5,1.5}}]

figE.pdf

Références

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