Solve ODEs with discontinuous input/forcing data

Question

I'm trying to solve a system of coupled, first-order ODEs in Python. I'm new to this, but the Zombie Apocalypse example from SciPy.org has been a great help so far.

An important difference in my case is that the input data used to "drive" my system of ODEs changes abruptly at various time points and I'm not sure how best to deal with this. The code below is the simplest example I can think of to illustrate my problem. I appreciate this example has a straightforward analytical solution, but my actual system of ODEs is more complicated, which is why I'm trying to understand the basics of numerical methods.

Simplified example

Consider a bucket with a hole in the bottom (this kind of "linear reservoir" is the basic building block of many hydrological models). The input flow rate to the bucket is R and the output from the hole is Q. Q is assumed to be proportional to the volume of water in the bucket, V. The constant of proportionality is usually written as , where T is the "residence time" of the store. This gives a simple ODE of the form

In reality, R is an observed time series of daily rainfall totals. Within each day, the rainfall rate is assumed to be constant, but between days the rate changes abruptly (i.e. R is a discontinuous function of time). I'm trying to understand the implications of this for solving my ODEs.

Strategy 1

The most obvious strategy (to me at least) is to apply SciPy's odeint function separately within each rainfall time step. This means I can treat R as a constant. Something like this:

import numpy as np, pandas as pd, matplotlib.pyplot as plt, seaborn as sn
from scipy.integrate import odeint

np.random.seed(seed=17)

def f(y, t, R_t):
    """ Function to integrate.
    """
    # Unpack parameters
    Q_t = y[0]

    # ODE to solve
    dQ_dt = (R_t - Q_t)/T

    return dQ_dt

# #############################################################################
# User input   
T = 10      # Time constant (days)
Q0 = 0.     # Initial condition for outflow rate (mm/day)
days = 300  # Number of days to simulate
# #############################################################################

# Create a fake daily time series for R
# Generale random values from uniform dist
df = pd.DataFrame({'R':np.random.uniform(low=0, high=5, size=days+20)},
                  index=range(days+20))

# Smooth with a moving window to make more sensible                  
df['R'] = pd.rolling_mean(df['R'], window=20)

# Chop off the NoData at the start due to moving window
df = df[20:].reset_index(drop=True)

# List to store results
Q_vals = []

# Vector of initial conditions
y0 = [Q0, ]

# Loop over each day in the R dataset
for step in range(days):
    # We want to find the value of Q at the end of this time step
    t = [0, 1]

    # Get R for this step
    R_t = float(df.ix[step])   

    # Solve the ODEs
    soln = odeint(f, y0, t, args=(R_t,))

    # Extract flow at end of step from soln
    Q = float(soln[1])

    # Append result
    Q_vals.append(Q)

    # Update initial condition for next step
    y0 = [Q, ]

# Add results to df
df['Q'] = Q_vals

Strategy 2

The second approach involves simply feeding everything to odeint and letting it deal with the discontinuities. Using the same parameters and R values as above:

def f(y, t):
    """ Function used integrate.
    """
    # Unpack incremental values for S and D
    Q_t = y[0]

    # Get the value for R at this t
    idx = df.index.get_loc(t, method='ffill') 
    R_t = float(df.ix[idx])       

    # ODE to solve
    dQ_dt = (R_t - Q_t)/T

    return dQ_dt

# Vector of initial parameter values
y0 = [Q0, ]

# Time grid
t = np.arange(0, days, 1)

# solve the ODEs
soln = odeint(f, y0, t)

# Add result to df
df['Q'] = soln[:, 0]

Both of these approaches give identical answers, which look like this:

However the second strategy, although more compact in terms of code, it much slower than the first. I guess this is something to do with the discontinuities in R causing problems for odeint?

My questions

1. Is strategy 1 the best approach here, or is there a better way?
2. Is strategy 2 a bad idea and why is it so slow?

Thank you!

Answer 1

1.) Yes

2.) Yes

Reason for both: Runge-Kutta solvers expect ODE functions that have an order of differentiability at least as high as the order of the solver. This is needed so that the Taylor expansion which gives the expected error term exists. Which means that even the order 1 Euler method expects a differentiable ODE function. Thus no jumps are allowed, kinks can be tolerated in order 1, but not in higher order solvers.

This is especially true for implementations with automatic step size adaptations. Whenever a point is approached where the differentiation order is not satisfied, the solver sees a stiff system and drives the step-size toward 0, which leads to a slowdown of the solver.

You can combine strategies 1 and 2 if you use a solver with fixed step size and a step size that is a fraction of 1 day. Then the sampling points at the day turns serve as (implicit) restart points with the new constant.