- Submitted by
- SciPy Central on 22 July 2011
- Download
- File
- Update history
- Revision 2 of 3: previous next
- Updated by
- SciPy Central on 22 July 2011
- Tags

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:

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.

Read the SciPy documentation for `ode`

.