Taylor Series and Euler methods

Suppose we wish to solve the initial value problem

       dy
       -- = f(x,y(x)),   with y(x0) = y0
       dx

There are two ways to derive Euler's method. First we apply the forward difference formula to dy/dx:

    y(x+h) - y(x)
    ------------- = f(x,y(x))
          h
which gives rise to
    y(x+h) = y(x) + h*f(x,y(x))
At x = x0, denoting y(x0+h) = y1, we obtain
    y1 = y0 + h*f(x0,y0),
for any index n we can write
    y    = y  + h*f(x , y )
     n+1    n        n   n
The second way to derive Euler's method is via Taylor series:
    y(x0+h) = y(x0) + h*y'(x0) + h^2/2*y"(x0) + O(h^3)
If we truncate after the term in h, and replace y'(x0) by f(x0,y0) -- we can do this because of the equation dy/dx = f(x,y(x)) -- we also obtain the formula for Euler's method.

For example, consider the very simple initial value problem:

   dy
   -- = y,   with y(0) = 1
   dx
Then the solution is y(x) = e^x. Euler's method is then very simple:
    y    = y  + h*y  = (1+h)*y
     n+1    n      n          n
In the following table we list values for h = 0.1:
+---+-----+--------------+--------------+---------+
| n |  x  | y = 1.1*y    | y(x ) = e^x  | error   |
|   |   n |  n       n-1 |    n       n |         |
+---+-----+--------------+--------------+---------+
| 0 | 0.0 |   1          |  1.000 000   | 0       |
| 1 | 0.1 |   1.1        |  1.105 171   | 5.2 E-3 |
| 2 | 0.2 |   1.21       |  1.221 403   | 1.1 E-2 |
| 3 | 0.3 |   1.331      |  1.349 859   | 1.9 E-2 |
| 4 | 0.4 |   1.464 1    |  1.491 825   | 2.8 E-2 |
+---+-----+--------------+--------------+---------+
We see the accuracy deteriorating.

To obtain a more accurate approximations, we go one step further in the Taylor series and apply the forward difference formula to y"(x0):

             y'(x0+h) - y'(x0)
   y"(x0) = -------------------
                    h
so that we obtain:
   y(x0+h) = y(x0) + h*y'(x0) + h/2*(y'(x0+h) - y'(x0)) + O(h^3)
           = y(x0) + h/2(y'(x0) + y'(x0+h)) + O(h^3)
For a general index n, we have the formula for the modified Euler method:
                 h  (                          )
    y    = y  + ---*( f(x , y ) + f(x   ,y   ) )
     n+1    n    2  (    n   n       n+1  n+1  )

Observe that the unknown value y_(n+1) for y(x0+(n+1)*h) appears on both sides of the equation. We apply the modified Euler method as a predictor-corrector method in two stages:

  1. Apply Euler's method to predict y_(n+1):
        yy = y  + h*f(x ,y )
              n        n  n
    

  2. Use the predicted value to correct the value:
                     h  (                        )
        y    = y  + ---*( f(x , y ) + f(x   ,yy) )
         n+1    n    2  (    n   n       n+1     )
    
    
The equation for the modified Euler method in general a nonlinear equation in y_(n+1) and the corrector method we applied in step two above can be interpreted as one step of a fixed-point iteration method. For increased accuracy we can apply more than one corrector step. To accelerate the convergence, Newton's method is recommended.

As numerical example, We consider again the problem dy/dx = y, y(0) = 1 and use h = 0.1 as same step size. In the table below we list the values obtained from the predictor (yy = 1.1*y_(n-1)) and corrector (1 step):

+---+-----+-----------+-----------+--------------+---------+
| n |  x  |    yy     |     y     | y(x ) = e^x  |  error  |
|   |   n |           |      n    |    n       n |         |
+---+-----+-----------+-----------+--------------+---------+
| 0 | 0.0 |    --     | 1         |  1.000 000   | 0       |
| 1 | 0.1 | 1.1       | 1.105     |  1.105 171   | 1.7 E-4 |
| 2 | 0.2 | 1.215 5   | 1.221 025 |  1.221 403   | 3.8 E-4 |
| 3 | 0.3 | 1.343 128 | 1.349 233 |  1.349 859   | 6.3 E-4 |
| 4 | 0.4 | 1.484 156 | 1.490 902 |  1.491 825   | 9.2 E-4 |
+---+-----+-----------+-----------+--------------+---------+
If we compare the values for y_n and the error with the values of the previous table, then we see that we have one more decimal place correct with the modified Euler method.