# Algorithms for (adaptive?) function plotting

I am looking for algorithms to draw standard 2d-graphs for functions that may or may not have singularities. The purpose is to write a "Mini-CAS", so I have no a priori knowledge of the types of functions, the users want to graph.

This problem is very old, so I imagine there must be some standard algorithms in the literature. For once I did not have much succes finding references via Google.

I did find one interesting algorithm, namely this one from "YACAS - Book of algorithms" named "Adaptive function plotting".

So in short:

• Is there standard algorithms?
• Is there a test suite for known difficult-to-plot functions?
• What are the interesting papers to read?
-
Maybe the question would be better understood with "function plotting" instead of "graph drawing"? I misinterpreted the title at first (graph theory). – astrojuanlu May 30 '12 at 7:36
@Juanlu001 Thanks for the suggestion. I changed the title. – soegaard May 30 '12 at 10:49
When you say 2D, do you mean plotting a one-variable function such as $f(x)$, or are you also interested in a two-variable function ($f(x,y)$) shown in 2D with e.g. different colours/shades representing different values? – Szabolcs May 31 '12 at 8:58
Well, I meant plotting a function of one variable. However, I'd also like to hear about algorithms for choosing which points to evaluate in a two-variable setting too. I am not that interested in hearing about colours and shading. – soegaard May 31 '12 at 9:22
For 2D functions, see my question and answer here. What I did there was quite limited, and won't work that well for arbitrary function. Also, there are some essential steps missing from the description without which the method won't converge properly: I needed to insert a new sampling point in the middle of each edge of the mesh which would disappear on the next retriangulation. (contd.) – Szabolcs May 31 '12 at 9:27

I have implemented Mathematica's adaptive sampling routine here on GitHub (It's a single C file, go up to the source tree for the header file). I found a description of the routine in a big book on Mathematica a long time ago, and I've been using variations on this implementation for some time now. It basically does a rough linear sample over the domain of interest, then goes back through to refine regions of high curvature. It is possible that some very sharp features are missed, but in practice I find this exceedingly rare. This file also contains the parallel version.

-
Which book is that? Is it the one I linked? Do you know what exactly changes in their implementation between versions 5 and 6? – Szabolcs May 31 '12 at 9:51
@Szabolcs: No, I believe it was in this book section 4.1.3. The description applies a very old version of Mathematica. The newer versions (maybe starting from v6) detect vertical asymptotes and remove the spurious vertical lines from the plots. The new versions definitely do a lot of sophisticated symbolic preprocessing to deal with discontinuities, undefined regions, and branch cuts. – Victor Liu May 31 '12 at 9:55
The symbolic preprocessing you are talking about is called "exclusion detection" in the documentation. It can be turned off either by Exclusions -> None or by hiding the structure of your function from Plot by defining it as f[x_?NumericQ] := .... This is not what I was referring to when I asked about the changes. I believe there were some changes to the algorithm, as v5 and v6 sampled at different points. Right now I can't test on v5 though to compare again. – Szabolcs May 31 '12 at 9:58
Thanks for the code sample and the reference. – soegaard Jun 2 '12 at 11:07
"The Mathematica Graphics Guidebook" contained a very good discussion of the problem. I especially liked that the short comings of the algorithm was also described. – soegaard Jun 5 '12 at 20:55

To my knowledge, Mathematica uses a variation on the following basic algorithm for plotting a one-variable function $f(x)$ or a parametric curve $(x(t), y(t))$ (I'm going to assume $f(x)$ for this description).

1. Start with a regularly spaced grid of points on the plotting domain. （In Mathematica, there's a parameter to control how many to take, called PlotPoints.)

2. Look at each pair of successive line segments (defined by three successive points, $(x_1, f(x_1)),\; (x_2, f(x_2)),\; (x_3, f(x_3))$) and insert a new sampling point in the middle of both segments ($\frac{x_1+x_2}{2}$ and $\frac{x_2+x_3}{2}$) if their angle is larger than a threshold.

3. If we haven't reached the iteration limit yet (set by MaxRecursion in Mathematica), repeat from step 2.

Some of this is discussed in the book Mathematica in Action by Stan Wagon, which you can see here on Google Books.

I implemented this algorithm before to have better control over how many times my expensive to calculate function was evaluated. Here's the Mathematica code for step 2:

nd[{points_, values_}] :=
Transpose@{(Drop[points, 1] + Drop[points, -1])/2,
Differences[values]/Differences[points]}

subdivide1d[result_, resolution_, maxAngle_: 10] :=
Module[
{deriv, angle, dangle, pos, nf},
deriv = nd[result\[Transpose]];
angle = ArcTan[#2] & @@@ deriv;
dangle = Differences[angle];
pos = Flatten@Position[dangle, d_ /; Abs[d] > maxAngle/180 Pi];
pos = Union[pos, pos + 1];
nf = Nearest[result[[All, 1]]];
Select[deriv[[pos, 1]], Abs[# - First@nf[#]] > resolution &]
]

-

The MathWorld web page on Function Graphs contains references to several papers which seem to be relevant on adaptative function plotting. Quoting the page:

Good routines for plotting graphs use adaptive algorithms which plot more points in regions where the function varies most rapidly (Wagon 1991, Math Works 1992, Heck 1993, Wickham-Jones 1994). Tupper (1996) has developed an algorithm [...]

On the other hand, on Google I stumbled upon a paper

www.cs.uic.edu/~wilkinson/Publications/plotfunc.pdf

that explains how to properly choose the domain and other things. I hope they are useful to you.

-
Thanks for the references. – soegaard Jun 2 '12 at 11:03