Mathematica - 寻找最大值 NDSolve 绘图

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]

- Carey

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

- Dr. belisarius

---

你不喜欢 FindMaximum[] 的原因是什么呢？ :) - user414706

@J. M. 因为NMaximize试图找到全局最大值，而FindMaximum适用于局部极值。对于振荡函数，NMaximize对我来说效果更好。例如，请参见此处http://reference.wolfram.com/mathematica/ref/FindMaximum.html下的**基本示例**中的第一个示例。 - Dr. belisarius

@J. M. 第二次可以使用FindMaximum[]。 - Dr. belisarius

好的，这是一个单变量函数，在尝试查找根或极值之前，肯定会绘制该函数图形，对吧？我这么说是因为在单变量情况下，引入NMaximize[]（默认情况下为Nelder-Mead）的机制有点浪费，当然可以使用FindMaximum[]的方法。当然，在优化之前绘图的原因是这些迭代方法倾向于偏离良好的起始点（即GIGO）。 - user414706

@J. M. 随着快速振荡函数的出现，有时我发现Plot[]是具有欺骗性的。这就是为什么我通常会尝试NMaximize。至于“浪费”的部分，我不同意...因为当它在运行时，我可以享受我的咖啡时间 :) - Dr. belisarius

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

< p > list 的第一个元素是图形本身，第二个元素是 Mathematica 在绘制图形时使用的 x 值列表。要获取最大值：

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

- Timo

---

我猜测在Plot（）中使用一些选项，可以深入扫描冲突的脉冲信号。 - Dr. belisarius

此外，Sowed值位于真实最大值左侧的相邻峰值处... - Dr. belisarius