# FindRoot Domain restriction

I am solving a system of four non-linear equations in four variables using FindRoot. I have some sense of the relationship between the variables so I don't want Mathematica to do its computations in certain funky regions of the domain. For example, I want to solve:

\begin{align} g_1(w,x,y,z) &=0\\ g_2(w,x,y,z) &=0\\ g_3(w,x,y,z) &=0\\ g_4(w,x,y,z) &=0\\ \end{align}

such that

\begin{align} x &\in& \{25,75\} \\ y &\in& \{0,80-f_1(x)\} \\ z &\in& \{x+f_2(x),800\} \\ w &\in& \{y,z\} \end{align}

and I would like my FindRoot to operate only in the above domain. How can I get it to do that?

Simply writing this doesn't seem to work:

FindRoot[{g1==0,..,..},{x,50,25,75},{y,a,0,80-f1(x)},{z,b,x+f2(x),800},{w,c,y,z}]

-
FindRoot[] has no explicit support for restrictions like that, I'm afraid. What do your $g_k(w,x,y,z)$ and $f_k(x)$ look like? Maybe something could be tailored to your problem... –  Ｊ. Ｍ. Sep 6 '12 at 2:10
Can you redefine your functions such that, for instance, they take some improbable (or complex) value outside the ranges you are interested in ? –  b.gatessucks Sep 6 '12 at 5:54
oh that sucks. My $g_k$ are very very large expressions involving $w,x,y,z$ raised to the powers of $k_i$ where $k_i$ are roots of a quartic equation and derivatives of such terms. They're huge and messy. Some of my domain restrictions, $f_i$ are linear and others I guess bound by linear relations. –  Amatya Sep 6 '12 at 9:05

Lets use FindMinimum

Another way to attack this problem is to transform this as a constrained nonlinear optimization problem. Given the equations $$f_1(x_1,..,x_n)=0,\cdots,f_k(x_1,..,x_n)=0\\$$ and constraints $$c_1(x_1,..,x_n)\bowtie0,\cdots,c_m(x_1,..,x_n)\bowtie0\\$$ where $\bowtie \in \{<,>,\leq,\geq\}$. We form the following objective function $$\Gamma(x_1,\cdots,x_n)=\sum_{i=1}^{k}{f_i^2(x_1,\cdots,x_n)}$$ Now it is simple to feed this problem to MMA superfunction FindMinimum. We use the examples of kguler

{g1, g2, g3} = {x y - z^3 + x, x y z - 2, x^2 + y^2 + z^3 - 5};
subdomain = -1 < x < 1 && -3 < y < x - x^3 && y < z < x;
res = FindMinimum[{Total[{g1, g2, g3}^2], subdomain}, {x, y, z},AccuracyGoal -> 11]


{5.08374*10^-26, {x -> 0.832027, y -> -2.32619, z -> -1.03335}}

This are indeed very good solutions as the residuals for the equations $g_1=0,\cdots,g_3=0$ are seen in the following to be practically close to zero.

{g1, g2, g3} /. res[[2]]


{1.92957*10^-13, 1.02585*10^-13, -5.50671*10^-14}

The constraints are fulfilled too

Boole[subdomain] /. res[[2]]


1

How good is FindRoot

On the other other hand if we use FindRoot the solution is not enough robust. It's convergence depends solely on the random choice of the initial guess. You can see the residual error for each equations for a run of $300$ times here. Black dots are the places where each of the three equations are satisfied with residual norm ($<10^{-8}$).

diffList = Transpose@{{0, 0, 0}};
neweqns = 1 - Boole[subdomain] + Boole[subdomain] {g1, g2, g3};
Monitor[For[i = 1, i <= 300, i++,
diff = (Transpose@{{g1, g2, g3}}) /.
Quiet@FindRoot[neweqns,
Transpose[{{x, y, z}, {RandomReal[{-1, 1}], RandomReal[{-3, 2}],
RandomReal[{-2, 0}]}}]];
val = MapThread[Join[#1, #2] &, {diffList, diff}];
diffList = val;
],
Quiet@ListLinePlot[diffList, Mesh -> All,
MeshStyle-> Directive[PointSize[.007], Red]]]


For a run of $300$ times we got a favorable solution using FindRoot $66$ times. For this example nonlinear system with a run of $10000$ the probability of finding a solution using FindRoot turned out to be a unimpressive $18.22\%$. So when it comes to robustness above constrained optimization trick using FindMinimum may not be a bad option.

-
PlatoManiac!!!! Thank you! I had thought about going the route of Lagrange multipliers and Kuhn-Tucker but I feared that I would be over complicating things but your fabulous and comprehensive post has changed my mind. I had not thought about the robustness issue at all. hmm...I think I'll try both FindRoot and FindMinimum for starters to compare what I get. Btw is it possible to accept two answers for a question here? Yours and Kgulers? –  Amatya Sep 8 '12 at 21:33
Yes, very interesting! I'm actually enrolled in a nonlinear optimization course this fall so it's nice to see an example. Where did your AccuracyGoal come from and why was the value of 11 selected? Does this type of optimization have a specific name I could reference or do you have a link I could look over? Thanks! –  john Sep 10 '12 at 14:58
@john I guess when we call FindMinimum without any explicit Method option it chooses some kind of a Newton solver in the background. In this case as we have a symbolic objective so MMA probably takes the advantage of easily constructing a closed form Jacobian and if needed the Hessian for its internal iterative method. Given that I used the AccuracyGoal option to force the algorithm a precision of 11 digits. This helps often the algorithm to converge to the real optimum which is zero. Yo can check more on this option here (reference.wolfram.com/mathematica/ref/AccuracyGoal.html). –  PlatoManiac Sep 10 '12 at 16:57

Using a simple example with three equations in three unknowns,

{g1, g2, g3} = {x y - z^3 + x, x y z - 2, x^2 + y^2 + z^3 - 5};
subdomain = -1 < x < 1 && -3 < y < x - x^3 && y < z < x;


one can modify the functions in the equations so that the modified functions coincide with the original functions when the constraints are satisfied. One of many possible ways doing this is to use Boole (or some other Piecewise function) as follows:

neweqns = 1 - Boole[subdomain] + Boole[subdomain] {g1, g2, g3};


Without the domain restrictions,

FindRoot[{g1, g2, g3},
Transpose[{{x, y, z},{RandomReal[{-1, 1}], RandomReal[{-3, 2}], RandomReal[{-2, 0}]}}]]


possibly after several trials perturbing the initial points, gives one of the four roots:

(* {x->-0.5950097551413702,y->2.6060720302882157,z->-1.2897914523282767} *)
subdomain /. %
(* False *)


The roots of neweqns are the roots of {g1, g2, g3} that also satisfy the constraints in subdomain. Again after possibly several trials with different initial points,

FindRoot[neweqns,
Transpose[{{x, y, z}, {RandomReal[{-1, 1}], RandomReal[{-3, 2}], RandomReal[{-2, 0}]}}]]


gives

(* {x -> 0.832027, y -> -2.32619, z -> -1.03335} *)
subdomain /. %
(* True *)


Update: If you have to use FindRoot, it is better to re-initialize automatically when a message is issued. For example:

 init = {RandomReal[{-1, 1}], RandomReal[{-3, 2}],RandomReal[{-2, 0}]};
While[err == Quiet@Check[sol = FindRoot[neweqns, Transpose[{{x, y, z}, init}]], err],
init = {RandomReal[{-1, 1}], RandomReal[{-3, 2}], RandomReal[{-2, 0}]};];
sol


As always, there are alternative approaches like reformulating the problem as constrained optimization and using FindMinimum as suggested by PlatoManiac.

NSolve is another alternative that allows use of constraints directly:

NSolve[{g1, g2, g3} == {0, 0, 0}, {x, y, z}, Reals]
(* {{x -> 0.832027, y -> -2.32619, z -> -1.03335},
{x -> 1.98191, y -> -1.25951, z -> -0.801208},
{x -> -3.01335, y -> 0.409805,  z -> -1.61958},
{x -> -0.59501, y -> 2.60607, z -> -1.28979}} *)
subdomain /. %
(* {True, False, False, False} *)

NSolve[subdomain && {g1, g2, g3} == {0, 0, 0}, {x, y, z}, Reals]
(* {{x -> 0.832027, y -> -2.32619, z -> -1.03335}} *)
`
-
damn! that's so cool! Thanks you so much! I'll comment again after I implement it. –  Amatya Sep 6 '12 at 18:19