AMTH247 Lecture 20

Differential Equations II

Reading: Heath §9.3.2-9.3.4

Contents

• Stability and Accuracy
  • Stability of Differential Equations
  • Stability of Numerical Methods
  • Accuracy
  • Local Truncation Error
  
  
• Backward Euler Method
  • Stability
  • Stiff Differential Equations
  
  
• About this document ...

Stability and Accuracy

Stability of Differential Equations

Consider a solution $y = y(t)$ of an initial value problem

$$\frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0$$ (1)

This solution $y(t)$ is said to be stable if solutions with nearby initial values remain close to $y(t)$. The solution $y(t)$ is unstable if solutions with nearby initial values diverge away from $y(t)$.

Example:

The initial value problem

$$\frac{dy}{dt} = y - 1, \quad y(0) = 1$$

has the solution

$$y(t) = 1 .$$

The general solution of the differential equation is

$$y(t) = 1 + A e^t .$$

The solution of the initial value problem has $A = 0$, but solutions with other initial values will have $A \neq 0$ and these will diverge away from the solution $y(t) = 1$. Here is a graph of some solutions:

-->a = [-0.0001 -0.001 -0.01 -0.1 0 0.1 0.01 0.001 0.0001]; 

-->x = 0:0.01:10;
 
-->for i = 1:9                                             
--> plot2d(x,1+a(i)*exp(x),1,"011",rect=[0 -5 10 5])       
-->end

Figure 1: An unstable solution.

It is extremely difficult to obtain accurate numerical solutions to unstable initial value problems. As soon as even the smallest numerical error occurs, the computed solution will move to a neighbouring solution curve and thus diverge from the true solution.

Stability of Numerical Methods

Instability of numerical solutions can be caused by the numerical method itself even when the solutions to the differential equation are stable. Consider the initial value problem

$$\frac{dy}{dt} = -\lambda y, \quad y(0) = y_0 .$$

The solution of this initial value problem is

$$y(t) = y_0 e^{-\lambda t} .$$

These solutions will be stable when $\lambda > 0$, since then all solutions are exponentially decaying to $y = 0$. Euler's method for the differential equation $dy/dt = f(t, y)$ is

$$y(t_0 + \Delta t) \approx y_0 + f(t_0, y_0) \Delta t .$$

Writing $y_k$ for the computed solution at $t_k$, that is after $k$ steps of size $h$, we can write Euler's method

$$y_{k+1} = y_k + h f(t_k, y_k) .$$

Applying Euler's method to the differential equation

$$\frac{dy}{dt} = -\lambda y$$

gives

$$y_{k+1} = y_k - h \lambda y_k = (1-h \lambda) y_k .$$

This is a difference equation which can be easily be seen to have the solution

$$y_k = (1 - \lambda h)^k y_0 .$$

If $\vert 1 - \lambda h \vert < 1$ then the computed solution will decay to zero. On the other hand, if $\vert 1 - \lambda h \vert > 1$ then computed solution will grow without bounds. For $\lambda > 0$, the differential equation itself is stable. For Euler's method to be stable we require, in addition, that

$$\vert 1 - \lambda h \vert < 1 .$$

That is

$$h < \frac{2}{\lambda} .$$

Example:

We will use as our example from the previous lecture

$$\frac{dy}{dt} = -2 t y, \quad y(0) = 1$$

which has the solution

$$y(t) = e^{-t^2} .$$

We saw above that Euler method for

$$\frac{dy}{dt} = -\lambda y$$

is unstable for

$$h > \frac{2}{\lambda} .$$

By analogy we might expect that Euler's method for

$$\frac{dy}{dt} = -2 t y$$

will be unstable for

$$h > \frac{2}{2t} = \frac{1}{t} .$$

Another way to think of this is that for fixed step size $h$, Euler's method for

$$\frac{dy}{dt} = -2 t y$$

should become unstable at

$$t = \frac{1}{h} .$$

Let us see what happens in practice. We will use the function euler from the previous lecture. The following example shows what happens:

-->function y = func(t, y)
-->  y = - 2*t*y          
-->endfunction  

-->[t,y] = euler(func, 0, 1, 0.1, 251);
 
-->plot2d(t, y)

Zooming in on the end of the graph shows clearly the onset of instability

This shows the instability apparently starting at $t \approx 25$, whereas according to the argument above it should start at $t \approx 1/h = 10$. Closer examination of the output shows that the instability does start at $t \approx 10$. The smallest value of the computed solution

-->[y1,k]=mini(abs(y))
 k  =
 
    102.  
 y1  =
 
    2.576E-57

occurs at k = 102, that is at $t = 10.1$. Plotting the solution near $t = 10$ clearly shows the switch from stability to instability:

-->plot2d(y(90:115))

Accuracy

It turns out that rounding errors are much less important than truncation errors in solving differential equations. There is still the problem that rounding error increases as the step size decreases, however the step sizes used in practice are usually such that truncation error is still the main source of error. There are two ways to look at the truncation error in a numerical method for solving differential equations. The global error is the difference between the exact and computed solutions. The local error is the error in one step of the method. The global error is not simply the sum of the local errors. The way the local errors compound together to form the global error is rather complicated. This is because the local error measures the error in one step assuming the starting point of the step is accurate. However the starting point of the step is itself subject to error compounded from errors at previous steps. There are no general reliable methods for estimating the global error, in general it is assumed that it is about the size the sum of the local errors.

Local Truncation Error

A numerical method for differential equations is said to be of order $p$ if the local truncation error is proportional to $h^{p+1}$. When the global error is the sum of the local errors, we then expect the global error to be proportional to $h^p$ since the number of steps taken to reach any particular value of $t$ is proportional to $1/h$. Euler's method is of order 1 since it is based on a Taylor series approximation and the first neglected term is proportional to $h^2$.

Example:

We will see how the error depends the step size in our example. We will take step sizes $h = 1/2, \dots, 1/2^{10}$ and for each $h$ compute the error in the computed solution at $t = 1$. The exact value is $e^{-1}$.

-->err = zeros(1,10);
 
-->for k = 1:10
-->  n = 2^k;
-->  [t,y] = euler(func, 0, 1, 1/n, n);
-->  err(k) = abs(y(n+1) - exp(-1));

-->err
 err  =
 
 
         column 1 to 5
 
!   0.1321206    0.0422768    0.0178353    0.0082540    0.0039755 !
 
         column  6 to 10
 
!   0.0019514    0.0009668    0.0004812    0.0002401    0.0001199 !
 
-->plot2d(log10(err))

It looks like the error is decreasing by a factor of about two when the step size halves. This is consistent with the Euler's method being of order 1.

Backward Euler Method

Euler's method is

$$y_{k+1} = y_k + h f(t_k, y_k) .$$

The backward Euler method is

$$y_{k+1} = y_k + h f(t_{k+1}, y_{k+1}) .$$

Given $(t_k, y_k)$, we need to solve the equation

$$y_{k+1} = y_k + h f(t_{k+1}, y_{k+1}) .$$

to compute $y_{k+1}$. Because $y_{k+1}$ is determined only implicitly by $y_k$, methods such as this are called implicit methods. Methods like Euler's method where $y_{k+1}$ can be determined directly from $y_k$ are called explicit methods. Each step of an implicit method requires the solution of a nonlinear equation.

Stability

Applying the backward Euler method to the differential equation

$$\frac{dy}{dt} = -\lambda y$$

gives

$$y_{k+1} = y_k - h \lambda y_{k+1} .$$

Rearranging gives

$$y_{k+1} = \frac{y_k}{1 + h \lambda} .$$

So for the backward Euler method to be stable we require

$$\left\vert \frac{1}{1 + h \lambda} \right\vert < 1$$

which, when $\lambda > 0$, holds for all step sizes $h > 0$. Like Euler's method, the backward Euler method is a first order method.

Stiff Differential Equations

The big disadvantage of implicit methods is that we need to solve a nonlinear equation at each step. However implicit methods, like the backward Euler method, have the property of being stable for all step sizes. Explicit methods, like Euler's method, are may require extremely small step sizes to maintain stability. Stiff differential equations may be loosely defined as those for which explicit numerical methods require very small step sizes to maintain stability. Stiff differential equations are common in applications. They typically arise where there is a slowly varying equilibrium solution with rapidly decaying transients. Stiff differential equations are very stable mathematically, but cause stability problems for explicit methods. Implicit methods are the most practical method for solving stiff systems; the disadvantage of having to solve a nonlinear system at each step is outweighed by the advantage of being able to take reasonably sized steps.

About this document ...

AMTH247 Lecture 20

Differential Equations II

This document was generated using the LaTeX2HTML translator Version 2K.1beta (1.61)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -split +0 lecture.tex

The translation was initiated by amth247 on 2003-06-02