8

I want to use PGF/Tikz to plot the points that satisfy the following recurrences:

\begin{align}
  q_{n+1}&=q_n+p_{n+1} \\
  p_{n+1}&=p_n+K\sin(q_n)
\end{align}

$p$ and $q$ are always taken modulo $2\pi$.

I want to plot it in 2D ($p\times q$) for certain values of K and plot it in 3D for all values of K between $0$ and $1$.

I know it has to look a little like it appears here: http://www.scholarpedia.org/article/Chirikov_standard_map but I have no idea on how to do it.

4
  • 2
    Please specify p[0], q[0] and K so that we get an idea what kind of picture you expect. It would be even better, if you could draw the plot with a computer algebra system, and than edit your question with and include a screenshot.
    – Matsmath
    May 16 '16 at 3:21
  • K=0.5 would be a good start. I have no idea on what should be p[0] and q[0] to get a plot like the one on the link on my post. I don't even know how to do this on a computer algebra system. I spend like an hour trying to do this on my nspire.
    – Gabriel
    May 16 '16 at 3:26
  • 1
    Put it on simulink with random initial conditions and get the data. Its pointless to generate this on TeX.
    – percusse
    May 16 '16 at 8:39
  • This TUGBoat article might give you a good starting point. May 16 '16 at 11:02
10

run with xelatex or latex->dvips->ps2pdf

\documentclass[border=5mm]{standalone}
\usepackage{pstricks}
\makeatletter
\define@key[psset]{}{K}[0.5]{\def\psk@K{#1 }}
\define@key[psset]{}{rDot}[0.5pt]{\pst@getlength{#1}\psk@rDot}
\psset{K=0.5,rDot=0.1pt}
\def\psChirikov{\pst@object{psChirikov}}
\def\psChirikov@i{%
  \begin@OpenObj
  \addto@pscode{
    /K \psk@K def
    /h 10 def
    0 1 h {
      /P ED
      0 1 h {
        /Q ED
        P h div TwoPi mul /p0 ED
        Q h div TwoPi mul /q0 ED
        0 1 100 {
          /i ED
          q0 TwoPi div 360 mul sin K mul p0 add /p1 ED
          p1 TwoPi gt { /p1 p1 TwoPi sub def } if
          p1 0 lt { /p1 p1 TwoPi add def } if         
          q0 p1 add /q1 ED
          q1 TwoPi gt { /q1 q1 TwoPi sub def } if
          q1 0 lt { /q1 q1 TwoPi add def } if
          q0 TwoPi div \pst@number\psxunit mul
          p0 TwoPi div \pst@number\psyunit mul
          \psk@rDot\space 0 360 arc \pst@usecolor\pslinecolor fill stroke
          /p0 p1 def /q0 q1 def
        } for
      } for
    } for
  }
  \end@OpenObj
  \ignorespaces
}
\makeatother
\begin{document}
\psset{unit=5,rDot=0.3pt}
\begin{pspicture}(0,0)(1,1)
\psChirikov
\end{pspicture}

\begin{pspicture}(0,0)(1,1)
\psChirikov[K=0.1,linecolor=red]
\end{pspicture}

\begin{pspicture}(0,0)(1,1)
\psChirikov[K=0.9,linecolor=blue]
\end{pspicture}

\begin{pspicture}(0,0)(1,1)
\psChirikov[K=5,linecolor=cyan]
\end{pspicture}
\end{document}

enter image description here

enter image description here

running it with pdflatex is also possible:

\documentclass[border=5mm]{standalone}
\usepackage{pstricks}
\usepackage{auto-pst-pdf}
\makeatletter
...

and then pdflatex --shell-escape <file>

A 3D image width K as z coordinate is also possible but I suppose it makes no real sense to use the above calculation. The next image has dK=0.1 and kMax=1, viewed by 30° horizontally and 10° vertically

enter image description here

2
  • Seems wonderful! I just have the same things to say as I said to Mark. 1- There is anyway to do this with PdfTex? My project has a lot of trouble with any other compiler. (I can always do the plot in other compiler and then import to my project, of course. But it would be nice to do it all in one file.) 2- Can you comment a little bit on the code (like: "this line does this, and this line does that...") so I can do similar stuff in other situations? Thanks
    – Gabriel
    May 16 '16 at 23:44
  • see edited answer for running with pdflatex
    – user2478
    May 17 '16 at 7:16
18

Metapost will be much much faster than TikZ or pgfplots. Previous answers (see the edit history) used lualatex for simple inclusion in to a document using luamplib (as well as lualatex with TikZ for extra slowness). The following example shows how metapost can be used with pdflatex using the gmp package, although the file must be compiled with --shell-escape.

\documentclass[border=5]{standalone}
\usepackage[shellescape,latex]{gmp}
\def\ChirikovMap#1#2#3#4#5{%
\begin{mpost}
beginfig(0)
  K := #1; h := #2; r := #3; u := #4;
  Z := 2 * 3.141592;
  pickup pencircle scaled r;
    for P = 0 upto h:
    for Q = 0 upto h:
      p0 := P / h * Z; q0 := Q / h * Z;
      for i = 0 upto 100:
        p1 := p0 + K * sind(q0 / Z * 360);
        if p1 > Z: p1 := p1 - floor(p1 / Z) * Z; fi
        if p1 < 0: p1 := (floor(abs(p1) / Z) + 1) * Z + p1; fi
        q1 := q0 + p1;
        if q1 > Z: q1 := q1 - floor(q1 / Z) * Z; fi
        if q1 < 0: q1 := (floor(abs(q1) / Z) + 1) * Z + q1; fi
        drawdot (q0 / Z * u, p0 / Z * u) withcolor #5;
        p0 := p1; q0 := q1;
      endfor
    endfor
  endfor
endfig;
\end{mpost}}
\begin{document}
\ChirikovMap{0.5}{10}{.25pt}{2cm}{red}
\ChirikovMap{0.971635}{10}{.25pt}{2cm}{0.5[green,black]}
\ChirikovMap{5}{10}{.25pt}{2cm}{blue}
\end{document}

enter image description here

5
  • Indeed. The Metapost version takes about 0.8s on my machine, the TikZ version about 10s. May 16 '16 at 18:27
  • @MarkWibrow The Metapost is incredible! I have just 2 things to say: 1- There is anyway to do this with PdfTex? My project has a lot of trouble with any other compiler. (I can always do the plot in other compiler and then import to my project, of course. But it would be nice to do it all in one file.) 2- Can you comment a little bit on the code (like: "this line does this, and this line does that...") so I can do similar stuff in other situations? Thanks
    – Gabriel
    May 16 '16 at 23:43
  • @GabrielRibeiro See the update for pdflatex. May 17 '16 at 8:24
  • It says "! Undefined control sequence.\gmp@latexpreamble ...g|documentclass[1\@ptsizept]{article}\ifx \empty \... \ChirikovMap{0.5}{10}{.25pt}{2cm}{red}" And says "Tex capacity exceeded". What is wrong?
    – Gabriel
    May 17 '16 at 20:56
  • @GabrielRibeiro No idea. It complies fine for me. May 18 '16 at 9:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.