程序多次运行结果不同

Question

Clear["Global`*"]
pw[t_] := Piecewise[{{1.0, 0 <= Mod[t, 5.1129] < 1.085},
    {1.5, 1.085 <= Mod[t, 5.1129] < 4.0279}},
   1.0];

sl = NDSolve[{x'[t] == y[t], 
   y'[t] == pw[t]*x[t] - (-2 + 2 x[t]^4)/x[t]^3,
   x[0] == 1.0, y[0] == 0}, {x, y}, {t, 0, 5.1128}]
x = x /. sl[[1, 1]];
y = y /. sl[[1, 2]];
ParametricPlot[Evaluate[{x[t], y[t]}], {t, 0, 5.1128}, 
 PlotRange -> All]

第一种情况

第二种情况

Comments

已经解决了，是不连续的原因，下边有解释，O(∩_∩)O谢谢！

---

基本确定是版本9引入的Bug，刚在版本8.0.4重复执行了50次，没有问题。

---

难道是9版本之后加入了不连续处理吗，不过你调试能力真不错，`Iterate，`StateData和NDSolve`ProcessEquations我都还没接触呢

---

哈……我也是最近才开始深入的。"DiscontinuityProcessing"应该是版本9引入的，版本8文档里查不到，硬加上就会冒出NDSolve没这个Method的警告。

---

在SE发现以前有人发现过这个问题，链接见我下面的答案。

Answer 1

可能是分段函数不连续导致的。

其实我也不太了解，帮助上说“NDSolve自动对诸如Sign的不连续函数进行处理”，“在一些时间积分方法下，解可能是相当不准确的”。所以关掉不连续处理方法之后就好了，但在不连续点可能失败。

详细的，还是自己看一下NDSolve帮助吧（Option > Method > DiscontinuityProcessing ）

sl = NDSolve[{x'[t] == y[t], 
   y'[t] == pw[t]*x[t] - (-2 + 2 x[t]^4)/x[t]^3, x[0] == 1.0, 
   y[0] == 0}, {x, y}, {t, 0, 5.1128}, 
  Method -> {"DiscontinuityProcessing" -> False}]

 

Comments

你的mma哪个版本的，上个朋友说是可能版本问题

---

mathematica 11.0.1

---

可能还是函数不连续引起的，我再看看吧

---

确实是，加入Method -> {"DiscontinuityProcessing" -> False}就好了，谢谢了

Answer 2

更新：

在Stackexchange翻了半天，发现以前有人发现过这个问题：

http://mathematica.stackexchange.com/q/21262/1871

---

不完整回答，算是对 @野鹤 答案的补充。可以确定：

1.问题和浮点数无关。

2.问题出在前处理上，但是不确定是前处理出错还是某些情况下前处理产生的数据不适合在后续步骤使用。

证据：

Clear[x, y]
pw[t_] = Rationalize[
   Piecewise[{{1.0, 0 <= Mod[t, 5.1129] < 1.085}, {1.5, 
      1.085 <= Mod[t, 5.1129] < 4.0279}}, 1.0], 0];
{state} = NDSolve`ProcessEquations[{x'[t] == y[t], 
    y'[t] == pw[t]*x[t] - (-2 + 2 x[t]^4)/x[t]^3, x[0] == 1, y[0] == 0}, {x, y}, {t, 0, 
    5}(*,Method\[Rule]{"DiscontinuityProcessing"\[Rule]False}*), WorkingPrecision -> 16];
state // InputForm
NDSolve`Iterate[state, 5];
{xsol, ysol} = {x, y} /. NDSolve`ProcessSolutions@state;
Plot[{xsol@t, ysol@t}, {t, 0, 5}]

反复执行直到结果开始出错，然后Ctrl+Shift+D从NDSolve`Iterate这行断开，再执行，会发现解的图形不再变化，也就是说“不同”是在NDSolve`ProcessEquations这步引入的，再对比两种情况下的state的InputForm，确实可以看到某些地方有不同，只可惜NDSolve`StateData的具体结构文档中没讲解，没法进一步分析。

版本8无此问题，因为那时候还不存在"DiscontinuityProcessing"这东西，看样子是这新加的步骤不稳定。

Comments

Stackexchange也是像这样的一个交流论坛吗？

---

这边应该算是参考那里建的吧……