Skip to content

Experiments with Schemes for Exponential Decay

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.

Table of contents

Mathematical problem
Numerical solution method
Implementation
Numerical experiments
      The Backward Euler method
      The Crank-Nicolson method
      The Forward Euler method
      Error vs \Delta t

Mathematical problem

We address the initial-value problem

u'(t) = -au(t), \quad t \in (0,T],

u(0)  = I,

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.

Numerical solution method

We introduce a mesh in time with points 0= t_0< t_1 \cdots < t_N=T.
For simplicity, we assume constant spacing \Delta t between the
mesh points: \Delta t = t_{n}-t_{n-1}, n=1,\ldots,N. Let
u^n be the numerical approximation to the exact solution at t_n.

The \theta-rule is used to solve (REF to equation ode not supported) numerically:

u^{n+1} = \frac{1 - (1-\theta) a\Delta t}{1 + \theta a\Delta t}u^n,

for n=0,1,\ldots,N-1. This scheme corresponds to

Implementation

The numerical method is implemented in a Python function
solver (found in the decay_mod module):

def solver(I, a, T, dt, theta):
    """Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""
    dt = float(dt)           # avoid integer division
    N = int(round(T/dt))     # no of time intervals
    T = N*dt                 # adjust T to fit time step dt
    u = zeros(N+1)           # array of u[n] values
    t = linspace(0, T, N+1)  # time mesh

    u[0] = I                 # assign initial condition
    for n in range(0, N):    # n=0,1,...,N-1
        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 vs \Delta t

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


Figure 1: Error versus time step.

Advertisements