# Numerical root-finding: Find leftmost root in a domain

This is surely a duplicate, but I cannot seem to find a similar question.

I use NDSolve to find an InterpolatingFunction y[t] which solves a differential equation with some initial conditions.

I have a function f[t] defined in terms of y[t] and some other stuff.

I consider the function f[t] over a particular domain {t, 6000000, 30000000}. The function always exhibits the same qualitative behavior over this domain: it decays toward 0 to the left, and oscillates between 0 and 3 to the right. However the location in the domain where it transitions from one behavior to the other is unknown.

I would like to identify the first time f[x] == 1 in my fixed domain. Below I show some examples indicating what the function might look like and which solution I wish to identify.

Obviously a naive use of FindRoot is not going to find the particular solution I am interested in, but yet it doesn't seem to me as though this ought to be a very hard problem. I only require 3 significant digits of precision or so.

I use Mathematica 9.0.1.

Note: Edited to include slightly more detail about the NDSolve portion of the problem, since it turned out to be relevant.

• You could record that root by using WhenEvent[] in the original NDSolve[] like here Mar 2, 2015 at 19:59
• @belisarius Wow, that's a really neat tool that I didn't know anything about! I think I have a working solution using that tip. Mar 2, 2015 at 20:18
• Yup, neat isn't it? Mar 2, 2015 at 20:47
• I suggest posting your working solution as an answer (once you have one :) ). It is indeed a nice trick that hasn't been explored too much in the site. Mar 2, 2015 at 20:57

@belisarius recommended I make use of the WhenEvent[] function in the original NDSolve[] that produced my interpolating function.

Since my NDSolve command finds y[t], but I am looking for solutions to an equation involving f[t], which is not defined until later, I cannot implement this suggestion in a completely straightforward fashion.

However, it happens to turn out that the maxima in f[t] always occur at the zeroes of y[t].

So I use WhenEvent to make NDSolve record the zeroes of y[t] as it solves:

zeroes = {};
NDSolve[{y'[t]==stuff, y[0]==a, y'[0]==b, WhenEvent[y[t]==0, AppendTo[zeroes, t]]}, y, {t, 0, 60000000}]


Now I know that the solution I want — and no others! — occurs between the left edge of my domain and zeroes[[1]].

So I use FindRoot with two initial times, one at the left edge of the domain, and the other just a smidgen to the left of the first maximum in f:

solution = t /. FindRoot[f[t] - 1, {t, 6000000, zeroes[[1]]-1000}]


Since f[t] is monotonic between those two points, FindRoot will stay in-bounds and Newton's method is sure to converge. Hooray!

• Why not use WhenEvent[f[t] == 0,...]? Mar 3, 2015 at 1:49
• Because f[t] is written in terms of some things that are not yet defined when I run the NDSolve, and it would require significant refactoring to change that. Which I could have done, but the solution I posted here worked perfectly well. Mar 3, 2015 at 16:48

Use a starting point to the left of the desired solution:

FindRoot[f[t], {t,0}]


or

FindRoot[f[t], {t,10^7}]


of course whatever is the leftmost point defined for your interpolating function.

• This was my initial thought as well, but if the starting point is indicated too far to the left of the desired solution, FindRoot seems to take steps which are too big and can wind up somewhere in the oscillating region, or even overshooting the right side of the domain entirely. Mar 2, 2015 at 19:52
• InterpolatingFunction::dmval: "Input value {6.6*10^7} lies outside the range of data in the interpolating function. Extrapolation will be used." Mar 2, 2015 at 19:52
• right as you see no guarantee FindRoot will find the closest root to where it starts. The {x,xstart,xmin,xmax} form might be useful, but since you are using NDsolve anyway WhenEvent is clearly the best approach Mar 2, 2015 at 23:21