Consider the IVP

^{(1)}(

*t*) = -21 y(t) +

*e*

^{-t}

y(0) = 0

The solution to this IVP is the function

A plot of this solution is shown in Figure 1.

Figure 1. The solution to the given initial value problem.

If we use Euler's method to approximate a solution to this IVP using a step
size of *h* = 0.1 for a few steps, we get the following values:

*y*

_{1}= 0.1

y(0.2) ≈

*y*

_{2}= -0.0195162582

y(0.3) ≈

*y*

_{3}= 0.1033409593

y(0.4) ≈

*y*

_{4}= -0.0395932331

We immediately note a mistake with *y*_{2} which is negative, and rather
than converging on the solution, the subsequent approximations appear to diverge. Figure 2 shows the solution
and the first twenty iterations using Euler's method.

Figure 2. The first 20 approximations using Euler's method.

To demonstrate this more clearly, we can draw in the slopes followed by Euler's method at each step (drawing in the slopes), as is demonstrated in Figure 3.

Figure 3. Euler's method approximations together with slopes.

To understand this better, Figure 4 focuses on the solution and the directional field in the close vicinity of the solution.

Figure 4. The direction field plot for the given differential equation with the solution for the initial value y(0) = 0.

We note that even close to the solution, the slope of the direction field are very large and
positive below the solution and negative above the solution. Such a differential
equation is termed *stiff*. One method of solving such problems is to use a smaller
value of *h*, however this may be unnecessarily expensive.

A better technique (rather than making *h* smaller) we will look at is a modification of Euler's method
which is called *backward-Euler's method*.

**Note:** this problem was specifically tailored to fail for Euler's method with a step size
of *h* = 0.1. While better techniques like 4th-order Runge Kutta will solve this problem
satisfactorily, the problem can be made more stiff by substituting -21 with a
larger negative number.

# Backward-Euler's Method

For Euler's method, we approximate y(*t*_{1}) with

*y*

_{1}=

*y*

_{0}+

*h*f(

*t*

_{0},

*y*

_{0})

where *h* = *t*_{1} − *t*_{0}.

Instead, evaluate the differential equation at *t* = *t*_{1}:

Replace the left-hand side with the backward divided-difference formula:

Next, replace y(*t*_{i}) ≈ *y*_{i}:

This gives us an *implicit* equation for *y*_{1}:

We can now solve for *y*_{1} using any of our standard root-finding
techniques. For example, substitute *y*_{1} = *υ*
and rewrite the equation as

and let *υ*_{0} = *y*_{0} and
*υ*_{1} = *y*_{0} + *h* f(*t*_{0}, *y*_{0}),
that is, y(*t*_{0}) and Euler's approximation to y(*t*_{1}), respectively. We
can now use the secant method to find solve for *υ* and we let *y*_{1} be
this approximation.

We can then repeat this process arbitrarily often to approximate *y*_{2}, *y*_{3},
etc.

To demonstrate this in our example, we will work out two steps:

### Finding *y*_{1}

Given *t*_{0} = 0, *t*_{1} = 0.1, and
*y*_{0} = 0, if we substitute these into:

we get the (rather simple) equation

In this case, we can trivially solve for *υ* = 0.02918830381, and therefore
y_{1} = 0.02918830381.

### Finding *y*_{2}

To find *y*_{2}, we use the equation

and substitute *t*_{2} = 0.2 and *y*_{1} = 0.02918830381 to get:

Again, this may be trivially solved for *υ* = 0.03582625132, and therefore
we set *y*_{2} = 0.03582625132.

### Further Solutions

If we continue this process, we get the approximations shown in Figure 5.

Figure 5. Approximations using the backward-Euler's method.

# Higher-order Backward-Euler Methods

We have used the first-order backward divided-difference formula. There is no reason we could not have substituted the second-order backward divided-difference such as:

which, given two initial values *y*_{0} and *y*_{1}, we could approximate
*y*_{2}, and then using *y*_{1} and *y*_{2}, we could approximate
*y*_{3}, and so on.

# Why does this work?

The critical question which you are probably asking is: why does this technique work when Euler's method does not? The following is a heuristic: consider the plot in Figure 6 which is a zoom of the solution shown in Figure 4.

Figure 6. A zoom of Figure 4.

Around the solution, all the direction arrows point in approximately the same direction of the solution, while further away, they tend to be more vertical.

If we apply Euler's method, we are following the slope, for example, in Figure 7 we follow the slope at the point (0.1, y(0.1)). This very quickly leads off into regions where the slope is large and negative, and therefore the next approximation with Euler's method will not point anywhere near the actual solution.

Figure 7. The tangent line at the point (0.1, y(0.1)).

If we apply Euler's method at (0.1, y(0.1)) for various values of *h*, Figure 8 shows
the next iteration. Very quickly, the second iteration begins to diverge very quickly from
the actual solution.

Figure 8. Two steps of Euler's method for values of *h* between 0.001 and 0.15.

Alternatively, the backward-Euler's method works as follows: going forward a distance *h*,
what is the eligible point which could have had (0.1, y(0.1)) as a point on the solution? Figure
9 chooses a number of points around 0.2 and the lines passing through those points.

Figure 9. 21 various (0.2, *y*_{1}) values and the backward paths.

Notice how only one of the points could have had (0.1, y(0.1)) as its predecessor. That
point is approximately (0.2, 0.39) (the actual value being (0.2, 0.03902971461)). Figure 10
shows both Euler's method and the backward-Euler's method with *h* = 0.1.

Figure 10. Euler's and the backward-Euler's method at (0.1, y(0.1)) with *h* = 0.1.

The second-order backward-Euler's method uses a second-order backward divided difference
approximation of the derivative. Given (0, 0) and (1, y(1)), Figure 11 shows the approximation
of y(0.2) using this technique. In this case, the 2nd-order technique is not as good as the
1st-order technique demonstrated in Figure 10, however, for smaller values of *h* it
becomes more accurate.

Figure 11. The second-order backward-Euler's method using (0, 0) and (0.1, y(0.1)).

Copyright ©2005 by Douglas Wilhelm Harder. All rights reserved.