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.

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

2 Answers 2

up vote 12 down vote accepted

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}]
share|improve this answer
3  
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 at 19:40
1  
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 at 19:40
2  
@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 at 19:50
1  
@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 at 20:14
1  
@Karsten7. I think that's why we speak of numerical recipes! It's really like cooking :) –  anderstood Aug 26 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.

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 ploting 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 them might be returned instead of the global one. But they would have to be very close, within error of the NDSolve computation. If one wanted to be sure, 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}} *)
share|improve this answer
    
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 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 at 1:35
    
+1 for the WhenEvent[] idea. Never thought about that use before. Great. –  belisarius Aug 27 at 5:49

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.