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$.