【发布时间】：2011-09-05 15:07:40
【问题描述】：
我想问一下我管理模拟绘制结果的以下方式是否是对 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 中执行上述操作时，它有一个很好的绘图工具，称为“hold on”。这样人们就可以绘制一个点，然后说“坚持”，这意味着下一个情节不会删除情节上已经存在的内容，而是会添加它。
我在 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是否发生变化。这可以做到，但需要额外的逻辑，我以前做过很多次，但它确实使事情复杂了一点。我在这个here.写了一个小笔记
另外，我注意到，随着时间的推移，我可以通过解决系统更小的时间段来获得更快的速度，而不是一次解决整个问题。 NDSolve 调用开销非常小。但是当 NDsolve 的持续时间很长时，当人们要求 NDSolve 提供更高的精度时，可能会出现问题，例如选项AccuracyGoal ->, PrecisionGoal ->，当时间间隔非常大时我无法做到这一点。
总体而言，与它在简化逻辑和速度方面的优势相比，为较小的段调用 NDSolve 的开销似乎要少得多（可能更准确，但我没有对此进行更多检查）。我知道继续调用 NDSolve 似乎有点奇怪，但是在尝试了这两种方法（一次全部，但添加逻辑以检查其他控制变量）与此方法之后，我现在倾向于此方法。
更新 2
我针对 2 个测试用例比较了以下 4 种方法：
tangle[j][j] 方法（贝利撒留）
AppendTo（由 Sjoerd 建议）
动态链表（Leonid）（带和不带SetAttributes[linkedList, HoldAllComplete]）
预分配缓冲区（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 和预分配数据的 buffer 方法一样好。
这是我用来生成上述结果的稍微修改的代码。
缓冲方法
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}
]
]
]
]
感谢大家的帮助。
【问题讨论】：
-
你有什么理由不只保留整个系列的数据吗？它真的不会占用那么多空间。您可以求解整个 t 值范围，它会生成一个插值函数，该函数存储最少量的数据来生成该曲线。
-
我现在正在比较这里讨论的方法对于相同测试用例的性能，稍后会发布结果。