Implicit Plotting code for Mathematica
I finally got a working implicit curve tracer running in Mathematica. It has severe limitations: the user must supply a seed point, tolerance, and step size, it only traces one connected curve from a given seed point, and if the curve escapes to infinity, so will the tracer.
On the upside, the algorithm does not require any analytic information– a density field would work as well as an analytic expression; all that’s necessary is that the surface is specified by the zero-level set of some function continuously differentiable on a compact neighborhood of the surface (in particular, normals must exist at each point in this set).
Overall, it’s a good starting point for a more robust algorithm.I’ve uploaded the Mathematica notebook– it’s uncommented, but is short enough to be self-explanatory, given the overview in this post. I’ll post any updates in the same place. Here’s a sample of its work for the surface 
, with starting point 
How it works: Given a point on the surface, the outer normal is calculated using a 2nd FDD scheme (the stepsize for this FDD is one of those annoyingly necessary parameters), and the counterclockwise tangent vector is found by rotating the normal. Then a step is taken along this tangent vector (size determined by another parameter) to start looking for another point on the surface. The next surface point is determined by finding a root of the function 
starting from this candidate point, using Newton’s method. Proceeding in this manner until you return close to the original point (another parameter determined how close that it), you get samples of the curve.
Essentially the same algorithm can be used in 3D to polygonize surfaces, except instead of sampling points, you’re sampling triangles in the tangent plane, and instead of having a 1-dimensional front, the algorithm proceeds along multiple 2-dimensional fronts. I could probably adapt some of the techniques used to resolve the more complicated problems in the 3D regime to ease some of the restrictions of the code.
My reference was A marching method for the triangulation of surfaces by Hartmann. BTW, how can I get my hands on a copy of Hartmann’s paper “Numerical implicitization for intersection and Gn-continuous blending of surfaces”, without actually paying for it– or even a paper with equivalent ideas? Apparently Caltech doesn’t have it available electronically or in a paper copy.
Update While considering how to analyze the error involved with approximating an implicit curve using this algorithm, I realized that this curve tracer is pretty much looking at an implicit curve as the solution to an ODE and solving it using Euler’s method. The projection step is like a reinitialization process that helps keep down cumulative error.
August 30th, 2007 at 7:03 pm
Lack of implicit plotting in mathematical applications is something that drives me to distraction. For 2d implicit plots, i.e. to plot f(x,y)=0, you can try a hack by plotting z=f(x,y) in 3D, and then restricting the z values of the plot to a small set around zero ([-0.01,0.01] say). You should get a little strip that corresponds to the curve f(x,y)=0.
Maple has an implicitplot3d function, which is fairly impressive as long as you give it a fine enough mesh and a lot of time. Unfortunately like any other piece of proprietary software, Maple does not play well with Linux, though I have gotten it to work. Changing or upgrading distros always seems to break it.
August 31st, 2007 at 6:26 pm
I’ve always been annoyed by the lack of a 3D implicit plotter in Mathematica. But, I can see why it doesn’t have one: if Wolfram did supply a built-in plotter, the implicit (pun intended) expectation would be that it’d be robust and reliable for a wide range of implicit surfaces, and that’s a tough hurdle to surmount. It’s easier on the company to not even bother trying, and let users develop their own algorithms, because their own code, being tailored to their specific circumstances, would be more robust– and if it fails, they can’t blame it on Wolfram. I suspect that’s pretty much the reason behind similar omissions in other CASes.
September 14th, 2007 at 8:39 pm
That sounds a lot like Euler-Newton continuation for finding the curve in x that satisfies F(x) = 0…? See Introduction to Numerical Continuation Methods by Allgower and Georg.
September 14th, 2007 at 8:41 pm
(Also in the appendix of Keller’s Numerical Methods for Two Point Boundary Value Problems, something about tracing bifurcations.)
September 16th, 2007 at 4:01 pm
Thanks justin. I see you’re a fellow techer