Solve the world's simplest differential equation

Mathematical problem

This exercise addresses the differential equation 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 constant 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

Derive the $ \theta $-rule scheme for solving \eqref{ode} numerically with time step $ \Delta t $: $$ u^{n+1} = \frac{1 - (1-\theta) a\Delta t}{1 + \theta a\Delta t}u^n, $$ Here, $ n=0,1,\ldots,N-1 $.

Hint.\n Set up the Forward Euler, Backward Euler, and the Crank-Nicolson (or Midpoint) schemes first. Then generalize to the $ \theta $-rule above.

Implementation

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

from numpy import linspace, zeros

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

Fix the values of where $ I $, $ a $, and $ T $. Then vary $ \Delta t $ for $ \theta=0,1/2,1 $. Illustrate that if $ \Delta t $ is not sufficiently small, $ \theta=0 $ and $ \theta=1/2 $ can give non-physical solutions (more precisely, oscillating solutions).

Perform experiments and determine empirically the convergence rate for $ \theta=0,1/2,1 $.