FindRoot Domain restriction

Question

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

Answer 1

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.

Answer 2

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}} *)