Mathematica - Find Maximum value NDSolve Plot

3
votes

After numerically solving a differential equation and plotting the results I would like to determine the single maximum value in the plotted range but do not know how.

The code below works for numerically solving the differential equation and plotting the results.

s = NDSolve[{x''[t] + x[t] - 0.167 x[t]^3 == 0.005 Cos[t + -0.0000977162*t^2/2], x[0] == 0, x'[0] == 0}, x, {t, 0, 1000}]

Plot[Evaluate[x[t] /. s], {t, 0, 1000}, 
Frame -> {True, True, False, False}, FrameLabel -> {"t", "x"}, FrameStyle -> Directive[FontSize -> 15], Axes -> False]

Mathematica graphics

2
wolfram-mathematica
Allow me to welcome you to Stack Overflow and remember you just three things we usually do here: 1) As you receive help, try to give it too answering questions in your area of expertise 2) Read the FAQ!! 3) When you see good questions and answers, upvote them using the gray triangles, as the credibility of the system is based on the reputation that users come to gain by sharing their knowledge. Also remember to accept the answer that better solves your problem, if any, by pressing the checkmark sign i.imgur.com/uqJeW.png - Dr. belisarius

2 Answers

4
votes

Use NMaximize

First approximation: 

s = NDSolve[{x''[t] + x[t] - 0.167 x[t]^3 ==  
            0.005 Cos[t + -0.0000977162*t^2/2], x[0] == 0, x'[0] == 0}, x[t], 
            {t, 0, 1000}]
NMaximize[{Evaluate[x[t] /. s[[1]]] , 100 < t < 1000}, t]  

{1.26625, {t -> 821.674}}

As your function is a rapid oscillation like this : alt text , it doesn't catch the real max value, as you may see below: 

Plot[{1.26625, Evaluate[x[t] /. s[[1]]]}, {t, 790, 830}, 
 Frame -> {True, True, False, False}, FrameLabel -> {"t", "x"}, 
 FrameStyle -> Directive[FontSize -> 15], Axes -> False, 
 PlotRange -> {{790, 830}, {1.25, 1.27}}]

alt text

So we refine our bounds, and tune a little the NMaximize function: 

NMaximize[{Evaluate[x[t] /. s[[1]]] , 814 < t < 816}, t, 
 AccuracyGoal -> 20, PrecisionGoal -> 18, MaxIterations -> 1000]  

NMaximize::cvmit: Failed to converge to the requested accuracy or 
                  precision within 1000 iterations. >>

{1.26753, {t -> 814.653}}

It failed to converge within the required accuracy, but now the result is good enough

Plot[{1.2675307922753962`, Evaluate[x[t] /. s[[1]]]}, {t, 790, 830}, 
 Frame -> {True, True, False, False}, FrameLabel -> {"t", "x"}, 
 FrameStyle -> Directive[FontSize -> 15], Axes -> False, 
 PlotRange -> {{790, 830}, {1.25, 1.27}}]

alt text

3
votes

You can use Reap and Sow to extract a list of values from any evaluation. For a simple Plot you would Sow the value of the function you are plotting and enclose the entire plot in a Reap:

list = Reap[
          Plot[Sow@Evaluate[x[t] /. s], {t, 0, 1000}, 
          Frame -> {True, True, False, False},
          FrameLabel -> {"t", "x"}, 
          FrameStyle -> Directive[FontSize -> 15],
          Axes -> False]];

The first element of list is the plot itself and the second element is the list of x-values Mathematica used in the plot. To get the Maximum:

In[1]  := Max[lst[[2, 1]]]
Out[1] := 1.26191