# SciPy Central

## 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)

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.

Read the SciPy documentation for ode.