# How to transform Do loop to more efficient codes?

I'm new to Mathematica. Here is my original program: How to rewrite the Do part?

Clear["Global*"];
da = {};
dt = 0.02;
endtime = 2;
Q = dQ = {0, 0, 0};
G = {0, -9.8, 0};
Do[
ddQ = G - Q;
dQ = dQ + ddQ*dt;
Q = Q + dQ*dt;
da = Append[da, dQ[[2]]]
, {t, 0, endtime, dt}
]
ListPlot[Flatten[da], DataRange -> {0, endtime}, PlotRange -> All]
Quit[];


Where Q denotes coordinates of a particle. ddQ denotes the acceleration which depends on the coordinates Q of the particle. In order to make the posted program readable and clean, I provided acceleration value in one element only.

-
Why bother with 3-dimensional vectors if you're only using the second component ? –  b.gatessucks Mar 12 '13 at 9:15

I had the same idea as Mr.Wizard, except that I thought it was even neater to treat both Q and dQ as part of the one state vector, and then select the column you actually want.

da2 = Rest@ NestList[#1 + dt*With[{ddq = (G - #1[[1 ;; 3]])},
Join[#1[[4 ;; 6]] + dt*ddq, ddq]] &, ConstantArray[0, {6}],
Floor[ endtime/dt] + 1][[All, 5]]


Notice the use of With to avoid calculating ddq twice. I have used ddq instead of your ddQ to avoid namespace clashes.

I note that you are updating dQ, which depends on the previous value of Q and then applying that to update the value of Q. So your problem isn't strictly iterating a single state vector. I wonder if that is really what you are intending to do, or if you actually mean to update Q with the previous value of dQ.

On my machine (a four-year-old MacBook Pro), my version is about 40 times faster than the version in your question using Do loops. Setting endtime=500 gives:

Do[ddQ = G - Q;
dQ = dQ + ddQ*dt;
Q = Q + dQ*dt;
da = Append[da, dQ[[2]]], {t, 0, endtime, dt}] // Timing

{2.624312, Null}


Versus

da2 = Rest@NestList[#1 +
dt*With[{ddq = (G - #1[[1 ;; 3]])},
Join[#1[[4 ;; 6]] + dt*ddq, ddq]] &, ConstantArray[0, {6}],
Floor[ endtime/dt] + 1][[All, 5]]; // Timing

{0.053916, Null}

-
Thanks Verbeia! The time difference is remarkable, on my laptop I got 10.14s, 0.78s and 0.094s runtime of 50000 iteration for Do[], Table[] and Rest@NestList[] respectively. Each about 10 times faster than the preceding. However, I'm worrying about the readablility of the program when the content in the loop becomes complex. –  novice Mar 13 '13 at 2:32

Many ways to do this in M. If you want to use a Do then you can use Reap and Sow to collect data. I changed all your UpperFirstCase letters with lowerFirstCase letter. Not good idea to use UPPERCASE in M.

## Using Do

Clear["Global*"];
dt = 0.02;
endtime = 2;
q = dQ = {0, 0, 0};
g = {0, -9.8, 0};

da = Rest@Reap[Do[ddQ = g - q;
dQ = dQ + ddQ*dt;
q = q + dQ*dt;
Sow[dQ[[2]]],
{t, 0, endtime, dt}
]
];


and now your data is in da which you can plot as you did

ListPlot[Flatten[da], DataRange -> {0, endtime}, PlotRange -> All]


## Another way using a Table

Clear["Global*"];
dt = 0.02;
endtime = 2;
q = dQ = {0, 0, 0};
g = {0, -9.8, 0};
nSteps = Length[Range[0, endtime, dt]];

makeStep[] := (ddQ = g - q; dQ = dQ + ddQ*dt; q = q + dQ*dt; dQ[[2]]);

da = Table[makeStep[], {nSteps}];

ListPlot[da, DataRange -> {0, endtime}, PlotRange -> All]

-

First, it is a good idea not to start any user symbols with capital letters, as these may conflict with internal functions. In the code below I will preserve your symbol names for easy reference.

For each method:

init[]:= (
dt = 0.02;
endtime = 2;
Q = dQ = {0, 0, 0};
G = {0, -9.8, 0};
)


Table:

init[]

Table[
dQ += (G - Q) dt;
Q  += dQ dt;
dQ[[2]],
{t, 0, endtime, dt}
] // ListPlot


NestList:

NestList[
With[{x = #[[1]] + (G - #[[2]]) dt}, {x, #[[2]] + x dt}] &,
{dQ, Q},
1 + Quotient[endtime, dt]
][[2 ;;, 1, 2]] // ListPlot


The use of (some) pure functions in NestList is particularly fast because it auto-compiles after (by default) 100 iterations (as controlled by SystemOptions["CompileOptions"]).

One may also observe that the outer two elements of your Q, dQ, and G lists are never used. Therefore, for the given output you could write this, which is quite fast:

dt = 0.02;
endtime = 20000;
Q = dQ = 0;
G = -9.8;

NestList[
With[{x = #[[1]] + (G - #[[2]]) dt}, {x, #[[2]] + x dt}] &,
{dQ, Q},
1 + Quotient[endtime, dt]
][[2 ;;, 1]] // Length // Timing

{0.187, 1000001}


Not bad for a million iterations, at least as Mathematica goes.

-

Yet another way :

da == Rest@NestList[{#[[1]] + G[[2]] dt - #[[2]] dt,
#[[2]] + (#[[1]] + G[[2]] dt - #[[2]] dt) dt} &,
{0, 0}, Ceiling[endtime/dt] + 1][[All, 1]]
(* True *)
`
-