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. – Rahul Feb 18 '15 at 20:14
• I think these plotting functions make use of Experimental`NumericalFunction for optimizations. ContourPlot3D may need millions of evaluations and still works reasonably fast when used with formulae (not numerical blackboxes). – Szabolcs Feb 18 '15 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. – xslittlegrass Feb 18 '15 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. – Szabolcs Feb 18 '15 at 23:06
• @xslittlegrass Search for "NumericalFunction" on this site to find out what has been discovered about them. – Michael E2 Feb 19 '15 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? – Hjulle Feb 19 '15 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. – george2079 Feb 19 '15 at 19:23