The simplest method of numerical integration is Euler's method, which is presented now. Consider the first order differential equation:
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)::
we will call this approximate value y*(t)
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.
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.
3) From this value find an approximate value for y*(to+h) using equation 2.
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:
which has an exact solution
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:
We then find the approximate solution after the time step:
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):
We then find the approximate solution after the second time step (again using our approximate value):
Our approximate solution is 1.08, the exact solution is 1.00.
Subsequent time steps
The procedure is repeated for subsequent time steps.
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
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
To next page
Comments or Questions?