# Plot the cusp catastrophe surface

I'm trying to draw (plot) the Cusp catastrophe surface into latex using any software logic. I am completely stuck because I'm definitely not a mathematician. The surface has the potential:

V(x) = F(x,u,v) = x^4 + ux^2 + vx


so it should be given by this implicit equation:

4x^3 + 2ux + v = 0


But I don't know how to implement this. Apparently gnuplot can't do this (even if I find something approaching) so I should go with a CAS like Octave or scripting something in Maxima first, then feed gnuplot or export the graphic. Ideally, it would be great if it was as easy as this Mathematica example:

F[x_, u_, v_] :=  x^4 + u x^2 + v*x

Needs["GraphicsContourPlot3D"]

ContourPlot3D[Evaluate[D[F[x, u, v], x]],
{u, -2.5, 2}, {v, -2.2, 3}, {x, -1.4, 1.3},
PlotPoints -> 7, ViewPoint -> {-1.5, 1.5, 1.4},
Axes -> True, ContourStyle -> {EdgeForm[]},
AxesLabel -> TraditionalForm /@ {u, v, x}] // Timing


So I dig I lot, apparently, I can't do it perfectly without the ability of a mathematician (that's for sure) si I decided I should stick with an implicit solving method. I can translate the Mathematica code into Sage:

u, v, x = var('u, v, x')
f(u,v,x) = x^4 + u*x^2 + v*x
g = f.diff(x)
implicit_plot3d(g, (u, -2, 2), (v, -2, 2), (x, -2, 2))


It produces neat results:

You can also check the result on cloud.sagemath (subscription required).

But unfortunately 3d plot can't be exported in vector format, maybe there's a method for raw (x,y,z) points export, dunno... I tried an other solution with a the experimental psplotImpIIID feature of ps-tricks, it's ugly and produce a lot of artefacts.

And to end up with pstricks, one could do something like this:

\begin{pspicture}(-6,-4)(6,5)
\pstThreeDCoor[xMin=-4,xMax=4,yMin=-4,yMax=4,zMin=-4,zMax=4,RotY=90,RotX=25]
\psplotThreeD[plotstyle=line,hiddenLine,% does not work in my case (xelatex?)
yPlotpoints=50,xPlotpoints=50,linewidth=.5pt,algebraic](-3,3)(-3,3){4*x^3 - 2*y*x}
\end{pspicture}


The advantage being that we can easily rotate the axis, the disadvantage (from the perspective of a pstricks noob) being styling, plus I need to clip the result image. There's also a solution that uses a R package which I tried and and it works great and outputs vector graphics. But honestly, I think I loop the web (I even scripted it as a Blender Python) and there is definitely no easy straight solution for this.

• I honestly didn't know that mathematicians were one of the few groups of people where a major catastrophe is just another day at the office ;) – Christian Apr 1 '14 at 18:00
• You might also look into SageTeX, if you have not already. Basically this just allows you to embed the sage code in your TeX file rather than having to make it separately. – Charles Staats Apr 3 '14 at 19:39

\documentclass[pstricks]{standalone}
\usepackage{pst-solides3d}
\begin{document}

\begin{pspicture}(-4,-6)(5,9)
\psset{viewpoint=100 30 40 rtp2xyz,lightsrc=viewpoint, Decran=120}
\psSurface[ngrid=.15 .15,incolor=yellow,hue=0 1,linewidth=0.1\pslinewidth,
algebraic,axesboxed](-1,-2)(1,2){ 4*x^3 - 2*y*x}
\end{pspicture}

\end{document}


The z max/min values can be handled only on PostScript level (will be changed in the future)

\documentclass[pstricks]{standalone}
\usepackage{pst-solides3d}
\begin{document}

\begin{pspicture}(-4,-6)(5,9)
\psset{viewpoint=100 30 40 rtp2xyz,lightsrc=viewpoint, Decran=120}
\psSurface[ngrid=0.1 0.15,incolor=yellow,hue=0 1,linewidth=0.1\pslinewidth,axesboxed,
% algebraic](-2,-2)(2,2){4*x^3 -2*y*x}
](-2,-2)(1,3){ 4 x 3 exp mul 2 y mul x mul sub
dup 4 gt { pop 4 } if
dup -4 lt { pop -4 } if
}
\end{pspicture}

\end{document}


However, the values are not ignored, they are only set to the max/min.

• Cool! Do you know if there is a method to clip the surface within a desired area (especially z like zmin,zmax in tikzpicture axis)? – nadous Apr 3 '14 at 18:31
• In general yes. but it doesn't work in all cases, like in this one ... But I'll have a look into the code. – user2478 Apr 3 '14 at 19:15
• Mr Voß? I'm reading the pst-solides3d doc right now, thanks for your work. – nadous Apr 3 '14 at 19:21
• @nadous In flesh and bone :) – percusse Apr 3 '14 at 20:35

You can actually do this via pgfplots but to make it a little nicer you have to switch to LuaLaTeX since PDFLaTeX chokes up on memory limitations (due to many samples per axis. You don't need that much but why not :P).

\documentclass{standalone}
\usepackage{pgfplots}
\pgfplotsset{compat=1.10}
\begin{document}
\begin{tikzpicture}
\begin{axis}[y domain=-4:4,
domain=-2:2,
restrict z to domain=-5:5,
samples=100,
view={43}{8},
mesh/interior colormap name=hot,
colormap/blackwhite,
xlabel=$x$,ylabel=$u$,zlabel=$v$,grid=both
]
\end{axis}
\end{tikzpicture}
\end{document}


If you decide the hide the rough edges you can get a lo-fi version which LaTeX can also handle

\documentclass[border=5mm]{standalone}
\usepackage{pgfplots}
\pgfplotsset{compat=1.10}
\begin{document}
\begin{tikzpicture}
\begin{axis}[y domain=-3:3,
domain=-2:2,zmin=-3,zmax=3,
restrict z to domain=-5:5,
samples=45,
view={48}{15},
mesh/interior colormap name=hot,
colormap/blackwhite,
xlabel=$x$,ylabel=$u$,zlabel=$v$,grid=both
]
\end{axis}
\end{tikzpicture}
\end{document}


• Thanks, it's a great answer. It would be awesome to have it flipped though (oriented like in the ePiX example). – nadous Apr 2 '14 at 7:15

I think what you are asking for it is precisely in ePiX's, (http://bay.uchicago.edu/tex-archive/graphics/epix/samples/butterfly.xp) examples.

If I understand the code, the strategy is to draw the surface in three patches. It would be fun to translate this into pgfplots and asymptote. Maybe asymptote can even define this curve implicitly (it is very good at that).

Here it is an image of the result (I didn't make it):

• Thanks, I will dig this. Still, it's basically what the 15 line of gnuplot code did in the "something approaching" example I signal in my post (see image) And by approaching, I mean that this solution does not plot the entire surface, neither the ePIX solution. I believe the only proper solution will be the derivation of the potential function as set in the Mathematica example. – nadous Apr 1 '14 at 20:50
• I don't understand what do you want exactly, to derive the equation? You can increase the range of the parameters to see the entire surface. But, do you understand that by doing so one branch will cover the whole thing and you won't be able to see the fold? – alfC Apr 1 '14 at 21:06
• Honestly, as I said, I'm not a mathematician, so I may say stuff I don't get either. For now, I will just try to install epix and get the example working. – nadous Apr 1 '14 at 22:40

Asymptote can plot an implicitly defined surface using the smoothcontour3 module. This module has been incorporated into the (just-released) version 2.33 of Asymptote, but if you don't have an earlier version, you can install the module using the link above. (Just copy the file into the directory containing your image.) Here's how the surface you describe might be drawn:

\documentclass{standalone}
\usepackage{asypictureB}
\begin{document}
\begin{asypicture}{name=catastrophe_surface}
settings.outformat = "png";
settings.render = 32;
size(8cm);
import smoothcontour3;
import graph3;

currentprojection=perspective(3*(4,-6,3));

real f(real v, real u, real x) { return 4x^3 + 2u*x + v; }
draw(implicitsurface(f,(-2,-2,-2),(2,2,2), n=15, overlapedges=true),
surfacepen=material(diffusepen=gray,
ambientpen=blue));

zaxis3("$x$", Bounds, InTicks, zmin=-2, zmax=2);
xaxis3("$v$", Bounds, InTicks);
yaxis3("$u$", Bounds, InTicks(beginlabel=false));
\end{asypicture}
\end{document}


It can even produce a vector graphics version if you change the first two lines inside the asypicture environment to

settings.outformat = "pdf";
settings.render = 0;


but the resulting pdf file will be enormous and will look bad in most viewers (in my experience, it looks okay in Adobe Reader, but not in any other I've tried).

Either way, for this particular problem, I think you're better off using the small bit of math

4x^3 + 2ux + v = 0
v = -4x^3 - 2ux


to turn this into the plot of a function rather than an implicitly defined surface. Implicitly defined surfaces generally don't look as nice.

Here's a non-implicit version with Asymptote. (Note that it requires my crop3D module.)

settings.outformat = "png";
settings.render = 16;
size(8cm);

import graph3;
import crop3D;
currentprojection=perspective(3*(4,-6,3));

triple F(pair ux) {
real u = ux.x;
real x = ux.y;
return (-4x^3 - 2*u*x, u, x);
}
surface s = surface(F, (-2,-2), (2,2), nu=20, usplinetype=Spline);
draw(crop(s,(-2,-2,-2),(2,2,2)), surfacepen=material(diffusepen=gray, ambientpen=blue));

zaxis3("$x$", Bounds, InTicks, zmin=-2, zmax=2);
xaxis3("$v$", Bounds, InTicks);
yaxis3("$u$", Bounds, InTicks(beginlabel=false));


It compiles more quickly than the implicit version and gives, in my opinion, a nicer result:

• Yeah, implicit was not the better way to go, but damn, math was so long ago that I absolutely forgot how to solve an equation! – nadous Apr 3 '14 at 19:24
• Thanks for your edition, that's a very very neat result! – nadous Apr 4 '14 at 6:48
• Sorry to undig this, but it appears I can't use your module. I have this error : /usr/local/texlive/2014/texmf-dist/asymptote/crop3D.asy: 45.67: cannot cast 'triple[][][](triple[][] P, real u=<default>)' to 'triple[][][](triple[][] P)' – nadous May 12 '15 at 14:05
• @nadous: That error has been fixed in the version at tex.stackexchange.com/a/159240/484 – Charles Staats May 12 '15 at 15:03