# Is Mathematica ContourPlot function really so efficient?

I wanted to find the roots of the function $f(x,y)=\sin(3.2x)\sin(1.3y)-2.1 \sin(1.3x)\sin(3.2y)$. This is what the function looks like:

f[x_, y_] = Sin[3.2 x]*Sin[1.3*y] - 2.1*Sin[1.3*x]*Sin[3.2*y]
Plot3D[f[x, y], {x, 0, 20}, {y, 0, 20}, PlotPoints -> 50, MeshFunctions -> {#3 &}, Mesh -> 1, MeshStyle -> {Red, Thick}, AxesLabel -> {"x", "y", "f(x,y)"}] Then, as I was looking for the roots (the $\color{red}{\text{red}}$ curves, actually), I started to do use some root-finding procedures, using FindRoot and perturbating manually the initial conditions. This worked OK, but I faced some problems (unequal roots density along a curve, missing some parts, etc.). Also, the computations were taking about 10-20 seconds with my procedure (surely not optimal).

Just to give an idea of what I'm talking about, this is an example of the result I had (with different parameters, but I does not matter) (the $\color{green}{\text{green}}$ points are the results of my computations and took 10 seconds to calculate): Then I switched my brain on and realized that Mathematica had already calculated many roots to plot the red curve above. So I tried:

plot = ContourPlot[f[x, y] == 0, {x, 0, 10}, {y, 0, 10}, PlotPoints -> 50, ContourStyle -> Red]


which yields: plot[[1,1]] contains more than 7000 points calculated in less than a second. The worst root Map[f, plot[[1, 1]]] // Max // Abs gives 0.01 corresponding to a "poor" accuracy, but using PlotPoints -> 50, MaxRecursion -> 7 lowered this upper bound to 0.0002 on 85000 points in 9 seconds, which remains very acceptable.

Question Is Mathematica really more efficient with ContourPlot than with other root-finding numerical functions (hard to believe)? Would it be possible to have the same efficiency using FindRoot or NSolve, etc.?

• You could do something like Table[{x, y} /. NSolve[f[x, y] == 0 && 0 <= y <= 10, y, Reals], {x, 0, 10, 0.1}]. The difference in speed is because NSolve tries to get all roots with machine-precision accuracy, about $10^{-16}$, whereas ContourPlot is content with much lower accuracy and possibly missing some roots.
– user484
Feb 18, 2015 at 20:14
• I think these plotting functions make use of ExperimentalNumericalFunction for optimizations. ContourPlot3D may need millions of evaluations and still works reasonably fast when used with formulae (not numerical blackboxes). Feb 18, 2015 at 20:18
• @Szabolcs where is the source code of this package located at? I looked up the folder AddOns-Packages-Experimental but the folder only contains a PacletInfo.m file. Feb 18, 2015 at 22:06
• @xslittlegrass It's not a package, it's an undocumented builtin function used for efficient computation of numerical functions. Or something like that. I'm not sure. Maybe Oleksandr knows more. Feb 18, 2015 at 23:06
• @xslittlegrass Search for "NumericalFunction" on this site to find out what has been discovered about them. Feb 19, 2015 at 0:18

maybe this will provide a little insight:

first look at the evaluation points used by ContourPlot:

 f[x_?NumericQ, y_?NumericQ] := (Sow[{x, y}];
Sin[3.2 x]*Sin[1.3*y] - 2.1*Sin[1.3*x]*Sin[3.2*y]);
{plot, dat} = Reap[ContourPlot[f[x, y] == 0, {x, 0, 2}, {y, 0, 2},
PlotPoints -> 50, ContourStyle -> Red]];
Row[{
Show[{Graphics@Point@dat , plot }],
Show[{Graphics@Point@dat , plot }, PlotRange -> {{1, 1.5}, {1, 1.25}}]}] what you see is ContourPlot recursively refines the plot near the contours, but only fairly coarsely, and then evidently does an interpolation to render the contour. ( on the right you see few if any eval points are actually on the curve )

Plot3D it seems is even cruder:

 {plot2, dat2} =
Reap[Plot3D[f[x, y], {x, 0, 2}, {y, 0, 2}, PlotPoints -> 50,
MeshFunctions -> {#3 &}, Mesh -> 1, MeshStyle -> {Red, Thick}]];

Show[{Graphics@Point@dat2 , plot }] here we see the evaluation grid is not refined at all in a search for the contours..

## edit

pulling the 2D contour lines out of the 3d plot is a bit of a trick:

 t2 = Cases[Normal[Cases[ plot2 , _GraphicsComplex  , Infinity][]] ,
Line[x_] :> Line[ x[[All, 1 ;; 2]] ], Infinity]
Show[{Graphics@Point@dat2 , plot, Graphics@{Blue, td} ,
Graphics[{PointSize[.01], Green,
Point[Table[ {x, y} /. FindRoot[ f[x, y]   , {y, 1} ]  ,
{x, .1, 2, .1}]]}] }, PlotRange -> {{.01, 1.99}, {.01, 1.99}}] you see a noticeable difference in the contours generated by the two methods. (green markers from FindRoot indicate the Contour result is quite good, while the Plot3D result is less accurate )

• Could you make the same plot for the NSolve/FindRoot variant as well, for comparison? Feb 19, 2015 at 14:13
• NSolve works very well for this example (see Rahul's comment ). It relies on analytic manipulation, so does not work with my _NumericQ` requirement, and so we cant make a direct comparison in terms of number of function evaluations. Feb 19, 2015 at 19:23