Mathematica - найти максимальное значение NDSolve Plot

Question

После численного решения дифференциального уравнения и построения графика результатов я хотел бы определить единственное максимальное значение в указанном диапазоне, но не знаю, как.

Приведенный ниже код работает для численного решения дифференциального уравнения и построения результатов.

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]

Answer 1

Использование NMaximize

Первое приближение:

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}}  

Поскольку ваша функция - это быстрое колебание, подобное этому: , оно не улавливает реальноемаксимальное значение, как вы можете видеть ниже:

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}}]

Итак, мы уточним наши границы и немного настроим функцию NMaximize:

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}}  

Не удалось сойтись с требуемой точностью, но теперь результат достаточно хороший

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}}]

Answer 2

Вы можете использовать Reap и Sow, чтобы извлечь список значений из любой оценки. Для простого Plot вы бы Sow указали значение функции, которую вы строите, и заключили бы весь график в 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]];

Первым элементом list является сам график, а вторым элементом является список значений x Mathematica, используемых на графике. Чтобы получить максимум:

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