# Updating Wagon's FindAllCrossings2D[] function

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? - Could someone come up with a Google Books link that works in Europe? – Szabolcs Jan 20 '12 at 13:04 @Szabolcs What about this here books.google.de/… – halirutan Jan 25 '12 at 2:29 As a tiny note: one unexpected benefit of the ContourPlot[] approach is that one can exploit the RegionFunction option if one is only interested in roots within a given region. – Guess who it is. Jan 25 '12 at 10:35 @J.M. I just posted a new version using ContourPlot[] - seems very short. – Vitaliy Kaurov Jan 26 '12 at 6:08 I noticed, @Vitaliy; sadly I can't upvote again... – Guess who it is. Jan 26 '12 at 6:09 ## 8 Answers 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 & ]]]  - Hi Dr. Wagon! Thanks for posting your updated code (and for the book in general; I hope I can get a copy of your new one some day); I'll do the tests later for this! I've also found that RegionFunction can be profitably used here to restrict the search to roots within a domain. – Guess who it is. Jan 26 '12 at 4:36 BTW: is there anything special about the increments chosen in {x, a-(b-a)/97, b+(b-a)/103}, {y, c-(d-c)/98, d+(d-c)/102}? – Guess who it is. Jan 26 '12 at 5:36 I'm afraid I cannot remember exactly why I chose 97 and 103, but the reason is related to the effect that when things are too symmetric then there can be problems at the centers of symmetry. So it is best to add in a little asymmetry. Also I note that someone commented that some new Mathematica functions do some of these sorts of things.. but only for nice functions. If one had the Zeta function involved, say, then I do not believe their new algebraic ideas would work. Of course, since typical usages might be in DEs, finding equil. pts., it might be worth adding some cases to detect niceness. – stan wagon Jan 26 '12 at 23:06 Oh. Would you happen to have a set of equations on hand where the asymmetry in ContourPlot[] is a necessity? – Guess who it is. Jan 26 '12 at 23:09 Not on hand. But this function is used many many times in the VisualDSolve project to find equilibrium points for two autonomous DEs and I must have found some difficult examples there. Now, it is possible as Mathematica changes internally that some of the problems disappear. Let me mention that this function was able to solve one of the "100 Digit Challenge" problems in optimizing a complicated function of 2 vbles by finding all of the critical points in the region, and there were well over 2000 of them. – stan wagon Jan 27 '12 at 15:21 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. - Sure; I was just looking for a test function for FindAllCrossings2D[] and I used Stan Wagon's example. I however use it for more complicated transcendental equations, and the proper analysis is just a liiitle bit out of reach... but thanks for giving the analytical solution; a shame I am out of votes for today. – Guess who it is. Jan 26 '12 at 16:54 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]] ] }]  - This is very nice (+1)! A side-question: what do you use for the zoomed area? Is it within Mathematica or externally? – gpap Apr 22 '14 at 14:07 @gpap Thanks. It's the crayon style in the Paint in Windows 8 :) But I believe it can be implemented in MMA. – Silvia Apr 22 '14 at 14:09 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 &] ] ]  - 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}]  - Nice work! I'm amazed by so short codes. – yulinlinyu May 17 '12 at 14:27 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]}] ]  - 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) &]

]

-

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.

-
Interesting approach. When I tried FindCrossings2D[{f, g}, {x, -7/2, 4}, {y, -9/5, 21/5}, 400, 0.8] on my system, it only got 59 out of the 67 roots within the region. (Increasing the value of the fourth argument gradually eventually yielded the roots, but takes quite a bit of time to execute.) I guess choosing appropriate parameters can be system-dependent... –  Guess who it is. Jan 25 '12 at 10:22
For the record, the new-in-8 functions ContourDetect (and CrossingDetect) would come handy with this approach. –  Matthias Odisio Jan 27 '12 at 13:02