How to reduce the residue in NSolve for a system of 3 non-linear equations

Question

I have written the following snippet in Mathematica to solve a system of 3 non-linear equations?

fr={-2.9*(10^(-15))/(x^4)+1.0/(x^2)-1.0/((y-x)^2)-1.0/((z-x)^2)==0,
-2.9*(10^(-15))/(y^4)+1.0/(y^2)+1.0/((y-x)^2)-1.0/((z-y)^2)==0,
-2.9*(10^(-15))/(z^4)+1.0/(z^2)+1.0/((z-x)^2)+1.0/((z-y)^2)};
fnprecise=SetPrecision[fr,32]
fnprecise[[1]]
sol=NSolve[fnprecise,{x,y,z},Reals,WorkingPrecision->2 $MachinePrecision]

But when I replace the solutions in the non-linear equations one by one; I get quite a bit large residues? Does any body know how can I increase the accuracy of my calculation?

{-2.9*(10^(-15))/(x^4)+1.0/(x^2)-1.0/((y-x)^2)-1.0/((z-x)^2)}/.sol[[1]]
{-2.9*(10^(-15))/(y^4)+1.0/(y^2)+1.0/((y-x)^2)-1.0/((z-y)^2)}/.sol[[1]]
{-2.9*(10^(-15))/(z^4)+1.0/(z^2)+1.0/((z-x)^2)+1.0/((z-y)^2)}/.sol[[1]]
{-2.9*(10^(-15))/(x^4)+1.0/(x^2)-1.0/((y-x)^2)-1.0/((z-x)^2)}/.sol[[2]]
{-2.9*(10^(-15))/(y^4)+1.0/(y^2)+1.0/((y-x)^2)-1.0/((z-y)^2)}/.sol[[2]]
{-2.9*(10^(-15))/(z^4)+1.0/(z^2)+1.0/((z-x)^2)+1.0/((z-y)^2)}/.sol[[2]]
{-2.9*(10^(-15))/(x^4)+1.0/(x^2)-1.0/((y-x)^2)-1.0/((z-x)^2)}/.sol[[3]]
{-2.9*(10^(-15))/(y^4)+1.0/(y^2)+1.0/((y-x)^2)-1.0/((z-y)^2)}/.sol[[3]]
{-2.9*(10^(-15))/(z^4)+1.0/(z^2)+1.0/((z-x)^2)+1.0/((z-y)^2)}/.sol[[3]]

the output of replacement is as follows:

{0.015625}
{0.00390625}
{0}
{0.015625}
{0.00390625}
{0}
{0}
{0}
{0.0136719}

As you can see, the residue of some replacements are still high, any hint on how to get lower residues? I mean increasing the accuracy of solutions. I want to make the residues extremely small. These are pretty big.

Thanks

Answer 1

Let's start with a system with exact coefficients:

frexact = fr /. Equal -> Subtract /. 
  x_Real :> Rationalize[10^16 (Rationalize@x)]/10^16
(*
{-(29/(10000000000000000 x^4)) + 1/x^2 - 1/(-x + y)^2 - 1/(-x + z)^2,
 -(29/(10000000000000000 y^4)) + 1/y^2 + 1/(-x + y)^2 - 1/(-y + z)^2,
 -(29/(10000000000000000 z^4)) + 1/z^2 + 1/(-x + z)^2 + 1/(-y + z)^2}
*)

The norm of the Jacobian at a point {x, y, z} gives the maximum change in the norm of the output per unit change in the norm of the input:

jac = D[frexact, {{x, y, z}}];
Norm /@ (jac /. sol) // N
(*  {2.4709*10^22, 2.4709*10^22, 1.43005*10^22, 1.43005*10^22}  *)

The norm of the solutions are around 2.5 * 10^-7 to 6 * 10^-7:

Norm /@ ({x, y, z} /. sol) // N
(*  {2.58198*10^-7, 2.58198*10^-7, 5.97*10^-7, 5.97*10^-7}  *)

Therefore, since the coefficients of fr are MachinePrecision, the smallest change the input is more than 2.5 * 10^-7 * $MachineEpsilon, and the error in the best computed solution should be less than half that. Hence, the best residual should be less than that times the norms of the Jacobian, which they are:

2.5*10^-7*$MachineEpsilon (Norm /@ (jac /. sol))/2
(*  {0.685812, 0.685812, 0.396919, 0.396919}  *)

Norm /@ (fr /. Equal -> Subtract /. sol)
(*  {0.0161059, 0.0161059, 0.0136719, 0.0136719}  *)

That does not mean that NSolve has found the best solutions -- the bounds are upper bounds. If they are not the best, it does show that NSolve has found ones that give residuals that are pretty close to the best possible.

So to get more accurate solutions, you need to use a system of equations with more accurate coefficients and use a correspondingly higher WorkingPrecision, as is done for fnprecise. You can also get the exact solutions:

solexact = Select[Solve[frexact == 0, {x, y, z}], FreeQ[N[#, 32], Complex] &]

However, N[solexact, 32] gives worse residuals to the machine-precision fr than the OP's sol. I'll leave pondering the complexities of the propagation of round-off error to the reader.

Answer 2

Try this. Use exact numbers in the equations

fr = {-(29/(10^16 x^4)) + 1/x^2 - 1/(y - x)^2 - 1/(z - x)^2 == 
   0, -(29/(10^16 y^4)) + 1/y^2 + 1/(y - x)^2 - 1/(z - y)^2 == 
   0, -(29/(10^16 z^4)) + 1/z^2 + 1/(z - x)^2 + 1/(z - y)^2 == 0}

Then solve

sol = Solve[fr, {x, y, z}] // Simplify;

very long answer. Check the first set

fr /. sol[[1]] // Simplify
(*{True, True, True}*)

So far precision is exact. Get a table of 60 precision results.

tab60 = Table[sol[[i]] // N[#, 60] &, {i, Length[sol]}];

The first 4 are real. look at the first one.

tab60[[1]]
{x -> 6.91952687538855886492684430491257331286052384388716806993577* 10^-8, 
    y -> 
-2.43974835105174540476188106423012828571010101879034468177723*10^-7, 
    z ->-4.85256661932407300745321948248699735260232051947286303702773*
        10^-8}

Check residues for the 4 real roots.

Table[Table[fr[[i, 1]] /. tab60[[n]], {i, Length[fr]}], {n, 4}]

(*{0.*10^-45    0.*10^-46   0.*10^-45
0.*10^-45   0.*10^-46   0.*10^-45
0.*10^-47   0.*10^-47   0.*10^-45
0.*10^-47   0.*10^-47   0.*10^-45
}*)

So it seems the residues are a good 15 orders of magnitude larger than the precision used, but they seem to be plenty small enough in this case. I'm not sure why NSolve doesn't do a better job.