Numerical investigations

Experiments with Schemes for Exponential Decay

Hans Petter Langtangen Simula Research Laboratory University of Oslo

August 20, 2012

Summary. This report investigates the accuracy of three finite difference schemes for the ordinary differential equation $ u'=-au $ with the aid of numerical experiments. Numerical artifacts are in particular demonstrated.

Mathematical problem

We address the initial-value problem $$ \begin{align} u'(t) &= -au(t), \quad t \in (0,T], \label{ode}\\ u(0) &= I, \label{initial:value} \end{align} $$ where $ a $, $ I $, and $ T $ are prescribed parameters, and $ u(t) $ is the unknown function to be estimated. This mathematical model is relevant for physical phenomena featuring exponential decay in time, e.g., vertical pressure variation in the atmosphere, cooling of an object, and radioactive decay.

Numerical solution method

We introduce a mesh in time with points $ 0= t_0< t_1 \cdots < t_{N_t}=T $. For simplicity, we assume constant spacing $ \Delta t $ between the mesh points: $ \Delta t = t_{n}-t_{n-1} $, $ n=1,\ldots,N_t $. Let $ u^n $ be the numerical approximation to the exact solution at $ t_n $. The $ \theta $-rule [1] is used to solve \eqref{ode} numerically: $$ u^{n+1} = \frac{1 - (1-\theta) a\Delta t}{1 + \theta a\Delta t}u^n, $$ for $ n=0,1,\ldots,N_t-1 $. This scheme corresponds to

The Forward Euler scheme when $ \theta=0 $
The Backward Euler scheme when $ \theta=1 $
The Crank-Nicolson scheme when $ \theta=1/2 $

Implementation

The numerical method is implemented in a Python function [2]:

def solver(I, a, T, dt, theta):
    """Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""
    Nt = int(round(T/float(dt)))  # no of intervals
    u = zeros(Nt+1)
    t = linspace(0, T, Nt+1)

    u[0] = I
    for n in range(0, Nt):
        u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*u[n]
    return u, t

Numerical experiments

We define a set of numerical experiments where $ I $, $ a $, and $ T $ are fixed, while $ \Delta t $ and $ \theta $ are varied. In particular, $ I=1 $, $ a=2 $, $ \Delta t= 1.25, 0.75, 0.5, 0.1 $.

The Backward Euler method

The Crank-Nicolson method

The Forward Euler method

Error versus time discretization

How $ E $ varies with $ \Delta t $ for $ \theta = 0, 0.5, 1 $ is shown below.

Observe: The data points for the three largest $ \Delta t $ values in the Forward Euler method are not relevant as the solution behaves non-physically.

Summary

$ \theta =1 $: $ E\sim \Delta t $ (first-order convergence).
$ \theta =0.5 $: $ E\sim \Delta t^2 $ (second-order convergence).
$ \theta =1 $ is always stable and gives qualitatively corrects results.
$ \theta =0.5 $ never blows up, but may give oscillating solutions if $ \Delta t $ is not sufficiently small.
$ \theta =0 $ suffers from fast-growing solution if $ \Delta t $ is not small enough, but even below this limit one can have oscillating solutions (unless $ \Delta t $ is sufficiently small).

Bibliography

A. Iserles. A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, 2009.
H. P. Langtangen. A Primer on Scientific Programming With Python, Springer, 2012.