How do I find all the solutions of three simultaneous equations within a given box?

Question

Sometimes, one needs to find all the solutions of three simultaneous nonlinear equations in three unknowns

$$\begin{align*}f(x,y,z)&=0\\g(x,y,z)&=0\\h(x,y,z)&=0\end{align*}$$

within a given cuboidal domain; that is, all triples $(x,y,z)$ satisfying the three equations given above, and within the region defined by $x_{\min}\leq x\leq x_{\max}$, $y_{\min}\leq y\leq y_{\max}$, $z_{\min}\leq z\leq z_{\max}$. (I restrict the discussion here to transcendental equations; algebraic equations are not too problematic for Mathematica (Solve[]/NSolve[], Resultant[], GroebnerBasis[]...))

How can I use Mathematica to find these solutions? FindRoot[] can only find one solution, and you still need an approximate location as a starting point for FindRoot[]. NSolve[] works (sometimes), but it takes long.

Answer 1

Stan Wagon's approach of using contour plotting to find solutions to simultaneous equations, as embodied in his function FindAllCrossings2D[], can be generalized to the three-dimensional case. The key is to use ContourPlot3D[] to generate space curves representing the intersection of two of the three functions, as embodied by Maxim Rytin's method (given here). Having done so, one can now use a strategy similar to what was done in FindAllCrossings2D[]; that is, evaluating the third function along the space curve(s), and locating sign changes. The locations of these sign changes can then be used as starting points for FindRoot[].

Here, then, is the routine FindAllCrossings3D[], whose use is completely analogous to its two-dimensional counterpart:

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

FindAllCrossings3D[funcs_?VectorQ,
                   {x_, xmin_, xmax_}, {y_, ymin_, ymax_}, {z_, zmin_, zmax_}, opts___] := 
 Module[{contourData, seeds, tt, fz = Compile[{x, y, z}, Evaluate[funcs[[3]]]]}, 
  contourData = Cases[Normal[ContourPlot3D[
      Evaluate[Most[funcs]], {x, xmin, xmax}, {y, ymin, ymax}, {z, zmin, zmax}, 
      BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> {}}, 
      ContourStyle -> None, Mesh -> None, Method -> Automatic, 
      Evaluate[Sequence @@ FilterRules[Join[{opts}, Options[FindAllCrossings3D]], 
         Options[ContourPlot3D]]]]], Line[l_] :> l, Infinity];
  seeds = Flatten[Pick[Rest[#], 
                  Most[#] Rest[#] &@Sign[Apply[fz, #, 2]], -1] & /@ contourData, 1];
  If[seeds === {}, seeds, 
   Select[Union[Map[{x, y, z} /. FindRoot[funcs, Transpose[{{x, y, z}, #}], 
         Evaluate[Sequence @@ FilterRules[Join[{opts}, Options[FindAllCrossings3D]], 
            Options[FindRoot]]]] &, seeds]],
          (xmin < #[[1]] < xmax && ymin < #[[2]] < ymax && zmin < #[[3]] < zmax) &]]]

As an example of how to use FindAllCrossings3D[]:

sols = FindAllCrossings3D[
        {Sin[x + y] Sin[y - z], Cos[x] Cos[y] - Sin[z], x^2 + y^2 + z^2 - 9},
        {x, -4, 4}, {y, -4, 4}, {z, -4, 4}]
     {{-2.80293, -0.756176, -0.756176}, {-2.78082, -0.360773, -1.06625},
      {-2.11276, 2.11276, 0.269309}, {-1.14056, -0.395145, 2.74645},
      {-1.14056, 2.74645, -0.395145}, {-0.883563, 0.883563, 2.72739},
      {-0.360773, -2.78082, -1.06625}, {0.360773, 2.78082, -1.06625},
      {0.883563, -0.883563, 2.72739}, {1.14056, -0.395145, 2.74645},
      {1.14056, 2.74645, -0.395145}, {2.11276, -2.11276, 0.269309},
      {2.78082, 0.360773, -1.06625}, {2.80293, -0.756176, -0.756176}}

The routine found $14$ solutions. To visualize the solutions, we can do the following:

l1 = Cases[Normal[ContourPlot3D[{Sin[x + y] Sin[y - z], Cos[x] Cos[y] - Sin[z]},
                                {x, -4, 4}, {y, -4, 4}, {z, -4, 4}, 
                                BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> {}}, 
                                ContourStyle -> None, Mesh -> None]],
           Line[l_] :> l, Infinity];

Graphics3D[{Line[l1], Sphere[{0, 0, 0}, 3], Sphere[sols, 1/10]}, 
 Axes -> Automatic]

where we used small spheres to mark the intersections of the space curves formed by the intersection of $\sin(x+y)\sin(y-z)=0$ and $\cos\,x\cos\,y=\sin\,z$, and the sphere $x^2+y^2+z^2=9$.