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 u'(t)=-au(t), for t in (0,T], and with initial condition u(0)=I. Here, I, a, 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.

Numerical solution method

We use a mesh in time with equally spaced points. The theta-rule is used to discretize the differential equation. 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:

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

    u[0] = I
    for n in range(0, N):
        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 the time step and theta are varied. In particular, I=1, a=2, dt=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 the time step, for the three discretization methods, is shown below.