\documentclass{article}
\usepackage[a4paper,margin=2cm]{geometry}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}%
\usepackage[garamond]{mathdesign}
\renewcommand{\ttdefault}{lmtt}
\usepackage{pst-node,pst-eqdf}
\usepackage{amsmath}
\usepackage{url}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Essai de Manuel Luque 19 février 2003
% transformé par Denis Girou le 25 février 2003
% révision le 14 septembre 2011
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\makeatletter
\pst@addfams{pst-plothsb}

\define@key[psset]{pst-plothsb}{HueBegin}{% Between 0 and 1
\edef\PstParametricplotHSB@HueBegin{#1}}

\define@key[psset]{pst-plothsb}{HueEnd}{% Between 0 and 1
\edef\PstParametricplotHSB@HueEnd{#1}}

\newif\ifPst@HSB
\define@key[psset]{pst-plothsb}{HSB}[true]{\@nameuse{Pst@HSB#1}}

% Default values
\psset{HueBegin=0,HueEnd=1,HSB=true}

\def\parametricplotHSB{\pst@object{parametricplotHSB}}
\def\parametricplotHSB@i{\@ifnextchar[{\parametricplotHSB@do}{\parametricplotHSB@do[]}}
\def\parametricplotHSB@do[#1]#2#3#4{{%
\psset{#1}%
\begin@ClosedObj
\addto@pscode{%
/t  #2 def
/dt #3 t sub \psk@plotpoints\space div def
/t t dt sub def
/Counter 0 def
\psk@plotpoints {
  /t t dt add def
  /Counter Counter 1 add def
/F@pstplot \ifPst@algebraic (#4)
                    tx@AlgToPs begin AlgToPs end cvx
                 \else { #4 }
          \fi
   def
\ifPst@algebraic
   F@pstplot
    \else  #4
\fi
  \pst@number\psyunit mul exch
  \pst@number\psxunit mul exch
  1 Counter eq
    {moveto}                    % First point
    {\ifPst@HSB                 % Other points than the first one
       /PointY exch def
       /PointX exch def
       Counter \psk@plotpoints\space div
           \PstParametricplotHSB@HueEnd\space
             \PstParametricplotHSB@HueBegin\space sub mul
           \PstParametricplotHSB@HueBegin\space add
         1 1 sethsbcolor
       PointX PointY lineto
       stroke
       PointX PointY moveto
     \else
       lineto
     \fi} ifelse
    } repeat}% fin du code ps
\end@ClosedObj}} % fin de la commande PSTricks

\makeatother

\makeatletter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% adaptation de \parametricplotHSB
% à \lisplotHSB
% le 10 juin 2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\listplotHSB{\pst@object{listplotHSB}}
\def\listplotHSB@i#1{%
\begin@ClosedObj
  \addto@pscode{%
  /Counter 0 def
     /XY [ #1 ] def
     XY 0 get
     XY 1 get
  \pst@number\psyunit mul exch
  \pst@number\psxunit mul exch
       moveto
     /L XY length def
     2 2 L 2 sub {/i exch def
   /x1  XY i get \pst@number\psxunit mul def
   /y1  XY i 1 add get \pst@number\psyunit mul def
    x1 y1
    \ifPst@HSB
    Counter L 2 div div
    \PstParametricplotHSB@HueEnd\space
    \PstParametricplotHSB@HueBegin\space sub mul
    \PstParametricplotHSB@HueBegin\space add
    1 1 sethsbcolor
    /Counter Counter 1 add def
    \fi
     lineto
     stroke
  x1 y1 moveto
  } for
  }%
\end@ClosedObj}
\makeatother
\def\formules{%
    /xObservateur abscisse Sin1 mul neg ordonnee Cos1 mul add def
    /yObservateur abscisse Cos1Sin2 mul neg ordonnee Sin1Sin2 mul sub cote Cos2 mul add def
    /zObservateur abscisse neg Cos1Cos2 mul ordonnee Sin1Cos2 mul sub cote Sin2 mul sub Dobs add def
    /xScreen DScreen xObservateur mul zObservateur div def
    /yScreen DScreen yObservateur mul zObservateur div def
  xScreen yScreen}

\makeatletter
\pst@addfams{pst-pnodeTroisD}

\define@key[psset]{pst-pnodeTroisD}{THETA}{\edef\psk@TroisD@Theta{#1 }}
\define@key[psset]{pst-pnodeTroisD}{PHI}{\edef\psk@TroisD@Phi{#1 }}
\define@key[psset]{pst-pnodeTroisD}{Dobs}{\edef\psk@TroisD@Dobs{#1 }}
\define@key[psset]{pst-pnodeTroisD}{Decran}{\edef\psk@TroisD@Ecran{#1 }}
\psset{THETA=300,PHI=20,Dobs=200,Decran=25}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\pnodeTroisD{\@ifnextchar[{\pst@pnodeTroisD}{\pst@pnodeTroisD[]}}

\def\pst@pnodeTroisD[#1](#2,#3,#4)#5{{%
\pst@killglue
\psset{#1}%
\pnode(!
 /THETA \psk@TroisD@Theta\space def
 /PHI \psk@TroisD@Phi\space def
 /Dobs \psk@TroisD@Dobs\space def
 /DScreen \psk@TroisD@Ecran\space def
 /Sin1 THETA sin def
 /Sin2 PHI sin def
 /Cos1 THETA cos def
 /Cos2 PHI cos def
 /Cos1Sin2 Cos1 Sin2 mul def
 /Sin1Sin2 Sin1 Sin2 mul def
 /Cos1Cos2 Cos1 Cos2 mul def
 /Sin1Cos2 Sin1 Cos2 mul def
 /xObservateur #2 Sin1 mul neg #3 Cos1 mul add def
 /yObservateur #2 Cos1Sin2 mul neg #3 Sin1Sin2 mul sub #4 Cos2 mul add def
 /zObservateur #2 neg Cos1Cos2 mul #3 Sin1Cos2 mul sub #4 Sin2 mul sub Dobs add def
 /xScreen DScreen xObservateur mul zObservateur div def
 /yScreen DScreen yObservateur mul zObservateur div def
  xScreen yScreen){#5}}}
\makeatother
\title{L'attracteur de Lorenz avec PSTricks}
\date{12 juin 2012}
\begin{document}
\maketitle
\section{Présentation}
Ceci n'est pas un cours sur l'attracteur de Lorenz, il existe sur internet un grand nombre de sites qui lui sont consacrés très complets et d'une grande qualité pédagogique. Il s'agit, plus simplement, de montrer que l'utilisation de PSTricks permet d'illustrer ce modèle très convenablement.

Pour résoudre numériquement le système de Lorenz :
\[
\left\{
\begin{array}{rcl}
\dfrac{\mathrm{d}x}{\mathrm{d}t}&=&\sigma(x-y)\\
\dfrac{\mathrm{d}y}{\mathrm{d}t}&=&rx-y-z\\
\dfrac{\mathrm{d}z}{\mathrm{d}t}&=&xy-\beta z
\end{array}
\right.
\]
j'utilise une version très simplifiée de la macro \verb+\psplotDiffEqn+ de Dominique Rodriguez qui est incluse dans le package \verb+pstricks-add+, mais qui était, à l'origine, un package `\verb+pst-eqdf+' à elle seule. Dominique m'a donné la permission de modifier sa macro afin de pouvoir sauvegarder dans un tableau \verb+postscript+, et/ou dans un fichier, toutes les valeurs nécessaires au dessin des figures et en particulier à la représentation en 3D du \textit{papillon }de Lorenz. Cette macro simplifiée, qui utilise la méthode Runge-Kutta~4 est incluse dans un tout petit package dont j'ai gardé le même nom \verb+pst-eqdf+, mais la macro a changé de nom \verb+\psequadiff+ car elle ne trace plus les fonctions directement.

Pour les couleurs, j'utilise une macro \verb+\listplotHSB+ qui est un complément à celle-ci : \verb+\parametricplotHSB+\footnote{\url{http://melusine.eu.org/syracuse/G/pstricks/courbes-arc-en-ciel/}}, que nous avions mise au point avec Denis Girou.

\section{Les étapes de la construction des figures}
\subsection{Équations, paramètres et conditions initiales}
Il faut d'abord définir les équations en notation algébrique pour la macro \verb+\psequadiff+ :
\begin{verbatim}
% y[0] y[1] y[2]
% x    y    z
\def\Lorenz{
Pr*(y[1]-y[0])|%
Rh*y[0]-y[1]-y[0]*y[2]|%
y[0]*y[1]-b*y[2]%
}
\end{verbatim}
ainsi que les paramètres :
\begin{verbatim}
\def\parametres{
    /Pr 10 def
    /Rh 28 def
    /b 8 3 div def
}%
\end{verbatim}
puis les conditions initiales :
\begin{verbatim}
                         x0 y0 z0
\def\conditionsInitiales{10 10 30}
\end{verbatim}
Ce sont ces deux dernières définitions qu'il faudra modifier si on veut changer les conditions initiales et les paramètres.
\subsection{L'utilisation de la commande \textbackslash{}psequadiff}
\begin{verbatim}
\pstVerb{\parametres}%
  \psequadiff[% x,y
  tabname=tabXY,
  method=rk4,
  whichabs=0,whichord=1,
  saveData,filename=LorenzXY.dat,
  plotpoints=2501,algebraic]{0}{25}{\conditionsInitiales}{\Lorenz}
\end{verbatim}
Dans cet exemple, l'intervalle de temps choisi est de 0 à 25 s. Les grandeurs prises sur la pile sont $(x,y)$. Le nombre de points calculés vaut 2501.
Les valeurs, sur cet intervalle de temps, sont stockées dans le tableau \verb+postscript+ :  \textsf{[\textbf{tabname}=tabXY]} avec \textsf{tabXY} qui est le nom que l'on choisit soi-même. 

Pour ne pas refaire le calcul, on peut choisir d'enregistrer la liste de valeurs du couple $(x,y)$, en positionnant le booléen \textsf{[saveData]} à true, on choisit un nom pour le fichier avec : \textsf{[\textbf{filename}=LorenzXY.dat]}. Par défaut, aucun fichier n'est enregistré.
\begin{verbatim}
 saveData,filename=LorenzXY.dat,
\end{verbatim}
Pour utiliser ce fichier, il suffit de faire :
\begin{verbatim}
    \pstVerb{/XY {(LorenzXY.dat) run } def}%
    \listplotHSB{XY}    
\end{verbatim}
Vous pouvez utiliser la commande \textbf{\textsf{lisplot}} du package \textsf{pst-plot} qui offre davantage d'options que \textsf{\textbf{listplotHSB}}.
\subsection{Le tracé de la courbe}
\begin{verbatim}
\listplotHSB[unit=0.25]{tabXY aload pop}%
\end{verbatim}
\verb+\listplotHSB+ permet de faire varier la couleur tout au long du tracé. Par défaut le paramètre \textsf{hue} de la définition d'une couleur dans le système \textsf{hsb} varie de 0 à 1. On peut fixer la plage de valeurs avec les 2 options suivantes :
\begin{itemize}
  \item \textsf{\textbf{HueBegin}} avec une valeur comprise entre 0 et 1.
  \item \textsf{\textbf{HueEnd}} avec une valeur comprise entre 0 et 1.
\end{itemize}
Pour avoir une couleur uniforme tout au long du tracé, il faut désactiver l'option \textsf{[\textbf{HSB}=false]} ainsi le tracé prendra la couleur que vous avez choisie avec l'option habituelle \textsf{linecolor}.
\section{Le tracé des courbes}
\subsection{Le tracé de XY}
% y[0] y[1] y[2]
 % x    y     z
\def\Lorenz{
Pr*(y[1]-y[0])|%
Rh*y[0]-y[1]-y[0]*y[2]|%
y[0]*y[1]-beta*y[2]%
}
\def\conditionsInitiales{10 10 30}
\def\parametres{
    /Pr 10 def
    /Rh 28 def
    /beta 8 3 div def
}%
\begin{center}
\begin{pspicture}(-8,-6)(8,8)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
\pstVerb{\parametres}%
  \psequadiff[% x,y
  tabname=tabXY,
  method=rk4,
  whichabs=0,whichord=1,
  saveData,filename=LorenzXY.dat,
  plotpoints=2501,algebraic]{0}{25}{\conditionsInitiales}{\Lorenz}
\listplotHSB[unit=0.25]{tabXY aload pop}%
\uput[r](8,0){$x$}
\uput[u](0,8){$y$}
\end{pspicture}
\end{center}
\begin{verbatim}
\begin{pspicture}(-8,-6)(8,8)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
\pstVerb{\parametres}%
  \psequadiff[% x,y
  tabname=tabXY,
  method=rk4,
  whichabs=0,whichord=1,
  saveData,filename=LorenzXY.dat,
  plotpoints=2501,algebraic]{0}{25}{\conditionsInitiales}{\Lorenz}
\listplotHSB[unit=0.25]{tabXY aload pop}%
\uput[r](8,0){$x$}
\uput[u](0,8){$y$}
\end{pspicture}
\end{verbatim}
\subsection{Le tracé de X(t)}
\begin{center}
\begin{pspicture}(0,-5)(13,5.5)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt](0,-5)(13,5)%
\pstVerb{\parametres}%
  \psequadiff[% t,x
  tabname=tabtX,
  method=rk4,
  whichord=0,
  saveData,filename=LorenztX.dat,
  plotpoints=2501,algebraic]{0}{25}{\conditionsInitiales}{\Lorenz}
\listplotHSB[yunit=0.2,xunit=0.5,HSB=false]{tabtX aload pop}%
\uput[r](13,0){$t$}
\uput[u](0,5){$x$}
\psline{<->}(13,0)(0,0)(0,5)
\end{pspicture}
\end{center}
\begin{verbatim}
\begin{pspicture}(0,-5)(13,5)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
\pstVerb{\parametres}%
  \psequadiff[% t,x
  tabname=tabtX,
  method=rk4,
  whichord=0,
  saveData,filename=LorenztX.dat,
  plotpoints=2501,algebraic]{0}{25}{\conditionsInitiales}{\Lorenz}
\listplotHSB[yunit=0.2,xunit=0.5,HSB=false]{tabtX aload pop}%
\uput[r](12.5,0){$t$}
\uput[u](0,5){$x$}
\end{pspicture}
\end{verbatim}
\subsection{Le tracé de XZ}
\begin{center}
\begin{pspicture}(-8,0)(8,12)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
\pstVerb{\parametres}%
  \psequadiff[% x,z
  tabname=tabXZ,
  method=rk4,
  whichabs=0,
  whichord=2,
  saveData,filename=LorenzXZ.dat,
  plotpoints=2501,algebraic]{0}{25}{\conditionsInitiales}{\Lorenz}
\listplotHSB[linecolor=green,unit=0.25,HueBegin=.25,HueEnd=0.65]{tabXZ aload pop}%
\uput[r](8,0){$x$}
\uput[u](0,12){$z$}
\end{pspicture}
\end{center}
\begin{verbatim}
\begin{pspicture}(-8,0)(8,12)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
\pstVerb{\parametres}%
  \psequadiff[% x,z
  tabname=tabXZ,
  method=rk4,
  whichabs=0,
  whichord=2,
  saveData,filename=LorenzXZ.dat,
  plotpoints=2501,algebraic]{0}{25}{\conditionsInitiales}{\Lorenz}
\listplotHSB[linecolor=green,unit=0.25,HueBegin=.25,HueEnd=0.65]{tabXZ aload pop}%
\uput[r](8,0){$x$}
\uput[u](0,12){$z$}
\end{pspicture}
\end{verbatim}
\subsection{Le tracé de YZ}
\begin{center}
\begin{pspicture}(-8,0)(8,12)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
\pstVerb{\parametres}%
  \psequadiff[% y,z
  tabname=tabYZ,
  method=rk4,
  whichabs=1,
  whichord=2,
  saveData,filename=LorenzYZ.dat,
  plotpoints=2501,algebraic]{0}{25}{\conditionsInitiales}{\Lorenz}
\listplotHSB[linecolor=magenta,unit=0.25,HueBegin=0.65,HueEnd=1]{tabYZ aload pop}%
\uput[r](8,0){$y$}
\uput[u](0,12){$z$}
\end{pspicture}
\end{center}
\begin{verbatim}
\begin{pspicture}(-8,0)(8,12)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
\pstVerb{\parametres}%
  \psequadiff[% y,z
  tabname=tabYZ,
  method=rk4,
  whichabs=1,
  whichord=2,
  saveData,filename=LorenzYZ.dat,
  plotpoints=2501,algebraic]{0}{25}{\conditionsInitiales}{\Lorenz}
\listplotHSB[linecolor=magenta,unit=0.25,HueBegin=0.65,HueEnd=1]{tabYZ aload pop}%
\uput[r](8,0){$y$}
\uput[u](0,12){$z$}
\end{pspicture}
\end{verbatim}
\subsection{Le tracé du papillon de Lorenz en 3D}
Le tracé en 3D est l'étape la plus délicate, puisqu'il faut réunir les fichiers $(x,y)$, $(x,z)$ et $(y,z)$ en un seul fichier $(x,y,z)$. Ce qu'on réalise ainsi :
\begin{verbatim}
\pstVerb{/XY [(LorenzXY.dat) run ] def
         /XZ [(LorenzXZ.dat) run ] def
         /YZ [(LorenzYZ.dat) run ] def
/tabXYZ [
    0 2 XY length 2 sub {/i exch def
    [ XY i get
      XY i 1 add get
      YZ i 1 add get
    ]
    } for
   ]def
/fichierpoints (XYZLorenz.dat) (w) file def
0 1 tabXYZ length 1 sub  {/i exch def
    tabXYZ i get aload pop
    /zi exch def
    /yi exch def
    /xi exch def
      fichierpoints xi 15 string cvs writestring
      fichierpoints 32 write %% espace
      fichierpoints yi 15 string cvs writestring
      fichierpoints 32 write %% espace
      fichierpoints zi 15 string cvs writestring
      fichierpoints 10 write %% CR
} for
fichierpoints closefile
\end{verbatim}
Les dernières lignes permettent d'enregistrer, si on le souhaite, le fichier de points.
\begin{center}
\begin{pspicture}(-6,-4)(6,12)
\psset{THETA=-70,PHI=20,Dobs=2000,Decran=400}
\pstVerb{/XY [(LorenzXY.dat) run ] def
         /XZ [(LorenzXZ.dat) run ] def
         /YZ [(LorenzYZ.dat) run ] def
/tabXYZ [
    0 2 XY length 2 sub {/i exch def
    [ XY i get
      XY i 1 add get
      YZ i 1 add get
    ]
    } for
   ]def
/fichierpoints (XYZLorenz.dat) (w) file def
0 1 tabXYZ length 1 sub  {/i exch def
    tabXYZ i get aload pop
    /zi exch def
    /yi exch def
    /xi exch def
      fichierpoints xi 15 string cvs writestring
      fichierpoints 32 write %% espace
      fichierpoints yi 15 string cvs writestring
      fichierpoints 32 write %% espace
      fichierpoints zi 15 string cvs writestring
      fichierpoints 10 write %% CR
} for
fichierpoints closefile
}%
\makeatletter
\pstVerb{
 /THETA \psk@TroisD@Theta def 
 /PHI \psk@TroisD@Phi def 
 /Dobs \psk@TroisD@Dobs def 
 /DScreen \psk@TroisD@Ecran def 
 /Sin1 THETA sin def
 /Sin2 PHI sin def
 /Cos1 THETA cos def
 /Cos2 PHI cos def
 /Cos1Sin2 Cos1 Sin2 mul def
 /Sin1Sin2 Sin1 Sin2 mul def
 /Cos1Cos2 Cos1 Cos2 mul def
 /Sin1Cos2 Sin1 Cos2 mul def
 /IIID {
 0 1 tabXYZ length 1 sub {/i exch def
 /TAB tabXYZ i get def
 /abscisse TAB 0 get def
 /ordonnee TAB 1 get def
 /cote TAB 2 get def
 \formules} for
 } def
 }%
 \makeatother
\listplotHSB[unit=1]{IIID}
\psset{unit=2.5,linestyle=dashed}
\pnodeTroisD(0,0,0){O}
\pnodeTroisD(0,0,5){Z}
\pnodeTroisD(5,0,0){X}
\pnodeTroisD(0,5,0){Y}
\pnodeTroisD(-10,-10,0){A}
\pnodeTroisD(-10,-10,20){B}
\pnodeTroisD(-10,10,20){C}
\pnodeTroisD(-10,10,0){D}
\pnodeTroisD(10,-10,0){E}
\pnodeTroisD(10,-10,20){F}
\pnodeTroisD(10,10,20){G}
\pnodeTroisD(10,10,0){H}
\pspolygon(A)(B)(C)(D)
\pspolygon(E)(F)(G)(H)
\psline(A)(E)
\psline(B)(F)
\psline(D)(H)
\psline(C)(G)
\psset{linestyle=solid,linecolor=red}
\psline{->}(O)(X)
\psline{->}(O)(Y)
\psline{->}(O)(Z)
\end{pspicture}
\end{center}
\newpage
Le tracé s'effectue grâce à la commande \verb+\pnodeTroisD+ et ses diverses options.
\begin{verbatim}
\makeatletter
\pstVerb{
 /THETA \psk@TroisD@Theta def
 /PHI \psk@TroisD@Phi def
 /Dobs \psk@TroisD@Dobs def
 /DScreen \psk@TroisD@Ecran def
 /Sin1 THETA sin def
 /Sin2 PHI sin def
 /Cos1 THETA cos def
 /Cos2 PHI cos def
 /Cos1Sin2 Cos1 Sin2 mul def
 /Sin1Sin2 Sin1 Sin2 mul def
 /Cos1Cos2 Cos1 Cos2 mul def
 /Sin1Cos2 Sin1 Cos2 mul def
 /IIID {
 0 1 tabXYZ length 1 sub {/i exch def
 /TAB tabXYZ i get def
 /abscisse TAB 0 get def
 /ordonnee TAB 1 get def
 /cote TAB 2 get def
 \formules} for
 } def
 }%
 \makeatother
\listplotHSB[unit=1]{IIID}
\psset{unit=2.5,linestyle=dashed}
\pnodeTroisD(0,0,0){O}
\pnodeTroisD(0,0,5){Z}
\pnodeTroisD(5,0,0){X}
\pnodeTroisD(0,5,0){Y}
\pnodeTroisD(-10,-10,0){A}
\pnodeTroisD(-10,-10,20){B}
\pnodeTroisD(-10,10,20){C}
\pnodeTroisD(-10,10,0){D}
\pnodeTroisD(10,-10,0){E}
\pnodeTroisD(10,-10,20){F}
\pnodeTroisD(10,10,20){G}
\pnodeTroisD(10,10,0){H}
\pspolygon(A)(B)(C)(D)
\pspolygon(E)(F)(G)(H)
\psline(A)(E)
\psline(B)(F)
\psline(D)(H)
\psline(C)(G)
\psset{linestyle=solid,linecolor=red}
\psline{->}(O)(X)
\psline{->}(O)(Y)
\psline{->}(O)(Z)
\end{verbatim}
\section{Animations avec le package animate}
En préparation.
\end{document}

