23
$\begingroup$

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)"}]

enter image description here

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):

enter image description here

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:

enter image description here

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

$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – Rahul Feb 18 '15 at 20:14
  • 2
    $\begingroup$ 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). $\endgroup$ – Szabolcs Feb 18 '15 at 20:18
  • $\begingroup$ @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. $\endgroup$ – xslittlegrass Feb 18 '15 at 22:06
  • $\begingroup$ @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. $\endgroup$ – Szabolcs Feb 18 '15 at 23:06
  • $\begingroup$ @xslittlegrass Search for "NumericalFunction" on this site to find out what has been discovered about them. $\endgroup$ – Michael E2 Feb 19 '15 at 0:18
28
$\begingroup$

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}}]}]

enter image description here

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 }]

enter image description here

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}}]

enter image description here

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 )

$\endgroup$
  • 1
    $\begingroup$ Could you make the same plot for the NSolve/FindRoot variant as well, for comparison? $\endgroup$ – Hjulle Feb 19 '15 at 14:13
  • $\begingroup$ 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. $\endgroup$ – george2079 Feb 19 '15 at 19:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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