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 am new to mathematica and so just experimenting with various programming constructs. Recently have been looking at NestList and how I could use this to implement Euler's method.

Euler[a0_, b0_, steps0_, init0_] := 
 Module[{a = a0, b = b0, steps = steps0, init = init0, t, y},
  dt = (b - a)/steps;
  f[{t0_, y0_}] := y0 Sin[t0];
  euler[{t_, y_}] := {t + dt, y + dt f[{t, y}]};
  NestList[euler, {a, init}, steps]
  ]
approx = ListPlot[Euler[0, 2 \[Pi], 30, 1]]

I came up with this implementation but when I set the number of steps to > 20 the performance is very poor and the function never returns. I've tried it for a number of different functions and the results agrees with DSolve. Is it better to use For, Do, etc.?


I thought i would experiment a little more with NestList & Euler so i implemented a small example that shows how to solve a 2nd order ODE for a damped system. As noted above in this example if i set c,k,m to be an integer value vs a decimal then mathematica evaluates them differently giving very different performance characteristics.

Euler[a0_,b0_,steps0_,x0_,v0_]:=Module[{a=a0,b=b0,
steps=steps0,xinit=x0,vinit=v0},
dt= (b-a)/steps;
k=5.0;
m=2.0;
c=0.5;
f[{t_,x_,v_}]:=v;
g[{t_,x_,v_}]:= -x (k/m) -c v; 
euler[{t_,x_,v_}]:={t+dt,
x+dt f[{t,x,v}],
v+dt g[{t,x,v}]};
NestList[euler,{a,xinit,vinit},steps]
]
result =Euler[0,20,100000,1,1];
ListPlot [ {result[[All,2]],result[[All,3]]},PlotRange->All]

Damped Oscillations

share|improve this question
1  
If you wrap your definition of f[{t0_, y0_}] := y0 Sin[t0]; around a N it should improve speed: f[{t0_, y0_}] := N[y0 Sin[t0]]; –  gpap Mar 25 '13 at 15:30
    
Thanks. Tried that and it seems to work but why do i need the N[] ? I understand what N[] does just not sure why its needed in this case ? –  David McHarg Mar 25 '13 at 15:33
    
The way you had it, Euler[a0, b0, n, init] was treating the expression symbolically. Try doing NestList[Sqrt[1 + #] &, 1, 10] as opposed to NestList[N@Sqrt[1 + #] &, 1, 10] to understand the difference. –  gpap Mar 25 '13 at 15:43
    
And NO, it isn't better to use For[]/Do[] at all. You're in the right path –  belisarius Mar 25 '13 at 15:53
    
Thanks for the info, much appreciated. NestList is actually a very neat solution to this problem. –  David McHarg Mar 25 '13 at 16:21

1 Answer 1

The good thing with Mathematica is that when you supply exact numbers to your computation, it will never make any evaluation which looses information (find more about this in this answer). If you want to see this in action, you could look at the expressions building up inside your f call during the computation. You can wrap Trace around your Euler call to see this

Trace[Euler[0, 2 Pi , 7, 1], f[__]]

Mathematica graphics

Consider that this is only for 7 steps! If you are using 20+ steps, the expression will fill books which makes it kind of hard for Mathematica to handle them fast.

Let's make one little adjustment and give one argument as non-exact number. This tells Mathematica that it can calculate all expression with finite precision

Trace[Euler[0,2.Pi ,7,1],f[__]]
(*
{{{{{{{f[{0,1}]}}}},{{{{f[{0.897598,1.}]}}}},
{{{{f[{1.7952,1.70177}]}}}},{{{{f[{2.69279,3.19098}]}}}},
{{{{f[{3.59039,4.43371}]}}}},{{{{f[{4.48799,2.70699}]}}}},
{{{{f[{5.38559,0.338121}]}}}}}}}
*)

It's hard to see, but I inserted a . after the 2 ;-).

Now you can check yourself

ListPlot[Euler[0, 2. \[Pi], 30, 1]] // Timing

Mathematica graphics

share|improve this answer
    
Thank you very much for such a detailed response. I now see exactly what was going on. –  David McHarg Mar 25 '13 at 16:22

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.