Integrating an initial value problem (single ODE)

Submitted by: SciPy Central, 22 July 2011
Update history: Revision 1 of 3: next
Updated by: SciPy Central, 22 July 2011
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
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]')
Software license: Creative Commons Zero. No rights reserved.
Users have permission to do anything with the code and other material on this page. (More details)
More information:

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.

Read the SciPy documentation for ode.