Topic 14.6: Stiff Differential Equations

There are a certain class of differential equations which the four numerical
solvers we have looked at (Euler, Heun, RK4 and RKF45) are numerically
unstable. Unfortunately, this class of problems is difficult to
define, though they appear quite regularly within the field of electrical
engineering. This topic is mostly here to make you aware of such problems
and to help you possibly recognize such situations. We will look at one
solver in particular — a modification of Euler's method — to solve
such problems.

We immediately note a mistake with y2 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(t1) with

y1 = y0 + h f(t0, y0)

where h = t1 − t0.

Instead, evaluate the differential equation at t = t1:

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

Next, replace y(ti) ≈ yi:

This gives us an implicit equation for y1:

We can now solve for y1 using any of our standard root-finding
techniques. For example, substitute y1 = υ
and rewrite the equation as

and let υ0 = y0 and
υ1 = y0 + h f(t0, y0),
that is, y(t0) and Euler's approximation to y(t1), respectively. We
can now use the secant method to find solve for υ and we let y1 be
this approximation.

We can then repeat this process arbitrarily often to approximate y2, y3,
etc.

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

Finding y1

Given t0 = 0, t1 = 0.1, and
y0 = 0, if we substitute these into:

we get the (rather simple) equation

In this case, we can trivially solve for υ = 0.02918830381, and therefore
y1 = 0.02918830381.

Finding y2

To find y2, we use the equation

and substitute t2 = 0.2 and y1 = 0.02918830381 to get:

Again, this may be trivially solved for υ = 0.03582625132, and therefore
we set y2 = 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 y0 and y1, we could approximate
y2, and then using y1 and y2, we could approximate
y3, 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, y1) 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.

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.

If we remove the h2 term, we get the equation which we
are solving for the backward-Euler's method. With further manipulation, it
can be argued that this implies that the error also drops with respect to
h2.

Questions

The following are not stiff differential equations, however, the
techniques may still be applied.

Question 1

Given the IVP

y(1)(t) = 1 - 0.25 y(t) + 0.2 t
y(0) = 1

approximate y(0.5), y(1), and y(1.5) using the backward-Euler's method.

Answer: 1.377777778, 1.76, 2.145454545.

Question 2

Given the same ODE as in Question 1, but with the initial condition
y(1) = 2, approximate y(1.5) and y(2.0).