Is Mathematica `ContourPlot` function really so efficient?

Question

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.?

Answer 1

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][[1]]] ,
       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 )