NDSolve with Euler method

Had the Euler method not been built-in, one could still use NDSolve[]'s method plug-in framework, which enables NDSolve[] to "know" how to use Euler's method.

Here's how to "teach" NDSolve[] the Euler method:

Euler[]["Step"[rhs_, t_, h_, y_, yp_]] := {h, h yp};
Euler[___]["DifferenceOrder"] := 1;
Euler[___]["StepMode"] := Fixed;

Plugging in the "new" method into NDSolve[] is a snap:

xa = x /. First @ NDSolve[{x'[t] == 0.5*x[t] - 0.04*(x[t])^2, x[0] == 1}, 
                          x, {t, 0, 10}, Method -> Euler, StartingStepSize -> 1];

Getting the corresponding table is easily done, thanks to the special methods for accessing the internals of an InterpolatingFunction[]:

pts = Transpose[Append[xa["Coordinates"], xa["ValuesOnGrid"]]]
   {{0., 1.}, {1., 1.46}, {2., 2.10474}, {3., 2.97991}, {4., 4.11467}, {5., 5.49478},
    {6., 7.03447}, {7., 8.57235}, {8., 9.91912}, {9., 10.9431}, {10., 11.6246}}

Showing the InterpolatingFunction[] and the points together in one plot is also easily done:

Plot[xa[t], {t, 0, 10}, Epilog -> {AbsolutePointSize[4], Red, Point[pts]}, Frame -> True]

solution of differential equation via Euler

If you wanted the derivatives as well in your table, it's easy to modify pts:

phs = Append[#, xa'[#[[1]]]] & /@ pts;

You can use these for the "phase plot" of your differential equation:

ParametricPlot[{xa[t], xa'[t]}, {t, 0, 10}, AspectRatio -> 1/GoldenRatio, 
               Epilog -> {AbsolutePointSize[4], Red, Point[Rest /@ phs]}, 
               Frame -> True]

phase plot

The phase plot looks ugly here, since the step size was not very small and the interval of integration was not sampled well. Due to this, the InterpolatingFunction[] is unable to make a smooth-looking derivative. A smaller value of StartingStepSize (say, $1/20$) would have resulted in something that looks smoother:

smoother phase plot

but of course at the expense of more evaluations of the right-hand side.


NDSolve has a slew of options that allow you to control the method. You can find the standard reference here.

There, we learn how to access Euler's method using NDSolve:

Clear[x];
x = x /. First[
  NDSolve[{x'[t] == 0.5*x[t] - 0.04*(x[t])^2, x[0] == 1}, x, {t, 0, 10},
  StartingStepSize -> 1, Method -> {"FixedStep", Method -> "ExplicitEuler"}]
];
grid = Table[{t, x[t]}, {t, 0, 10, 1}]
ListLinePlot[grid]

plot of differential equation solution

It is also quite easy to program this from scratch, if you prefer.


As Mark wrote "It is also quite easy to program this from scratch, if you prefer."

For a solution without NDSolve try..

MyEuler[start_, end_, initialvalue_, nrOfsteps_] :=Module[
       {a = start, b = end, j, m = nrOfsteps},
       h = (b - a)/m; (* fixed step-size *)
       T = Table[ a + (j - 1) h, {j, 1, m + 1}]; 
       Y = Table[ initialvalue, {j, 1, m + 1}];  
       For[ j = 1, j <= m, j++,
            Y[[j + 1]] = Y[[j]] + h  f[T[[j]], Y[[j]]];
           ]; 
       Transpose@{T, Y}] 

Testing

f[t_, x_] = 0.5*x - 0.04*x^2;(* rhs of your ODE *) 
pts = MyEuler[0.0, 10.0, 1.0, 10]; 
ListLinePlot[pts, Mesh -> All, MeshStyle -> Red, Frame -> True]

enter image description here