I am having issues finding the root of an interpolating function obtained from NDSolve. For example:

sol = NDSolve[{y''[t] == -y[t], y[0] == 1, y'[0] == 0}, y, {t, 0, 10}]

(*  {{y -> InterpolatingFunction[{{0., 10.}}, <>]}}  *)

Of course, the solution is simply $y=\cos(t)$ (this is a simpler example, my actual function doesn't have an analytic solution), but absurd things happen when I try to find roots:

FindRoot[y[t] == 0 /. sol, {t, 0}]

(*  {t -> -1.77636*10^-15}  *)

It gives me a number close to 0, which is not a root of $\cos(t)$. If I start elsewhere I get other stuff that is also wrong. For example, starting near 0.1, I get 10.99 which is not a zero either...

FindRoot[y[t] == 0 /. sol, {t, 0.1}]

(*  {t -> 10.9973}  *)

If I start the search near $\pi/2$ then I do get that answer, but I want to find the first positive zero and my actual function has a parameter that can push the zero arbitrarily close to $t=0$, so I need to start the search near 0 but then it just always returns something close to zero which is not the solution... is there any way to make this work?

My second attempt was to use FindInstance, but with interpolating functions it doesn't evaluate at all, it just returns the same thing I wrote...

FindInstance[y[t] == 0 /. sol, t]

Out[36]= FindInstance[{InterpolatingFunction[{{0., 10.}}, <>][t] == 0}, t]

And it just does that and won't give me a number no matter what I do... am I doing something wrong?

Thanks for the help

  • 1
    $\begingroup$ Show us your real function. In your simple case providing two starting values to avoid using derivatives (e.g. FindRoot[y[t] == 0 /. sol, {t, 0.5, 0.7}]) leads to the $\pi/2$ solution. On the other hand, if you start with $0$ but restrict the search domain to that over which you solved the differential equation (e.g. FindRoot[y[t] == 0 /. sol, {t, 0, 0, 10}]) you get an error explaining the weird behavior you observed ("The point {0.} is at the edge of the search region {0.,10.} in coordinate 1 and the computed search direction points outside the region."). $\endgroup$ – MarcoB Jan 27 '16 at 23:19
  • $\begingroup$ I did not get any error messages, but that makes sense. I tried changing the domain of the differential equation's solution to {t,-10,10}, however, and it is still returning 0 as a root when I search near 0. Giving it two starting points does give me a root, but since I want it to work for arbitrary parameters, there is no way to warranty that that root will be the smallest one. $\endgroup$ – Jorge Jan 28 '16 at 0:08
  • $\begingroup$ My actual problem is a system of equations: eq = {phib''[ t] == (m2*l2*ls)/(m1*l1^2 + m2*l2^2 + Ib)*(phis''[t]* Cos[phib[t] - phis[t]] + (phis'[t])^2* Sin[phib[t] - phis[t]]) - g*(m1*l1 - m2*l2 + mb*(l2 - l1)/2)/(m1*l1^2 + m2*l2^2 + Ib)* Sin[phib[t]], phis''[t] == l2/ls*(phib''[t]*Cos[phis[t] - phib[t]] + (phib'[t])^2* Sin[phis[t] - phib[t]]) - g/ls*Sin[phis[t]], phib[0] == 3*Pi/4, phis[0] == Pi/2, phib'[0] == 0, phis'[0] == 0} and I need a function that returns the value of t for which phib[t]=0 as a function of l1,l2,ls. $\endgroup$ – Jorge Jan 28 '16 at 0:11
  • $\begingroup$ Could you put the actual equations back into the main question instead of as a comment? I tried running them from the comment but got back "NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0..", so I'm guessing the cut and paste didn't work. $\endgroup$ – barrycarter Jan 28 '16 at 1:11

You're experiencing the typical and, in the simple example, expected limitations of searching for roots. The two FindRoot results are easily understood in terms of Newton's method. The best way to proceed, assuming given the example is typical, is to use WhenEvent.

sol = NDSolve[{y''[t] == -y[t], y[0] == 1, y'[0] == 0, 
   WhenEvent[y[t] == 0, firstzero = t; "RemoveEvent"]}, y, {t, 0, 10}];
(*  1.5708  *)

If firstzero is not accurate enough, use it as a starting point for FindRoot.

FindRoot[y[t] == 0 /. sol, {t, firstzero, 0, 10}]
(*  {t -> 1.5708}  *)

Addendum -- What's happening in Newton's method

Basically the first step in both of the OP's examples goes outside the domain of the solution. Stepping outside the domain of an interpolating function usually leads to unpredictable results.

FindRoot[y[t] == 0 /. sol, {t, 0, 0, 10}]

FindRoot::reged: The point {0.} is at the edge of the search region {0.,10.} in coordinate 1 and the computed search direction points outside the region. >>

  {t -> 0.}

Note the derivative is slightly positive at the initial point t == 0, indicating that the next step should be in the negative direction.

y'[0] /. sol

In the second case, the starting is just not close enough to the root.

FindRoot[y[t] == 0 /. sol, {t, 0.1, 0, 10}]

FindRoot::reged: The point {10.} is at the edge of the search region {0.,10.} in coordinate 1 and the computed search direction points outside the region. >>

  {t -> 10.}

The derivative is so small that the first step goes just beyond the domain again:

y'[0.1] /. sol
0.1 - y[0.1]/y'[0.1] /. sol

If you do not have any information about the roots, then a direct search, which is what WhenEvent is essentially doing above, seems necessary. In terms of built-in functions, NDSolve seems hard to beat when you are already using it to construct the function. One could use Plot, but it seems more likely to miss two close zero-crossings. Even if a function f[t] is obtained by other means, if it is differentiable, this should be efficient and robust:

NDSolve[{y'[t] == f'[t], y[0] == f[0], 
  WhenEvent[y[t] == 0, firstzero = t; "StopIntegration"]}, y, {t, 0, 10}];

If the derivative f'[t] cannot be calculated, then Plot is probably the best bet. Then one could use Brent's method.

Related: How to find numerically all roots of a function in a given range?


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.