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

I would want to have a general-purpose, reasonably robust method of finding the next numerical root above a specific value of x. I'm stumped by the fact NArgMin fails spectacularly even on this very simple example:

NArgMin[{x, x^2 - 1 == 0 && x > -1/2}, x, Reals]

During evaluation of NArgMin::nosat: Obtained solution does not satisfy
the following constraints within Tolerance -> 0.001: {1-x^2==0}. >>

(* -0.5 *)


Correct answer here is, of course, 1. All Methods seem to produce the same result. (I am actually a bit concerned it even produces a result, since it's very clearly wrong.) FindRoot isn't of help either, at least I can't understand how it would be.

Symbolic methods are, of course, "the solution" in this case. Problems with them arise when there are infinite amount of non-trivial roots, say in the case of Sin[1 / (1 - x)] == x && x > 0. Yet, on most of these cases location of the "next" root is almost glaringly obvious. How to solve this, preferably by not explicitly implementing a root-searching algorithm in own code?

• Try NMinimize[{x, x^2 - 1 == 0 && x > -1/2}, x, Method->"NonlinearInteriorPoint"]. I realize this might not be documented. Alternatively use FindMinimum or FindArgMin. – Daniel Lichtblau Mar 3 '16 at 21:12
• @DanielLichtblau Please make this an answer! Why FindMinimum and FindArgMin work when NArgMin doesn't? Also, should it actually be considered a bug that NArgMin returns a value that clearly contradicts with specified constraints? – kirma Mar 3 '16 at 21:20
• NArgMin should work with that Method option setting. As for the failure, I don't know the answer but will ask. – Daniel Lichtblau Mar 3 '16 at 21:25
• @DanielLichtblau Method -> "NonlinearInteriorPoint" seems to be the most robust approach to get an answer, and the intended root with couple other problematic cases I'm playing with. – kirma Mar 3 '16 at 21:27
• Its important to realize none of these methods operate by first identifying the discrete space defined by your constraints. The "NonlinearInteriorPoint" method re-evaluates the objective function over 5000 times to get your solution. – george2079 Mar 3 '16 at 22:44

### 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];

• This is cool and probably closest to using mechanics which make most sense for the problem at hand, although I wished a bit less code would be sufficient. It takes a while to understand! :) – kirma Mar 7 '16 at 7:11
• @kirma Thanks. I'll add some explanation later. I just ran out of time. getSysOpts returns pseudo-code that is converted into actual code via localize. I used some of the same string names as NDSolveSolutionDataComponent (see reference.wolfram.com/language/tutorial/NDSolveStateData.html, reference.wolfram.com/language/tutorial/NDSolvePlugIns.html) in the hopes that it would be familiar to some. Of course I had to add some variables. Another approach would be to share variables between the function in a package. – Michael E2 Mar 7 '16 at 12:39
• @kirma Updated with explanation and slight improvements to code. – Michael E2 Mar 8 '16 at 15:59

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 didn't remember NSolve can return a partial set like that! Now I wonder how it judges how many solutions to list on a non-complete set... – kirma Mar 4 '16 at 5:22
• Amusingly this approach breaks down with roots like Sin[x] == 1/2. Quick workaround for this case is Min @@ ((x /. NSolve[{Sin[x] == 1/2, x > 0}, x, Reals]) /. ConditionalExpression[v_, c_] :> MinValue[{v, c}, Union@Cases[c, C[_], \[Infinity]]]) - are there more robust methods? – kirma Mar 4 '16 at 5:43
• ... and now that I think of it, periodic roots were the reason why I didn't use NSolve in the first place. Maybe it's still the least problematic option. – kirma Mar 4 '16 at 6:05

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.}


Inspired by @BobHanlon's answer and FindRoot, I developed my own findNextRoot:

ClearAll[findNextRoot];

Options[findNextRoot] = {Direction -> 1};

findNextRoot[eqn_, {var_, val_}, OptionsPattern[]] :=
Module[{select, compare, findvalue, sol},
{select, compare, findvalue} =
If[OptionValue[Direction] >= 0,
{Min, Greater, NMinValue}, {Max, Less, NMaxValue}];
sol = NSolve[{eqn, compare[var, val]}, var, Reals];
(((var /. # & /@ sol) /.
ConditionalExpression[v_, c_] :>
Root[{Function[var, eqn /. Equal -> Subtract],
findvalue[{v, c},
Union@Cases[c,
C[_], \[Infinity]]]}]) /. {nonempty__} :> {var ->
select[nonempty]}) /; Head@sol === List]


It works for (rational) polynomials with starting points of "wrong" slope:

With[{eq1 = x^2, eq2 = 1},
Plot[{eq1, eq2}, {x, -2, 2},
Epilog -> {PointSize@Large,
Point[{x, eq1}] /. findNextRoot[eq1 == eq2, {x, -1/2}]}]]


With[{eq1 = (-15 + 2 x (-104 + x (-93 + 22 x)))/
(10 + 20 x), eq2 = 0},
Plot[{eq1, eq2}, {x, -2, 6},
Epilog -> {PointSize@Large,
Point[{x, eq1}] /. findNextRoot[eq1 == eq2, {x, 1}]}]]


With Direction, one can also search roots backwards:

With[{eq1 = (-15 + 2 x (-104 + x (-93 + 22 x)))/
(10 + 20 x), eq2 = 0},
Plot[{eq1, eq2}, {x, -2, 6},
Epilog -> {PointSize@Large,
Point[{x, eq1}] /.
findNextRoot[eq1 == eq2, {x, 1}, Direction -> -1]}]]


Cyclic roots are also handled gracefully:

With[{eq1 = Sin[x], eq2 = 1/2},
Plot[{eq1, eq2}, {x, 0, 2},
Epilog -> {PointSize@Large,
Point[{x, eq1}] /. findNextRoot[eq1 == eq2, {x, 0}]}]]


... and an infinite amount of untrivial roots:

With[{eq1 = Sin[1/(1 - x)], eq2 = x},
Plot[{eq1, eq2}, {x, 0, 1},
Epilog -> {PointSize@Large,
Point[{x, eq1}] /. findNextRoot[eq1 == eq2, {x, 0}]}]]


With[{eq1 = x Sin[1/x], eq2 = x/2},
Plot[{eq1, eq2}, {x, -1, 1},
Epilog -> {PointSize@Large,
Point[{x, eq1}] /.
findNextRoot[eq1 == eq2, {x, 1/3}, Direction -> -1]}]]


When no root is found, empty set is returned like in the case of FindRoot:

findNextRoot[x == 1, {x, 2}]

{}


This version seems to also work equally well using either NSolve or Solve. With Solve, one can get nice symbolic Roots.

• Counterexamples to these succeeding roots are welcome. – kirma Mar 4 '16 at 18:28
• One counterexample is essentially the same as Solve[Sin[1/(1 - E^(1 - x))] == 1 - x && x > 0, x, Reals], for which (N)Solve seems to fail to find a root. – kirma Mar 4 '16 at 20:25
• use an interval: NSolve[Sin[1/(1 - E^(1 - x))] == 1 - x && 0 < x < 0.99, x][[1]] – Bob Hanlon Mar 17 '16 at 17:59