High-Precision NSolve

Question

I need to calculate the intersection of two curves

f1[x_] := ((Zl ρ ) Exp[-x])
f2[x_] := (α k e^2 /x^2)

Where

Zl = 2.05*10^-15 
α = 1.6381
ρ = 0.326*10^-10
k = 8.9875517873681764*10^9
e = 1.602176565*10^\[Minus]19

But NSolve[f1[x]==f2[x],x,WorkingPrecision->100] returned with a warning

Solve was unable to solve the system with inexact coefficients. The \ answer was obtained by solving a corresponding exact system and \ numericizing the result.

and the result is not in the desired precision. As the precisions of k and e are too high, Rationalize can not work properly.But as f1 and f2 are two curves that are very close to each other, without a high working precision the error could be huge. How can I improve the precision in this case? (Here I can just lower the precision of k and e and let NSolve do its job, but what if I have to use all the ultra-precise numbers?)

Answer 1

Here is a case where you should take a close look at the magnitude of your quantities and do some manipulation to normalize things before throwing the system at the computer:

  Zl = 2.05*10^-15
  α = 1.6381
  ρ = 0.326*10^-10
  k = 8.9875517873681764*10^9
  e = 1.602176565*10^\[Minus]19

divide both of your expressions by ( ρ e ) :

  f1[x_] = Zl  Exp[-x] / e 
  f2[x_] = α k e/x^2/ρ 

note by the way there is no need for delayed definitions here. note also I did this by hand, you do not want an expression with e^2 / e

the expression you are solving is now fairly tame: f1[x]==f2[x]

12795.1 E^-x == 72.356/x^2

 NSolve[ f1[x] == f2[x], x]

{{x -> -0.0725216}, {x -> 0.0781981}, {x -> 9.72452}}

with a warning that there could be other solutions, but you can easily convince yourself that there are not.

 Plot[{ f1[x] , f2[x] } , {x, 5, 12}]

restoring the original definitions we get:

this result is machine precision

9.724519615727882`

the result I get rationalizing everything (ie. rationalizing each input quantity and never doing any floating operations ) and using FindRoot[ WorkingPrecision->100 ] is

9.72451961572788227332283528658...

you would be hard pressed to explain how the high precision result is somehow meaningful.

Answer 2

The code

NSolve[Rationalize[f1[x] == f2[x], 0], x, Reals], 100]

yields three solutions (with or without N), which is the minimum number of solutions if Zl ρ is positive.

The following, which sets the precision of the input to match the working precision,

NSolve[SetPrecision[f1[x] == f2[x], 100], x, Reals, WorkingPrecision -> 100]

also yields three solutions, after a while, which agree to 26 digits or so with the first code. (Using SetPrecision[f1[x] == f2[x], Infinity] gives an equivalent result much faster.)

Precision: The input is basically MachinePrecision. The OP's code, NSolve[f1[x] == f2[x], x, Reals, WorkingPrecision -> 100] (with Reals added per OP's comment), is asking for a working precision that is higher than the input precision. The output consists of three accurate solutions at MachinePrecision, which is what I would expect. But I also expected a warning about the precision mis-match, like the one you get with FindRoot. The solutions are exactly the same as those obtained without WorkingPrecision -> 100. It seems the option is ignored, but without a warning.

I do not see anything wrong with the results. Generally in Mathematica to get a high-precision output, you have to tell it explicitly that the input is high-precision. I think this is a good approach, since if the error in the input is not as small as the precision claims, misleading results occur.

One might wonder which of the two results above is more accurate. That depends on which of the codes

Rationalize[f1[x] == f2[x], 0]
SetPrecision[f1[x] == f2[x], 100]

have coefficients that are closer to the true values of the coefficients. There is no way to know that from what is given in the question. If the digits of each parameter reflect the precision to which they are known, then the machine precision results are already more precise than the input and an error analysis would be appropriate.

Answer 3

f1[x_] = ((Zl ρ) Exp[-x]);
f2[x_] = (α k e^2/x^2);

Zl = 2.05*10^-15 // Rationalize[#, 0] &;
α = 1.6381 // Rationalize[#, 0] &;
ρ = 0.326*10^-10 // Rationalize[#, 0] &;
k = 8.9875517873681764*10^9 // Rationalize[#, 0] &;
e = 1.602176565*10^\[Minus]19 // Rationalize[#, 0] &;

sol = NSolve[f1[x] == f2[x], x, WorkingPrecision -> 100];

NSolve::ifun: Inverse functions are being used by NSolve, so some solutions may not be found; use Reduce for complete solution information. >>

f1[x] - f2[x] /. sol

{0.*10^-125, 0.*10^-125, 0.*10^-129}