0
$\begingroup$

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))?

$\endgroup$
2
$\begingroup$

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}

$\endgroup$
1
$\begingroup$

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}  *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.