Updating Wagon's FindAllCrossings2D[] function

Question

Stan Wagon's Mathematica in Action (second edition; I haven't read the third edition and I'm hoping to eventually see it), demonstrates a nifty function called FindAllCrossings2D[]. What the function basically does is to augment FindRoot[] by using ContourPlot[] to find crossings that FindRoot[] can subsequently polish. Here, Wagon uses the function to assist in solving one of the questions of the SIAM hundred-digit challenge.

ContourPlot[] changed quite a bit starting from version 6 (e.g., it now outputs GraphicsComplex[] objects), and FilterRules[] has superseded the old standby FilterOptions[] With these in mind, I set out to update FindAllCrossings2D[]:

Options[FindAllCrossings2D] = 
  Sort[Join[Options[FindRoot], {MaxRecursion -> Automatic, 
     PerformanceGoal :> $PerformanceGoal, PlotPoints -> Automatic}]];

FindAllCrossings2D[funcs_, {x_, xmin_, xmax_}, {y_, ymin_, ymax_}, opts___] := 
 Module[{contourData, seeds, tt, fy = Compile[{x, y}, Evaluate[funcs[[2]]]]},

  contourData = Map[First, Cases[
     Normal[
      ContourPlot[funcs[[1]], {x, xmin, xmax}, {y, ymin, ymax},
       Contours -> {0}, ContourShading -> False, 
       PlotRange -> {Full, Full, Automatic}, 
       Evaluate[
        Sequence @@ 
         FilterRules[Join[{opts}, Options[FindAllCrossings2D]], 
          DeleteCases[Options[ContourPlot], Method -> _]]]
       ]], _Line, Infinity]];

  seeds = Flatten[Map[#[[
       1 + Flatten[Position[Rest[tt = Sign[Apply[fy, #, 2]]] Most[tt], -1]]
                ]] &, contourData], 1];

  If[seeds == {}, seeds,
   Select[
    Union[Map[{x, y} /. 
        FindRoot[{funcs[[1]] == 0, 
          funcs[[2]] == 0}, {x, #[[1]]}, {y, #[[2]]}, 
         Evaluate[
          Sequence @@ 
           FilterRules[Join[{opts}, Options[FindAllCrossings2D]], 
            Options[FindRoot]]]] &, 
      seeds]], (xmin < #[[1]] < xmax && ymin < #[[2]] < ymax) &]]]

The function works splendidly, it seems. I tried out the same example Wagon used in his book:

f[x_, y_] := -Cos[y] + 2 y Cos[y^2] Cos[2 x];
g[x_, y_] := -Sin[x] + 2 Sin[y^2] Sin[2 x];

pts = FindAllCrossings2D[{f[x, y], g[x, y]}, {x, -7/2, 4}, {y, -9/5, 21/5},
                         Method -> {"Newton", "StepControl" -> "LineSearch"}, 
                         PlotPoints -> 85, WorkingPrecision -> 20] // Chop;

ContourPlot[{f[x, y], g[x, y]}, {x, -7/2, 4}, {y, -9/5, 21/5},
        Contours -> {0}, ContourShading -> False, 
        Epilog -> {AbsolutePointSize[6], Red, Point /@ pts}]

Whew, that preamble was quite long. Here's my question, then:

Are there "neater" (for some definition of "neater") ways to update/reimplement FindAllCrossings2D[] than my attempt?

Answer 1

Here is my latest code for this function, from Chapter 12 of the third edition of "Mathematica in Action". It is pretty short, but I will let you work out if it is faster or more robust than yours. Note the PlotPoints option for difficult cases.

FindRoots2D::usage = 
  "FindRoots2D[funcs,{x,a,b},{y,c,d}] finds all nontangential solutions to
   {f=0, g=0} in the given rectangle."; 

Options[FindRoots2D] = {PlotPoints -> Automatic, MaxRecursion -> Automatic}; 

FindRoots2D[funcs_, {x_, a_, b_}, {y_, c_, d_}, opts___] := Module[
  {fZero, seeds, signs, fy}, 
  fy = Compile[{x, y}, Evaluate[funcs[[2]]]]; 

  fZero = Cases[Normal[
     ContourPlot[
        funcs[[1]] == 0, 
        {x, a-(b-a)/97, b+(b-a)/103}, {y, c-(d-c)/98, d+(d-c)/102}, 
        Evaluate[FilterRules[{opts}, Options[ContourPlot]]]]], 
     Line[z_] :> z, Infinity]; 

  seeds = Flatten[(
     (signs = Sign[Apply[fy, #1, {1}]]; 
      #1[[1 + Flatten[Position[Rest[signs*RotateRight[signs]], -1]]]]) &
     ) /@ fZero, 1];
  If[seeds == {}, {}, 
     Select[
        Union[({x, y} /.
           FindRoot[{funcs[[1]], funcs[[2]]}, {x, #1[[1]]}, {y, #1[[2]]}, 
              Evaluate[FilterRules[{opts}, Options[FindRoot]]]] & ) /@ seeds, 
           SameTest -> (Norm[#1 - #2] < 10^(-6) & )], 
        a <= #1[[1]] <= b && c <= #1[[2]] <= d & ]]]

Answer 2

This is ContourPlot based but seems much shorter:

FindCrossings2D[{f_, g_}, {x_, xmin_, xmax_}, {y_, ymin_, ymax_}] := 
  {x, y} /. (FindRoot[{f[x, y] == 0, g[x, y] == 0}, {{x, #[[1]]}, 
  {y, #[[2]]}}] & /@ (ContourPlot[{f[x, y] == 0, g[x, y] == 0}, 
  {x, xmin, xmax}, {y, ymin, ymax}][[1, 1]]))

It works:

f[x_, y_] := -Cos[y] + 2 y Cos[y^2] Cos[2 x];
g[x_, y_] := -Sin[x] + 2 Sin[y^2] Sin[2 x];

pts = FindCrossings2D[{f, g}, {x, -7/2, 4}, {y, -9/5, 21/5}];

ContourPlot[{f[x, y] == 0, g[x, y] == 0}, {x, -7/2, 4}, {y, -9/5, 
  21/5}, Epilog -> {AbsolutePointSize[6], Red, Point /@ pts}]

Answer 3

Let me give a different approach. FindRoot does a good job, but maybe we can calculate the seed-points in a different way. When you want to find the common roots of $f(x,y)$ and $g(x,y)$ you can transform the problem into one equation which has the same roots $$0=f(x,y)^2+g(x,y)^2$$ The nice property here, which I will use is that the right hand side of this equation is always greater than zero. The bad thing is, that functions that don't cross zero are kind of difficult for some numerical methods, but that will be no concern here.

When we look at our function and think of it as a kind of water-pool with a very low level of water, you would get something like this

Now the idea is to go around each of those small pools which have maybe one, maybe more local zeroes in it, and start from every coast point a root-search.

That's the time where the image-processing kicks in. Our function is always positive which gives a really nice image (I inverted gray-levels):

Cutting off the image at sea-level is just a binarization of the image.

Finding the coast-line of each pool is simply implemented by an image subtraction and a dilation of the binarized image.

The complete method is therefore to raster the above function, extract all coast-pixel with image processing and run FindRoot for each coast-point.

FindCrossings2D[{f_, g_}, {x_, xmin_, xmax_}, {y_, ymin_, ymax_}, 
                 n_, threshold_] := Module[
{seeds = ImageData[ImageSubtract[Dilation[#, 1], #] &@
  Binarize[ColorNegate[Image[
         Table[f[x, y]^2 + g[x, y]^2, 
          {y, ymin, ymax, (ymax - ymin)/(n - 1.0)}, 
          {x, xmin, xmax, (xmax - xmin)/(n - 1.0)}]]], 
       threshold], "Bit"]},
  DeleteDuplicates[Last@Last@Reap[MapIndexed[
    If[#1 === 1, Sow[{x, y} /. FindRoot[{f[x, y] == 0, g[x, y] == 0}, 
     {{x, Rescale[#2[[2]], {1, n}, {xmin, xmax}]},
      {y, Rescale[#2[[1]], {1, n}, {ymin, ymax}]}}]]] &, seeds, 2]
   ], (Norm[#1 - #2] < 10.^(-6)) &]
] 

Here n is the raster-size and thresh is the binarization threshold which should be a bit smaller than 1.

f[x_, y_] := -Cos[y] + 2 y Cos[y^2] Cos[2 x];
g[x_, y_] := -Sin[x] + 2 Sin[y^2] Sin[2 x];

roots = FindCrossings2D[{f, g}, {x, -7/2, 4}, {y, -9/5, 21/5}, 400, 0.8];

This approach has clearly the disadvantage of having a fixed raster size, while ContourPlot uses adaptive sampling. Nevertheless, for raster-sizes from 200-500 and thresholds from ?-0.95 the method finds all or at least many roots.

Answer 4

You could use MorphologicalBranchPoints in combination with ContourPlot to find the seeds for the intersection points. Consider for example

f[x_, y_] := -Cos[y] + 2 y Cos[y^2] Cos[2 x];
g[x_, y_] := -Sin[x] + 2 Sin[y^2] Sin[2 x];

range = {{-7/2, 4}, {-9/5, 21/5}};

First we create a binarized, thinned image of the contour plot.

binPlot = Thinning@Binarize[Image[ContourPlot[{f[x, y], g[x, y]},
   {x, range[[1, 1]], range[[1, 2]]},
   {y, range[[2, 1]], range[[2, 2]]},
   Contours -> {0}, PlotPoints -> 30,
   ContourStyle -> White, Background -> Black, 
   PlotRangePadding -> 0, Frame -> False]]];

After applying MorphologicalBranchPoints on this image we find

intersPlot = MorphologicalBranchPoints[binPlot];
GraphicsGrid[{{binPlot, Dilation[intersPlot, 1]}}]

Then the seeds are just the rescaled positions of 1 in the ImageData of intersPlot.

seeds = DeleteDuplicates[
  Position[Reverse[ImageData[intersPlot]], 2].{{0, 1}, {1, 0}}, 
  (ChessboardDistance[#1, #2] <= 1) &];
seeds = Transpose[N@MapThread[
   Rescale[#, {1, #2}, #3] &, 
   {Transpose[seeds], ImageDimensions[intersPlot], range}, 1]];

The intersection points can then be found using FindRoot as before

crossp = {x, y} /. Quiet@FindRoot[{f[x, y] == 0, g[x, y] == 0}, 
  {x, #1}, {y, #2}] & @@@ seeds;
DeleteDuplicates[crossp, (Norm[#1 - #2] < .0001 &)]

Show[ContourPlot[{f[x, y], g[x, y]}, {x, -7/2, 4}, {y, -9/5, 21/5}, 
    Contours -> {0}, PlotPoints -> 30],
  Graphics[{PointSize[Large], Red, Point[crossp]}]
]

Answer 5

This method is based on the MeshFunctions. Detailed description is in this post:

Clear[f, g]
f[x_, y_] := -Cos[y] + 2 y Cos[y^2] Cos[2 x];
g[x_, y_] := -Sin[x] + 2 Sin[y^2] Sin[2 x];

Show[{
        ContourPlot[{f[x, y] == 0, g[x, y] == 0},
            {x, -7/2, 4}, {y, -9/5, 21/5},
            ContourStyle -> {Lighter[Brown, .7], GrayLevel[.7]}],
        ContourPlot[f[x, y] == 0,
            {x, -7/2, 4}, {y, -9/5, 21/5},
            ContourStyle -> None,
            MeshFunctions -> Function[{x, y, z}, g[x, y]],
            Mesh -> {{0}},
            MeshStyle -> Directive[Red, AbsolutePointSize[4]]
            ]
        }]

Note the amplified area, the cross point there is actually out of the specified range. Expanding the x range and using option PlotPoints -> 300, we can obtain this point:

Show[{
        ContourPlot[{f[x, y] == 0, g[x, y] == 0},
            {x, -7/2, 4.0002}, {y, -9/5, 21/5},
            ContourStyle -> {Lighter[Brown, .7], GrayLevel[.7]}],
        ContourPlot[f[x, y] == 0,
            {x, -7/2, 4.0002}, {y, -9/5, 21/5},
            PlotPoints -> 300,
            ContourStyle -> None,
            MeshFunctions -> Function[{x, y, z}, g[x, y]],
            Mesh -> {{0}},
            MeshStyle -> Directive[Red, AbsolutePointSize[4]]
            ]
        }]

Answer 6

Here is my revision of Stan Wagon's 3rd Edition function.

It is again faster and IMHO cleaner.

FindRoots2D::usage = 
  "FindRoots2D[funcs,{x,a,b},{y,c,d}] finds all nontangential solutions to
   {f=0, g=0} in the given rectangle."; 

Options[FindRoots2D] = {PlotPoints -> Automatic, 
   MaxRecursion -> Automatic};

FindRoots2D[
  funcs : {f1_, f2_}, {x_, a_, b_}, {y_, c_, d_}, opts : OptionsPattern[]
] :=
 Module[{fZero, seeds, fy = Compile[{x, y}, f2]},

  fZero =
    Cases[
      Normal @ ContourPlot[
        f1 == 0,
        {x, a - (b-a)/97, b + (b-a)/103},
        {y, c - (d-c)/98, d + (d-c)/102},
        Evaluate @ FilterRules[{opts}, Options @ ContourPlot] ],
      Line[z_] :> z,
      Infinity
   ];

  seeds = 
    Pick[Rest@#, Rest[#]Most[#]& @ Sign @ Apply[fy, #, 2], -1] & /@ fZero;

  With[{seq = FilterRules[{opts}, Options @ FindRoot]},
    Select[
      Union[
        {x, y} /. FindRoot[funcs, {x, #}, {y, #2}, seq] & @@@ Join @@ seeds,
        SameTest -> (Norm[# - #2] < 1*^-6 &)],
      a <= #[[1]] <= b && c <= #[[2]] <= d &] ]
 ]

Answer 7

Here's a way that's somewhat more efficient than FindRoots2D. I use Plot3D instead of ContourPlot, but more importantly for speed, I use ListPlot to approximate the zeros of the second function along the first. These are then polished with FindRoot as in other answers.

ClearAll[findAllRoots2D];
Options[findAllRoots2D] = Join[Options[FindRoot], Options[Plot3D]];

findAllRoots2D[{f1_, f2_}, {x_, a_, b_}, {y_, c_, d_}, opts___] :=
 Module[{f1plot, f2plot},
  f1plot = Plot3D[f1, {x, a, b}, {y, c, d},
    MeshFunctions -> {Function @@ {{x, y}, f1}},
    Mesh -> {{0}}, PlotStyle -> None,
    PlotRange -> All, BoundaryStyle -> None, Method -> Automatic,
    Evaluate@FilterRules[{opts}, Options[Plot3D]]];
  f2plot = ListLinePlot[
    Cases[Normal@f1plot, Line[pts_] :> pts[[All, {1, 2}]], Infinity],
    MeshFunctions -> {Function @@ {{x, y}, f2}},
    Mesh -> {{0}}
    ];
  Quiet[Check[
      FindRoot[{f1 == 0, f2 == 0}, {x, #[[1]], a, b}, {y, #[[2]], c, d},
        Evaluate@FilterRules[{opts}, Options[FindRoot]]], 
      Unevaluated@Sequence[], FindRoot::reged], FindRoot::reged] & /@ 
   Cases[Normal@f2plot, Point[p_] :> p, Infinity]
  ];

On the example:

f[x_, y_] := -Cos[y] + 2 y Cos[y^2] Cos[2 x];
g[x_, y_] := -Sin[x] + 2 Sin[y^2] Sin[2 x];

pts = {x, y} /. 
    findAllRoots2D[{f[x, y], g[x, y]}, {x, -7/2, 4}, {y, -9/5, 21/5}, 
     Method -> {"Newton", "StepControl" -> "LineSearch"}, 
     PlotPoints -> 85, WorkingPrecision -> 20]; // AbsoluteTiming
(*  {11.8102, Null}  *)

ContourPlot[{f[x, y], g[x, y]}, {x, -7/2, 4}, {y, -9/5, 21/5}, 
 Contours -> {0}, ContourShading -> False, 
 Epilog -> {AbsolutePointSize[6], Red, Point@pts}]

FindRoots2D takes about 7 seconds longer:

FindRoots2D[{f[x, y], g[x, y]}, {x, -7/2, 4}, {y, -9/5, 21/5},
     Method -> {"Newton", "StepControl" -> "LineSearch"}, 
    PlotPoints -> 85, WorkingPrecision -> 20]; // AbsoluteTiming
(*  {18.6097, Null}  *)

Mr. Wizard's improvement was 1 sec. faster than FindRoots2D.

Answer 8

Here is your own code cleaned up a bit. It runs about a third faster and most of the time is spent on ContourPlot.

Options[FindAllCrossings2D] = 
  Sort[Join[Options[FindRoot], {MaxRecursion -> Automatic, 
     PerformanceGoal :> $PerformanceGoal, PlotPoints -> Automatic}]];

FindAllCrossings2D[
 {func1_, func2_},
 {x_, xmin_, xmax_},
 {y_, ymin_, ymax_},
 opts___
] :=
 Module[{contourData, seeds, optsflt, fy = Compile[{x, y}, func2]},

  optsflt[fname_] := Sequence @@
    FilterRules[{opts} ~Join~ Options@FindAllCrossings2D, Options@fname];

  contourData =
    Cases[ Normal @ ContourPlot[
       func1, {x, xmin, xmax}, {y, ymin, ymax}, Contours -> {0}, 
       ContourShading -> False, PlotRange -> {Full, Full, Automatic}, 
       Method -> Automatic, Evaluate[optsflt @ ContourPlot] ],
     L_Line :> L[[1]],
     Infinity
    ];

  seeds = 
    Pick[Rest@#, Rest[#]Most[#]& @ Sign @ Apply[fy, #, 2], -1] & /@ contourData;

  Select[
   Union @ With[{seq = optsflt @ FindRoot},
     {x, y} /. FindRoot[{func1 == 0, func2 == 0}, {x, #1}, {y, #2}, seq] &
       @@@ Join @@ seeds],
   (xmin < #[[1]] < xmax && ymin < #[[2]] < ymax) &]

 ]

Answer 9

I am not sure if it is worth pointing out (it might even be mentioned in Stan Wagon’s book - I have the second edition, but as I am away from home I can’t check it) that these particular system of equations can be solved by Mathematica (or, more precisely with the help of Mathematica) exactly (by means of Reduce) so that all these initial points etc., are quite unnecessary. Here is how you do it:

eq1 = TrigExpand[g[x, y]];

eq2 = TrigExpand[f[x, y]] /. Sin[x]^2 -> 1 - Cos[x]^2;

eq = Eliminate[{eq1 == 0, eq2 == 0}, Cos[x]];

solsNonzero = Reduce[Sin[x] != 0 && eq && -5 <= y <= 5, y];

solsZero = Reduce[Sin[x] == 0 && eq && -5 <= y <= 5, y];

sols1 = {x, y} /. 
   N[{ToRules[
      Reduce[solsNonzero && 
        eq1 == 0 && -5 <= x <= 5 && -5 <= y <= 5, {y, x}]]}];

sols2 = {x, y} /. 
   N[{ToRules[
      Reduce[solsZero && f[x, y] == 0 && 
        g[x, y] == 0 && -5 <= x <= 5 && -5 <= y <= 5, {y, x}]]}];

sols = Join[sols1, sols2];

ContourPlot[{f[x, y], g[x, y]}, {x, -5, 5}, {y, -5, 5}, 
 Contours -> {0}, ContourShading -> False, 
 Epilog -> {AbsolutePointSize[6], Red, Point /@ sols}]

If you evaluate all the above you should see the already all too familiar picture.