Integrating an initial value problem (multiple ODEs)

Submitted by: SciPy Central, 22 July 2011
Update history: Revision 2 of 2: previous 
Updated by: SciPy Central, 06 August 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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
# This code is an extension of http://scpyce.org/12/

import numpy as np
from scipy import integrate
from matplotlib.pylab import *

def tank(t, y):
    """
    Dynamic balance for a CSTR

    C_A = y[0] = the concentration of A in the tank, [mol/L]
    T   = y[1] = the tank temperature, [K]

    Returns dy/dt = [F/V*(C_{A,in} - C_A) - k*C_A^2       ]
                    [F/V*(T_in - T) - k*C_A^2*HR/(rho*Cp) ]
    """
    F = 20.1     # L/min
    CA_in = 2.5  # mol/L
    V = 100.0    # L
    k0 = 0.15    # L/(mol.min)
    Ea = 5000    # J/mol
    R = 8.314    # J/(mol.K)
    Hr = -590    # J/mol
    T_in = 288   # K
    rho = 1.050  # kg/L

    # Assign some variables for convenience of notation
    CA = y[0]
    T  = y[1]

    # Algebraic equations
    k = k0 * np.exp(-Ea/(R*T))  # L/(mol.min)
    Cp = 4.184 - 0.002*(T-273)  # J/(kg.K)

    # Output from ODE function must be a COLUMN vector, with n rows
    n = len(y)      # 2: implies we have two ODEs
    dydt = np.zeros((n,1))
    dydt[0] = F/V*(CA_in - CA) - k*CA**2
    dydt[1] = F/V*(T_in - T) - (Hr*k*CA**2)/(rho*Cp)
    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 = 45.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
T_t_zero = 295.0
r.set_initial_value([CA_t_zero, T_t_zero], t_start)

# Additional Python step: create vectors to store trajectories
t = np.zeros((num_steps, 1))
CA = np.zeros((num_steps, 1))
temp = np.zeros((num_steps, 1))
t[0] = t_start
CA[0] = CA_t_zero
temp[0] = T_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]
    temp[k] = r.y[1]
    k += 1

# All done!  Plot the trajectories in two separate plots:
fig = figure()
ax1 = subplot(211)
ax1.plot(t, CA)
ax1.set_xlim(t_start, t_final)
ax1.set_xlabel('Time [minutes]')
ax1.set_ylabel('Concentration [mol/L]')
ax1.grid('on')

ax2 = plt.subplot(212)
ax2.plot(t, temp, 'r')
ax2.set_xlim(t_start, t_final)
ax2.set_xlabel('Time [minutes]')
ax2.set_ylabel('Temperaturere [K]')
ax2.grid('on')

fig.savefig('coupled-ode-plot.png')
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 extend the example from http://scpyce.org/12/ here to integrate 2 coupled ODEs and include several algebraic equations.

Consider a chemical reactor with a second order reaction \({\sf A} \rightarrow {\sf B}\). Our aim is calculate the concentration, \(C_{\sf A}(t)\), and tank temperature, \(T(t)\), as functions of time using an ODE integrator.

We know the reaction rate is \(r = k C_{\sf A}^2\), and is an algebraic function of temperature: \(k = 0.15 e^{- E_a/(RT)}\) L/(mol.min). Furthermore, the reaction is known to release heat, with \(\Delta H_r = -590\) J/mol. The time-dependent mass and energy balance for this system can be derived:

\[\begin{split}\frac{dC_{\sf A}(t)}{dt} &= \frac{F(t)}{V} \bigg( C_{{\sf A},\text{in}} - C_{\sf A} \bigg) - 0.15 e^{- E_a/(RT)} C_{\sf A}^2 \\ \frac{dT(t)}{dt} &= \frac{F(t)\rho C_p(T)}{V \rho C_p(T)}\bigg(T_\text{in} - T(t) \bigg) - \frac{0.15 e^{- E_a/(RT)} C_{\sf A}^2 V (\Delta H_r)}{V \rho C_p}\end{split}\]

Notice that these are coupled ODEs, as the first ODE is a function of \(T(t)\), while the second ODE is a function of \(C_{\sf A}(t)\). We also know that:

  • \(C_{{\sf A},\text{in}} = 2.5\) mol/L,
  • \(T_\text{in} = 288\) K
  • a constant volume of \(V = 100\) L
  • constant inlet flow of \(F(t) = 20.1\) L/min, though we could easily make it a function of time and adjust the function tank to use the time variable, t
  • molar heat capacity is a weak function of the system temperature: \(C_p(T) = 4.184 - 0.002(T-273)\) J/(kg.K),
  • liquid phase density is \(\rho=1.05\) kg/L.
  • \(E_a = 5000\) J/mol
  • \(R = 8.314\) J/(mol.K),

We need two initial conditions, one for each of the ODEs:

  • \(C_{\sf A}(t=0) = 0.5\) mol/L
  • \(T(t=0) = 295\) K

In the code we will integrate the ODE from \(t_\text{start} = 0.0\) minutes up to \(t_\text{final} = 45.0\) minutes and plot the time-varying trajectories of temperature and concentration in the tank using matplotlib.

http://scipy-central.org/media/scipy_central/images/201108/coupled-ode-plot.jpg

Read the SciPy documentation for ode.