# Root finding in a non-rectangular solution region

If I want to find a root $$(x^*,y^*)$$ of functions $$f(x,y),g(x,y)$$ in the rectangular region $${\cal R}=[x_i,x_f]\times [y_i,y_f]$$ I can write

FindRoot[{f[x,y],g(x,y)},{x,x0,xi,xf},{y,y0,yi,yf}]

which starts searching for the solution near $$(x_0,y_0)$$. Is there any simple way to do the same thing in a non-rectangular region? I am thinking in particular about a triangular region $$y>x,x\in[0,x_f]$$.

I am relatively confident there is only one root in the region I care about, and that there are many roots outside of it, furthermore the function evaluation of $$f$$ is much more costly outside of it - I therefore don't want to just try lots of initial guesses $$(x_0,y_0)$$ and throw away the ones outside of $${\cal R}$$.

• NSolve[] has improved a lot that it should now be possible to make region-based constraints; otherwise, have you already seen this? Dec 10, 2021 at 16:53
• f is a function of 2 variables. Therefore it will probably be zero on some curve not only in one point. Dec 10, 2021 at 16:57
• @DanielHuber my mistake, I have two functions (I am minimising something, so the functions are partial derivatives of a different function). Dec 10, 2021 at 17:04
• Is NMinimize[F[x, y], {x, y} \[Element] Disk[], Method -> {"RandomSearch", "SearchPoints" -> 1}] helpful? Dec 10, 2021 at 17:33
• Multiply f[[x,y] and g[x,y] by Piecewise[{{1, x >= xi && x <= xj && y - yi >= (yj - yi)*(x - xi)/(xj - xi)}, {0, True}}]. I am not sure whether FindRoot handles discontinuous functions. Dec 10, 2021 at 20:42

If you want to use FindRoot  in any case, you can implement restrictions y > x or y < x with an additional variable a. Choose random starting values for x,y,a until all solutions are found. ( Excluded th solution at x==0)

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]; plot =
ContourPlot[{f[x, y] == 0, g[x, y] == 0, x == y}, {x, 10^-3, 2}, {y,
10^-3, 2}, PlotPoints -> 50, GridLines -> Automatic]


fr1 := (ff =
FindRoot[{f[x, y] == 0, g[x, y] == 0,
y == x - a}, {{x, RandomReal[{10^-3, 2}], 10^-3, 2}, {y,
RandomReal[{10^-3, 2}], 10^-3, 2}, {a, RandomReal[{10^-3, 2}],
10^-3, 2}}, Method -> "Secant", WorkingPrecision -> 15]; ff)

fr2 := (ff =
FindRoot[{f[x, y] == 0, g[x, y] == 0,
y == x + a}, {{x, RandomReal[{10^-3, 2}], 10^-3, 2}, {y,
RandomReal[{10^-3, 2}], 10^-3, 2}, {a, RandomReal[{10^-3, 2}],
10^-3, 2}}, Method -> "Secant", WorkingPrecision -> 15]; ff)

w1 := (While[Check[fr1, True], fr1] // Quiet; ff)

w2 := (While[Check[fr2, True], fr2] // Quiet; ff)

({x, y, a} /. Table[w1, {30}]) //
Union[#, SameTest -> (Rationalize[#1, 10^-5] ==
Rationalize[#2, 10^-5] &)] &

(*   {{1.3163, 1.29964, 0.0166549}}   *)

({x, y, a} /. Table[w2, {30}]) //
Union[#, SameTest -> (Rationalize[#1, 10^-5] ==
Rationalize[#2, 10^-5] &)] &

(*   {{0.24248, 0.510362, 0.267882}, {0.769709, 1.66914, 0.899426}}   *)


Using the system defined by cvgmt but also including the boundaries

Clear["Global*"]

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

eqns = {f[x, y] == 0, g[x, y] == 0, 0 <= x <= 2, x <= y <= 2};


The solutions are

sol1 = NSolve[eqns, {x, y}, WorkingPrecision -> 10]

(* {{x -> 0, y -> 0.4584575234}, {x -> 0, y -> 1.188862545}, {x -> 0.2424802952,
y -> 0.5103619318}, {x -> 0.7697090252, y -> 1.669135177}} *)


Verifying the solutions

(And @@ eqns) /. sol1

(* {True, True, True, True} *)


Graphically,

rgn = ImplicitRegion[{y >= x, 0 <= x <= 2}, {x, y}];

ContourPlot[{f[x, y], g[x, y]}, {x, -0.025, 2}, {y, 0, 2},
Contours -> {{0}},
RegionFunction -> Function[{x, y}, {x, y} ∈ rgn],
PlotPoints -> 200,
MaxRecursion -> 5,
Epilog -> {Red, AbsolutePointSize[4],
Point[{x, y} /. sol1]},
PlotLegends -> Placed[
Thread[{f[x, y], g[x, y]} == 0], {2/3, 1/5}]]


sol2 = FindRoot[{f[x, y] == 0, g[x, y] == 0},
{{x, #[[1]]}, {y, #[[2]]}}, WorkingPrecision -> 10] & /@
{{0, 1/2}, {0, 6/5}, {1/4, 1/2}, {3/4, 5/3}}

(* {{x -> 0, y -> 0.4584575234}, {x -> 0, y -> 1.188862545}, {x -> 0.2424802952,
y -> 0.5103619318}, {x -> 0.7697090252, y -> 1.669135177}} *)


Verifying,

(And @@ eqns) /. sol2

(* {True, True, True, True} *)


Edit

FindInstance work fine when the region doesn't contain it's boundary.

f[x_, y_] = x - y^2 Cos[y];
g[x_, y_] = -y + x*Sin[x];
reg = ImplicitRegion[{-10 < x < 10, x < y < 10}, {x, y}];
pts1 = {x, y} /.
FindInstance[{f[x, y] == 0,
g[x, y] == 0, {x, y} ∈ reg}, {x, y}, Reals, 20];
ContourPlot[{f[x, y], g[x, y]}, {x, y} ∈ reg,
PlotPoints -> 100, MaxRecursion -> 0,
Epilog -> {PointSize -> Large, Red, Point /@ pts1}]



f[x_, y_] = x - y^2 Cos[y];
g[x_, y_] = -y + x*Sin[x];
reg = ImplicitRegion[{-10 < x < 10, x < y < 10}, {x, y}];
plot = ContourPlot[{f[x, y] == 0, g[x, y] == 0}, {x, y} ∈
reg, PlotPoints -> 100, MaxRecursion -> 0, PlotRange -> All,
AspectRatio -> Automatic];
intersections =
GraphicsMeshFindIntersections[plot,
GraphicsMeshAllPoints -> False];
roots = {x, y} /.
FindRoot[{f[x, y] == 0, g[x, y] == 0}, {{x, #1}, {y, #2}}] & @@@
intersections;
Show[plot, Graphics[{PointSize[Large], Red, Point /@ roots}]]


Original

Since FindInstance or NMinimize can not work for the case as below, we have to try to use ContourPlot to draw such plot and locate the initial point.

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]; plot =
ContourPlot[{f[x, y] == 0, g[x, y] == 0}, {x, y} ∈
ImplicitRegion[{y > x, 0 < x < 2}, {x, y}], PlotPoints -> 50,
MaxRecursion -> 2, PlotRange -> All, AspectRatio -> Automatic];
pts = GraphicsMeshFindIntersections[plot,
GraphicsMeshAllPoints -> False]
FindRoot[{f[x, y] == 0, g[x, y] == 0}, {{x, #1}, {y, #2}}] & @@@ pts

(* FindInstance[{f[x,y]==0,g[x,y]==0,y>x,0<x<2},{x,y},Method-> Automatic]//N *)
(*  NMinimize[{f[x,y]^2+g[x,y]^2,y>x,0<x<2},{x,y}] *)


{x -> 0.24248, y -> 0.510362}

• With tighter bounds from the ContourPlot you can use NSolve[{f[x, y] == 0, g[x, y] == 0, 0 < x < 1/2, 0 < y < 1}, {x, y}]` Dec 12, 2021 at 1:39