3
$\begingroup$

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.

Three example plots of function

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.

$\endgroup$
  • 3
    $\begingroup$ You could record that root by using WhenEvent[] in the original NDSolve[] like here $\endgroup$ – Dr. belisarius Mar 2 '15 at 19:59
  • $\begingroup$ @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. $\endgroup$ – thecommexokid Mar 2 '15 at 20:18
  • $\begingroup$ Yup, neat isn't it? $\endgroup$ – Dr. belisarius Mar 2 '15 at 20:47
  • $\begingroup$ 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. $\endgroup$ – Dr. belisarius Mar 2 '15 at 20:57
2
$\begingroup$

@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!

$\endgroup$
  • 1
    $\begingroup$ Why not use WhenEvent[f[t] == 0,...]? $\endgroup$ – Michael E2 Mar 3 '15 at 1:49
  • $\begingroup$ 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. $\endgroup$ – thecommexokid Mar 3 '15 at 16:48
0
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – thecommexokid Mar 2 '15 at 19:52
  • $\begingroup$ InterpolatingFunction::dmval: "Input value {6.6*10^7} lies outside the range of data in the interpolating function. Extrapolation will be used." $\endgroup$ – thecommexokid Mar 2 '15 at 19:52
  • $\begingroup$ 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 $\endgroup$ – george2079 Mar 2 '15 at 23:21

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.