Intro

Many dynamical systems involve processes with time delays. Such a time-delay system can be described by rate equations containing functions of previous time points. This type of differential equations is called delay differential equations (DDEs).

We consider again the exponential growth model, but now the rate equation depends on \(N\) an earlier time point \(t - \tau\):

\[\frac{dN}{dt} = r N(t-\tau) \tag {1}\] Here, \(N\) is the population size, \(t\) is time, \(\tau\) is a constant time delay. We consider \(\tau = 1\) here. Unlike ODEs where an initial condition (i.e., \(N(t_0)\) here) needs to be specified, the DDE in Equation (1) requires \(N(t)\) for \(t \in [t_0-\tau, t_0]\). If a time-delay system has an steady state \(N_s\), \(N_s\) is also the steady state of the corresponding ODE with time delays. (This is because a steady state is time independent) Thus, same method for ODEs can be used to identify steady state solutions. However, the stability of the system may change because of time delays.

Euler method for DDE

We consider the DDE in Equation (1) and \(N(t - \tau) = N_0\) for \(t \leq \tau\). (Be careful, this seems intuitive but may cause issues as the derivatives of the solutions are discontinuous. Because of the discontinuity of the time trajectories, DDE often is ill defined. One needs to be cautious when interpreting the simulation results of DDE.) The most common way to simulate a DDE is to generate time steps and \(N\) solutions for the whole interval \(t \in [-\tau, t_{max}]\). \(N(t-\tau)\) can be retrieve anytime during DDE integration. This is not the most efficient numerical method, due to high storage requirement for saving the whole trajectory (Q: think how to make it more efficient, see the section Method of steps).

The function dN_delays specifies the derivative function for a system with a constant time delay. Note that we supply \(N\) for both current time N and previous time N_old.

# derivative function for delay exponential growth model
# t: time
# N: population size of the current time N(t)
# N_old: population size of the previous time N(t-tau). 
#        We will use the DDE integrator to specify the constant time delay.
# r: parameter, growth rate
dN_delays <- function(t, N, N_old, r) {
   return(r*N_old)
}

The DDE integrator can be modified from an ODE integrator to incorporate data from previous time points. Below, we present an implementation using the Euler method.

\[N_{n+1} = N_n + h*f(t_n, N_n, N_{n-d}) \tag {2}\]

Here, \(n - d\) represents the index of the time step for delays.

# Euler method for DDEs with a constant time delay
DDE_Euler <- function(derivs, t0, X0, t.total, dt, tau, ...){
  # derivs: the function of the derivatives 
  # t0: initial time
  # X0: initial condition
  # t.total: total simulation time
  # dt: time step size 
  # tau: a constant time delay 
  # (this illustrate will not work for cases using multiple time delays)
  t_all = seq(t0-tau, t.total, by=dt)   # save data for all previous time points too
  n_all = length(t_all)
  X_all = numeric(n_all)
  
  n_delay = as.integer(tau/dt) # number of delay time steps
  X_all[1:(n_delay+1)] = X0   # a constant initial condition between t in [-tau, t.total]
  
  for (i in (n_delay+1):(n_all-1)) {
    X_all[i+1] = X_all[i] + dt * derivs(t_all[i],X_all[i],X_all[i-n_delay],...)
  }
  return(cbind(t_all, X_all))   # the output is a matrix of t & N(t) for all time steps
}

We set \(r = -0.3\). It’s a typical exponential decay.

r1 = - 0.3
results_1 = DDE_Euler(derivs = dN_delays, t0 = 0, X0 = 1, t.total = 40, dt = 0.01, tau = 1, r = r1)

plot(results_1[,1], results_1[,2], xlab="t", ylab="N", xlim=c(-1,40), ylim=c(0,1),
     type = "l", col=2)
legend("topright", inset=0.02, legend = "r = -0.3", col=2, lty=1, cex=0.8)

When \(r\) is slightly larger than a critical point \(\pi/2\), the system exhibits a damping oscillation. Note that \(N\) would be negative in this model, which is unrealistic to represent the population size.

r2 = - 1.4
results_2 = DDE_Euler(derivs = dN_delays, t0 = 0, X0 = 1, t.total = 40, dt = 0.01, tau = 1, r = r2)

plot(results_2[,1], results_2[,2], xlab="t", ylab="N", xlim=c(-1,40), ylim=c(-1,1),
     type = "l", col=2)
legend("topright", inset=0.02, legend = "r = -1.4", col=2, lty=1, cex=0.8)

When \(r\) is slightly smaller than a critical point \(\pi/2\), the system exhibits an increasing oscillation.

r3 = - 1.7
results_3 = DDE_Euler(derivs = dN_delays, t0 = 0, X0 = 1, t.total = 40, dt = 0.01, tau = 1, r = r3)

plot(results_3[,1], results_3[,2], xlab="t", ylab="N", xlim=c(-1,40), ylim=c(-10,11),
     type = "l", col=2)
legend("topleft", inset=0.02, legend = "r = -1.7", col=2, lty=1, cex=0.8)

A 2nd order DDE integrator (Heun’s method)

For a better performance, we convert a second order ODE integrator, the Heun’s method, for modeling DDE. The Heun’s method is very similar to 2nd order Runge-Kutta but uses only integer time steps. So it’s much convenient to convert the ODE integrator to the DDE version using Heun’s method.

\[\begin{equation} \begin{aligned} k_1 &= h*f(t_n, N_n, N_{n-d}) \\ N_{n+1} &= N_n + k_1 \\ k_2 &= h*f(t_{n+1}, N_{n+1}, N_{n+1-d}) \\ N_{n+1} &= N_n + (k_1+k_2)/2 \end{aligned}\tag{3} \end{equation}\]

# Heun's method for DDEs with a constant time delay
DDE_Heun <- function(derivs, t0, X0, t.total, dt, tau, ...){
  # derivs: the function of the derivatives 
  # t0: initial time
  # X0: initial condition
  # t.total: total simulation time
  # dt: time step size 
  # tau: a constant time delay 
  # (this illustrate will not work for cases using multiple time delays)
  t_all = seq(t0-tau, t.total, by=dt)   # save data for all previous time points too
  n_all = length(t_all)
  X_all = numeric(n_all)
  
  n_delay = as.integer(tau/dt) # number of delay time steps
  X_all[1:(n_delay+1)] = X0   # a constant initial condition between t in [-tau, 0]
  
  for (i in (n_delay+1):(n_all-1)) {
    k1 = dt * derivs(t_all[i],X_all[i],X_all[i-n_delay],...)
    X_all[i+1] = X_all[i] + k1   # X_all_i+1 is first assigned here, so that the function works for tau = 0
    k2 = dt * derivs(t_all[i+1],X_all[i+1], X_all[i+1-n_delay],...)
    X_all[i+1] = X_all[i] + (k1+k2)/2
  }
  return(cbind(t_all, X_all))   # the output is a matrix of t & N(t) for all time steps
}

results_3_2 = DDE_Heun(derivs = dN_delays, t0 = 0, X0 = 1, t.total = 40, dt = 0.01, tau = 1, r = r3)

plot(results_3[,1], results_3[,2], xlab="t", ylab="N", xlim=c(-1,40), ylim=c(-11,12),
     type = "l", col=2)
lines(results_3_2[,1], results_3_2[,2], col = 3)
legend("topleft", inset=0.02, legend = c("Euler", "Heun"), col=c(2,3),lty=1, cex=0.8)

The time trajectories from the two methods become more different in later time points, with the Euler trajectory slightly overshooting.