7
$\begingroup$

Probably a fairly easy one; but I'm drawing a blank I would like to find the maximum of an interpolated data set. I've tried using

maxCurrent = 
FindMaximum[
Interpolation[dataNi59s, InterpolationOrder -> 2], {y, 
dataNi59s[[1]]}]

Just for clarification - I would like to return the maximum value of the interpolation of the data, rather than just the largest value in the data set.


Edit 1: A lot of the answers have (I presume correctly) outputted the answer as '{28.9179, {x -> 177.957}}', however the intended result is to output ONLY the largest value of y that lies on the blue interpolated line - by eye this looks to be about 23.9uA. (The x-co-ordinate is at about 400, and if anyone knows how to return that separately that would also be very helpful)

I've included this image to help clarify - I would like to return the y co-ordinate of the point labelled with a red arrow

graph - showing global maxima


Edit 2:

For the record - following Aardvark's advice that finding the peak of the ListLinePlot interpolation line was physically not meaningful, i decided to just output the highest current value of one of the data points. Aardvark's answer does show how to do this, but technically it outputs all points that share that maximum and the x and y co-ordinate, whereas i just wanted the y-co-ordinate once. So for posterity, here is the code I used to do that, based on Aardvark's answer

maxCurrent = MaximalBy[dataNi59s, Last][[1, 2]]

This returns (for my data) just: 23.4, which is the new intended result.

For any future reader, if you'd like to know how to find the maximum of the ListLinePlot as I clarified in Edit 1 with the red arrow on the graph, then it can be found in update 2 of Aardvark's answer.


Data:

dataNi59s = {{26, 8.07}, {28, 9.12}, {36, 11.9}, {40, 13.4}, {66, 
17.5}, {69, 17.9}, {163, 18.9}, {165, 21.4}, {191, 21.3}, {200, 
21.6}, {214, 21.5}, {227, 21.2}, {260, 20.9}, {342, 20.8}, {362, 
23.4}, {497, 23.4}, {510, 23.3}, {535, 22.9}, {566, 22.4}, {577, 
22.2}, {596, 21.9}, {607, 21.7}, {671, 21.0}, {682, 20.7}, {689, 
20.5}, {700, 20.9}, {704, 20.8}, {799, 19.6}, {809, 19.5}, {823, 
19.6}, {831, 19.5}, {841, 19.3}, {852, 19.3}, {875, 19.1}, {881, 
19.1}, {901, 19.0}, {911, 19.0}, {928, 19.0}, {955, 18.8}, {973, 
18.7}, {978, 18.6}, {1048, 18.5}, {1067, 18.4}, {1090, 
18.2}, {1111, 18.1}, {1126, 18.0}, {1134, 17.9}, {1196, 
17.7}, {1200, 17.6}, {1215, 17.6}, {1227, 17.5}, {1270, 
17.4}, {1367, 16.6}, {1377, 16.6}, {1381, 16.5}, {1494, 
16.3}, {1534, 16.3}, {1584, 16.0}, {1686, 15.4}, {1700, 
15.3}, {1707, 15.3}, {1805, 14.5}, {1869, 21.6}, {1878, 
21.3}, {1905, 21.1}, {1968, 21.0}, {1980, 21.0}, {2010, 
20.9}, {2025, 20.8}, {2063, 20.5}, {2070, 20.5} , {2109, 
20.1}, {2146, 19.8}, {2232, 19.3}, {2249, 19.1}, {2280, 
18.8}, {2347, 18.1}, {2370, 17.9}, {2394, 17.8}, {2400, 
17.6}, {2416, 17.5}, {2430, 17.4}, {2452, 17.2}, {2462, 
17.1}, {2488, 17.0}, {2508, 16.8}, {2528, 16.5}, {2599, 
15.7}, {2604, 15.2}, {2612, 15.0}, {2637, 14.7}, {2651, 
14.6}, {2660, 14.5}, {2760, 13.6}, {2769, 13.6}, {2805, 
13.2}, {2844, 12.9}, {2880, 12.6}, {2904, 12.4}, {3040, 
12.4}, {3081, 12.1}, {3189, 11.6}, {3221, 11.5}, {3314, 
10.7}, {3602, 12.5}, {3610, 12.1}, {3624, 12.0}, {3635, 
12.0}, {3640, 11.7}, {3735, 11.1}, {3765, 11.1}, {3797, 
11.0}, {4066, 7.71}, {4120, 5.91}, {4123, 5.78}, {4129, 
5.61}, {4169, 4.77}, {4174, 4.69}, {4183, 4.56}, {4188, 
4.45}, {4194, 4.40}, {4200, 4.36}, {4221, 4.13}, {4345, 
4.47}, {4360, 4.39}, {4492, 3.83}, {4550, 3.99}, {4567, 
3.94}, {4574, 3.93}, {4580, 3.92}, {4600, 3.95}, {4652, 
3.79}, {4679, 3.56}, {4794, 3.49}, {4842, 3.59}, {5049, 
4.37}, {5094, 4.47}, {5124, 3.94}, {5167, 3.11}, {5172, 
3.07}, {5210, 2.92}, {5278, 2.78}, {5314, 2.76}, {5327, 
2.75}, {5442, 2.81}, {5479, 2.79}, {5565, 2.72}, {5580, 
2.71}, {5768, 2.26}, {5929, 1.03}, {5941, 1.03}, {5948, 
1.03}, {6000, 1.01}, {6058, 0.992}, {6059, 1.14}, {6210, 
1.04}, {6328, 0.977}, {6330, 0.998}, {6333, 1.01}, {6351, 
1.01}, {6420, 1.00}, {6434, 0.964}, {6441, 1.02}, {6611, 
1.01}, {6661, 1.00}, {6696, 1.03}, {6908, 0.968}};

Graph Code:

{Ni59s[1], Ni59s[2]} = SplitBy[dataNi59s, #[[1]] < 1850 &];
avgCurrent = 
Integrate[
Interpolation[dataNi59s, InterpolationOrder -> 2][x], {x, 
 dataNi59s[[1, 1]], dataNi59s[[-1, 1]]}]/(Subtract @@ 
 dataNi59s[[{-1, 1}, 1]]);
maxCurrent = 
FindMaximum[
Interpolation[dataNi59s, InterpolationOrder -> 2], {y, 
dataNi59s[[1]]}]
noPoints = Length[dataNi59s];
lifeTime = dataNi59s[[-1, 1]];
ListLinePlot[{Ni59s[1], Ni59s[2]}, PlotRange -> {{0, 7000}, {0, 25}}, 
PlotStyle -> {Directive[RGBColor[0., 0.75, 0.85], 
AbsoluteThickness[2.5]]}, PlotTheme -> "Detailed", 
Epilog -> {{Dashed, Line[{{1805, 14.5`}, {1869, 21.6`}}]}, {Red, 
AbsolutePointSize[4], Point[dataNi59s]}, {GrayLevel[.5], 
InfiniteLine[{{0, 1}, {1, 1}}]}, 
Inset[Framed[
 Grid[{{"Average Current: ", avgCurrent, 
    " \[Mu]A"}, {"Cathode Lifetime: ", lifeTime, 
    "minutes"}, {"No. of Data Points: ", noPoints}}], 
 Background -> White, RoundingRadius -> 5], Scaled[{.85, .9}]]}, 
PlotLabel -> 
"Comparison of Current Yields from Nickel Cathode 59 (Interpolated) \
(Anomolous Points Removed)", 
PlotLegends -> 
Placed[LineLegend[{"Cathode 59 - Nickel"}, 
LegendFunction -> "Panel"], {0.85, 0.5}], 
LabelStyle -> {Black, Bold}, Frame -> {{True, False}, {True, False}},
FrameLabel -> {"Cathode Runtime (minutes)", "Current (\[Mu]A)"}, 
ImageSize -> {850, 550}, InterpolationOrder -> 2]

Graph Code (following this question's answers):

{Ni59s[1], Ni59s[2]} = SplitBy[dataNi59s, #[[1]] < 1850 &];
avgCurrent = 
Integrate[
Interpolation[dataNi59s, InterpolationOrder -> 2][x], {x, 
 dataNi59s[[1, 1]], dataNi59s[[-1, 1]]}]/(Subtract @@ 
 dataNi59s[[{-1, 1}, 1]]);
maxCurrent = MaximalBy[dataNi59s, Last][[1, 2]];
noPoints = Length[dataNi59s];
lifeTime = dataNi59s[[-1, 1]];
ListLinePlot[{Ni59s[1], Ni59s[2]}, PlotRange -> {{0, 7000}, {0, 25}}, 
PlotStyle -> {Directive[RGBColor[0., 0.75, 0.85], 
AbsoluteThickness[2.5]]}, PlotTheme -> "Detailed", 
Epilog -> {{Dashed, Line[{{1805, 14.5`}, {1869, 21.6`}}]}, {Red, 
AbsolutePointSize[4], Point[dataNi59s]}, {GrayLevel[.5], 
InfiniteLine[{{0, 1}, {1, 1}}]}, 
Inset[Framed[
 Grid[{{"Average Current: ", avgCurrent, 
    " \[Mu]A"}, {"Maximum Current: ", maxCurrent, 
    " \[Mu]A"}, {"Cathode Lifetime: ", lifeTime, 
    "minutes"}, {"No. of Data Points: ", noPoints}}], 
 Background -> White, RoundingRadius -> 5], Scaled[{.85, .9}]]}, 
 PlotLabel -> 
 "Comparison of Current Yields from Nickel Cathode 59 (Interpolated) \
 (Anomolous Points Removed)", 
 PlotLegends -> 
 Placed[LineLegend[{"Cathode 59 - Nickel"}, 
 LegendFunction -> "Panel"], {0.85, 0.5}], 
 LabelStyle -> {Black, Bold}, Frame -> {{True, False}, {True, False}},
 FrameLabel -> {"Cathode Runtime (minutes)", "Current (\[Mu]A)"}, 
 ImageSize -> {850, 550}, InterpolationOrder -> 2]
$\endgroup$
  • $\begingroup$ Question editted to provide clarification on intended result. Read: Edit. $\endgroup$ – Epideme Dec 14 '17 at 11:09
13
$\begingroup$

NMaximize is good for finding global maxima:

NMaximize[{f[x], 26 < x < 6908}, x]

(* {28.9179, {x -> 177.957}} *)

For FindMaximum (and FindMaxValue, and other numerical searches for that matter) it's often a good idea to help Mathematica along with a good starting value, to avoid getting trapped in local maxima. In this case, plotting your function first can help with this:

f = Interpolation[dataNi59s, InterpolationOrder -> 2]
Plot[f[x], {x, 26, 6910}]

enter image description here

Clearly, f has plenty of local maxima that you don't want to get stuck in, and the global maximum occurs for somewhere for x < 1000. Manipulate can be helpful here:

Manipulate[
 Plot[f[x], {x, xmin, xmax}],
 {{xmin, 26}, 26, 999}, {{xmax, 1000}, xmin + 1, 1000}]

enter image description here

So we're looking for the global maximum somewhere around x == 180 (of course, we already know this from NMaximize, but just play along).

With InterpolatingFunctions, you also need to be careful that the search doesn't go out of range. So we've got to put some constraints on x, as well as specifying a starting value:

FindMaximum[{f[x], 26 < x < 6908}, {x, 180}]

(* {28.9179, {x -> 177.957}} *)

UPDATE1: The Truth is Subject to Interpolation

The OP has expressed concern that the maximum {28.9179, {x -> 177.957}} is not what it should be. This is a valid point, and the discrepancy between the given and expected answers comes down to two things:

  1. Interpolation[data, InterpolationOrder -> 2] and ListLinePlot[data, InterpolationOrder -> 2] do not give the same interpolated function.

  2. InterpolationOrder -> 2 can give weird artifacts in either context.

Consider the following two plots:

Show[Plot[
  Interpolation[dataNi59s, InterpolationOrder -> 2][x], {x, 130, 240}], 
  ListPlot[dataNi59s, PlotStyle -> Directive[PointSize[Medium], Red]]
]
Show[ListLinePlot[
  dataNi59s, InterpolationOrder -> 2, PlotRange -> {{130, 240}, All}], 
  ListPlot[dataNi59s, PlotStyle -> Directive[PointSize[Medium], Red]]
]

enter image description here

enter image description here

Clearly, there are problems with both these plots. In the first, the maximum (which turns out to be the global maximum) is an artefact of InterpolationOrder -> 2. The second isn't even a function, so I wouldn't put much stock in what that says.

The underlying questions, then, are:

What precisely do you want to find the maximum of?

Why? What will that specific thing (e.g. that specific interpolation of your data points) tell you about your Nickel cathode that the data by itself doesn't.

  1. If, ultimately, you only want the maximum of the data, then, of course, that's easily done: MaximalBy[dataNi59s, Last] returns {{362, 23.4}, {497, 23.4}}.
  2. If you actually want the maximum of an interpolated function with second order interpolation, then you want {28.9179, {x -> 177.957}}. But you must be aware that artefacts may creep in. In this case, the artefact creates a global maximum that's nowhere near the maximum of the data.
  3. If it's the ListLinePlot interpolation, then you need to consider what that would mean in the context of your system, because it's not what I would think of as a meaningful interpolation of your data. The interpolation in ListLinePlot isn't intended for this.

My intuition would be to use InterpolationOrder -> 1 unless you have very solid reasons for doing otherwise. That would mean that the maximum of the interpolated function would just be the maximum of the data. If you want more precision than that, you will have to think a bit more carefully about how to do the interpolation. Or, ideally, get more data points in the range x ∈ [362, 497].

For completeness, here are plots of the range between those two points using the two interpolations.

Show[Plot[
  Interpolation[dataNi59s, InterpolationOrder -> 2][x], {x, 300, 550},
   PlotRange -> {{300, 550}, {19, 25}}], 
 ListPlot[dataNi59s, PlotStyle -> Directive[PointSize[Medium], Red]]
]
Show[ListLinePlot[
  dataNi59s, InterpolationOrder -> 2, PlotRange -> {{300, 550}, {19, 25}}], 
 ListPlot[dataNi59s, PlotStyle -> Directive[PointSize[Medium], Red]]
]

enter image description here

enter image description here

Ask yourself if either of them are more trustworthy than the data points themselves.

UPDATE 2: ListLinePlot Interpolation Maximum

Just for completeness, if you really want to find the maximum of the ListLinePlot interpolation, you can extract the data from the plot using @Mr.Wizard's answer here.

plot = ListLinePlot[dataNi59s, InterpolationOrder -> 2, PlotRange -> {All, {0, 25}}];
data = Cases[plot, Line[data_] :> data, -4, 1][[1]];
maxpoint = MaximalBy[data, Last]
Show[plot, Graphics[{PointSize[Large], Red, Point[maxpoint]}]]

(* {{404.7, 23.6727}} *)

enter image description here

Which is probably closer to what you're looking for, but I hope you will question its validity.

$\endgroup$
  • $\begingroup$ I managed to get a result to output using maxCurrent = NMaximize[{Interpolation[dataNi59s, InterpolationOrder -> 2][x], Min[dataNi59s[[All, 1]]] <= x <= Max[dataNi59s[[All, 1]]]}, {x, 0, 7000}], however this is the incorrect result. It outputs {28.9179, {x -> 177.957}}, as you suggested it should. However the intended result it just to output the largest yalue of y that lies on the interpolated line (the blue one) - so the desired output would be about 23.9uA. What am i missing? $\endgroup$ – Epideme Dec 14 '17 at 11:00
  • $\begingroup$ @Epideme See my update. Probably not quite what you were looking for, but I hope it's illuminating. $\endgroup$ – aardvark2012 Dec 15 '17 at 1:48
  • $\begingroup$ Thank you for the thorough explanation Aardvark - both scientifically and on the inner workings of Mathematica's interpolation function. A well earned best answer from me. In the end, I heeded your advice that the peak of the interpolated line was ultimately a fairly arbitrary artifact of the algorithm used, and used 23.4uA as the maximum.` InterpolationOrder->2` comes from an old model of how this sputtering process works physically, which is related to surface area. And thank you for adding in the way to find the listlineplot maximum for future reference. $\endgroup$ – Epideme Dec 19 '17 at 15:49
8
$\begingroup$

You want to restrict FindMaximum to the region

region = Min[First /@ dataNi59s] < time < Max[First /@ dataNi59s] && Min[Last /@ dataNi59s] < current < Max[Last /@ dataNi59s]
(* 26 < time < 6908 && 0.964 < current < 23.4 *)

With this,

maxCurrent = FindMaximum[{Interpolation[dataNi59s, InterpolationOrder -> 2][{time, current}][[2]], region}, {{time, dataNi59s[[1, 1]]}, {current, dataNi59s[[1, 2]]}}] //Quiet
(* {6.49271, {time -> 820.675, current -> 23.4}} *)

Visually speaking, this is the expected result:

enter image description here

$\endgroup$
4
$\begingroup$

Use NMaximize:

NMaximize[{Interpolation[dataNi59s, InterpolationOrder -> 2][x],Min[dataNi59s[[All, 1]]] <= x <= Max[dataNi59s[[All, 1]]]}, {x,3000, 4000}] 

In the list {x,3000, 4000} I gave a search range for NMaximize to get the right maximum. In your problem you can get the 'right' interval in advance by analysing your data

$\endgroup$
  • $\begingroup$ So, I tried this line... maxCurrent = NMaximize[{Interpolation[dataNi59s, InterpolationOrder -> 2][x], Min[dataNi59s[[All, 1]]] <= x <= Max[dataNi59s[[All, 1]]]}, {x, 0, 7000}] It returns: {28.9179, {x -> 177.957}} I used 0-7000 as the interval (essentially the x range of the whole function, as i'd like this code to work for several data sets). The intended result is just the highest y co-ordinate on the blue interpolated line - which by eye, looks about 23.8uA. How would I go about correcting this so mathematica returns the intended answer? $\endgroup$ – Epideme Dec 12 '17 at 16:40
  • $\begingroup$ If you analyse your data you'll find two maxima {{362, 23.4}, {497, 23.4}}, which you could use as startvalues for NMaximize... $\endgroup$ – Ulrich Neumann Dec 14 '17 at 14:08
  • $\begingroup$ I see your point now - aardvark's explanation helped with this particular issue, Thank you Ulrich. $\endgroup$ – Epideme Dec 19 '17 at 15:50
3
$\begingroup$

First, the syntax of FindMaximum is (only list relevant parts)

FindMaximum[f, x]
searches for a local maximum in f, starting from an automatically selected point.
FindMaximum[f, {x, x0}]
searches for a local maximum in f, starting from the point x = x0. 

The f is supposed to be an expression with variable x, so in your case you need to write

FindMaximum[
 Interpolation[dataNi59s, InterpolationOrder -> 2][y], {y, 
  dataNi59s[[1]]}]

because Interpolation[...] returns a InterpolatingFunction object.

Moreover, the second argument should be either y or {y, y0}, and dataNi59s[[1]] in your case is {26, 8.07}, which is not what you want. Just let Mathematica automatically choose the point to find would be good enough.

So you need

FindMaximum[Interpolation[dataNi59s, InterpolationOrder -> 2][y], y]

When the second argument (x0) is a list, Mathematica will understand that the variable is a list of that format. For example:

FindMaximum[-(x[[1]] - 2)^2 - (x[[2]] + 3)^2, {x, {1, 1}}]

Mathematica will understand that x is a list of 2 numbers, and starts from point {1, 1}, and gives the result {0., {x -> {2., -3.}}} (this is not a perfect example because Mathematica throws a lot of warnings). That's similar to:

FindMaximum[-(a - 2)^2 - (b + 3)^2, a, b]

which outputs

{0., {a -> 2., b -> -3.}}
$\endgroup$
  • $\begingroup$ Hi, I tried this line as you suggested... maxCurrent = FindMaximum[Interpolation[dataNi59s, InterpolationOrder -> 2][y], y]; It returns the following: InterpolatingFunction: Input value {1.} lies outside the range of data in the interpolating function. Extrapolation will be used. General: Further output of InterpolatingFunction will be suppressed during this calculation. $\endgroup$ – Epideme Dec 12 '17 at 16:28
  • $\begingroup$ Also: FindMaximum: The line search decreased the step size to within the tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient increase in the function. You may need more than MachinePrecision digits of working precision to meet these tolerances. $\endgroup$ – Epideme Dec 12 '17 at 16:33
  • $\begingroup$ This solution as you explained it to me, definitely makes sense to me, but i don't see why Mathematica doesn't much like it. $\endgroup$ – Epideme Dec 12 '17 at 16:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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