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?
  • 2
    $\begingroup$ Maybe the question would be better understood with "function plotting" instead of "graph drawing"? I misinterpreted the title at first (graph theory). $\endgroup$
    – user782
    Commented May 30, 2012 at 7:36
  • $\begingroup$ @Juanlu001 Thanks for the suggestion. I changed the title. $\endgroup$
    – soegaard
    Commented May 30, 2012 at 10:49
  • $\begingroup$ 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? $\endgroup$
    – Szabolcs
    Commented May 31, 2012 at 8:58
  • $\begingroup$ 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. $\endgroup$
    – soegaard
    Commented May 31, 2012 at 9:22
  • 1
    $\begingroup$ 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.) $\endgroup$
    – Szabolcs
    Commented May 31, 2012 at 9:27

4 Answers 4


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.

  • 1
    $\begingroup$ Which book is that? Is it the one I linked? Do you know what exactly changes in their implementation between versions 5 and 6? $\endgroup$
    – Szabolcs
    Commented May 31, 2012 at 9:51
  • 1
    $\begingroup$ @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. $\endgroup$
    – Victor Liu
    Commented May 31, 2012 at 9:55
  • $\begingroup$ 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. $\endgroup$
    – Szabolcs
    Commented May 31, 2012 at 9:58
  • $\begingroup$ "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. $\endgroup$
    – soegaard
    Commented Jun 5, 2012 at 20:55
  • $\begingroup$ I can no longer find the GitHub file, did it move? $\endgroup$
    – Andrei
    Commented Jun 21, 2014 at 11:51

Knowing how other CAS do this might help you.

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,

subdivide1d[result_, resolution_, maxAngle_: 10] :=
    {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


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


I found this topic and thought I should share the developer issue page for adding this to the Julia library Plots.jl. We tried a bunch of techniques to see what would give good results, starting from the notes on Mathematica's implementation. Adding some pruning, a small perturbation to not start exactly at the interval endpoints, a recursion limit, and a dual mesh error estimator were all necessary to "get it just right". The thread also points you to the open source code for the implementation. So it did take a bit of tweaking, but adding those features made it pretty robust (according to the tests, as shown in the thread).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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