title: Solve ODEs with discontinuous input/forcing data
body: <p>I'm trying to solve a system of coupled, first-order ODEs in Python. I'm new to this, but the <a href="http://wiki.scipy.org/Cookbook/Zombie_Apocalypse_ODEINT" rel="nofollow noreferrer">Zombie Apocalypse example</a> from SciPy.org has been a great help so far.</p>
  
  <p>An important difference in my case is that the input data used to "drive" my system of ODEs changes <strong>abruptly</strong> 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.</p>
  
  <p><strong>Simplified example</strong></p>
  
  <p>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 <strong>R</strong> and the output from the hole is <strong>Q</strong>. <strong>Q</strong> is assumed to be proportional to the volume of water in the bucket, <strong>V</strong>. The constant of proportionality is usually written as <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=%5C%5B%5Cfrac%7B1%7D%7BT%7D%5C%5D%0A" alt="formula1">, where <strong>T</strong> is the "residence time" of the store. This gives a simple ODE of the form</p>
  
  <p><img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=%5C%5B%5Cfrac%7BdQ%7D%7Bdt%7D%20%3D%20%5Cfrac%7B(R%20-%20Q)%7D%7BT%7D%5C%5D%0A%0A" alt="fromula2"></p>
  
  <p><a href="https://i.stack.imgur.com/Tv33Hm.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Tv33Hm.png" alt="Simple linear reservoir"></a></p>
  
  <p>In reality, <strong>R</strong> is an observed time series of daily rainfall totals. <em>Within</em> each day, the rainfall rate is assumed to be constant, but <em>between</em> days the rate changes abruptly (i.e. <strong>R</strong> is a <strong>discontinuous</strong> function of time). I'm trying to understand the implications of this for solving my ODEs.</p>
  
  <p><strong>Strategy 1</strong></p>
  
  <p>The most obvious strategy (to me at least) is to apply SciPy's <code>odeint</code> function separately within each rainfall time step. This means I can treat <strong>R</strong> as a constant. Something like this:</p>
  
  <pre><code>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
  </code></pre>
  
  <p><strong>Strategy 2</strong></p>
  
  <p>The second approach involves simply feeding everything to <code>odeint</code> and letting it deal with the discontinuities. Using the same parameters and <strong>R</strong> values as above:</p>
  
  <pre><code>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]
  </code></pre>
  
  <p>Both of these approaches give identical answers, which look like this:</p>
  
  <p><a href="https://i.stack.imgur.com/tyYQM.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/tyYQM.png" alt="Result"></a></p>
  
  <p>However the second strategy, although more compact in terms of code, it <strong>much</strong> slower than the first. I guess this is something to do with the discontinuities in <strong>R</strong> causing problems for <code>odeint</code>?</p>
  
  <p><strong>My questions</strong></p>
  
  <ol>
  <li>Is <strong>strategy 1</strong> the best approach here, or is there a <strong>better way</strong>?</li>
  <li>Is <strong>strategy 2</strong> a bad idea and why is it so slow?</li>
  </ol>
  
  <p>Thank you!</p>