An efficient data structure or method to manage plotting data that grow with time

Question

I'd like to ask if the following way I manage plotting result of simulation is efficient use of Mathematica and if there is a more 'functional' way to do it. (may be using Sow, Reap and such).

The problem is basic one. Suppose you want to simulate a physical process, say a pendulum, and want to plot the time-series of the solution (i.e. time vs. angle) as it runs (or any other type of result).

To be able to show the plot, one needs to keep the data points as it runs.

The following is a simple example, that plots the solution, but only the current point, and not the full 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 not interesting, as one only sees a point moving, and not the full solution path.

The way currently I handle this, is allocate, using Table[], a buffer large enough to hold the largest possible time-series size that can be generated.

The issue 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 run (which is 10 seconds here), then I know how much to allocate.

I also need an 'index' to keep track of the buffer. Using this method, here is a way to do it:

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 facility for plotting, called 'hold on'. So that one can plot a point, then say 'hold on' which means that the next plot will not erase what is already on the plot, but will add it.

I did not find something like this in Mathematica, i.e. update a current plot on the fly.

I also did not want to use Append[] and AppendTo[] to build the buffer as it runs, as that will be slow and not efficient.

My question: Is there a more efficient, Mathematica way (which can be faster and more elegent) to do a typical task such as the above, other than what I am doing?

thanks,

UPDATE:

On the question on why not solving the ODE all at once. Yes, it is possible, but it simplifies things alot 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 were to solve the system once before, then they need to watch out if the IC changes. This can be done, but need extra logic and I have done this before many times, but it does complicate things a bit. I wrote a small note on this here.

Also, I noticed that I can get much better speed by solving the system for smaller time segments as time marches on, than the whole thing at once. NDSolve call overhead is very small. But when the time duration to NDsolve for is large, problems can result when one ask for higher accuracy from NDSolve, as in options AccuracyGoal ->, PrecisionGoal ->, which I could not when time interval is very large.

Overall, the overhead of calling NDSolve for smaller segments seems to much less compare to the advantages it makes in simplifing the logic, and speed (may be more accurate, but I have not checked on this more). I know it seems a bit strange to keep calling NDSolve, but after trying both methods (all at once, but add logic to check for other control variables) vs. this method, I am now leaning towards this one.

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 this, is by running it over 2 cases, one for 10,000 points, and the second for 20,000 points. I did leave the Plot[[] command there, but do not display it on the screen, this is to eliminate any overhead of the actual rendering.

I used Timing[] around a Do loop which iterate over the core logic which called NDSolve and iterate over the time span using delT increments as above. No Manipulate was used.

I used Quit[] before each run.

For Leonid method, I changed the Column[] he had by the Do loop. I verified at the end, but plotting the data using his getData[] method, that the result is ok.

All the code I used is below. I made a table which shows the results for the 10,000 points and 20,000. Timing is per 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}}
  ]

Clearly, unless I did something wrong, but code is below for anyone to verify, Leonid method wins easily here. I was also surprised that AppendTo did just as well as the buffer method which pre-allocated data.

Here are the slightly modified code I used to generate 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 everyone help.

Answer 1

I don't know how to get what you want with Manipulate, but I seem to have managed getting something close with a custom Dynamic. The following code will: use linked lists to be reasonably efficient, stop / resume your plot with a button, and have the data collected so far available on demand at any given time:

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]
  ]
]

Linked lists here replace your buffer and you don't need to pre-allocate and to know in advance how many data points you will have. The plot[] is a custom low-level plotting function, although we probably could just as well use ListPlot. You use the "Play" button to both stop and resume plotting, and you use the custom "Restart" button to reset the parameters.

You can call getData[] at any given time to get a list of data accumulated so far, like so:

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