Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

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

4 Answers 4

up vote 4 down vote accepted

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}
share|improve this answer
    
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
add comment

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]

Mathematica graphics

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]
share|improve this answer
add comment

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

Mathematica graphics

NestList:

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

Mathematica graphics

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.

share|improve this answer
add comment

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 *)
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.