Euler's Method for First Order Differential Equations

The simplest method of numerical integration is Euler's method, which is presented now. Consider the first order differential equation:

wpe3D.gif (1777 bytes)                (eq.1)

Starting at some time, to, he value of y(to+h) can then be approximated by the value of y(to) plus the time step multiplied by the slope of the function, which is the derivative of y(t) (Note: this is simply a first order Taylor expansion)::

wpe4C.gif (1550 bytes)

we will call this approximate value y*(t)

wpe4F.gif (1560 bytes)            (eq. 2)

So, if we can calculate the value of dy/dt at time to (using equation 1) then we can generate an approximation for the value of y at time  to+h using equation 2.  We can then use this new value of y (at to) to find dy/dt (at to) and repeat.  Although this seems circular, if properly used it can generat an approximate solution.  This is refered to as Euler's method.

This process is shown graphically below for a sinlge interation of the process.

wpe6C.gif (3907 bytes)

In this diagram you can see how the time step, h, and the slope of the function, k1, are used to generate y*(to+h) from y(to).

With this simple background the Euler method for first order differential equations can be stated as follows:

1) Starting at time to, choose a value for h, and find initial condition y(to).

2) From y(to) calculate the derivative of y(t) at t=to using equation 1. Call this k1.  The slope is shown as the red line in the drawing above.

wpe57.gif (1247 bytes)

3) From this value find an approximate value for y*(to+h) using equation 2. 

wpe64.gif (1274 bytes)

This is labeled in red in the diagram above.  The exact value is shown in blue.  Note that the difference between approximate and exact solutions is quite large; in practice h is chosen small enough to make this error small.

4) Let to=to+h, and y(to)=y*(to+h).

5) Repeat steps 2 through 4 until the solution is finished.

 

Numerical Example (1st order differential equation - Euler)

Now lets go through a numerical example with to=0 for the equation:

wpe6D.gif (1561 bytes)

which has an exact solution

wpe6E.gif (1231 bytes)

First Time Step

First, let's use h=0.1   To find y*(to+h)=y*(h) (since to=0) we first evaluate the derivative:

wpe71.gif (2004 bytes)

We then find the approximate solution after the time step:

wpe72.gif (1776 bytes)

Our approximate solution is 1.1, the exact solution is 1.04.

Second Time Step

To find the next approximation we take  to=0.1 and repeat.   First we evaluate the derivative (using our approximation for the value of y(to):

wpe75.gif (2156 bytes)

 

We then find the approximate solution after the second time step (again using our approximate value):

wpe78.gif (2030 bytes)

Our approximate solution is 1.08, the exact solution is 1.00.

Subsequent time steps

The procedure is repeated for subsequent time steps.

Results

The table below show the results of the Euler approximations, as well as the exact solution up to 4 seconds.  Here is a link to Matlab code, and C code, that performs the approximation.  The Matlab code printed the results and made the plot (I added results for h=0.01 -- see below)..

t=0.000, y(t)=1.0000, y*(t)=1.0000
t=0.100, y(t)=1.0413, y*(t)=1.1000
t=0.200, y(t)=1.0018, y*(t)=1.0811
t=0.300, y(t)=0.9202, y*(t)=0.9997
t=0.400, y(t)=0.8205, y*(t)=0.8901
t=0.500, y(t)=0.7167, y*(t)=0.7726
.
.
.
t=3.900, y(t)=0.0010, y*(t)=0.0006
t=4.000, y(t)=0.0008, y*(t)=0.0004

wpe7C.gif (4803 bytes)

If h is decreased to 0.01 accuracy is improved, but the complexity (and time taken) of the calculation is greatly increased.  Compare resulsts with h=0.1 at  t=0.1, and t=4.0.  Note that on the graph above, the exact solution, and the approximate solution with h=0.01, are almost identical.

Note: if h is decreased too much, accuracy can actually suffer because there are so many calculations that cummulative error to the roundoffs inherent in computer calculations becomes significant.

t=0.000, y(t)=1.0000, y*(t)=1.0000
t=0.010, y(t)=1.0093, y*(t)=1.0100
t=0.020, y(t)=1.0173, y*(t)=1.0186
t=0.030, y(t)=1.0240, y*(t)=1.0259
t=0.040, y(t)=1.0296, y*(t)=1.0320
t=0.050, y(t)=1.0340, y*(t)=1.0370
t=0.060, y(t)=1.0374, y*(t)=1.0408
t=0.070, y(t)=1.0397, y*(t)=1.0436
t=0.080, y(t)=1.0411, y*(t)=1.0454
t=0.090, y(t)=1.0417, y*(t)=1.0462
t=0.100, y(t)=1.0413, y*(t)=1.0462
t=0.110, y(t)=1.0402, y*(t)=1.0454
.
.
.
t=3.890, y(t)=0.0010, y*(t)=0.0010
t=3.900, y(t)=0.0010, y*(t)=0.0010
t=3.910, y(t)=0.0010, y*(t)=0.0010
t=3.920, y(t)=0.0010, y*(t)=0.0009
t=3.930, y(t)=0.0010, y*(t)=0.0009
t=3.940, y(t)=0.0009, y*(t)=0.0009
t=3.950, y(t)=0.0009, y*(t)=0.0009
t=3.960, y(t)=0.0009, y*(t)=0.0009
t=3.970, y(t)=0.0009, y*(t)=0.0008
t=3.980, y(t)=0.0009, y*(t)=0.0008
t=3.990, y(t)=0.0009, y*(t)=0.0008
t=4.000, y(t)=0.0008, y*(t)=0.0008

AG00051_.gif (1652 bytes)To next page



Comments or Questions?