7
$\begingroup$

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

$\endgroup$
6
  • 1
    $\begingroup$ 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]]; $\endgroup$
    – gpap
    Mar 25, 2013 at 15:30
  • $\begingroup$ 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 ? $\endgroup$ Mar 25, 2013 at 15:33
  • $\begingroup$ 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. $\endgroup$
    – gpap
    Mar 25, 2013 at 15:43
  • $\begingroup$ And NO, it isn't better to use For[]/Do[] at all. You're in the right path $\endgroup$ Mar 25, 2013 at 15:53
  • $\begingroup$ Thanks for the info, much appreciated. NestList is actually a very neat solution to this problem. $\endgroup$ Mar 25, 2013 at 16:21

1 Answer 1

6
$\begingroup$

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

$\endgroup$
1
  • $\begingroup$ Thank you very much for such a detailed response. I now see exactly what was going on. $\endgroup$ Mar 25, 2013 at 16:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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