How can I avoid a singular Jacobian when solving a system of equations with FindRoot?

Question

I am trying to solve to following equation with findRoot.

rho0 = 1.0; rho1 = 1.0; rho2 = 10^(-10); 
fr = 
  {-(rho0 rho1)/(x^6) + (rho0 rho1)/(x) - (rho1 rho2)/(y - x) == 0,
   -(rho0 rho2)/(y^6) + (rho0 rho2)/y + (rho1 rho2)/(y - x) == 0};
frprecise = SetPrecision[fr, 32]
sol = 
  FindRoot[frprecise,{{x,0.015}, {y,0.8}},
    AccuracyGoal -> 20,
    PrecisionGoal -> 20, 
    WorkingPrecision -> 2 $MachinePrecision, 
    MaxIterations -> 2000000]
{-(rho0 rho1)/(x^6) + (rho0 rho1)/(x) - (rho1 rho2)/(y - x)} /. sol

No matter what values I take for the initial points, I still get a singular Jacobian. Would somebody kindly tell me how can I avoid the singular Jacobian and actually solve tis system to low residues (less than 10^(-30))?

Answer 1

I believe that you cannot avoid a singular Jacobian because the Jacobian is singular at the solution you find.

However, you can obtain more precision by using exact numbers for the inputs and by setting a higher working precision.

rho0 = 1; rho1 = 1; rho2 = 10^(-10);
fr = {-(rho0 rho1)/(x^6) + (rho0 rho1)/(x) - (rho1 rho2)/(y - x) == 
    0, -(rho0 rho2)/(y^6) + (rho0 rho2)/y + (rho1 rho2)/(y - x) == 
    0};
sol = FindRoot[fr, {{x, 15/1000}, {y, 8/10}}, WorkingPrecision -> 50, PrecisionGoal -> 30, MaxIterations -> 1000]
{-(rho0 rho1)/(x^6) + (rho0 rho1)/(x) - (rho1 rho2)/(y - x)} /. sol

{-2.3492174793046040979*10^-30}

Answer 2

Both NSolve[N@fr, {x, y}] and NSolve[frprecise, {x, y}, WorkingPrecision -> 32] give 30 complex solutions and no real solutions. To use FindRoot to search for a real solution is highly unlikely to be fruitful.

The error (singular Jacobian) in FindRoot arises from FindRoot pursuing a zero by going out toward infinity. The magnitude of the functions decrease roughly as 1/Norm[{x, y}]. Eventually a point is reached at which numerical error diminishes the rank of the Jacobian from 2 to 1.

sol = FindRoot[frprecise, {{x, 0.015}, {y, 0.8}}, AccuracyGoal -> 20, 
  PrecisionGoal -> 20, WorkingPrecision -> 2 $MachinePrecision, 
  MaxIterations -> 2000000]
(*  {x -> 1., y -> 591.241}  (FindRoot::jsing error omitted) *)

jac = D[frprecise /. Equal -> Subtract, {{x, y}}];
MatrixRank[jac /. sol]
(*  1  *)

Note the solution is reported in MachinePrecision even though higher precision was requested (maybe it's a bug). One can get the requested high precision by giving high precision initial values. Here I up the precision to 64, to get a residual less the the OP's requested goal.

frprecise = SetPrecision[fr, 64];
sol = FindRoot[frprecise, {{x, 0.015`64}, {y, 0.8`64}}, 
  AccuracyGoal -> 20, PrecisionGoal -> 20, 
  WorkingPrecision -> 4 $MachinePrecision, MaxIterations -> 2000000]
(*  {x -> 0.9999999999999999..., y -> 3.5708094605541272691...*10^26}  (error omitted) *)

jac = D[frprecise /. Equal -> Subtract, {{x, y}}];
RowReduce[jac /. sol]
(*  {{1, 1.5685437702672624548...*10^-64}, {0, 0}}  *)

{-(rho0 rho1)/(x^6) + (rho0 rho1)/(x) - (rho1 rho2)/(y - x)} /. sol
(*  {-2.80049*10^-37}  *)

We can get a much smaller residual by hand by going out further toward infinity, which does not mean it's getting closer to a root:

{-(rho0 rho1)/(x^6) + (rho0 rho1)/(x) - (rho1 rho2)/(y - x)} /.
   {x -> 1, y -> 10^100}
(*  {-1.*10^-110}  *)