0
$\begingroup$

EDIT: I just realized this is closely related to a question I asked two years ago, Is FindRoot wrong about its WorkingPrecision? (which is still unanswered). I don't know if the current question counts as a duplicate, but it may ask the same thing more clearly. I hope you still read this one and see if the different phrasing triggers any insights.

I have a similar question as Precision of FindRoot but slightly different. Consider getting an InterpolatingFunction f[t] from NDSolve then finding a root with FindRoot:

Clear[f] fsol = NDSolve[{f''[t] - f'[t] + f[t] - 1 == 0, f[0] == 1, f'[0] == 1}, f, {t, 0, 20}, WorkingPrecision -> $MachinePrecision]; f[t_] = Evaluate[f[t] /. fsol]; Plot[f[SetPrecision[t, Infinity]], {t, 0, 20}, PlotRange -> All] t0 = t /. FindRoot[f[t] == 10, {t, 18}, WorkingPrecision -> 50] Precision@f[t0] Precision@t0 

Output:

f[t]

18.136956334574359755720315216764489747727661914414 11.6327 50. 

The precision of my root t0 is equal to the WorkingPrecision I gave in FindRoot, despite that that is higher than the Precision of f itself, which came from NDSolve. My intuition about precision as a concept may be failing me here, but this does not seem right. My hunch is that if you use FindRoot to solve f[t]==number, then your answer for t should be no more precise than the precision of f at that t.

Even if Mathematica is technically right (is it?), I feel it's failing my purposes. There's real error introduced by the imprecision of my f[t]. Once I get t0 and plug it in (in my actual code it won't be back into f[t]!), I want it to reflect that t0 is imprecise because it was solving for an f[t] that was imprecise.

So I'd like answers to (a) Is Mathematica misusing precision here? And if not, then (b) Is there a way to make it reflect precision in a way that seems honest to me?

$\endgroup$

1 Answer 1

1
$\begingroup$

This doesn't answer your question but rather notes that the differential equation can be solved exactly.

Clear[f] eqns = {f''[t] - f'[t] + f[t] - 1 == 0, f[0] == 1, f'[0] == 1}; fsol = DSolve[eqns, f, t][[1]]; 

Verifying that the solution satisfies the equations

eqns /. fsol // Simplify (* {True, True, True} *) f[t_] = f[t] /. fsol // Simplify (* 1 + (2*E^(t/2)*Sin[(Sqrt[3]*t)/2])/ Sqrt[3] *) 

There are four solutions for f[t] == 10. These can all be found using NSolve (or Solve or Reduce) rather than FindRoot

f10 = t /. NSolve[{f[t] == 10, 0 < t < 20}, t, Reals, WorkingPrecision -> 50] (* {7.4711859387099462515468807380311790899637056073604, \ 10.843003786966755891571562660309141693335994867762, \ 14.516733262382665885039289324062217368693394531490, \ 18.136956466215141134109455148326120335487744137331} *) Plot[f[t], {t, 0, 20}, PlotRange -> {-800, 150}, Epilog -> {Red, AbsolutePointSize[6], Point[{#, 10} & /@ f10]}] 

enter image description here

Precision /@ f10 {50., 50., 50., 50.} 

However, some precision is lost in calculating the function

f /@ f10 (* {10.00000000000000000000000000000000000000000000000, \ 10.0000000000000000000000000000000000000000000000, \ 10.000000000000000000000000000000000000000000000, \ 10.00000000000000000000000000000000000000000000} *) Precision /@ % (* {48.4669, 47.602, 46.6845, 45.8028} *) 

The precision decreases with increasing magnitude of the first derivative of f

Abs[f'[#]] & /@ f10 // N (* {45.6818, 221.584, 1424.41, 8672.91} *) 
$\endgroup$
1
  • $\begingroup$ Apologies for not specifying, but the equation in my question is just a random example. My actual equation is a set of four nonlinear equations in four variables: x, px, y, py. After solving these with NDSolve, I use FindRoot to find the t0 when y[t0]=1. From this I get my desired solution, which is the pair {x[t0], px[t0]}. As far as I know FindRoot is the best tool to get t0 is this situation, but I'd be happy to be contradicted on that opinion. $\endgroup$ Commented Jul 15, 2017 at 4:33

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.