# How to find all the local minima/maxima in a range

I want to find :

• all local maxima in range
• all local minima in range

From those points I can interpolate and combine functions upper and lower boundary. What I am really interested in, is the mean function of those boundaries.

Data model for this plot:

GetRLine3[MMStdata_, IO_: 1][x_: x] := ListInterpolation[#, InterpolationOrder -> IO,
Method -> "Spline"][x] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@ MMStdata);
data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4],
Tanh[#] + (Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@ Range[-4, 4, 0.08]}];
xLimits = {Min@#1, Max@#1} & @@ Transpose[data];
f = D[GetRLine3[{data}, 3][x], x];


Edit: As my effort:

minimums = DeleteDuplicates[Round[x /. Last[FindMinimum[f, {x, #}]] & /@ Transpose[data][[1]], 0.0001]]
minimumvalues = (f /. x -> #)[[1]] & /@ minimums;
minimumData := Transpose[{minimums, minimumvalues}];
maximums = DeleteDuplicates[Round[x /. Last[FindMaximum[f, {x, #}]] & /@ Transpose[data][[1]], 0.0001]];
maximumsvalues = (f /. x -> #)[[1]] & /@ maximums;
maximumsData := Transpose[{maximums, maximumsvalues}];

maxf = Max[{GetRLine3[{maximumsData}, 3][x], f}]
minf = Min[{GetRLine3[{minimumData}, 3][x], f}]
mf = Mean[{maxf, minf}]


This was what I was trying to make:

I still get quite few warnings and I'm sure it's not the best solution. I don't like the DeleteDuplicates@Round@ part, but it was necessarily to get the interpolation function working.

• There's always a maximum in between two minima. So if you have found all minima (I don't know how to be sure of that, though), you can just look for a maximum in each of the intervals between subsequent minima (and for possibly two more on the left/right of the first/last minimum). Commented May 15, 2012 at 10:50

This can be done using event location within NDSolve. I start off as below (note f is slightly modified from what you have, mostly to rescale it).

GetRLine3[MMStdata_, IO_: 1][x_: x] :=
ListInterpolation[#, InterpolationOrder -> IO, Method -> "Spline"][
x] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@ MMStdata);
data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4],
Tanh[#] + (Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@
Range[-4, 4, 0.08]}];

xLimits = {Min@#1, Max@#1} & @@ Transpose[data];
f = First[100*D[GetRLine3[{data}, 3][x], x]];


We'll recapture f using NDSolve, and locate the points where the derivative vanishes in the process.

vals = Reap[
soln = y[x] /.
First[NDSolve[{y'[x] == Evaluate[D[f, x]],
y[-9.9] == (f /. x -> -9.9)}, y[x], {x, -9.9, 30},
Method -> {"EventLocator", "Event" -> y'[x],
"EventAction" :> Sow[{x, y[x]}]}]]][[2, 1]];


Visual check:

Plot[f, {x, -9.9, 30},
Epilog -> {PointSize[Medium], Red, Point[vals]}]


• Unbeatable simplicity (+1).
– Jens
Commented May 15, 2012 at 16:21
• I assume you forgot {Take[val, {2, -1, 2}],Drop[val, {2, -1, 2}]} and to differentiate maximum list from minimum list, we just need to compare first element. Not bad at all actually +1. (Not sure however, if this would work, if function has a segment with a constant growth of 0) Commented May 15, 2012 at 19:48
• Very very cool! Commented Mar 19, 2015 at 13:37

Sounds like a job for Ted's RootSearch package.

Clear[f];
SeedRandom[1];
GetRLine3[MMStdata_, IO_: 1][x_: x] := ListInterpolation[#, InterpolationOrder -> IO,
Method -> "Spline"][x] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@ MMStdata);
data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4],
Tanh[#] + (Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@ Range[-4, 4, 0.08]}];
xLimits = {Min@#1, Max@#1} & @@ Transpose[data];
f[x_] = First[D[GetRLine3[{data}, 3][x], x]];


Note that I changed the way yo defined f slightly.

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];
Needs["ErsekRootSearch"]
criticalPoints = x /. RootSearch[f'[x] == 0, {x, a, b}];
mins = Select[criticalPoints, f''[#] > 0 &];
maxs = Select[criticalPoints, f''[#] < 0 &];


Well, I get an error so I'm not certain how well it worked, but it does look good.

Plot[f[x], {x, -9, 30}, Epilog ->{PointSize[Medium],
Blue, Point[{#, f[#]} & /@ mins],
Darker[Red], Point[{#, f[#]} & /@ maxs]
}]]


• I didn't actually notice the importance of separating the maxima from the minima right away and have edited accordingly
• Of course, it's feasible that $f'(x)=0$ without a max or min occurring at $x$, although I think that the random nature of your function makes the probability of this occurring zero.
• Not sure we can really prove we've got all the extremes.

Edit in Response to Artes

If you use FindMaximum from two different starting values you will likely get two slightly different approximations to the same root. Here's an example with exactly one max.

g[x_] = Exp[-(Cos[x] - x)^2];
x1 = x /. Last[FindMaximum[g[x], {x, 0.8}]];
x2 = x /. Last[FindMaximum[g[x], {x, 0.7}]];
x1 - x2


8.35447*10^-9

• Evaluating e.g. FindMaximum[f, {x, #}] & /@ {-5, -4.5} we got 2 maxima, where the difference between them : Subtract @@ (#2[[1, 2]] & @@@ %)is about 1.33227*10^-14. Does RootSearch treat them as two different ones or as the only one local maximum ? Commented May 15, 2012 at 12:08
• @Artes Due to the random nature of the function, I can't reproduce this behavior. I think it is more likely, though that your FindMaximum code is finding two slightly different approximations to the same max. See my edit for an example. Commented May 15, 2012 at 12:15
• Indeed, it seems to be the same maximum. In your example we can set appropriately WorkingPrecision etc. to ensure the issue, however for some experimental or interpolated data it could be difficult do decide. Commented May 15, 2012 at 16:37

Since you're fitting your data with a cubic spline, and defining $f(x)$ as the derivative of the interpolation function, then if we solve for $f'(x)=0$ for the critical points, the resulting function is a piecewise linear function. We can pull the interpolation x-grid data and use it to construct linear functions which are easy to solve for the roots without resorting to the FindMinimum, etc...

Clear[f];
SeedRandom[1];
GetRLine3[MMStdata_, IO_: 1][x_: x] :=
ListInterpolation[#, InterpolationOrder -> IO,Method -> "Spline"] & /@
(({{#[[1]]}, #[[2]]}) & /@ # & /@MMStdata);
data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4], Tanh[#] +
(Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@Range[-4, 4, 0.08]}];
xLimits = {Min@#1, Max@#1} & @@ Transpose[data];
f = First[GetRLine3[{data}, 3][x]]

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];
grid = First@InterpolatingFunctionCoordinates[f];

Show[{Plot[f''[x], {x, -9, 30}, PlotRange -> All],ListPlot[Thread[{grid,
f''[grid]}], PlotStyle -> Red],
[{grid, f''[grid]}], 2, 1]]]}]


The dashed curve is the reconstruction from the interpolation grid points shown in red. Next, $\it{eqn}$ defines the linear parameterizations where we have an intersection with the $x-$axis if $0<t<1$ when $y=0$. Since the interpolation function is a polynomial, the minima and maxima alternate and can be pulled from the zeros list in order.

NB: The dashed curve is not exactly $f''(x)$ here, as the red interpolation sample points don't lie exactly on the critical points in the first figure. According to the documentation, taking the derivative of an InterpolatingFunction returns a new InterpolatingFunction, however, it seems the grid values are the same for both. I guess you can always find $f'''(x)=0$ to get those.

pts = Thread[{grid, f''[grid]}];
eqn = Map[#[[1]] (1 - t) + #[[2]] (t) &, Partition[pts, 2, 1]];
zeros = Select[
Flatten[Table[{x, t} /. Solve[eqn[[j]] == {x, 0}, {x, t}], {j,
Length[eqn]}], 1], 0 <= #[[2]] <= 1 &][[All, 1]];
zerosM = zeros[[1 ;; -1 ;; 2]];
zerosN = zeros[[2 ;; -1 ;; 2]];

Plot[f'[x], {x, -9, 30},
Epilog -> {PointSize[Medium], Red, Point[{#, f'[#]} & /@ zerosM],
Blue, Point[{#, f'[#]} & /@ zerosN]}, PlotRange -> All]


• What a great answer. At first, it did not work with my real data and I assumed your code was improperly written - it was not the case. The bug was that I attached [x] at the end of a interpolation function and it messed up further computations. I can only imagine how much time you saved me, it's pity that I can +1 only once. Commented May 16, 2012 at 12:49

I just come to sell my new MeshFunctions powered "golden hammer" :)

GetRLine3[MMStdata_, IO_: 1][x_: x] :=
ListInterpolation[#, InterpolationOrder -> IO, Method -> "Spline"][
x] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@ MMStdata);
data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4],
Tanh[#] + (Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@
Range[-4, 4, 0.08]}];
xLimits = {Min@#1, Max@#1} & @@ Transpose[data];

f = Function[x, Evaluate[D[GetRLine3[{data}, 3][x], x][[1]]]];

funcline = Plot[f[x], {x, ##}, PlotStyle -> GrayLevel[.8]] & @@ xLimits;

maximapts = Plot[f[x], {x, ##},
PlotStyle -> None,
PlotRange -> All,
MeshFunctions -> Function[{x, y}, f'[x]],
Mesh -> {{0}},
MeshStyle -> Directive[AbsolutePointSize[4], Red],
RegionFunction -> Function[{x, y}, f''[x] < 0]
] & @@ xLimits;

minimapts = Plot[f[x], {x, ##},
PlotStyle -> None,
PlotRange -> All,
MeshFunctions -> Function[{x, y}, f'[x]],
Mesh -> {{0}},
MeshStyle -> Directive[AbsolutePointSize[4], Blue],
RegionFunction -> Function[{x, y}, f''[x] > 0]
] & @@ xLimits;

Show[{funcline, maximapts, minimapts}]


# Update

To extract the found extreme points, we can do it like this:

Cases[maximapts,
GraphicsComplex[ pts_, body__ ]:>(
pts[[#]]& /@
Cases[{body}, Point[idxes_]:>idxes, ∞]
),
∞]

{{{{1.43947,0.041135},{4.75921,0.130593},<<30>>,{21.4528,-0.000886152}}}}

• your method works very well for me. Can you please tell me how to obtain a list of the coordinates of the maxima and minima? Commented Nov 27, 2014 at 17:56
• @Pincopallino Thanks! To extract the coordinates, you might want to check the InputForm of the graphics for clues. (e.g. Cases[maximapts,GraphicsComplex[pts_,body__]:>Column[Short/@{pts,Cases[{body},Point[idxes_]:>idxes,∞]},Frame->All],∞][[1]] -- you might want to re-enter the Point because of SE's formatting on linebreaks in comment. ) Commented Nov 27, 2014 at 18:14
• @Pincopallino Please see my update in the answer. Commented Nov 27, 2014 at 18:32
• Thank you very much for your prompt response! Commented Nov 29, 2014 at 17:23
• @Pincopallino No problem! But I should add a warning that this method will work only where the points are differentiable. Commented Nov 29, 2014 at 18:14

Just for the record, if you are willing to (re)sample your function to have {x, y} pairs instead of a regular function or InterpolatingFunction, there are multiple methods available (see for example here and here):

Original function:

getRLine3[MMStdata_, IO_: 1][x_: x] :=
ListInterpolation[#, InterpolationOrder -> IO, Method -> "Spline"][
x] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@ MMStdata);
data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4],
Tanh[#] + (Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@
Range[-4, 4, 0.08]}];

(* sample function at regular times *)
times = Range[0, 30, .01];
vals = f[times];


Using FindPeaks:

{maxP, max} = Transpose@FindPeaks[vals];
{minP, min} = Transpose@FindPeaks[-vals];

Plot[f[x], {x, 0, 30}, Epilog -> {AbsolutePointSize[5],
Green, Point@Transpose@{times[[maxP]], max},
Red, Point@Transpose@{times[[minP]], -min}}]


Using MaxDetect and MinDetect:

{minD, minD} = {MinDetect@vals, MaxDetect@vals};

Plot[f[x], {x, 0, 30}, Epilog -> {AbsolutePointSize[5],
Green, Point@Pick[Transpose@{times, vals}, maxD, 1],
Red, Point@Pick[Transpose@{times, vals}, minD, 1]
}]


MaxDetect+MinDetect is at least an order of magnitude faster than FindPeaks but of course the latter offers smoothing capabilities and further methods to finetune detection.