I would like to ask whether the construction of the simulation result can be controlled in the following way - this is the effective use of Mathematica, and if there is a more “functional” way to do it. (Sow, Reap, etc. can be used).
The problem is basic. Suppose you want to simulate a physical process, say, a pendulum, and you want to build time series of a solution (that is, time against an angle) as it starts (or any other type of result).
To be able to display a graph, you need to save the data points as they are launched.
Below is a simple example that displays a solution, but only the current point, and not a complete time series:
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) ]
The above is uninteresting, since it sees only the movement of the point, and not the full path of the solution.
The way I'm currently processing this allocates, using Table[] , a buffer large enough to store the maximum possible size of the time series that can be generated.
The problem is that the time step can change, and the smaller it is, the more data will be generated.
But since I know the smallest possible time step (which is 0.1 seconds in this example), and I know the total time to start (which is 10 seconds here), then I know how much to allocate.
I also need an “index” to track the buffer. Using this method, this is the way to do this:
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) ]
For reference, when I do something like the above in Matlab, it has a nice build object called “hold”. So you can build a point, and then say “hold”, which means that the next plot will not erase what is already on the plot, but will add it.
I did not find something similar in Mathematica, i.e. updated the current chart on the fly.
I also did not want to use Append [] and AppendTo [] to create a buffer during its launch, as this would be slow and inefficient.
My question is: is there a more efficient Mathematica way (which can be faster and more elegant) to perform a typical task like the one above, besides what I do?
thanks,
UPDATE:
On the question of why not solve the ODE all at once. Yes, it is possible, but it simplifies things to do it in pieces, also for performance reasons. Here is an example with ode with initial conditions:
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) ]
Now in one it was necessary to solve the system once earlier, then they need to monitor whether the IC is changing. It can be done, but you need additional logic, and I have done it many times, but it complicates the situation a bit. I wrote a short note about this here.
In addition, I noticed that I can get much better speed by solving the system for smaller segments of time as time march than all at once. The overhead of NDSolve is very small. But when the length of time for NDsolve is long, problems can arise when you request higher accuracy from NDSolve, as in the AccuracyGoal ->, PrecisionGoal -> parameters, which I could not when the time interval is very large.
In general, the overhead of calling NDSolve for smaller segments seems to be much less comparable to the advantages it makes in simplifying logic and speed (maybe more accurate, but I haven't tested it more yet). I know that it’s a little strange to maintain an NDSolve call, but after trying both methods (right away, but add logic to check for other control variables) compared to this method, I am now inclined to this.
UPDATE 2
I compared the following 4 methods for 2 test cases:
tangle [j] [j] method (Belisarius) AppendTo (suggested by Sjoerd) Dynamic linked list (Leonid) (with and without SetAttributes[linkedList, HoldAllComplete] ) preallocate buffer (Nasser)
The way I did it was to launch it in 2 cases, one for 10,000 points, and the second for 20,000 points. I left the Plot [[] command there, but not displaying it on the screen, this eliminates any overhead of the actual rendering.
I used Timing [] around the Do loop, which iterates through the core logic called NDSolve, and iterates over the time interval using delT increments, as mentioned above. Manipulate was not used.
I used Quit [] before each run.
For Leonid’s method, I changed the [] column that he had in the Do loop. I confirmed at the end, but built the data using its getData [] method so that the result is in order.
All the code I used is below. I made a table that shows the results for 10,000 points and 20,000. Time in seconds:
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}} ]
It is clear that if I did not do something wrong, but the code below is for anyone to check, the Leonid method easily wins here. I was also surprised that AppendTo did the same as the buffer method that pre-allocated the data.
Here is the slightly modified code that I used to create the above results.
buffer method
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 Method
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} ] ]
tangle [j] [j] method
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} ] ]
dynamic linked list method
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} ] ] ] ]
Thanks for helping everyone.