How to transform Do loop to more efficient codes?

Question

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.

Answer 1

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}

Answer 2

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]

Answer 3

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.

Answer 4

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 *)