Wie plottet man einen Lorenz-Attraktor?

convert -density 100 alpha off <mydoc>.pdf <mydoc>.gif
composite -stereo +0 <mydoc>.gif[0] <mydoc>.gif[36] Stereo.gif
i=1
while [ $i -lt 36 ]
do
composite -stereo +0 <mydoc>.gif[$i] <mydoc>.gif[$((i+36))] tmp.gif
convert -loop 0 Stereo.gif tmp.gif Stereo.gif
((i++))
done

% arara: latex
% arara: latex
% arara: dvips
% arara: ps2pdf
%%%  A C H T U N G    - 'kann' ein Weilchen dauern! %%%

\documentclass{article}
\usepackage{multido}
\usepackage{pst-3dplot}
\usepackage{pst-ode}
\usepackage{animate}

\begin{document}
%Lorenz’ set of differential equations
\def\lorenz{%
10*(x[1]-x[0]) | %dx/dt
x[0]*(28-x[2]) - x[1] | %dy/dt
x[0]*x[1] - 8/3*x[2] %dz/dt
}%
%
%write timeline file
\newwrite\OutFile%
\immediate\openout\OutFile=lorenz.tln%
\multido{\iLorenz=0+1}{101}{%
\immediate\write\OutFile{::\iLorenz x0}%
}%
\immediate\write\OutFile{::c,101}%
\multido{\iLorenz=102+1}{89}{%
\immediate\write\OutFile{::\iLorenz}%
}%
\immediate\closeout\OutFile%
%
\psset{unit=0.155,linewidth=0.5pt}%
\begin{animateinline}[
timeline=lorenz.tln,
controls,poster=last,
begin={\begin{pspicture}(-39,-13)(39,60)},
end={\end{pspicture}}
]{10}
%coordinate axes
\psset{Alpha=120,Beta=20}%
\pstThreeDCoor[xMax=33,yMax=33,zMax=55,linecolor=black]%
\newframe
\pstODEsetOrRestoreState{10 10 30} %initial condition
\multiframe{100}{i=0+1,rtZero=0+0.25,rtOne=0.25+0.25}{%
\pstODEsaveState{savedState}%
%compute attractor segments, stored in ‘lorenzXYZseg’
\pstODEsolve[algebraic]{lorenzXYZseg}{0 1 2}{\rtZero}{\rtOne}{26}{savedState}%
{\lorenz}%
%compute the whole attractor, stored in ‘lorenzXYZall’, with the same stepping
\ifnum\i=0\relax%
\pstODEsolve[algebraic]{lorenzXYZall}{0 1 2}{\rtZero}{\rtOne}{26}{savedState}%
{\lorenz}%
\else%
\pstODEsolve[algebraic,append]{lorenzXYZall}{0 1 2}{\rtZero}{\rtOne}{26}%
{savedState}{\lorenz}%
\fi%
\listplotThreeD[plotstyle=line]{lorenzXYZseg}%
}%
\newframe% required between two \multiframe
\multiframe{90}{rAlpha=116+-4}{% fly-around
\psset{Alpha=\rAlpha,Beta=20}%
\pstThreeDCoor[xMax=33,yMax=33,zMax=55,linecolor=black]%
\listplotThreeD[plotstyle=line]{lorenzXYZall}%
}%
\end{animateinline}
\end{document}

\documentclass[tikz]{standalone}
\usepackage{luacode,pgfplots}
\begin{luacode*}
-- Differential equation of Lorenz attractor
function f(x,y,z)
    local sigma = 3
    local rho = 26.5
    local beta = 1
    return {sigma*(y-x),-x*z + rho*x - y,x*y - beta*z}
end

-- Code to write PGFplots data as coordinates
function print_LorAttrWithEulerMethod(h,npoints,option)
    -- The initial point (x0,y0,z0)
    local x0 = 0.0
    local y0 = 1.0
    local z0 = 0.0
    -- add random number between -0.25 and 0.25
    local x = x0 + (math.random()-0.5)/2
    local y = y0 + (math.random()-0.5)/2
    local z = z0 + (math.random()-0.5)/2
    if option ~= [[]] then
        tex.sprint("\\addplot3[" .. option .. "] coordinates{")
    else
        tex.sprint("\\addplot3 coordinates{")
    end
    -- dismiss the first 100 points
    -- to go into the attractor
    for i=1, 100 do
        m = f(x,y,z)
        x = x + h * m[1]
        y = y + h * m[2]
        z = z + h * m[3]
    end
    for i=1, npoints do
        m = f(x,y,z)
        x = x + h * m[1]
        y = y + h * m[2]
        z = z + h * m[3]
        tex.sprint("("..x..","..y..","..z..")")
    end
    tex.sprint("}")
end
\end{luacode*}
\newcommand\addLUADEDplot[3][]{%
  \directlua{print_LorAttrWithEulerMethod(#2,#3,[[#1]])}%
}
\pgfplotsset{width=.9\hsize}
\begin{document}
\begin{tikzpicture}
  \begin{axis}
    % SYNTAX: Solution of the Lorenz system
    % with step h=0.02 sampled at 1000 points.
    \addLUADEDplot[color=red,smooth]{0.02}{1000};
    \addLUADEDplot[color=green,smooth]{0.02}{1000};
    \addLUADEDplot[color=blue,smooth]{0.02}{1000};
    \addLUADEDplot[color=cyan,smooth]{0.02}{1000};
    \addLUADEDplot[color=magenta,smooth]{0.02}{1000};
    \addLUADEDplot[color=yellow,smooth]{0.02}{1000};
  \end{axis}
\end{tikzpicture}
\end{document}