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

Question

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.

Answer 1

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

Answer 2

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["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Needs["Ersek`RootSearch`"]
{{a, b}} = InterpolatingFunctionDomain[Head[f[x]]];
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]
}]]

Some comments

• 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

Answer 3

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["DifferentialEquations`InterpolatingFunctionAnatomy`"];
grid = First@InterpolatingFunctionCoordinates[f];

Show[{Plot[f''[x], {x, -9, 30}, PlotRange -> All],ListPlot[Thread[{grid,  
f''[grid]}], PlotStyle -> Red], 
Graphics[Map[{Dashed, Line[#]} &,Partition[Thread
[{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]

Answer 4

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

Answer 5

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]}];
f = Head@First@D[getRLine3[{data}, 3][x], x];

(* 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.