Integrating an initial value problem (single ODE)
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]')
|
Users have permission to do anything with the code and other material on this page. (More details)
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.
Page views: 207 (past 60 days)
Identifier #: 12
Permalink to this revision: http://scpyce.org/12/3/
Permalink to latest revision: http://scpyce.org/12/