How to find the next root larger than a specified value, numerically?

NDSolve as an equation solver

NDSolve can search for the next root, using WhenEvent.

With[{f = x^2 - 2},
  NDSolve[{
    {y'[x] == D[f, x], y[x] == f /. x -> 1/2},
    WhenEvent[f == 0, x0 = x; "StopIntegration"]},
   {}, {x, 1/2, Infinity}]
  ];
x0
(*  1.41421  *)

By default WhenEvent does not find the root as accurately as possible, but it does a very good job:

x0^2 - 2
(*  -2.73115*10^-14  *)

Strengths and limitations

Accuracy. NDSolve, since it uses local error estimates to control step size, is unlikely to skip a root. It could happen in a DAE, if the underlying ODE allowed for large step sizes (it's not entirely clear to me how error estimation works in DAE, but I have observed step sizes too large for an accurate solution of the algebraic part of the system.) The accuracy of the solution x0 can be controlled through the "LocationMethod" option, or by post-processing with FindRoot.

Multiple roots. Another issue is that WhenEvent detects roots by sign changes, so double roots would be missed. These can be found by checking extrema or by checking when the value of the function is close to zero.

Smooth functions. Dealing with differentiable functions is easier than with non-differentiable ones. The approach above is unlikely to miss a zero-crossing. If the accumulation of numerical error is an issue, one can use the "Projection" Method for NDSolve to keep error drift in check. For an even-multiple roots, one can check for extrema that are close enough to be zero to be considered zero. The ODE tends to become stiff near a multiple zero, so Method -> "StiffnessSwitching" may appropriate in this case.

Non-differentiable and singular functions. Functions that don't have a derivative, such as Abs, cannot be integrated as above, because they cannot be evaluated numerically. One can sometimes replace them with equivalents, such Abs[x] with Sqrt[x^2] for real x. Other functions, such as the cube root, have singularities that can give NDSolve trouble. In some cases, NDSolve will either be insufficient or need extra help.

General function

The function findNextRoot (code below) returns a list of the next n roots (default 1).

findNextRoot[
 eqn,          (* equation to be solved, or function whose next root is sought *)
 {var,         (* the variable on which eqn depends *)
  val,         (* starting value for search *)
  v2},         (* optional: upper limit for search; default = Infinity *)
 n,            (* optional: max. number of roots to seek; default = 1 *)
 opts          (* options *)
 ]

Options and their defaults:

Direction -> 1,               (* direction of search *)
"ZeroTolerance" -> Automatic, (* if residual is less, then a possible zero *)
"PostProcess" -> None,        (* function to apply, e.g., FindRoot *)
"WhenEventOptions" ->         (* options to pass to WhenEvent *)
  {"LocationMethod" -> {"Brent", AccuracyGoal -> 4, PrecisionGoal -> 4}},
"NDSolveOptions" -> {}}       (* options to pass to NDSolve *)

For differentiable functions, "ZeroTolerance" is 10^-8; otherwise it is 10^-6. There are two events, one to detect simple crossings and one to detect multiple roots (or non-sign-change zeros). In my (perhaps limited) experience, it is better to post-process with FindRoot than to use it through WhenEvent.

Examples

findNextRoot[x^2 - 2 == 0, {x, 1/2}]
(*  {{x -> 1.41421}}  *)

Find a root in the other direction:

findNextRoot[x^2 - 2 == 0, {x, 1/2}, Direction -> -1]
(*  {{x -> -1.41421}}  *)

Find up to 15 roots on the interval 1/2 < x < 20 (there are only six):

findNextRoot[Sin[x] == 0, {x, 1/2, 20}, 15]
(*
  {{x -> 3.14159}, {x -> 6.28319}, {x -> 9.42478},
   {x -> 12.5664}, {x -> 15.708},  {x -> 18.8496}}
*)

Find the next 15 roots:

findNextRoot[Sin[x] == 0, {x, 1/2}, 15]
(*
  {{x -> 3.14159}, {x -> 6.28319}, {x -> 9.42478},
   {x -> 12.5664}, {x -> 15.708},  {x -> 18.8496}, {x -> 21.9911},
   {x -> 25.1327}, {x -> 28.2743}, {x -> 31.4159}, {x -> 34.5575},
   {x -> 37.6991}, {x -> 40.8407}, {x -> 43.9823}, {x -> 47.1239}}
*)

Find multiple roots of Cos[x]^6 == 1. This requires some tweaking. In particular "StiffnessSwitching" and post-processing by FindRoot. Note that this equation has the same solutions as the previous example Sin[x] == 0 and that NSolve cannot solve the multiple-root equation (although Solve can).

(sol = findNextRoot[Cos[x]^6 == 1, {x, 1.5, 50}, 20
   , "ZeroTolerance" -> 1*^-6
   , "PostProcess" -> FindRoot
   , "WhenEventOptions" -> {"DetectionMethod" -> "DerivativeSign"}
   , "NDSolveOptions" -> {Method -> "StiffnessSwitching", 
      PrecisionGoal -> 10, AccuracyGoal -> 10}]) // AbsoluteTiming
NSolve[Sin[x] == 0 && 1.5 < x < 50, x] // AbsoluteTiming
NSolve[Cos[x]^6 == 1 && 1.5 < x < 50, x] // AbsoluteTiming
(*
  {0.029014, {{x -> 3.14159}, ..., {x -> 47.1239}}}
  {0.048492, {{x -> 3.14159}, ..., {x -> 47.1239}}}
  {0.085889, {}}
*)

Cos[x]^6 - 1 /. sol
(*  {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}  *)

Five points of intersection:

sols = findNextRoot[Sin[x] == 0.4 Cos[10 Sqrt[2] x], {x, 2}, 5]

Plot[{Sin[x], 0.4 Cos[10 Sqrt[2] x]}, {x, 2., 8.}, 
 Epilog -> {Red, PointSize[Medium], Point[{x, Sin[x]} /. sols]}]

(*  {{x -> 3.02062}, {x -> 3.23834}, {x -> 3.39528}, {x -> 5.9535}, {x -> 6.06944}}  *)

Non-differentiable case. Here the problem FindRoot faces is that starting below Sqrt[2] leads to trouble while starting above does not.

findNextRoot[Abs[x - Sqrt[2]], {x, 1/2}(*,"ZeroTolerance"\[Rule]2*10^-8*)]
Abs[x - Sqrt[2]] /. %

findNextRoot[Abs[x - Sqrt[2]], {x, 1/2}, "ZeroTolerance" -> 2*10^-8]
Abs[x - Sqrt[2]] /. %

FindRoot::lstol: The line search...was unable to find a sufficient decrease in the merit function.

(*
{{x -> 1.41421}}
{9.32587*10^-15}
{{x -> 1.41421}}
{2.22045*10^-16}
*)

You can verify with

FindRoot[Abs[x - Sqrt[2]], {x, Sqrt[2] - 1.*^-14}]
FindRoot[Abs[x - Sqrt[2]], {x, Sqrt[2] + 1.*^-14}]

Another non-differential case:

findNextRoot[Abs[Cos[x]] == 0, {x, 1/2}, 3]
(*  {{x -> 1.5708}, {x -> 4.71239}, {x -> 7.85398}}  *)

Code dump

The approach in findNextRoot[] is to determine whether the derivative of the residual function can be evaluated and then to set up an appropriate system. I say the function is "differentiable" when it can be evaluated. The two principle workhorses are NDSolve and WhenEvent, so the ability to pass options to them is provided via the options "NDSolveOptions" and "WhenEventOptions". In the non-differentiable case, FindRoot is sometimes called; one might want to add a way to pass options to it, too.

The function that constructs the system for NDSolve is getSysOpts. The system it constructs needs to be in terms of the local Module variables used in the NDSolve call. These could be passed as arguments, which would add ten more arguments to getSysOpts; or they could be localized as private symbols in a package context. Instead I adapted the method in NDSolve Method Plugin Framework, in which string tokens represent the variables of the system. The tokens in the code returned by getSysOpts are replaced by local variables before passing it to NDSolve.

The string tokens represent the following variables:

"X"    (* dependent variable                                                *) 
"T"    (* independent variable (for time integration)                       *)
"F"    (* the function whose root is to be found                            *)
"FP"   (* the derivative of "F"                                             *)
"XDAE" (* dependent variable for DAE integration                            *) 
"T0"   (* initial value of "T"                                              *) 
"T1"   (* temporary variable for storing result of FindRoot (non-diff case) *)  
"T2"   (* dummy symbolic variable for FindRoot (non-diff case)              *)
"CNT"  (* counts number of roots found so far                               *)
"MAX"  (* maximum number of roots to find                                   *)

The basic systems for integration are

{"X"'["T"] == "FP", "X"["T0"] == "F" /. "T" -> "T0"}          (* diff case *)
{"X"["T"] == "F", "XDAE"'["T"] == "F", "XDAE"["T0"] == "T0"}  (* non-diff case *)

In both cases "X" is integrated to be the residual function "F". In the DAE system, "XDAE" is the integral of "F". The hope here is that NDSolve will better track the features of "F". The alternative "XDAE"'["T"} == 1 tends to be faster, but the larger step sizes sometimes miss roots.

The strategies for root-finding are the following events:

TYPE OF ROOT                 DIFF CASE       NON-DIFF CASE
Sign change (oddroot):       "F" == 0        "F" == 0
No sign change (evenroot):   "FP" == 0       Abs["X"["T"]] < tol

The no-sign-change case is tricky. In each case, further processing is required. In the differentiable case, one is detecting all the local extrema, so one has to check whether they are close enough to zero be considered zero. In the non-differentiable case, one cannot use the derivative, so one can only check when the function gets close to zero, and then search for a zero. What is considered "close" is given by the value of the option "ZeroTolerance". This has limitations, of course. Basic numerical difficulties aside, one particular one is that there seems to be no easy way to bracket the root, so that FindRoot does not go on a wild goose chase. The expression for "ZeroTolerance" may be given in terms of NDSolve variables "T", "X", "F" and so forth and the appropriate localized variables. For example:

findNextRoot[Abs[Exp[-x] Sin[x]] Cos[x], {x, Pi/3}, 15, 
 "ZeroTolerance" -> 10^-1 Exp[-"T"]]

The variable x may be used instead of "T" here because "T" represents the variable passed. Overall, however, one should expect difficult cases to be rare and to require special attention by the user.

Note: The event oddroot in the non-differentiable case is often superseded by the criterion in evenroot: the residual will be less than the tolerance before it crosses zero, unless (1) the step size is large enough that the residual is not less than the tolerance or (2) the function remains below the tolerance. There might be a better design.

Full code:

ClearAll[findNextRoot, iFindNextRoot, getSysOpts];
Options[findNextRoot] = Options[iFindNextRoot] =
   {Direction -> 1,
    "ZeroTolerance" -> Automatic, (* if residual is less, then a possible zero *)
    "PostProcess" -> None,  (* function to apply, e.g., FindRoot *)
    "WhenEventOptions" -> {"LocationMethod" -> {"Brent", 
        AccuracyGoal -> 4, PrecisionGoal -> 4}},
    "NDSolveOptions" -> {}};

getSysOpts[differentiable_, opts_] := getSysOpts[differentiable,
  OptionValue[findNextRoot, opts, "ZeroTolerance"],
  OptionValue[findNextRoot, opts, "WhenEventOptions"],
  OptionValue[findNextRoot, opts, "NDSolveOptions"]]

(* differentiable f *)
getSysOpts[True, tol0_, eventopts_, ndsolveopts_] :=
  Module[{oddroot, evenroot, ode, ndsolveopts2},
   With[{tol = tol0 /. Automatic -> 1*^-8},
    ode = Hold@{"X"'["T"] == "FP", "X"["T0"] == "F" /. "T" -> "T0"};
    oddroot = WhenEvent["F" == 0,
      "CNT"++; Sow["T", "iFindNextRoot"];
      If["CNT" >= "MAX", "StopIntegration", {}],    (* found all: done *)
      eventopts];
    evenroot = WhenEvent["FP" == 0,
      If[Abs["X"["T"]] < tol,                       (* close enough to be a zero? *)
       "CNT"++; Sow["T", "iFindNextRoot"];
       If["CNT" >= "MAX", "StopIntegration", {}],   (* found all: done *)
       {}],
      eventopts];
    ndsolveopts2 = Join[ndsolveopts,
      {Method -> {"Projection",
                  "Invariants" -> {"X"["T"] == "F"},
                  Method -> "ExplicitRungeKutta"}}]];
   {{ode, oddroot, evenroot}, ndsolveopts2}];

(* non-differentiable f *)
getSysOpts[False, tol0_, eventopts_, ndsolveopts_] :=
  Module[{oddroot, evenroot, ode},
   With[{tol = tol0 /. Automatic -> 10^-6},
    ode = Hold@{"X"["T"] == "F", "XDAE"'["T"] == "F", "XDAE"["T0"] == "T0"};
    oddroot = Hold@With[{f1 = "F" /. "T" -> "T2"},
       WhenEvent["F" == 0,
        "CNT"++; Sow["T", "iFindNextRoot"];
        If["CNT" >= "MAX", "StopIntegration", {}],   (* found all: done *)
        eventopts]];
    evenroot = Hold@With[{f1 = "F" /. "T" -> "T2"},
       WhenEvent[Abs["X"["T"]] < tol,                (* close to root? *)
        "T1" = "T2" /. FindRoot[f1, {"T2", "T"}];    (* search *)
        If[NumberQ["T1"],                            (* search succeeded? *)
         "CNT"++; Sow["T1", "iFindNextRoot"];
         If["CNT" >= "MAX", "StopIntegration", {}],  (* found all: done *)
         {}],
        eventopts]]];
   {{ode, oddroot, evenroot}, ndsolveopts}];

iFindNextRoot[f_, d_, n_, {var_, val_, v2_: Automatic}, opts : OptionsPattern[]] :=
  Module[{vtmp, var2, y, daeY, diffQ, nfound, roots,
    sys, ndsolveopts, localize, post, dir},
   dir = Sign[OptionValue[Direction]];
   diffQ = NumberQ[d /. var -> N@val]; (* differentiable Q *)
   localize = # /. {"X" -> y, "T" -> var, "F" -> f, "FP" -> d, "XDAE" -> daeY, 
       "T0" -> val, "T1" -> vtmp, "T2" -> var2, "CNT" -> nfound, "MAX" -> n} &;
   {sys, ndsolveopts} = ReleaseHold@ localize@ getSysOpts[diffQ, {opts}];
   nfound = 0;
   roots = Last@ Reap[
      NDSolve[sys, {}, {var, val, v2 /. Automatic :> dir*Infinity}, ndsolveopts],
      "iFindNextRoot"] /. {r_List} :> r;
   post = OptionValue["PostProcess"] /. 
     Automatic | None -> ({First[#2] -> Last[#2]} &);
   post[f, {var, #}] & /@ roots
   ];

findNextRoot[eqn_, {var_, val_, v2_: Automatic}, n_Integer: 1, 
   opts : OptionsPattern[]] :=
  With[{dir = Sign[OptionValue[Direction]],
    f = eqn /. Equal -> Subtract},
   iFindNextRoot[
     f,           (* objective function *)
     D[f, var],   (* Jacobian *)
     n,           (* (max) number of roots to find *)
     {var, val, v2}, opts] /; Abs@dir == 1];

NSolve is pretty robust

eqn1 = x^2 - 1 == 0;

soln1 = NSolve[{eqn1, x > -1/2}, x, Reals][[1]]

(*  {x -> 1.}  *)

eqn1 /. soln1

(*  True  *)

eqn2 = Sin[1/(1 - x)] == x;

OrderedQ[x /. NSolve[{eqn2, x > 0}, x, Reals]] // Quiet

(*  True  *)

Since the roots are ordered, just take the first

soln2 = NSolve[{eqn2, x > 0}, x, Reals][[1]] // Quiet

(*  {x -> 0.599745}  *)

eqn2 /. soln2

(*  True  *)

Plot[{Sin[1/(1 - x)], x}, {x, 0, 0.95},
 Epilog -> {Red, AbsolutePointSize[5],
   Point[{x, x} /. soln2]}]

I don't know if this is a bug or just goes with the territory when doing constrained optimization; my guess is the latter. Workarounds include using FindArgMin, which has fewer method coices and basically is forced to use one that works, or to specify a nondefault Method for NArgMin. Below we see that either of these will behave better on this problem.

In[196]:= NArgMin[{x, x^2 - 1 == 0 && x > -1/2}, x, 
 Method -> "NonlinearInteriorPoint"]

Out[196]= 0.999999992549

In[198]:= FindArgMin[{x, x^2 - 1 == 0 && x > -1/2}, x]

Out[198]= {1.}