Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

enter image description here

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: enter image description here

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.

share|improve this question
    
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). –  celtschk May 15 '12 at 10:50
add comment

4 Answers 4

up vote 40 down vote accepted

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

enter image description here

share|improve this answer
    
Unbeatable simplicity (+1). –  Jens May 15 '12 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) –  Margus May 15 '12 at 19:48
add comment

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

enter image description here

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

share|improve this answer
    
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 ? –  Artes May 15 '12 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. –  Mark McClure May 15 '12 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. –  Artes May 15 '12 at 16:37
add comment

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

enter image description here

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]

enter image description here

share|improve this answer
    
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. –  Margus May 16 '12 at 12:49
add comment

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

styled minima and maxima on line

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.