# FindRoot constrained to maximum

I am using FindRoot to find the local maximum op some function, by examining where its derivative is zero. However, keeping the starting value the same, sometimes FindRoot returns the local minimum of my function, as the derivative is, obviously, also zero there.

How do I require FindRoot to look whether the derivative is positive (local minimum) or negative (local maximum), before accepting a particular root?

For example:

Plot[{x^3 - 4 x, 3 x^2 - 4}, {x, -3, 3}]
FindRoot[3 x^2 - 4, {x, 1}]
FindRoot[3 x^2 - 4, {x, -1}]


The real problem:

Plot[{(-9-90 s-160 s^3+s^2 (-16+Sqrt[-(((9+16 s^2) (-9-180 s+320 s^3-16 s^2 (1+200 Log[5])+1600 s^2 Log[2 \[Pi] (1/4+9/(64 s^2))]))/s^4)]))/(320 s^2)},{s,0,3}]
FindRoot[(-9-90 s-160 s^3+s^2 (-16+Sqrt[-(((9+16 s^2) (-9-180 s+320 s^3-16 s^2 (1+200 Log[5])+1600 s^2 Log[2 \[Pi] (1/4+9/(64 s^2))]))/s^4)]))/(320 s^2),{s,1}]
FindRoot[(-9-90 s-160 s^3+s^2 (-16+Sqrt[-(((9+16 s^2) (-9-180 s+320 s^3-16 s^2 (1+200 Log[5])+1600 s^2 Log[2 \[Pi] (1/4+9/(64 s^2))]))/s^4)]))/(320 s^2),{s,1/10}]

• Yes, it has to be FindRoot. It is an equation which cannot be solved analytically, nor via NSolve, as it is not a polynomial. Commented Mar 19, 2014 at 23:11
• reference.wolfram.com/mathematica/tutorial/… The first paragraph suggested it couldn't and it doesn't work in my case. So, my question remains: Can FindRoot be asked to take the derivative at the root into account? Commented Mar 19, 2014 at 23:20
• Ok, just wanted to be sure that you've tried.
– Kuba
Commented Mar 19, 2014 at 23:27
• Do you mind posting the function? (not the derivative) Commented Mar 20, 2014 at 0:04
• Can you use FindMaximum instead of FindRoot? Or you only have the derivative, not the function itself? Commented Mar 20, 2014 at 0:44

You can use WhenEvent with NDSolve. It's a little slower to integrate up to the root you're looking for than to use FindRoot directly, but NDSolve finds the root with FindRoot (basically) once it detects a zero-crossing. You can see they agree on the answer with 30 digits of precision. The condition f'[s] < 0 in WhenEvent detects a crossing from positive to negative, so it will detect a local maximum. Use f'[s] > 0 to detect a local minimum. (You could even find both in the same command, if desired.) See the documentation for a similar example. The only trick is to find a starting point to the left of all the roots in the interval of integration.

df = (-9 - 90 s - 160 s^3 +
s^2 (-16 +
Sqrt[-(((9 + 16 s^2) (-9 - 180 s + 320 s^3 -
16 s^2 (1 + 200 Log[5]) +
1600 s^2 Log[2 π (1/4 + 9/(64 s^2))]))/
s^4)]))/(320 s^2);

{f0, {smax}} =
Reap@NDSolveValue[{f'[s] == df, f[1/10] == 0,
WhenEvent[f'[s] < 0, Sow[s]; "StopIntegration"]},
f, {s, 1*^-8, 3},
WorkingPrecision -> 30];

smax
FindRoot[df, {s, 1}, WorkingPrecision -> 30]
(*
{1.18047643840721719526394598817}
{s -> 1.18047643840721719526394598817}
*)

• This seems the easiest, and most theoretically satisfying method. Thank you. Commented May 3, 2014 at 16:35
• Hey Michael, could you help me generalise this to two dimensions? Is that even possible? For example: NDSolveValue[{D[f[sH,sL],sH]==3-sH,D[f[sH,sL],sL]==1-sL,f[1,sL]==0,f[sH,1]==0,WhenEvent[D[f[sH,sL],sH]<0&&D[f[sH,sL],sL]<0,mySow={sH,sL};"StopIntegration"],WhenEvent[sH>10&&sL>10,mySow=0;"StopIntegration"]},f,{sH,1,11},{sL,1,11},WorkingPrecision->6] Commented Jul 14, 2014 at 20:08
• @LBogaardt My answer basically integrates the gradient but we know the path (the real number line). In 2D+, we don't know the path. It's easy to construct a path that follows the gradient from any given point, but the path is not guaranteed to lead to the global maximum. If you're lucky enough to know there's only one maximum, then you might run into a saddle point, which is a singularity of the gradient field, but that is improbable.... Commented Jul 15, 2014 at 1:02
• @LBogaardt ...Example: grad = {3 - sH[t], 1 - sL[t]}; laststep = {1, 0}; sol = NDSolve[{{sH'[t], sL'[t]} == grad, {sH[0], sL[0]} == laststep, WhenEvent[Norm[{sH[t], sL[t]} - laststep] < 10^-6, "StopIntegration"]}, {sH, sL}, {t, 0, Infinity}, WorkingPrecision -> 6, StepMonitor :> (laststep = {sH[t], sL[t]})]. Commented Jul 15, 2014 at 1:03
• Cool! Thanks :) Commented Jul 15, 2014 at 13:38

For the equation you posted you can find a good starting point manually. I assume that is not your actual problem and that your actual problem might contain a parameter (so the good starting point may depend on the parameter),

Someone asked me to solve a similar problem recently. The problem was to find a root of -1 + x + a x Log[x] that is smaller than 1 (where a < -1). 1 is always a root. A random starting point in FindRoot may converge to either the desired root or to 1, depending on the value of a. By plotting the function we see that it has a single maximum:

Manipulate[Plot[-1 + x + a x Log[x], {x, 0, 2}], {a, -5, -1}]


If we choose the starting point to the left of the maximum, FindRoot will converge to the desired solution.

So what I did was find the maximum ...

Normal@Solve[D[-1 + x + a x Log[x], x] == 0, x]


... and use it for determining the starting point:

root[a_?NumericQ] := FindRoot[-1 + x + a x Log[x], {x, 0.9 E^((-1 - a)/a)}]


You might try something similar, but you'd need to find the maximum numerically, for example

f[s_] = (-9 - 90 s - 160 s^3 +
s^2 (-16 +
Sqrt[-(((9 + 16 s^2) (-9 - 180 s + 320 s^3 -
16 s^2 (1 + 200 Log[5]) +
1600 s^2 Log[2 \[Pi] (1/4 + 9/(64 s^2))]))/
s^4)]))/(320 s^2)

argmax = s /. Last@FindMaximum[{f[s], s > 0}, {s, 1}]

FindRoot[f[s], {s, 0.9 argmax}]

(* ==> {s -> 0.297118} *)

FindRoot[f[s], {s, 1.1 argmax}]

(* ==> {s -> 1.18048} *)


This exploits the fact that your function has a single maximum for s>0.

Now that I understood your starting point, I'll suggest another approach:

Since you only have the derivative, you can still work with one of the maximization functions by doing a numerical integral of the derivative first. This seems clumsy at first, but it turns out to be quite fast to evaluate, and easy to write:

df = (-9 - 90 s - 160 s^3 +
s^2 (-16 +
Sqrt[-(((9 + 16 s^2) (-9 - 180 s + 320 s^3 -
16 s^2 (1 + 200 Log[5]) +
1600 s^2 Log[2 Pi (1/4 + 9/(64 s^2))]))/
s^4)]))/(320 s^2);

ff[x_?NumericQ] := NIntegrate[df, {s, 1, x}]

FindMaximum[ff[x], {x, 1, 3}]

(* ==> {0.00710499, {x -> 1.18048}} *)


Here I did two things to avoid any attempts at symbolic processing: first, define the integral ff only for numeric arguments, and second use FindMaximum with two limits for the search variable x, which causes it to not do symbolic derivatives. You could also just specify one single starting value, as in FindMaximum[ff[x], {x, 1}], but that's slightly slower.

• Thank you so much. However, the maximization problem I used to illustrate my problem was really just an example. For the function I am working with, I only have the derivative in analytical form. So I still really need to use FindRoot. Commented Mar 19, 2014 at 23:29
• @LBogaardt OK, thanks for clarifying your problem. I changed my answer to account for this.
– Jens
Commented Mar 20, 2014 at 4:23