View source for Assignment 6 - 2010 - Solutions

{{OtherSidebar
| due_dates = 01 December 2010
| dates_alt_text = Due date(s)
| questions_PDF = Assignment-6-2010-Solutions.pdf
| questions_text_alt = Assignment questions and solutions
| solutions_PDF = 
| solutions_text_alt = Solutions, prepared by Ali, Elliot and Kevin
| other_instructions = Hand-in at class.
}}

<rst>
<rst-options: 'toc' = False/>
<rst-options: 'reset-figures' = False/>
.. |m3| replace:: m\ :sup:`3`

.. rubric:: To integrate ODEs - both by hand and with computer software.
 	
.. rubric:: Solutions by Ali, Elliot and Kevin.

.. Potential questions:

	Exothermic reaction in a PFR; does temperature exceed 300C over the distance of the pipe?  (See reactor design project)

	Epidemic example, p 557 in Friendly book
	
	Stiff system?
	
	
Question 1 [2]
===============

In assignment 1 you derived the dynamic balance for a waste water treatment system (aerobic heterotrophic growth using Monod kinetics); and in assignment 3 you used Newton's method to find the steady-state solution to the nutrient (:math:`N`) and biomass (:math:`B`) balances:

.. math::
		\frac{dN}{dt} &= N_\text{in}(t)\frac{F(t)}{V} - N(t)\frac{F(t)}{V} - \left(\frac{1}{Y_B}\right)\left(\mu_{\text{max}}\frac{N}{K+N(t)}\right)B(t) \\ 		\frac{dB}{dt} &= - B(t)\frac{F(t)}{V} + \left(\mu_{\text{max}}\frac{N(t)}{K+N(t)}\right)B(t)

*	The tank volume is assumed constant at :math:`V = 1600` |m3|
*	The conversion efficiency is constant at :math:`Y_B=0.80`
*	The two Monod kinetics parameters are: :math:`\mu_\text{max} = 5 \,\,\text{hr}^{-1}` and :math:`K = 20` g/|m3|
*	The inlet flow rate is :math:`F(t) = 5000` |m3|/hr
*	The nutrient inflow is :math:`N_\text{in}(t) = 100` g/|m3|

These dynamic balances describe how nutrients in the inlet (i.e. the :math:`N_\text{in} = 100` g/|m3| of waste in the wastewater) are converted to biomass, B. The stream leaving the well-mixed CSTR should have lower nutrient levels, :math:`N(t)`.

#.	Code these equations in a MATLAB or Python function and use the built-in integrators in either software package to find the steady state solutions when using starting values of :math:`N(t=0) = 100` g/|m3| and :math:`B(t) = 20` g/|m3| and show that you obtain the same steady state values as found in assignment 3.

#.	Most real processes are never at steady state.  In this reactor a (slightly) more realistic input, to account for the diurnal nature of the flow, is: 

	.. math::
	
		F(t) = 5000+900\sin\left(\displaystyle \frac{2\pi}{24}t -\frac{\pi}{5}\right)

	*	Plot 100 hours of process operation at these conditions; use the steady state values from part 1 as your initial conditions.  
	*	Environmental limits prohibit outlet nutrient concentration from exceeding 60 g/|m3|.  Does this system meet those standards?

#.	Using the simulation from part 2 (with the sinusoidal input) for the biomass and nutrient profiles, see what happens if your plant accidentally dumps nutrient into the inlet stream for 10 hours when an upstream unit fails.  Use an increased nutrient influx of :math:`N_\text{in}(t) = 150` g/|m3| between :math:`55 \leq t \leq 65` hours.

	*	Does the plant still manage to treat this accident without violating the 60 g/|m3| outlet limit on :math:`N(t)`?
	*	Would it handle this problem had it occurred between :math:`75 \leq t \leq 85` hours?

.. Solution:  see code/wastewater.py file

Solution
------------

#.	The code simulating the dynamical system is given below. Note that the same code was used to solve all parts of this problem. So, some of the input values or settings may need to change to suit the conditions given in the different parts of this question.

	Solving the dynamical system with the given initial conditions results in the following plots for :math:`N(t)` and :math:`B(t)`. The steady-state values of :math:`N_s=33.33` and :math:`B_s=53.33` are seen to match the results obtained in Assignment 3 using nonlinear solution of the steady-state equations.

	.. figure:: ../che3e4/Assignments/Assignment-6/images/Q1P1.png
		:scale: 50
		:align: center
		
#.	 The steady-state values from Part 1 are used as initial conditions, i.e. :math:`N(0)=33.33` and :math:`B(0)=53.33`. The input flow rate is now fluctuating rather than being constant. This can be seen in code shown earlier. The simulation results given the new conditions are plotted below.

	.. figure:: ../che3e4/Assignments/Assignment-6/images/Q1P2.png
		:scale: 50
		:align: center

	As seen from the plot, the outlet nutrient concentration :math:`N(t)` does not exceed the maximum allowable limit of 60 g/|m3| (It hits this point once though.)

3.	In case of an upstream failure where the nutrient influx is suddenly increased to :math:`N_\text{in}(t) = 150`, we can simulate the system by modelling the nutrient influx as a function of time. It is :math:`N_\text{in}(t) = 150` if :math:`55 \leq t \leq 65` and :math:`N_\text{in}(t) = 100` otherwise (see the code that shows the implementation). Re-integrating the dynamic equations results the following plots below.

	.. figure:: ../che3e4/Assignments/Assignment-6/images/Q1P3a.png
		:scale: 50
		:align: center

	It is seen that due to the sudden increase in the nutrient influx at :math:`t=55`, the nutrient out flow also spikes and exceeds the allowable limit of 60 g/|m3|. However, this situation does not last as the nutrient outlet decreases as a result of its conversion to :math:`B`.

	If the accident occurred between :math:`75 \leq t \leq 85`, the plant would behave as follows, which shows that the nutrient outlet would remain within its limit.

	.. figure:: ../che3e4/Assignments/Assignment-6/images/Q1P3b.png
		:scale: 50
		:align: center
		
The following code was used in all 3 parts.

.. twocolumncode::
    :code1: ../che3e4/Assignments/Assignment-6/code/wastemodel.m
    :language1: matlab
    :header1: MATLAB code
    :code2: ../che3e4/Assignments/Assignment-6/code/wastewater.py
    :language2: python
    :header2: Python code

Question 2 [3]
===============

.. Reactor system, with recycle loop; disturbance in middle reactor at t = 10; see the disturbance propagate through the system and the recycle loop

.. Series of CSTR's with reentry into middle CSTR

Consider three CSTR's in series, operating at constant volume of :math:`V = 100` L.  A first order reaction, :math:`{\sf A \rightarrow B + C}` takes place in the system, with reaction rate = :math:`r = 0.09 C_{\sf A}`.  The inlet conditions to the first reactor are constant with time: inlet concentration = :math:`C_{{\sf A},0} = 2.7` mol/L at a flow rate of :math:`F_\text{in}^{0}` = 15.8 L/min.  The outlet from the first reactor, :math:`C_{{\sf A},1}`, is the inlet to the second reactor, and the outlet from the second reactor, :math:`C_{{\sf A},2}`, is passed to third reactor.  The outlet from the third reactor, :math:`C_{{\sf A},3}`, can be used to calculate the overall conversion = :math:`\displaystyle \frac{C_{{\sf A},0} - C_{{\sf A},3}}{C_{{\sf A},0}}`.

The dynamic balances in each reactor are:

.. math::
	
	\frac{d C_{{\sf A},1}(t)}{dt} &= \frac{F}{V_1}\bigg(C_{{\sf A},0} - C_{{\sf A},1} \bigg) - kC_{{\sf A},1} \\ 	\frac{d C_{{\sf A},2}(t)}{dt} &= \frac{F}{V_2}\bigg(C_{{\sf A},1} - C_{{\sf A},2} \bigg) - kC_{{\sf A},2} \\	\frac{d C_{{\sf A},3}(t)}{dt} &= \frac{F}{V_3}\bigg(C_{{\sf A},2} - C_{{\sf A},3} \bigg) - kC_{{\sf A},3} \\
	
1.	Describe how you would ordinarily solve the system of equations :math:`f({\bf x}) = {\bf 0}` where :math:`{\bf x} = \big[C_{{\sf A},1}, C_{{\sf A},2}, C_{{\sf A},3} \big]` to find the steady-state values of :math:`\big[C_{{\sf A},1}, C_{{\sf A},2}, C_{{\sf A},3} \big]`.   Now use a computer-based ODE solver to integrate the equations to steady state and report the steady-state values of :math:`{\bf x} = \big[C_{{\sf A},1}, C_{{\sf A},2}, C_{{\sf A},3} \big]` and the overall conversion.

	.. [1.720, 1.095, 0.698], conversion = (2.7 - 0.698)/2.7 = 74.1%

You are tasked with increasing the conversion of the overall system.  You are investigating two possibilities: adding a recycle stream, and operating at larger tank volumes.

2.	What are the new steady state values of :math:`{\bf x} = \big[C_{{\sf A},1}, C_{{\sf A},2}, C_{{\sf A},3} \big]` and the overall conversion when :math:`V_1 = V_2 = V_3 = 150` L?  Why does operating at a larger tank volumes increase conversion?  You can answer this question even if you haven't taken the reactor design course yet - *hint*: interpret the above modelling equations.

	.. [1.456, 0.785, 0.423], conversion = (2.7 - 0.423)/2.7 = 84.3%

3.	Let's take a look at what happens when a disturbance enters the system.  Simulate the case when a step change of +1.2 mol/L in :math:`C_{{\sf A},0}` is made at time :math:`t = 50` (i.e. the inlet concentration is 2.7 mol/L for 50 minutes, then increased to 3.9 mol/L for the rest of the time).  Choose suitable initial conditions so that your system is roughly at steady state before :math:`t=50`, continue the simulation on until the new steady state is reached.

	What do you notice about the shape of the concentration *vs* time trajectories from each reactor (superimpose the 3 trajectories on the same plot).


Adding a recycle stream from the final reactor to the first or second reactor is an exercise that you can complete on your own.  It leads to a counter-intuitive phenomenon: you'd expect an increase in conversion, but ....

Solution
---------

#.	Ordinarily one would use Newton's method to solve a system of :math:`N` non-linear equations in :math:`N` unknowns.  In this case however, the system of equations are linear, so we could use any of the linear algebra methods, such as Gauss reduction, LU decomposition,  Jacobi method, *etc*.

	Using the ``ode`` solvers in computer software and integrating for a *long time* from any starting point will get us the steady state values of :math:`{\bf x} = \big[C_{{\sf A},1}, C_{{\sf A},2}, C_{{\sf A},3} \big] = \big[1.720, 1.096, 0.6982 \big]` mol/L and the overall conversion is 
	
	.. math::
		
		X_\text{overall} = \frac{C_{{\sf A},0} - C_{{\sf A},3}}{C_{{\sf A},0}} = \frac{2.7 - 0.6982}{2.7} = {\bf74.1\%}

#.	Changing the volumes to 150 L in the code and running the integrator to steady state gives :math:`\big[C_{{\sf A},1}, C_{{\sf A},2}, C_{{\sf A},3} \big] = \big[ 1.456, 0.785, 0.423\big]` mol/L, and an overall conversion of :math:`{\bf 84.3\%}`.

	Operating at larger tank volumes can be explained as follows.  Consider the first tank at steady state and calculate its conversion.  We  let :math:`\tau = \frac{F}{V}` for convenience (but it also has an interpretation, called *residence time*, which you will/have learn about in reactor design).
	
	.. math::
	
		0 &= \frac{F}{V}\bigg(C_{{\sf A},0} - C_{{\sf A},1} \bigg) - kC_{{\sf A},1} \\		\left(k+\tau\right)C_{{\sf A},1} &= \tau C_{{\sf A},0} \\		X_1 = \frac{C_{{\sf A},0} - C_{{\sf A},1}}{C_{{\sf A},0}} &= 1 - \frac{C_{{\sf A},1}}{C_{{\sf A},0}} = 1 - \frac{C_{{\sf A},0}}{C_{{\sf A},0}}\cdot \frac{\tau}{k+\tau}\\		X_1 &= \frac{\tau}{k + \tau}
		
	And you can show for tanks 2 and 3 that the same equation holds:
	
	.. math::
	
		X_2 &= \frac{\tau}{k + \tau} \\		X_3 &= \frac{\tau}{k + \tau} 
		
	So the outlet concentration is the product of all conversions:
	
	.. math::
	
		C_{{\sf A},3}  = X_3 C_{{\sf A},2}  = X_3 X_2 C_{{\sf A},1} = X_3 X_2 X_1 C_{{\sf A},0} = \left( \frac{\tau}{k + \tau} \right)^3 C_{{\sf A},0}
		
	Overall conversion can be increased in any way that decreases :math:`C_{{\sf A},3}`: 
	
		*	increasing the rate constant, :math:`k` (maybe with a catalyst?)
		*	making the value of :math:`\tau` smaller, which implies using a lower inlet flow rate, :math:`F` or using increased volumes, :math:`V`.
		
	However, notice the diminishing returns for the last bullet point above, because the term appears both in the numerator and denominator.  This can be seen in the plot here of reactor volume (for each CSTR) against overall conversion:
	
	.. figure:: ../che3e4/Assignments/Assignment-6/images/volume-vs-conversion.png
		:scale: 80
		:align: center
		
#.	The trajectories for this disturbance, which enters the first tank, are shown below.

	.. figure:: ../che3e4/Assignments/Assignment-6/images/disturbance-propogation.png
		:scale: 80
		:align: center

	The key point is that the first CSTR response has a sharp, rising edge (discontinuity) as the trajectory rises immediately when the disturbance hits.  The second reactor has a slightly smoother response, while the third reactor responds even slower still.  This is typical behaviour of all series reactors -- disturbances in early reactors get "washed out" in later reactors.
	
	
The following Python code was used to solve this question; MATLAB code (not available) would be similar to this.

.. literalinclude:: ../che3e4/Assignments/Assignment-6/code/series_of_cstrs.py
		
Question 3 [2]
================

.. note:: This question is adapted from the textbook by Bradie, listed under the `Suggested readings <http://modelling3e4.connectmv.com/wiki/Suggested_readings>`_ sections of the course website.

A genetic switch is a biological mechanism to activate whether or not a particular protein product from a cell is synthesized.  For example, under certain conditions a pigment protein can be synthesized, which indicates the presence/absence of another chemical (the switch).

This model has been proposed for a genetic switch:

.. math::

	\frac{dg}{dt} = 0.01s - 1.51g + 3.03 \frac{g^2}{1+g^2}\qquad g(0) = 0.0

In this notation, the pigment protein concentration is :math:`g(t)` and :math:`s(t)` is the concentration of the chemical (the switch) that activates the gene that produces the pigment.

The :math:`-1.51g` term indicates the protein's natural decay rate is proportional to its concentration, while the last term describes the positive feedback  the protein exerts on its own formation.

A switch has two important characteristics:

	*	Below a threshold value of :math:`s` the gene concentration must be close to zero (i.e. little/no pigment produced means the switch is off), while above a threshold, the gene concentration, :math:`g`, must be high (i.e. the switch has been turned on).
	
	*	Once turned on (i.e. high values of :math:`g`), the switch must remain on, even if :math:`s` is returned below the threshold value.  Similarly, if the switch was never on to begin with, it must still remain "off" if :math:`s` is below the threshold.  This phenomenon is known as the hysteresis effect.
	

#.	First we wish to find the threshold value of :math:`s` and interpret what :math:`s` means.  Plot and label four trajectories of :math:`g(t)` using :math:`s=20.0, 20.33, 20.67, 21.0` over a time range of 100 minutes, and describe what the plots are showing.  

	*	Be specific: how would you use these plots to design your genetic switch?  
	*	What does :math:`s` mean in your design? (Hint: it has two interpretations)

#.	Now we wish to verify if the hysteresis effect is present, because it is important the switch cannot be turned off.  For each of the four values of :math:`s` in part 1, run a simulation where :math:`s(t>150)=0` (i.e. at time :math:`t=150` the concentration of :math:`s` is set from its current value to zero).  Show your results over a time horizon of 300 minutes and **comment on the results**.

.. solution: code/genetic_switch.py

Solution
-----------

#.	Trajectories, at different levels of :math:`s`, are simulated over a time scale of 100 minutes in the plot.  


	:math:`s` must be at least greater than 20.0 for the switch to be activated; values at or below 20.0 never activate the switch, while values above this will activate it, though slowly at first (e.g 20.33 takes a long time to show).  So one interpretation of :math:`s` is that it is the threshold level above which we can detect it.  The other interpretation is that the higher the value of :math:`s`, the faster it will be detected by the gene.  So :math:`s` affects both the threshold of the switch, and the time taken to activate it.

	.. figure:: ../che3e4/Assignments/Assignment-6/images/Q3P1.png
		:scale: 50
		:align: center

#.	The hysteresis effect is present. The simulations show that for all four levels of :math:`s`, that if the switch is activated (:math:`s>20.0`) that the gene will remain present, even if :math:`s` is removed.

	We can set an arbitrary threshold on the gene concentration, e.g. 0.8, so that above this level we can consider the switch to be on, and below it the switch is off.

	.. figure:: ../che3e4/Assignments/Assignment-6/images/genetic-switch.png
		:scale: 50
		:align: center

The following code was used in both parts of the question.

.. twocolumncode::
    :code1: ../che3e4/Assignments/Assignment-6/code/bioswitch.m
    :language1: matlab
    :header1: MATLAB code
    :code2: ../che3e4/Assignments/Assignment-6/code/genetic_switch.py
    :language2: python
    :header2: Python code

Question 4 [2]
==============

.. note:: This question continues on from question 5, assignment 5.  It was adapted from a question in the final exam, 2009, worth 3+1 out of 50 marks.

Consider a well-mixed tank containing species :math:`{\sf A}` at concentration :math:`C_{\sf A}=1` mol/L. There is no inlet or outlet from the tank, i.e. it is a batch system. A fluid of constant density is considered throughout this question.  The reactor is filled to 5.2 m\ :sup:`3`.

A catalyst is introduced in the tank at time :math:`t=0` minutes to initiate the reaction :math:`{\sf A} \rightarrow {\sf B} + {\sf C}`. Due to quick catalyst deactivation, the dynamic equation describing the evolution of concentration :math:`C_{\sf A}` in the tank is given by:

.. math::

	\frac{{\rm d}C_{\sf A}(t)}{{\rm d}t} = - k\,C_{\sf A}(t)\,e^{-\frac{t}{2}}
	
with :math:`k=1\ \rm min^{-1}`.

#.	Determine the concentration of A in the reactor at time :math:`t=20` minutes using the classical Runge-Kutta method with step size :math:`h=10` minutes.

#.	Compare and discuss your answer with respect to the analytical solution: 

	.. math::
	
		C_{\sf A}(t) = C_{\sf A}(t=0)\ \exp\left[2\,k\left(e^{-\frac{t}{2}}-1\right)\right]
	

Solution
----------

.. Solution: see code/catalyst.py

For notation simplicity, let us denote :math:`y=C_{\sf A}` and :math:`f(t,y)=\displaystyle\frac{{\rm d}y(t)}{{\rm d}t}=- k\,y(t)\,e^{-\frac{t}{2}}`. :math:`y_0= y(t=0)=1` and the step size :math:`h=10`.

Classical 4th-order Runge-Kutta formula:

.. math::

	y_1 &= y_0 + \frac{1}{6}[K_1 + 2K_2 + 2K_3 + K_4 ] h \\ 	K_1 &= f(t_0, y_0)= f(0,1) = -1\times 1e^0 = -1\\		K_2 &= f(t_0+\frac{1}{2}h,y_0 + \frac{1}{2}hK_1)=f(5,-4)=0.328\\	K_3 &= f(t_0+\frac{1}{2}h,y_0 + \frac{1}{2}hK_2)=f(5,2.642)=-0.217\\	K_4 &= f(t_0+h, y_0 + hK_3)=f(10,-1.17)=0.008\\	y_1 &= 1 + \frac{1}{6}\bigg[-1 + 2\times 0.328 + 2\times -0.217 +0.008 \bigg] \times 10={\bf -0.283}

Does the result make sense? We ended up a negative concentration! Let us double-check with the analytical solution. The actual value of :math:`y(t=10)` using the analytical solution is  0.1372 [mol/L]. So why did we get a terribly wrong answer using the Runge-Kutta? The answer lies in the stability of ODE methods. Some ODE methods have certain stability restrictions that impose a maximum length on the integration step size. Here we chose :math:`h=10` which is too large for our system, thereby leading to an *unstable solution*. A workaround is to choose a very smaller step size and reach :math:`t=10` in multi-steps instead. This is not requested in this question though. So, let us continue with :math:`h=10` and see what we get as :math:`y(t_2=20)` using the obtained :math:`y(t_1=10)`

.. math::

	y_2 &= y_1 + \frac{1}{6}[K_1 + 2K_2 + 2K_3 + K_4 ] h \\	K_1 &=f(t_1, y_1)=f(10,-0.283)= 0.0019	\\	K_2 &=f(t_1+\frac{1}{2}h, y_1 + \frac{1}{2}hK_1)=f(15,-0.2735)= 1.5127\times 10^{-4}\\	K_3 &=f(t_1+\frac{1}{2}h, y_1 + \frac{1}{2}hK_2)=f(15,-0.2822)= 1.5608\times 10^{-4}\\	K_4 &=f(t_1+h, y_1 + hK_3)=f(20,-0.2814)= 1.2776\times 10^{-5} \\	y_2 &= -0.283 + \frac{1}{6}\bigg[0.0019 + 2\times 1.5127\times 10^{-4}  + 2\times 1.5608\times 10^{-4}  + 1.2776\times 10^{-5} \bigg] \times 10= {\bf -0.2788}

The actual value of :math:`X(t=20)` using the analytical solution is  0.1353 [mol/L] (compare to the above approximate solution!)


Question 5 [3]
==============

.. From Tutorial 9, 2008

Consider the CSTR shown below, where an irreversible reaction of the form :math:`{\sf A}\rightarrow {\sf B}` takes place.

.. figure:: ../che3e4/Assignments/Assignment-6/images/CSTR_1input.png
	:scale: 68
	:align: center
	
**Model parameters**

* :math:`C_{\sf A}^{\rm in} = 0.79 \rm kmol\ m^{-3}`
* :math:`T^{\rm in} = 352.6 \rm K`
* :math:`F = 0.2 \rm m^3\ min^{-1}` 
* :math:`V = 0.1 \rm m^3` 
* :math:`k_{0} = 7.2\times 10^{10} \rm min^{-1}` 
* :math:`E = 8.314\times 10^{4} \rm kJ\ kmol{-1}` 
* :math:`\Delta H = -4.78\times 10^4 \rm kJ\ kmol^{-1}` 
* :math:`C_{p} = 0.239 \rm kJ\ kg^{-1}\ K^{-1}` 
* :math:`\rho = 1000 \rm kg\ m^{-3}` 
* :math:`R = 8.314 \rm kJ\ kmol^{-1}\ K^{-1}` 

#.	Derive, by hand, the differential equations describing the evolution of the concentration of species :math:`{\sf A}`, :math:`C_{\sf A}(t)`, as well as of the temperature of the outlet stream, :math:`T(t)`, based on the following assumptions:

 	*	The concentrations of :math:`{\sf A}` in the inlet stream is :math:`C_{\sf A}^{\rm in}`, the inlet and outlet volumetric flow rates are constant and both equal to :math:`F`, and the temperature of the inlet stream is :math:`T^{\rm in}`.
	*	The reaction has first-order kinetics, with the dependence of the reaction rate :math:`r` (moles of :math:`{\sf A}` consumed per unit volume and per unit time) on the temperature given by :math:`r=k_0\, e^{-\frac{E}{RT}}\, C_{\sf A}`.
			
		where :math:`E` is the activation energy; :math:`R`, the universal gas constant; :math:`T`, the temperature at which the reaction is taking place; and :math:`k_0`, a constant.
		
	*	The heat of reaction is :math:`\Delta H` and no heat is lost or added to the reactor, other than through the liquid streams.
	
	*	The density and specific heat capacity of the fluid are :math:`\rho` and :math:`C_p` respectively.

#.	For the values of the process parameters listed in the table and the following initial conditions :math:`C_{\sf A}(0) = 0.7\ \rm kmol/m^3` and :math:`T(0)=355\ \rm K`, compute, by hand, the value of :math:`\big[C_{\sf A}(t),\; T(t)\big]` at :math:`t=0.2`, using:

	*	Euler's method with :math:`h=0.1`, and
	*	Heun's predictor-corrector method with :math:`h = 0.1`, and
	*	the fourth-order Runge-Kutta method with :math:`h=0.2`
	
#.	Implement these equations in computer software and compare the integration by hand with the computer version at :math:`t=0.2`.

Solution
--------

.. Solution: see code/cstr_integration.py


#.	First note that the rate of consumption of component :math:`{\sf A}` is :math:`-r_{\sf A} = kC_{\sf A} = k_0 \exp\left(-\frac{E}{RT}\right) C_{\sf A}`.

	A mole balance on :math:`{\sf A}` can be done to determine the evolution of the concentration of species :math:`{\sf A}`:

	.. math::

	    \dfrac{d V C_{\sf A}}{dt} &= F C_{\sf A}^{\rm in} - F C_{\sf A} - r_{\sf A} V \\	    V\dfrac{d C_{\sf A}}{dt}  &= F C_{\sf A}^{\rm in} - F C_{\sf A} - k_0 \exp\left(-\frac{E}{RT}\right) C_{\sf A} V \\

	On the other hand, an energy balance on the reactor can be done to determine the evolution of the temperature in the reactor:

	.. math::

		\dfrac{d (\rho V C_p (T-T_{\rm ref}))}{dt} &=  \rho C_p F (T_{\rm in} - T_{\rm ref}) - \rho C_p F (T - T_{\rm ref})  + (-\Delta H ) r_{\sf A} V  \\		\rho V C_p\dfrac{d T}{dt} &= \rho C_p F (T_{\rm in} - T) - \Delta H k_0 \exp\left(-\frac{E}{RT}\right) C_{\sf A} V

	
#.	Substituting the values from the table in the dynamic balances gives:

	.. math::
	
		\frac{dy_1}{dt} =& = f_1(t, {\bf y}) = 1.58 - 2y_1 - 7.2\times 10^{10} \; e^{-\frac{10^4}{y_2}}  y_1\\		\frac{dy_2}{dt} =& = f_2(t, {\bf y}) = 705.2 - 2y_2 + 1.44 \times 10^{13} \; e^{-\frac{10^4}{y_2}} y_1,
		
	with :math:`y_1(t)=C_{\sf A}(t)$ and $y_2(t)=T(t)`.
	
	*	**Explicit Euler method**
	
		From :math:`t` = 0 to 0.1
		
		At :math:`t=0`, :math:`\boldsymbol{\phi}={\bf f}(0,{\bf y}^{(0)}) = \begin{bmatrix} 0.1506 \\ 1.0859 \end{bmatrix}`. The estimated solution at :math:`t=0.1` is:
		
		.. math::
			
			{\bf y}^{(1)} &= {\bf y}^{(0)} + h \boldsymbol{\phi} = \begin{bmatrix} 0.7 \\ 355 \end{bmatrix} + 0.1 \begin{bmatrix} 0.1506 \\  1.0859 \end{bmatrix} = \begin{bmatrix} 0.7151 \\  355.1086 \end{bmatrix}

		From :math:`t` = 0.1 to 0.2:
		
		At :math:`t=0.1`, :math:`\boldsymbol{\phi}={\bf f}(0.1,{\bf y}^{(1)}) = \begin{bmatrix} 0.1196 \\ 1.0473 \end{bmatrix}`. The estimated solution at :math:`t=0.2` is:
		
		.. math::
		
			{\bf y}^{(2)} &= {\bf y}^{(1)} + h \boldsymbol{\phi} \\			&= \begin{bmatrix} 0.727 \\  355.2 \end{bmatrix}

		Overall, we thus obtain the following solution:
		
		:math:`C_A(0.2) = 0.727\ \rm kmol/m^{3} \quad\text{and}\quad T(0.2) = 355.2\ \rm K`
		
	*	**Heun's predictor-corrector method**
	
		From :math:`t` = 0 to 0.1
		
		The first step of Heun's method (predictor) uses the Euler slope, calculated above.  The corrector step calculates the slope at :math:`t=0.1` and the Heun slope is just the average of the two.
		
		.. math::
			\hat{\bf y}^{(1)} &= \begin{bmatrix} 0.7151 \\  355.1086 \end{bmatrix} \\			\boldsymbol{\phi}_\text{Heun} &= 0.5 \begin{bmatrix} 0.1506 \\ 1.0859 \end{bmatrix} + 0.5 {\bf f}(0.1, \hat{\bf y}^{(1)}) \\                                          &= 0.5 \begin{bmatrix} 0.1506 \\ 1.0859 \end{bmatrix} + 0.5 \begin{bmatrix} 0.1196 \\ 1.0473 \end{bmatrix}\\										  &= \begin{bmatrix} 0.1351 \\ 1.0667 \end{bmatrix} \\			{\bf y}^{(1)} &= {\bf y}^{(0)} + h \boldsymbol{\phi}_\text{Heun} = \begin{bmatrix} 0.7 \\ 355 \end{bmatrix} + 0.1 \begin{bmatrix} 0.1351 \\ 1.0667 \end{bmatrix} = \begin{bmatrix} 0.7135 \\  355.107 \end{bmatrix}

		From :math:`t` = 0.1 to 0.2:
		
		.. math::
			\hat{\bf y}^{(2)} &= \begin{bmatrix} 0.7135 \\  355.107 \end{bmatrix} + 0.1\begin{bmatrix} 0.1227 \\ 1.0365 \end{bmatrix} \\			\hat{\bf y}^{(2)} &= \begin{bmatrix} 0.72577 \\ 355.210 \end{bmatrix} \\			\boldsymbol{\phi}_\text{Heun} &= 0.5 \begin{bmatrix} 0.1227 \\ 1.0365 \end{bmatrix} + 0.5 {\bf f}(0.2, \hat{\bf y}^{(2)}) \\                                          &= 0.5 \begin{bmatrix} 0.1227 \\ 1.0365 \end{bmatrix} + 0.5 \begin{bmatrix} 0.097435 \\0.9851 \end{bmatrix}\\										  &= \begin{bmatrix} 0.1101 \\ 1.0108 \end{bmatrix} \\			{\bf y}^{(2)} &= {\bf y}^{(1)} + h \boldsymbol{\phi}_\text{Heun} = \begin{bmatrix} 0.7135 \\ 355.107 \end{bmatrix} + 0.1 \begin{bmatrix} 0.1101 \\ 1.0108 \end{bmatrix} = \begin{bmatrix} 0.7245 \\  355.21 \end{bmatrix}

		Overall, we thus obtain the following solution:
		
		:math:`C_A(0.2) = 0.7245 \ \rm kmol/m^{3} \quad\text{and}\quad T(0.2) = 355.21 \ \rm K`
		
	*	**4th-order Runge-Kutta method (known as the classical Runge-Kutta method)**
	
		From :math:`t` = 0 to 0.2:
		
		At :math:`t=0`, :math:`{\bf y}^{(0)} = \begin{bmatrix} 0.7 \\ 355 \end{bmatrix}`. The intermediate quantities :math:`{\bf K}_1,\ldots,{\bf K}_4` are now 2-dimensional vectors:
		
		 
		.. math::
		
			{\bf K}_1 &= \begin{bmatrix} 0.1506 \\ 1.0859 \end{bmatrix} \\			{\bf K}_2 &= \begin{bmatrix} 0.1196 \\ 1.0473 \end{bmatrix} \\			{\bf K}_3 &= \begin{bmatrix} 0.1259 \\ 1.0269 \end{bmatrix} \\			{\bf K}_4 &= \begin{bmatrix} 0.0986 \\ 0.9870 \end{bmatrix}    

			\boldsymbol{\phi} &= \dfrac{1}{6} {\bf K}_1 + \dfrac{1}{3} {\bf K}_2 +\dfrac{1}{3} {\bf K}_3 + \dfrac{1}{6} {\bf K}_4 = \begin{bmatrix} 0.1234 \\ 1.0369 \end{bmatrix}
			
		The estimated solution at :math:`t=0.2` is:
		
		:math:`{\bf y}^{(1)} = {\bf y}^{(0)} + h \boldsymbol{\phi} = \begin{bmatrix} 0.7247 \\ 355.2074 \end{bmatrix}`

		Overall, we thus obtain the following solution:
		
		:math:`C_A(0.2) = 0.725\ \rm kmol/m^{3} \quad\text{and}\quad T(0.2) = 355.2\ \rm K`, which is in good agreement with the solution found earlier with the explicit Euler method and Heun's method.	

#.	Integrating the ODE equations with MATLAB's ``ode45`` function, or Python's ``integrate`` function gave the following results at :math:`t=0.2`: :math:`{\bf y} = \left[0.72468, 355.2073\right]`.

	:math:`C_A(0.2) = 0.725\ \rm kmol/m^{3} \quad\text{and}\quad T(0.2) = 355.2\ \rm K`, which agrees very well with the Runge-Kutta solution, and less well with the Heun and Euler methods.

.. raw:: latex

	\vspace{0.25cm}
	\hrule
	\begin{center}END\end{center}
</rst>