2010-12-23 82 views
3

数值求解微分方程并绘制结果后,我想确定绘图范围内的单个最大值,但不知道如何。Mathematica - 查找最大值NDSolve曲线

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

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

+0

请允许我欢迎您来到Stack Overflow,并记住我们通常在这里做的三件事:1)当您收到帮助时,请尝试t o也给予答案**在您的专业领域回答问题** 2)**阅读常见问题解答!! ** 3)当您看到好的问题和答案时,请使用灰色三角形**作为可信度系统的基础是用户通过分享知识获得的声誉。还请记住接受更好地解决您的问题的答案,如果有的话,**按复选标记** http://i.imgur.com/uqJeW.png – 2010-12-27 01:54:20

回答

4

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

由于你的函数是一种快速振荡是这样的:alt text,它并没有赶上真正的最大值,如下所示:

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

所以我们完善我们的边界,并调整了一点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}}] 

alt text

3

您可以使用ReapSow从任何评价中提取值的列表。对于一个简单的Plot你会Sow你绘制函数的值和包围在一个Reap整个剧情:

list = Reap[ 
      Plot[[email protected][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 
+0

我想在剧情[]中使用一些选项是可能的扫描深度冲突尖峰 – 2010-12-23 21:37:31

+0

此外,Sowed值位于实际最大值左边的相邻峰值... – 2010-12-23 22:49:48