rk4-02.tex

\documentclass[svgnames,a4paper,10pt]{article}
\ifdefined\XeTeXinterchartoks
\usepackage{fontspec}
\usepackage[math-style=TeX]{unicode-math}
\usepackage{xltxtra}
\setmathfont{texgyrepagella-math.otf}
\else
\usepackage{amsmath}
\usepackage[T1]{fontenc}
\usepackage{bera}
\fi
 
\usepackage[french]{babel}
\usepackage[margin=2.5cm]{geometry}
\usepackage{shortvrb}
\usepackage[dvipsnames]{xcolor}
\usepackage{graphicx}
\usepackage{pstricks-add}
\usepackage{multicol}
\usepackage{xintexpr}
 
\ifdefined\XeTeXinterchartoks
% ici a priori c'est ok de rester à l'intérieur de la branche conditionnelle
\defaultfontfeatures{Mapping=tex-text} 
%\setromanfont [Ligatures={Common}, Numbers={OldStyle}]{Hoefler Text}
\setmonofont[
    Path=fontes/PFDINMono/,
    Scale=0.95,
    BoldFont = PFDINMono-Bold.otf,
    ItalicFont = PFDINMono-Italic.otf,
    BoldItalicFont = PFDINMono-BoldItalic.otf
  ]{PFDINMono.otf}
\setmainfont[
    Ligatures=Common,
    Scale=1.1,
    Path=fontes/Athelas/,
    BoldFont = Athelas-Bold.otf,
    ItalicFont = Athelas-Italic.otf,
    BoldItalicFont = Athelas-BoldItalic.otf
  ]{Athelas-Regular.otf}
\setsansfont[
    Ligatures=Common,
    Path=fontes/OptimaLTStd/,
    BoldFont = OptimaLTStd-Bold.otf,
    ItalicFont = OptimaLTStd-Italic.otf,
    BoldItalicFont = OptimaLTStd-BoldItalic.otf
  ]{OptimaLTStd.otf}
\fi
 
\MakeShortVerb{\|}
\newcommand\ftf[2]{\xintieval{[#1]\xintfloatexpr#2\relax}}
\newcommand\ft[1]{\xintfloatexpr #1\relax}
\input{sty-jfbu}
\input{xintf-rk4-4x4}
\input{xintf-io1}
 
%% ========================================================================== %%
%% Code déplacé pour éviter les problèmes avec ; et ?  (actifs avec french)   %%
%% ========================================================================== %%
\def\FigureIV#1{
 \psset{xunit=1.5cm,yunit=1.5cm}
 \begin{pspicture}(-0.1,-2)(6.5,2)
   \psaxes{->}(0,0)(-0.1,-2)(6.5,2)
   \multido{\i=1+1}{\ftf0{N}}{%
      \xintdeffloatvar ti,Ki := LT\i,K(pvc(\i));
      \psdot[linecolor=LimeGreen,dotsize=2pt](\ftf4{LT\i},\ftf4{#1(Ki-K0)})
      \psdot[linecolor=DodgerBlue,dotsize=2pt](\ftf4{LT\i},\ftf4{Ki})
    }%
  \end{pspicture}
}
 
 
\begin{document}
\parindent0pt
\setlength{\parskip}{1ex}
% pour les short verb |...|
\def\MicroFont
    {\ttfamily\color[named]{DarkCyan}\makestarlowast\makequotesstraight}
% pour les everbatim
\def\MacroFont {\ttfamily\color[named]{Olive}}
% Note 1 : everbatim fait \makestarlowast\makequotesstraight en interne
%          alors que pour les short verb il faut l'ajouter via \MicroFont
% Note 2 : everbatim fait un \small via \def\everbatimtop{\MacroFont\small}
 
{\Large Construction de l'orbite de \textbf{Hill}  avec \textbf{xint}}\\
\rule{\linewidth}{1pt}
\medskip
 
\textbf{George William Hill} (1838-1914) est un astronome et mathématicien
américain qui a beaucoup étudié le problème des trois ou quatre corps, en
particulier le \emph{problème des trois corps restreint} (corps de masse
nulle dans le champ de gravitation de deux corps massifs).
 
Une illustration de ce dernier problème peut-être vue dans la trajectoire
d'un satellite artificiel dans le champ de gravitation du couple
\emph{Terre-Lune}. Il y a une grande variété d'orbites possibles sous le
seul effet de l'attraction universelle, l'une d'entre elles est
fascinante, c'est justement celle qui a été mise à jour par \textbf{G.W.
Hill}.
 
Nous allons envisager le satellite dans le plan du mouvement de la Terre
et de la Lune autour de leur centre de gravité\footnote{Le mouvement de
la Lune et de la Terre est considéré comme circulaire, c'est déjà assez
compliqué comme ça!}. La masse de la Lune sera notée $\mu$ et celle de la
Terre $1-\mu$, de sorte que la masse totale soit $1$. Dans le repère
tournant autour du centre de gravité
\emph{Terre-Lune}, les coordonnées de la Terre seront $(-\mu,0)$ et celle
de la Lune $(1-\mu,0)$; ainsi l'unité de distance sera fixée.
 
Pour faire court, après changement de référentiel (passage dans le repère
tournant - vitesse de rotation: $1$) et transformation canonique (pour
disposer de variables conjuguées), nous avons les équations du mouvement:
 
$$\begin{matrix}
\dot{Q}_x&=&P_x+Q_y\hfill\\[2mm]
\dot{Q}_y&=&P_y-Q_x\hfill\\[2mm]
\dot{P}_x&=&P_y+(1-\mu)\frac{\partial}{\partial Q_x}\left(\frac1{R_T}\right)+\mu\frac{\partial}{\partial Q_y}\left(\frac1{R_L}\right)\\[2mm]
\dot{P}_y&=&-P_x+(1-\mu)\frac{\partial}{\partial Q_x}\left(\frac1{R_T}\right)+\mu\frac{\partial}{\partial Q_y}\left(\frac1{R_L}\right)
\end{matrix}$$
où $(Q_x,Q_y)$ est la position du satellite, $(P_x,P_y)$ les moments
conjugués, $R_T=\sqrt{(Qx+\mu)^2+Q_y^2}$,
$R_L=\sqrt{(Q_x-1+\mu)^2+Q_y^2}$.
 
\subsection*{Mise en place du second membre}
 
La transcription du système ci-dessus est un peu lourde mais il ne faut
reculer devant rien avec \textbf{xint}.
\begin{everbatim*}
 \xintdeffloatvar m := 1/82.45;% <-- mu
 \xintdeffloatvar M := 1-m;%
 %
 \xintdeffloatfunc r32(x,y,c) := ((x+c)^2+y^2)^(1.5);%
 \xintdeffloatfunc sm(t,qx,qy,px,py) := px+qy,py-qx,%
                      py-M(qx+m)/r32(qx,qy,m)-m(qx-M)/r32(qx,qy,-M),%
                      -px-qy*(M/r32(qx,qy,m)+m/r32(qx,qy,-M));%
\end{everbatim*}
 
\subsection*{Conditions initiales}
 
Pour obtenir l'orbite périodique, les conditions initiales sont:
 
\begin{everbatim*}
 \xintdeffloatvar t0  := 0;%
 \xintdeffloatvar tf  := 6.1927;%    Période de l'orbite de Hill
 \xintdeffloatvar hi  := 0.01;%      Pas de temps initial;
 \xintdeffloatvar qx0 := 1.2;%
 \xintdeffloatvar qy0 := 0;%
 \xintdeffloatvar px0 := 0;%
 % La vitesse initiale est -1.04935750483 selon l'axe des y
 \xintdeffloatvar py0 := 1.2-1.04935750483;%
\end{everbatim*}
 
\newpage
\subsection*{Go !}
 
\begin{everbatim*}
\rkSeuils{}
\rkSolIV(t0,tf,hi,qx0,qy0,px0,py0,hill01) % <- Les données seront sauvegardées dans
                                          % les fichiers hill01_*.xint
\xinteval{N}, \ftf4{hmin}, \ftf4{hmax}
\end{everbatim*}
 
Vu les pas minimal et maximal, l'attraction a dû être forte à certains
moments alors qu'à d'autres le mouvement était tranquille...
 
À ce niveau là de notre scénario, les points calculés de la trajectoire
ont été sauvegardés; lors de toute recompilation, ce sont ces données
sauvegardées qui seront utilisées, on gagne du temps et on peut se livrer
à toutes les analyses souhaitées ainsi que multiplier les
représentations. Pour recalculer une trajectoire, avec des paramètres
différents par exemple, il suffit de supprimer le fichier
|hill01_par.xint| ou de faire de nouveau appel à |\rkSolIV| avec un
\emph{préfixe} de fichier différent (le huitième argument).
 
\subsection*{Enfin, la trajectoire de Hill !}
 
\begin{everbatim*}
\begin{center}
 \psset{xunit=5cm,yunit=5cm}
 \begin{pspicture}(-1.3,-0.7)(1.3,0.7)%
  \psdot[linecolor=DarkGoldenrod,dotsize=4pt](\ftf4{-m},0)
  \psdot[linecolor=DimGray,dotsize=2.5pt](\ftf4{M},0)
  \multido{\i=1+1}{\ftf0{N}}{%
   \psdot[linecolor=DodgerBlue,dotsize=2pt](\ftf4{QX\i},\ftf4{QY\i})
  }%
  \uput{8pt}[-90](\ftf4{-m},0){\sffamily\textbf{T}}
  \uput{8pt}[-90](\ftf4{M},0){\sffamily\textbf{L}}
 \end{pspicture}
\end{center}
\end{everbatim*}
 
Nous allons interpoler la courbe avec |\listplot|.
\begin{everbatim*}
\xintNewFunction{point}[1]{%
  \ftf4{QX\xinteval{#1}},\ftf4{QY\xinteval{#1}}}%
\edef\Points{%
  \xintthespaceseparated \xintfloatexpr seq(point(i),i=1..N)\relax
}%
\end{everbatim*}
\newpage
 
\begin{everbatim*}
\begin{center}
  \psset{xunit=5cm,yunit=5cm}
  \begin{pspicture}(-1.3,-0.7)(1.3,0.7)%
    \psdot[linecolor=DarkGoldenrod,dotsize=4pt](\ftf4{-m},0)
    \psdot[linecolor=DimGray,dotsize=2.5pt](\ftf4{M},0)
    \listplot[plotstyle=curve,linecolor=LightSteelBlue,linewidth=2pt]{\Points}%
    \multido{\i=1+1}{\ftf0{N}}{%
      \psdot[linecolor=DodgerBlue,dotsize=2pt](\ftf4{QX\i},\ftf4{QY\i})
    }%
    \uput{8pt}[-90](\ftf4{-m},0){\sffamily\textbf{T}}
    \uput{8pt}[-90](\ftf4{M},0){\sffamily\textbf{L}}
  \end{pspicture}
\end{center}
\end{everbatim*}
 
\subsection*{Et dans le repère fixe ?}
 
Dans la normalisation du mouvement, les choix ont conduit à ce que la
vitesse de rotation ($\omega$) du repère tournant soit égale à $1$.
 
Préparons la mise en scène de la rotation.
\begin{everbatim*}
\xintdeffloatfunc rotation(x,y,t) := cos(t)*x-sin(t)*y,sin(t)*x+cos(t)*y;%
\xintNewFunction{pointt}[1]{%
  \ftf4{QX\xinteval{#1}},\ftf4{QY\xinteval{#1}},\ftf4{LT\xinteval{#1}}}%
\edef\PointsR{%
  \xintthespaceseparated \xintfloatexpr seq(rotation(pointt(i)),i=1..N)\relax
}%
\end{everbatim*}
 
Regardons! Les trajectoires de la Terre et de la Lune sont représentées
sur un cycle complet, ce qui correspond à peu près, à la période de
l'orbite de \textbf{Hill}.
 
\begin{everbatim*}
\begin{center}
  \psset{xunit=5cm,yunit=5cm}
  \begin{pspicture}(-1.6,-1)(1.6,1)%
    \pscircle[linecolor=DarkGoldenrod](0,0){\ftf4{5m}}
    \pscircle[linecolor=DimGray](0,0){\ftf4{5M}}
    \listplot[plotstyle=curve,linecolor=LightSteelBlue,linewidth=2pt]{\PointsR}%
    \multido{\i=1+1}{\ftf0{N}}{%
      \psdot[linecolor=DodgerBlue,dotsize=2pt](\ftf4{rotation(pointt(\i))})
    }%
    % \uput{8pt}[-90](\ftf4{-m},0){\sffamily\textbf{S}}
    % \uput{8pt}[-90](\ftf4{M},0){\sffamily\textbf{L}}
  \end{pspicture}
\end{center}
\end{everbatim*}
\medskip
 
C'est moins spectaculaire, néanmoins amusant: nos amis sélènes peuvent
envisager, en partant très tôt, de venir faire leurs courses chez nous et
de retourner pour la collation de minuit lunaire; à soumettre à \emph{Elon Musk}! 
 
 
\subsection*{Mais quelle précision ?}
 
Ici, on ne dispose pas de la solution exacte du problème pour faire une
comparaison effective avec la solution calculée.
 
Tout n'est pas perdu, il y une \emph{intégrale première},
l'\textbf{hamiltonien} $K$, dont nous allons pouvoir vérifier la stabilité
tout au long des calculs.
 
Voici l'expression de $K$:
$$K=\frac12\left(P_x^2+P_y^2\right)+P_xQ_y-P_yQ_x-\left(\frac{1-\mu}{R_T}+\frac{\mu}{R_L}\right)$$
 
Et sa transcription avec \textbf{xint}:
\begin{everbatim*}
\xintdeffloatfunc R(x,y,c) := sqrt((x-c)^2+y^2);%
\xintdeffloatfunc K(qx,qy,px,py) := 0.5(px^2+py^2)+px*qy-py*qx-(M/R(qx,qy,-m)+m/R(qx,qy,M));%
\xintdeffloatvar K0 := K(qx0,qy0,px0,py0);%
\ftf9{K0}
\end{everbatim*}
 
Nous aurons besoin des paires de variables conjuguées:
\begin{everbatim*}
\xintNewFunction{pvc}[1]{%
 QX\xinteval{#1},QY\xinteval{#1},PX\xinteval{#1},PY\xinteval{#1}%
}%
\end{everbatim*}
 
\begin{everbatim}
\begin{center}
 \psset{xunit=1.5cm,yunit=1.5cm}
 \begin{pspicture}(-0.1,-2)(6.5,2)
   \psaxes{->}(0,0)(-0.1,-2)(6.5,2)
   \multido{\i=1+1}{\ftf0{N}}{%
      \xintdeffloatvar ti,Ki := LT\i,K(pvc(\i));
      \psdot[linecolor=Crimson,dotsize=2pt](\ftf4{LT\i},\ftf4{10000(Ki-K0)})
      \psdot[linecolor=DodgerBlue,dotsize=2pt](\ftf4{LT\i},\ftf4{Ki})
    }%
  \end{pspicture}
\end{center}
\end{everbatim}
 
\begin{center}
\FigureIV{10000}
\end{center}
 
En bleu, l'hamiltonien, en vert, la différence entre l'hamiltonien de
l'origine et celui de l'instant multipliée par 10000. On constate les
dégâts (\emph{relatifs}) à proximité de la Terre.
 
Essayons d'être plus contraignant sur la conservation du pas.
\begin{everbatim*}
\xintdeffloatvar epsilon := 1e-6;
\rkSolIV(t0,tf,hi,qx0,qy0,px0,py0,hill02)
\xinteval{N}, \ftf6{hmin}, \ftf6{hmax}
\end{everbatim*}
\medskip
 
\begin{center}
\FigureIV{1000000}
\end{center}
 
Là, la différence entre l'hamiltonien de l'origine et celui de l'instant a
été multipliée par $10^6$. On peut estimer que l'on dispose, avec cette
trajectoire, d'une précision de $10^{-6}$; ce dont je vais me satisfaire
pour l'instant!
 
\textbf{Idée pour la suite}: Mettre en