Mathematica - 查找最大值 NDSolve 图答案

【问题标题】：Mathematica - Find Maximum value NDSolve PlotMathematica - 查找最大值 NDSolve 图
【发布时间】：2010-12-23 17:13:18
【问题描述】：

在数值求解微分方程并绘制结果后，我想确定绘制范围内的单个最大值，但不知道如何。

下面的代码用于数值求解微分方程并绘制结果。

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]

【问题讨论】：

请允许我欢迎您来到 Stack Overflow，并记住我们通常在这里做的三件事：1) 当您获得帮助时，请尝试在您的专业领域 2) 阅读常见问题解答！！ 3) 当您看到好的问题和答案时，使用灰色三角形为它们投票，因为系统的可信度基于用户通过分享他们的知识而获得的声誉。还记得接受更好地解决您的问题的答案，如果有的话，按复选标记 i.imgur.com/uqJeW.png

标签： wolfram-mathematica

【解决方案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}}]

【讨论】：

有什么理由不喜欢FindMaximum[]？ :)
@J. M. 因为 NMaximize 试图找到一个 GLOBAL 最大值，而 FindMaximum 对 LOCAL 极端值很有效。对于振荡函数，NMaximize 对我来说效果更好。请参阅此处的示例reference.wolfram.com/mathematica/ref/FindMaximum.html 基本示例下的第一个示例
@J.不过，M. 本可以第二次使用 FindMaximum[]。
嗯，它是一个 单变量 函数，肯定会在尝试找到根或最优值之前绘制它，对吧？我这样说是因为在单变量情况下引入NMaximize[]（默认为 Nelder-Mead）的机制有点浪费，而FindMaximum[] 的方法足够了。当然，在优化之前进行绘图的原因是，如果没有给定一个好的起点（即 GIGO），这些迭代方法往往会误入歧途。
@J. M. 由于具有快速振荡功能，有时我发现 Plot[] 具有欺骗性。这就是为什么我通常会尝试 NMaximize。至于浪费部分，我不同意......因为它正在运行时我得到了我的咖啡:)

【解决方案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 的第一个元素是绘图本身，第二个元素是绘图中使用的 Mathematica 的 x 值列表。获得最大值：

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

【讨论】：

我猜使用 Plot[] 中的一些选项可以深入扫描冲突的尖峰
此外，Sowed 值在实际最大值左侧的相邻峰值中...