Integrating an initial value problem (single ODE)

Submitted by
SciPy Central on 22 July 2011
Download
File
Update history
Revision 3 of 3: previous
Updated by
SciPy Central on 06 August 2011
Tags
import numpy as np from scipy import integrate from matplotlib.pylab import * def tank(t, y): """ Dynamic balance for the chemical reactor (stirred tank) C_A = y[0] = the concentration of A in the tank, mol/L Returns dy/dt = F/V*(C_{A,in} - C_A) - k*C_A^2 """ F = 20.1 # L/min CA_in = 2.5 # mol/L V = 100 # L k = 0.15 # L/(mol.min) # Assign some variables for convenience of notation CA = y[0] # Output from this ODE function must be a COLUMN vector, with n rows n = len(y) dydt = np.zeros((n,1)) dydt[0] = F/V*(CA_in - CA) - k*CA**2 return dydt # The "driver" that will integrate the ODE(s): # ----------- # Start by specifying the integrator: # use ``vode`` with "backward differentiation formula" r = integrate.ode(tank).set_integrator('vode', method='bdf') # Set the time range t_start = 0.0 t_final = 10.0 delta_t = 0.1 # Number of time steps: 1 extra for initial condition num_steps = np.floor((t_final - t_start)/delta_t) + 1 # Set initial condition(s): for integrating variable and time! CA_t_zero = 0.5 r.set_initial_value([CA_t_zero], t_start) # Create vectors to store trajectories t = np.zeros((num_steps, 1)) CA = np.zeros((num_steps, 1)) t[0] = t_start CA[0] = CA_t_zero # Integrate the ODE(s) across each delta_t timestep k = 1 while r.successful() and k < num_steps: r.integrate(r.t + delta_t) # Store the results to plot later t[k] = r.t CA[k] = r.y[0] k += 1 # All done! Plot the trajectories: plot(t, CA) grid('on') xlabel('Time [minutes]') ylabel('Concentration [mol/L]')

We want to integrate a single equation \(\displaystyle \frac{dy(t)}{dt} = f(t, y)\) with a given initial condition \(y(t=0)=y_0\).

There are several integrators available in SciPy. This tutorial uses the VODE integrator. It is a good general-purpose solver for both stiff and non-stiff systems.

The example is a common modelling reaction: a liquid-based stirred tank reactor, with (initially) constant physical properties, a second order chemical reaction where species A is converted to B according to \({\sf A} \rightarrow {\sf B}\), with reaction rate \(r = k C_{\sf A}^2\). One can find the time-dependent mass balance for this system to be:

\[\frac{dC_{\sf A}(t)}{dt} = \frac{F(t)}{V} \bigg( C_{{\sf A},\text{in}} - C_{\sf A} \bigg) - k C_{\sf A}^2\]

where

  • \(C_{{\sf A},\text{in}} = 5.5\) mol/L,
  • we will initially assume constant volume of \(V = 100\) L
  • constant inlet flow of \(F(t) = 20.1\) L/min, and
  • a reaction rate constant of \(k = 0.15 \frac{\text{L}}{\text{mol}.\text{min}}\).

We must specify an initial condition for each differential equation: we will assume \(C_{\sf A}(t=0) = 0.5\) mol/L.

In the code we integrate the ODE from \(t_\text{start} = 0.0\) minutes up to \(t_\text{final} = 10.0\) minutes and plot the time-varying trajectory of concentration in the tank using matplotlib. The plot shows the reactor starts with a concentration of \(C_{\sf A}(t=0) = 0.5\) mol/L and reaches a steady state value of about \(C_{\sf A}(t=\infty) = 1.28\) mol/L.

Take a look at the Python output in printed form:

>>> t

array([[ 0. ],

       [ 0.1],

       [ 0.2],

       [ 0.3],

        ....

       [ 9.7],

       [ 9.8],

       [ 9.9],

       [10. ]])



>>> CA

array([[ 0.5       ],

       [ 0.53581076],

       [ 0.57035065],

       [ 0.60362937],

        ....

       [ 1.27573058],

       [ 1.27592064],

       [ 1.27609987],

       [ 1.27626888]])

The evenly spaced periods of time were by design. Internally, the code calculates the ODE value at many points within each of the delta_t time steps. But it only reports the values at the boundaries of each delta_t, not at these intermediate points.

/media/images/201108/single-ode-plot.png

Read the SciPy documentation for ode.

Please Sign in or Register to leave a comment

  • Creative Commons Zero. No rights reserved.
    Users have permission to do anything with the code and other material on this page. (More details)
  • Trademarks are property of their respective owners. Code and comments are owned by their respective posters. © 2013 All Rights Reserved