NDSolve with Euler method

Question

I want to solve this equation with NDSolve[] using the Euler method:

x'[t] == 0.5*x[t]-0.04*(x[t])^2

with initial condition and step size

x[0] == 1, h == 1

and final time t = 10.

I also don't know how to make a table including x'[t] and x[t] and time.

score 20 · Accepted Answer · answered Oct 11 '12 at 23:26

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.

+1. This is neat, I did not know about NDSolve's plug-in framework, very useful. — RunnyKine, May 28 '15 at 04:25

score 17 · Answer 2 · edited Oct 12 '12 at 11:27

17

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.

edited Oct 12 '12 at 11:27

J. M.'s missing motivation

124,525
11
401
574

answered Oct 11 '12 at 20:18

Mark McClure

32,469
3
103
161

6

No need for "FixedStep"; since Mathematica already "knows" that the Euler method is a fixed step method, one can do this instead: NDSolve[{x'[t] == 0.5*x[t] - 0.04*(x[t])^2, x[0] == 1}, x, {t, 0, 10}, StartingStepSize -> 1, Method -> "ExplicitEuler"] – J. M.'s missing motivation Oct 11 '12 at 23:05

PlatoManiac · Answer 3 · 2012-10-11T20:47:25.980

11

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

edited Oct 11 '12 at 20:47

answered Oct 11 '12 at 20:24

PlatoManiac

14,723
2
42
74

2

something like myEuler[t0_, t1_, initial_, steps_] := With[{h = (t1 - t0)/steps}, FoldList[#1 + h f[#2, #1] &, initial, Range[t0 + h, t1, h]]] could do – Dr. belisarius Oct 11 '12 at 20:53
@belisarius Sure! But the OP wanted some Table. I guess by table he meant the array to store the discrete states as commonly done in most textbooks. I felt that mundane For will also better suit a classic textbook Euler .. – PlatoManiac Oct 11 '12 at 21:09
You could also go one step further and show how to exten NDSolve with your method. – Oct 11 '12 at 23:24
@ruebenko I think J.M already has done exactly that while we were sleeping in Europe..right? – PlatoManiac Oct 12 '12 at 07:04
yes, I think he must have posted this a few minutes after the comment... – Oct 12 '12 at 10:19
FWIW, I didn't see ruebenko's comment before I wrote my answer. After writing it, I didn't think to check for new comments, so... – J. M.'s missing motivation Oct 12 '12 at 11:19

NDSolve with Euler method

3 Answers3

Linked