Find maximum value of interpolation function - obviously wrong result

Yet, I observe a strange behaviour. Well, I can understand that the wrong guess of Mathematica is due to several maxima, but still I'm surprised I don't manage to get a better result.

Generation of the considered function

Here is how I get the interpolation function Theta[t]:

OmegaS = 1.32; OmegaP = 1; l0 = 1; tmax = 100;

eqn1[t_] = r''[t]/l0 - r[t]*theta'[t]^2/l0 + OmegaS^2*(r[t]/l0 - 1) - OmegaP^2*Cos[theta[t]];
eqn2[t_] = r[t]*theta''[t]/l0 + 2 r'[t]*theta'[t]/l0 + OmegaP^2*Sin[theta[t]];

sol = NDSolve[{eqn1[t] == 0, eqn2[t] == 0, r == .9, r' == 0,
theta == 0.001, theta' == 0}, {r, theta}, {t, 0, tmax}][];

Theta[t_] = theta[t] /. sol


This is what Theta looks like:

Plot[Theta[t], {t, 0, 100}, PlotRange -> All,  AxesLabel -> {"t", "Theta(t)"}] Attempts to identify the maximum

I now want to find the maximum of Theta[t] for $t\in[0,100]$. I underline that the peak could be anywhere in $[0,100]$ so I cannot indicate any relevant initial guess.

I tried all the following solutions, which obviously give an incorrect result:

MaxValue[{Theta[t], 0 <= t <= tmax}, t]
(* 0.243314 *)
NMaxValue[{Theta[t], 0 <= t <= tmax}, t]
(* 0.243314 *)
FindMaximum[{Theta[t], 0 <= t <= tmax}, t]
(* {0.000999881, {t -> 0.0146385}} *)
Maximize[{Theta[t], 0 <= t <= tmax}, t]
(* {0.243314, {t -> 92.1206}} *)
NMaximize[{Theta[t], 0 <= t <= tmax}, t]
(* {0.243314, {t -> 92.1206}} *)


The best (inelegant) way I found to approximate the maximum is brute-force:

Max[Table[Theta[i], {i, 0, tmax, .001}]]
(* 0.687071 *)


Question

Do you know:

• why my attempts fail?
• a good way to find the (global!) maximum?
• You use FindMaximum, which only looks for a local maximum. It would do what you want if you gave it a starting point near the peak you were interested in: FindMaximum[{Theta[t], 0 <= t <= tmax}, {t, 50}]. You could try NMaximize, which tries to find a global maximum, but I doubt that will work very well for arbitrary interpolating functions. – KAI Aug 26 '14 at 19:08
• @KAI Unfortunately the peak can be anywhere in $[0,100]$. I hadn't put NMaximize in the list but I guess Maximize automatically switched to NMaximize. Anyway, the give the same (wrong) answer, unfortunately. [I'll add NMaximize in my question]. – anderstood Aug 26 '14 at 19:13
• Probably the only foolproof way would be to peer into the internal structure of the InterpolatingFunction: it does a piecewise polynomial interpolation, and one might be able to reliably find the maximum of each piece. Unfortunately I'm afraid I don't know how to do that. – Rahul Aug 26 '14 at 19:19
• @Pickett: True! But Minimize[{-Theta[t], 0 <= t <= tmax}, t] does not work. – anderstood Aug 26 '14 at 19:25
• @anderstood Ah yes, I deleted my comment when I realized. Thanks for asking this question, I learned something from it :) – C. E. Aug 26 '14 at 19:53

As described e.g. in the tutorial Numerical Nonlinear Global Optimization there are different optimization methods available.
For your problem "SimulatedAnnealing" seems to work:

NMaximize[{Theta[t], 0 <= t <= tmax}, t, Method -> "SimulatedAnnealing"]

{0.687071, {t -> 48.2449}}


"DifferentialEvolution" will work, if the population is of sufficient size:

NMaximize[{Theta[t], 0 <= t <= tmax}, t, Method -> {"DifferentialEvolution", "SearchPoints" -> 30}]


"RandomSearch" will also work, if the number of points used to start local searches is big enough, e.g. bigger than 13 for your example.

NMaximize[{Theta[t], 0 <= t <= tmax}, t, Method -> {"RandomSearch", "SearchPoints" -> 14}]

• Map onto all methods and select one that works best: {#, NMaximize[{Theta[t], 0 <= t <= tmax}, t, Method -> #]} & /@ { Automatic, "DifferentialEvolution", "NelderMead", "RandomSearch", "SimulatedAnnealing"} // Grid – Bob Hanlon Aug 26 '14 at 19:40
• I unaccepted you question (before your edit) because the solution with SimulatedAnnealing is not very robust: it does not work if you change OmegaS to 1.96. I'm going to try the other suggestion right now, thank you. – anderstood Aug 26 '14 at 19:40
• @anderstood Robustness will increase, if you increase the number of initial points. It can also be increased for "SimulatedAnnealing", e.g. Method -> {"SimulatedAnnealing", "SearchPoints" -> 50} – Karsten 7. Aug 26 '14 at 19:50
• @anderstood I tried different number of points. But now I see that the speed depends on your OmegaS. For OmegaS = 1.32 "SimulatedAnnealing" is faster, for OmegaS = 1.96 "RandomSearch" is faster. So, you are right, it is a bit more complicated. – Karsten 7. Aug 26 '14 at 20:14
• @Karsten7. I think that's why we speak of numerical recipes! It's really like cooking :) – anderstood Aug 26 '14 at 20:25

Two approaches: Finding maxima (1) of an InterpolatingFunction and (2) via NDSolve.

InterpolatingFunction

To find the extrema of an InterpolatingFunction, one should start with the data contained in the function.

ifn = theta /. sol


Then ifn["Grid"] and ifn["VaiuesOnGrid"] contain the abscissas (t) and ordinates (theta[t]) computed by NDSolve. We should start our search for a maximum near the maximum value stored in ifn. The position of the maximum can be calculated with Ordering. So the best starting point is found with

{{t0}} = ifn["Grid"] ~Part~ Ordering[ifn["ValuesOnGrid"], -1]
(* {{48.2466}} *)


Then one can use FindMaximum:

FindMaximum[ifn[t], {t, t0}]

(* FindMaximum::lstol complaint... *)

(* {0.687071, {t -> 48.24493512854019}} *)


If the complaint is worrisome, then one can check the result by plotting the result, increasing working precision or trying, say, NMaximize. These return the same solution with no warnings:

FindMaximum[ifn[t], {t, t0}, WorkingPrecision -> $MachinePrecision] NMaximize[{ifn[t], 0 < t < 100}, t, Method -> {"RandomSearch", "InitialPoints" -> {{t0}}}]  NDSolve An easier approach, if the InterpolatingFunction was created with NDSolve, is to use WhenEvent. One can detect a local maximum with the event theta'[t] < 0, where the derivative becomes negative, and Sow the corresponding values of t and theta[t]. Then one selects the absolute maximum among them with Ordering. {sol, {localmax}} = Reap[NDSolveValue[{eqn1[t] == 0, eqn2[t] == 0, r == .9, r' == 0, theta == 0.001, theta' == 0, WhenEvent[theta'[t] < 0, Sow[{theta[t], {t1 -> t}}]]}, {r, theta}, {t, 0, tmax}] ]; # ~Part~ First@Ordering[#[[All, 1]], -1] &@ localmax (* {0.687071, {t1 -> 48.2449}} *)  Index for the last line of code: localmax : a list of local maxima in the form {{y.yyy, {t1 -> x.xxx}},...} #[[All, 1]] : {y.yyy,...} -- list of local maxima First@ Ordering[#[[All, 1]], -1] : position of the greatest local maximum = global maximum # ~Part~ "" : extracts the part of localmax that is the global maximum  Re global maximization -- updated In general, global optimization is difficult, especially when there are many local extrema. Finding all of them can be uncertain and time consuming. In turn one cannot be certain about the global extremum found. Sometimes personal insight has to be applied to optimizing a given objective function. In this case there seems to be an issue with precision in the InterpolatingFunction. Starting with the critical points found above, if we examine the values of ifn, we see that they constant in a neighborhood of the critical point. This probably has something to do with the optimization problem, and almost certainly with the FindMaximum::lstol message. tfindroot = {t -> 48.24493512854019}; (* FindRoot *) tndsolve = {t -> 48.24493514514225}; (* NDSolve *) Table[ifn[t (1 + 2^20$MachineEpsilon dt)] - ifn[t] /. tfindroot, {dt, -2, 2}]
Table[ifn[t (1 + 2^20 $MachineEpsilon dt)] - ifn[t] /. tndsolve, {dt, -2, 2}] (* {-3.33067*10^-16, -2.22045*10^-16, 0., 0., 0.} {-1.11022*10^-16, 0., 0., 0., -1.11022*10^-16} *) Table[ifn[t (1 + 2^21$MachineEpsilon dt)] - ifn[t] /. tfindroot, {dt, -2, 2}]
Table[ifn[t (1 + 2^21 $MachineEpsilon dt)] - ifn[t] /. tndsolve, {dt, -2, 2}] (* {-9.99201*10^-16, -3.33067*10^-16, 0., 0., -2.22045*10^-16} {-5.55112*10^-16, -1.11022*10^-16, 0., -1.11022*10^-16, -5.55112*10^-16} *)  Interestingly, using arbitrary precision helps NMaximize on the posted problem, but it does not work for all values of the parameter OmegaS. Sample code: sol = First@NDSolve[ Rationalize@ {eqn1[t] == 0, eqn2[t] == 0, r == 0.9, r' == 0, theta == 0.001, theta' == 0}, {r, theta}, {t, 0, tmax}, WorkingPrecision ->$MachinePrecision];
ifn = theta /. sol;
NMaximize[{ifn[t], 0 <= t <= tmax}, t, WorkingPrecision -> $MachinePrecision] (* {0.6870714027852798, {t -> 48.24493489787723}} *)  The NDSolve method is probably the best way to proceed. Theoretically, if there are local maxima that are close in value to the global maximum, then one of them might be returned instead of the global one. But they would have to be very close, within the error of the NDSolve computation. If one wanted further evidence, one could use the "RandomSearch" method of NMaximize with the local maxima (localmax from the NDSolve approach above) as initial points: NMaximize[ifn[t], t, Method -> {"RandomSearch", "InitialPoints" -> ({t1} /. localmax[[All, -1]])}, WorkingPrecision ->$MachinePrecision]
(* {0.6870714027850805, {t1 -> 48.24493514549923}} *)

• Thank you for this very instructive answer, which I have not tried yet. I guess you meant ifn["..."] instead of if["..."] in the first paragraph but I prefer to let you change it, if need be. – anderstood Aug 27 '14 at 1:32
• @anderstood You're welome, and thanks for pointing out if. I changed from if to ifn because if looked too much like an English word. Guess I missed some. – Michael E2 Aug 27 '14 at 1:35
• +1 for the WhenEvent[] idea. Never thought about that use before. Great. – Dr. belisarius Aug 27 '14 at 5:49
• Nowadays, PeakDetect[]` is a usable alternative for the first approach in this answer. – J. M. is away Dec 9 '16 at 20:02