一种有效的数据结构或方法,用于管理随时间增长的绘图数据

Question

我想问一下,我管理绘图结果的以下方式是否有效地使用Mathematica,以及是否有更"实用"的方法.(可能正在使用Sow,Reap等).

问题是基本问题.假设您想要模拟一个物理过程,例如钟摆,并希望在运行时(或任何其他类型的结果)绘制解决方案的时间序列(即时间与角度).

为了能够显示绘图,需要在运行时保留数据点.

以下是一个简单的示例,它绘制了解决方案,但只绘制了当前点,而不是完整的时间序列:

Manipulate[
 sol = First@NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == Pi/4, y'[0] == 0}, 
                      y, {t, time, time + 1}];

 With[{angle = y /. sol},
  (
   ListPlot[{{time, angle[time]}}, AxesLabel -> {"time", "angle"}, 
            PlotRange -> {{0, max}, {-Pi, Pi}}]
   )
  ],

 {{time, 0, "run"}, 0, max, Dynamic@delT, ControlType -> Trigger},
 {{delT, 0.1, "delT"}, 0.1, 1, 0.1, Appearance -> "Labeled"},

 TrackedSymbols :> {time},
 Initialization :> (max = 10)
 ]

上面的内容并不有意思,因为只看到一个点移动,而不是完整的解决方案路径.

目前我处理这个问题的方法是使用Table[]一个足够大的缓冲区来分配可以生成的最大可能时间序列大小.

问题是时间步长可以改变,而且越小,产生的数据就越多.

但是因为我知道可能的最小时间步长(在这个例子中是0.1秒),并且我知道运行的总时间(这里是10秒),然后我知道要分配多少.

我还需要一个'索引'来跟踪缓冲区.使用此方法,这是一种方法:

Manipulate[

 If[time == 0, index = 0];

 sol = First@NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == Pi/4,y'[0] == 0}, 
                      y, {t, time, time + 1}];

 With[{angle = y /. sol},
  (
   index += 1;
   buffer[[index]] = {time, angle[time]};

   ListPlot[buffer[[1 ;; index]], Joined -> True, AxesLabel -> {"time", "angle"}, 
            PlotRange -> {{0, 10}, {-Pi, Pi}}]
   )
  ],

 {{time, 0, "run"}, 0, 10, Dynamic@delT, AnimationRate -> 1, ControlType -> Trigger},
 {{delT, 0.1, "delT"}, 0.1, 1, 0.1, Appearance -> "Labeled"},

 {{buffer, Table[{0, 0}, {(max + 1)*10}]}, None},
 {{index, 0}, None},

 TrackedSymbols :> {time},
 Initialization :> (max = 10)
 ]

作为参考,当我在Matlab中执行上述操作时,它具有很好的绘图功能,称为"保持".因此,人们可以绘制一个点,然后说'保持',这意味着下一个绘图将不会删除已经在绘图上的内容,但会添加它.

我没有在Mathematica中找到类似的东西,即即时更新当前的情节.

我也不想在运行时使用Append []和AppendTo []来构建缓冲区,因为缓慢而且效率不高.

我的问题:除了我正在做的事情之外,是否有更高效的Mathematica方式(可以更快,更优雅)来完成上述典型任务？

谢谢,

更新:

关于为什么不立刻解决ODE的问题.是的,这是可能的,但是出于性能原因,它简化了很多事情.以下是带有初始条件的ode的示例:

Manipulate[
 If[time == 0, index = 0];
 sol = First@
   NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == y0, 
     y'[0] == yder0}, y, {t, time, time + 1}];

 With[{angle = (y /. sol)[time]},
  (
   index += 1;
   buffer[[index]] = {time, angle};

   ListPlot[buffer[[1 ;; index]], Joined -> True, 
    AxesLabel -> {"time", "angle"}, 
    PlotRange -> {{0, 10}, {-Pi, Pi}}])],

 {{time, 0, "run"}, 0, 10, Dynamic@delT, AnimationRate -> 1, 
  ControlType -> Trigger}, {{delT, 0.1, "delT"}, 0.1, 1, 0.1, 
  Appearance -> "Labeled"},
 {{y0, Pi/4, "y(0)"}, -Pi, Pi, Pi/100, Appearance -> "Labeled"},
 {{yder0, 0, "y'(0)"}, -1, 1, .1, Appearance -> "Labeled"},

 {{buffer, Table[{0, 0}, {(max + 1)*10}]}, None},
 {{index, 0}, None},

 TrackedSymbols :> {time},
 Initialization :> (max = 10)
 ]

现在,其中一个是先解决系统,然后他们需要注意IC是否发生变化.这可以做到,但需要额外的逻辑,我已经多次这样做了,但它确实使事情复杂化了一些.我在这里写了一个小小的说明.

此外,我注意到,随着时间的推移,我可以通过在较小的时间段内解决系统而获得更好的速度,而不是一次性完成.NDSolve调用开销非常小.但是当NDsolve的持续时间很长时,如果从NDSolve请求更高的精度,就会出现问题,就像在选项中那样AccuracyGoal ->, PrecisionGoal ->,当时间间隔非常大时,我无法做到这一点.

总体而言,与简化逻辑和速度(可能更准确,但我没有更多地检查)相比,为较小的段调用NDSolve的开销似乎要少得多.我知道继续调用NDSolve似乎有点奇怪,但在尝试了两种方法(一次性完成,但添加逻辑以检查其他控制变量)与此方法之后,我现在倾向于这一方法.

更新2

我将2个测试用例的以下4种方法进行了比较:

纠结[j] [j]方法(Belisarius)AppendTo(由Sjoerd建议)动态链表(Leonid)(有和没有SetAttributes[linkedList, HoldAllComplete])preallocate buffer(Nasser)

我这样做的方法是运行2个案例,一个运行10,000个点,第二个运行20,000个点.我确实在那里留下了Plot [[]命令,但是没有在屏幕上显示它,这是为了消除实际渲染的任何开销.

我在一个Do循环周围使用了Timing [],循环遍历调用NDSolve的核心逻辑,并使用上面的delT增量迭代时间跨度.没有使用Manipulate.

我在每次运行前都使用了Quit [].

对于Leonid方法,我通过Do循环改变了他的Column [].我最后验证了,但是使用他的getData []方法绘制数据,结果还可以.

我使用的所有代码如下.我制作了一张表格,显示10,000点和20,000点的结果.时间是每秒:

 result = Grid[{
   {Text[Style["method", Bold]], 
    Text[Style["number of elements", Bold]], SpanFromLeft},
   {"", 10000, 20000},
   {"", SpanFromLeft},
   {"buffer", 129, 571},
   {"AppendTo", 128, 574},
   {"tangle[j][j]", 612, 2459},
   {"linkedList with SetAttribute", 25, 81},
   {"linkedList w/o SetAttribute", 27, 90}}
  ]

显然,除非我做错了什么,但代码低于任何人都要验证,Leonid方法在这里很容易获胜.我也很惊讶AppendTo和预先分配数据的缓冲方法一样好.

以下是我用于生成上述结果的略微修改的代码.

缓冲方法

delT = 0.01; max = 100; index = 0;
buffer = Table[{0, 0}, {(max + 1)*1/delT}];

Timing[
 Do[
  sol = First@
    NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == Pi/4, 
      y'[0] == 0}, y, {t, time, time + 1}];

  With[{angle = y /. sol},
   (index += 1;
    buffer[[index]] = {time, angle[time]};
    foo = 
     ListPlot[buffer[[1 ;; index]], Joined -> True, 
      AxesLabel -> {"time", "angle"}, 
      PlotRange -> {{0, 10}, {-Pi, Pi}}]
    )
   ], {time, 0, max, delT}
  ]
 ]

AppendTo方法

Clear[y, t];
delT = 0.01; max = 200;
buffer = {{0, 0}};  (*just a hack to get ball rolling, would not do this in real code*)

Timing[
 Do[
  sol = First@
    NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == Pi/4, 
      y'[0] == 0}, y, {t, time, time + 1}];

  With[{angle = y /. sol},
   (AppendTo[buffer, {time, angle[time]}];
    foo = 
     ListPlot[buffer, Joined -> True, AxesLabel -> {"time", "angle"}, 
      PlotRange -> {{0, 10}, {-Pi, Pi}}]
    )
   ], {time, 0, max, delT}
  ]
 ]

缠结[j] [j]方法

Clear[y, t];
delT = 0.01; max = 200;
Timing[
 Do[
  sol = First@
    NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == Pi/4, 
      y'[0] == 0}, y, {t, time, time + 1}];
  tangle[time] = y /. sol;
  foo = ListPlot[
    Table[{j, tangle[j][j]}, {j, .1, max, delT}],
    AxesLabel -> {"time", "angle"},
    PlotRange -> {{0, max}, {-Pi, Pi}}
    ]
  , {time, 0, max, delT}
  ]
 ]

动态链表方法

Timing[
 max = 200;

 ClearAll[linkedList, toLinkedList, fromLinkedList, addToList, pop, 
  emptyList];
 SetAttributes[linkedList, HoldAllComplete];
 toLinkedList[data_List] := Fold[linkedList, linkedList[], data];
 fromLinkedList[ll_linkedList] := 
  List @@ Flatten[ll, Infinity, linkedList];
 addToList[ll_, value_] := linkedList[ll, value];
 pop[ll_] := Last@ll;
 emptyList[] := linkedList[];

 Clear[getData];

 Module[{ll = emptyList[], time = 0, restart, plot, y},

  getData[] := fromLinkedList[ll];

  plot[] := Graphics[
    {
     Hue[0.67`, 0.6`, 0.6`],
     Line[fromLinkedList[ll]]
     },
    AspectRatio -> 1/GoldenRatio,
    Axes -> True,
    AxesLabel -> {"time", "angle"},
    PlotRange -> {{0, 10}, {-Pi, Pi}},
    PlotRangeClipping -> True
    ];

  DynamicModule[{sol, angle, llaux, delT = 0.01},

   restart[] := (time = 0; llaux = emptyList[]);
   llaux = ll;

   sol := 
    First@NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == Pi/4, 
       y'[0] == 0}, y, {t, time, time + 1}];
   angle := y /. sol;

   ll := With[{res = 
       If[llaux === emptyList[] || pop[llaux][[1]] != time,
        addToList[llaux, {time, angle[time]}],
        (*else*)llaux]
      },
     llaux = res
     ];

   Do[
    time += delT;
    plot[]
    , {i, 0, max, delT}
    ]
   ]
  ]
 ]

谢谢大家的帮助.

Answer 1

我不知道如何得到你想要的东西Manipulate，但我似乎已经设法得到一些接近定制的东西Dynamic。以下代码将：使用链表来提高效率，使用按钮停止/恢复绘图，并在任何给定时间按需提供迄今为止收集的数据：

ClearAll[linkedList, toLinkedList, fromLinkedList, addToList, pop, emptyList];
SetAttributes[linkedList, HoldAllComplete]; 
toLinkedList[data_List] := Fold[linkedList, linkedList[], data];
fromLinkedList[ll_linkedList] := List @@ Flatten[ll, Infinity, linkedList];
addToList[ll_, value_] := linkedList[ll, value];
pop[ll_] := Last@ll;
emptyList[] := linkedList[];


Clear[getData];
Module[{ll = emptyList[], time = 0, restart, plot, y},
   getData[] := fromLinkedList[ll];
   plot[] := 
      Graphics[{Hue[0.67`, 0.6`, 0.6`], Line[fromLinkedList[ll]]}, 
        AspectRatio -> 1/GoldenRatio, Axes -> True, 
        AxesLabel -> {"time", "angle"}, PlotRange -> {{0, 10}, {-Pi, Pi}}, 
        PlotRangeClipping -> True];
   DynamicModule[{sol, angle, llaux, delT = 0.1},
     restart[] := (time = 0; llaux = emptyList[]);
     llaux = ll;
     sol := First@
        NDSolve[{y''[t] + 0.1 y'[t] + Sin[y[t]] == 0, y[0] == Pi/4, y'[0] == 0}, 
             y, {t, time, time + 1}];
     angle := y /. sol;
     ll := With[{res =
              If[llaux === emptyList[] || pop[llaux][[1]] != time, 
                 addToList[llaux, {time, angle[time]}],
                 (* else  *)
                 llaux]},
              llaux = res];
     Column[{
        Row[{Dynamic@delT, Slider[Dynamic[delT], {0.1, 1., 0.1}]}],
        Dynamic[time, {None, Automatic, None}],
        Row[{
          Trigger[Dynamic[time], {0, 10, Dynamic@delT}, 
               AppearanceElements -> { "PlayPauseButton"}], 
          Button[Style["Restart", Small], restart[]]
        }],
        Dynamic[plot[]]
      }, Frame -> True]
  ]
]

这里的链接列表取代了您的列表buffer，您不需要预先分配并提前知道您将拥有多少个数据点。这plot[]是一个自定义的低级绘图函数，尽管我们也可以使用ListPlot. 您可以使用“播放”按钮来停止和恢复绘图，并使用自定义的“重新启动”按钮来重置参数。

您可以在任何给定时间调用getData[]以获取迄今为止累积的数据列表，如下所示：

In[218]:= getData[]
Out[218]= {{0,0.785398},{0.2,0.771383},{0.3,0.754062},{0.4,0.730105},{0.5,0.699755},
{0.6,0.663304},{0.7,0.621093},{0.8,0.573517},{0.9,0.521021},{1.,0.464099},
{1.1,0.403294},{1.2,0.339193},{1.3,0.272424}}