20
$\begingroup$

Minimizing and maximizing interpolation function has already been asked and answered, see here for example.

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[0] == .9, r'[0] == 0, 
theta[0] == 0.001, theta'[0] == 0}, {r, theta}, {t, 0, tmax}][[1]];

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

This is what Theta looks like:

Plot[Theta[t], {t, 0, 100}, PlotRange -> All,  AxesLabel -> {"t", "Theta(t)"}]

plot Theta

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?
$\endgroup$
5
  • 1
    $\begingroup$ 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. $\endgroup$
    – KAI
    Aug 26, 2014 at 19:08
  • 1
    $\begingroup$ @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]. $\endgroup$
    – anderstood
    Aug 26, 2014 at 19:13
  • $\begingroup$ 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. $\endgroup$
    – user484
    Aug 26, 2014 at 19:19
  • $\begingroup$ @Pickett: True! But Minimize[{-Theta[t], 0 <= t <= tmax}, t] does not work. $\endgroup$
    – anderstood
    Aug 26, 2014 at 19:25
  • 1
    $\begingroup$ @anderstood Ah yes, I deleted my comment when I realized. Thanks for asking this question, I learned something from it :) $\endgroup$
    – C. E.
    Aug 26, 2014 at 19:53

2 Answers 2

18
$\begingroup$

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}]
$\endgroup$
9
  • 5
    $\begingroup$ Map onto all methods and select one that works best: {#, NMaximize[{Theta[t], 0 <= t <= tmax}, t, Method -> #]} & /@ { Automatic, "DifferentialEvolution", "NelderMead", "RandomSearch", "SimulatedAnnealing"} // Grid $\endgroup$
    – Bob Hanlon
    Aug 26, 2014 at 19:40
  • 1
    $\begingroup$ 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. $\endgroup$
    – anderstood
    Aug 26, 2014 at 19:40
  • 2
    $\begingroup$ @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} $\endgroup$
    – Karsten7
    Aug 26, 2014 at 19:50
  • 1
    $\begingroup$ @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. $\endgroup$
    – Karsten7
    Aug 26, 2014 at 20:14
  • 1
    $\begingroup$ @Karsten7. I think that's why we speak of numerical recipes! It's really like cooking :) $\endgroup$
    – anderstood
    Aug 26, 2014 at 20:25
15
$\begingroup$

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.

Let's start with the function itself;

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[0] == .9, 
     r'[0] == 0, theta[0] == 0.001, theta'[0] == 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] == 0.9, 
     r'[0] == 0, theta[0] == 0.001, theta'[0] == 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}} *)
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – anderstood
    Aug 27, 2014 at 1:32
  • $\begingroup$ @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. $\endgroup$
    – Michael E2
    Aug 27, 2014 at 1:35
  • $\begingroup$ +1 for the WhenEvent[] idea. Never thought about that use before. Great. $\endgroup$ Aug 27, 2014 at 5:49
  • $\begingroup$ Nowadays, PeakDetect[] is a usable alternative for the first approach in this answer. $\endgroup$ Dec 9, 2016 at 20:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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