Tutorial 5 - 2010 - Solution/Question 1

From Process Model Formulation and Solution: 3E4
Jump to navigation Jump to search
Back to all questions Next step: Question 2 →

<rst> <rst-options: 'toc' = False/> <rst-options: 'reset-figures' = False/>

Question 1 [1]

===

Show the first 4 iterations of the bisection method to solve for :math:`T`, justifying your choice for the initial bracket.

For iteration 2, 3, and 4, please also report these two relative errors:

.. math::

\varepsilon_\text{rel,x}^{(k)} = \displaystyle \left|\frac{x_m^{(k)} - x_m^{(k-1)}}{x_m^{(k)}} \right| \qquad\qquad \varepsilon_\text{rel,f}^{(k)} = \displaystyle \left|\frac{f(x_m^{(k)}) - f(x_m^{(k-1)})}{f(x_m^{(k)})} \right|

Solution

Since we are attempting to solve for the temperature associated with :math:`\Delta H_{r}^{0}(T) = -23505` cal/mol we are essentially asked to solve for the zero of the following function:

.. math::

f(T) = - 24097 - 0.26T + 1.69\text{x}10^{-3}T^{2} + \frac{1.5\text{x}10^{5}}{T} - (-23505)

where

.. math::

\Delta H_{r}^{0}(T) = - 24097 - 0.26T + 1.69\text{x}10^{-3}T^{2} + \frac{1.5\text{x}10^{5}}{T}

To apply the bisection method we first need to select a starting interval that the root is known to fall within. More specifically, we need to also be sure that the root is a "true root" in the sense that it crosses the y-axis (and therefore the function changes sign on either side of the root). **The bisection method is unable to locate roots that only touch the x-axis but do not cross it, e.g. a double root**. Given that this is a physically perceivable system of one variable we can select the lower and upper bounds graphically or based on prior engineering knowledge.

If we choose to attack the problem based on our previous knowledge then selecting the lower bound is easy since we are dealing with the Kelvin temperature scale. Therefore we know the lowest allowable temperature is 0 K. Selecting the upper bound based on prior knowledge is slightly trickier. If we knew what reaction this equation is describing then we might have some basic knowledge of how the heat of reaction changes with T. In this case we would be able to guess a suitable upper temperature. Other than this we are essentially stuck trying to overguess to ensure we have a root contained in our bounds. This technique, of course, carries with it the critical issue that if the function contains multiple roots (as it will be shown this one does) an over guess could easily contain more than one root and so would either confuse the solver or cause the solver to think that no roots exist (if the function is a quadratic shape, for example, then it is possible that the function could have the same sign on either end of a bounding interval despite containing both roots). For this reason I would recommend the graphical approach.

Plotting the function in MATLAB/Python is accomplished using the following scripts.

.. twocolumncode:: :code1: ../che3e4/Tutorials/Tutorial-5/code/Q1_Plot.m :language1: matlab :header1: MATLAB code :code2: ../che3e4/Tutorials/Tutorial-5/code/Q1_Plot.py :language2: python :header2: Python code

The result of this script is shown below. From the plot below it is clear that this function has two roots. One between the bounds (200,400) and one between the bounds (400,600).

.. figure:: ../che3e4/Tutorials/Tutorial-5/images/Q1_Figure.png :scale: 70 :align: center

As a final note before getting into the bisection solution, it is believed that some confusion has likely arisen over the function error tolerance equation. In this case f(x) refers to the original :math:`\Delta H_{r}^{0}(T)` function and not the vertically shifted function you are trying to find the roots of (i.e. :math:`\Delta H_{r}^{0}(T) + 23505`). Since in your adjusted function you are trying to solve for when the function is 0, it does not make sense to calculate a relative comparison to 0 (the numbers will explode and contain essentially no useful information). This is actually an error that can crop up when a root also happens to be zero, but naturally this is less easy to predict. So to quickly wrap up, when root finding, your best bets for stopping criteria are either the relative error on x from step to step (unless the root is zero), the value of the root finding function (this would be the most versatile, as you are no longer dividing by any value and you know your function will be going to zero), the relative error of the original function from step to step (unless this also is equal to zero at the root...).

For the bisection test we will consider the lower of the two bounds. Utilizing the procedure given on slide 8 of the *Section C: Nonlinear Algebraic Equations* slide set.

    • ITERATION 1**
  • STEP 0: INITIALIZATION*
  • :math:`T_{U}^{(0)} = 400`
  • :math:`T_{L}^{(0)} = 200`
  • :math:`\epsilon_{tol} = 10^{-6}` (for the sake of this question the selection of :math:`\epsilon_{tol}` was arbitrary)

Calculating the function value at upper and lower bounds.

.. math::

\begin{array}{rl} f\left(T_{U}^{(0)}\right) &= - 24097 - 0.26\left(T_{U}^{(0)}\right) + 1.69\text{x}10^{-3}\left(T_{U}^{(0)}\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(T_{U}^{(0)}\right)} - \left(-23505\right)\\ f\left(T_{U}^{(0)}\right) &= - 24097 - 0.26\left(400\right) + 1.69\text{x}10^{-3}\left(400\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(400\right)} - \left(-23505\right)\\ f\left(T_{U}^{(0)}\right) &= -50.6000\\ & \\ f\left(T_{L}^{(0)}\right) &= - 24097 - 0.26\left(T_{L}^{(0)}\right) + 1.69\text{x}10^{-3}\left(T_{L}^{(0)}\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(T_{L}^{(0)}\right)} - \left(-23505\right)\\ f\left(T_{L}^{(0)}\right) &= - 24097 - 0.26\left(200\right) + 1.69\text{x}10^{-3}\left(200\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(200\right)} - \left(-23505\right)\\ f\left(T_{L}^{(0)}\right) &= 173.6000\\ & \\ \Delta H_{r}^{0}\left(T_{U}^{(0)}\right) &= - 24097 - 0.26\left(T_{U}^{(0)}\right) + 1.69\text{x}10^{-3}\left(T_{U}^{(0)}\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(T_{U}^{(0)}\right)}\\ \Delta H_{r}^{0}\left(T_{U}^{(0)}\right) &= - 24097 - 0.26\left(400\right) + 1.69\text{x}10^{-3}\left(400\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(400\right)}\\ \Delta H_{r}^{0}\left(T_{U}^{(0)}\right) &= -23555.6\\ & \\ \Delta H_{r}^{0}\left(T_{L}^{(0)}\right) &= - 24097 - 0.26\left(T_{L}^{(0)}\right) + 1.69\text{x}10^{-3}\left(T_{L}^{(0)}\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(T_{L}^{(0)}\right)}\\ \Delta H_{r}^{0}\left(T_{L}^{(0)}\right) &= - 24097 - 0.26\left(200\right) + 1.69\text{x}10^{-3}\left(200\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(200\right)}\\ \Delta H_{r}^{0}\left(T_{L}^{(0)}\right) &= -23331.4 \end{array}

Testing if the initial bounds contain a root.

.. math::

\begin{array}{rl} Check &= f\left(T_{U}^{(0)}\right)\times f\left(T_{L}^{(0)}\right) \geq 0 ?\\ Check &= \left(-50.6000\right)\times\left(173.6000\right) \geq 0 ?\\ Check &= \left(-8784.16\right) \geq 0 ? \;\; \rightarrow \;\; \text{No, therefore a root is bounded} \end{array}

  • STEP 1: CALCULATE MID-POINT* :math:`T_{M}^{(1)}`

.. math::

\begin{array}{rl} T_{M}^{(1)} &= \frac{T_{U}^{(0)} + T_{L}^{(0)}}{2}\\ T_{M}^{(1)} &= \frac{\left(400\right) + \left(200\right)}{2}\\ T_{M}^{(1)} &= 300 \end{array}

.. math::

\begin{array}{rl} f\left(T_{M}^{(1)}\right) &= - 24097 - 0.26\left(T_{M}^{(1)}\right) + 1.69\text{x}10^{-3}\left(T_{M}^{(1)}\right)^{2} + \frac{1.5\text{x}10^{5}}{T_{M}^{(1)}} - \left(-23505\right)\\ f\left(T_{M}^{(1)}\right) &= - 24097 - 0.26\left(300\right) + 1.69\text{x}10^{-3}\left(300\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(300\right)} - \left(-23505\right)\\ f\left(T_{M}^{(1)}\right) &= -17.9000\\ & \\ \Delta H_{r}^{0}\left(T_{M}^{(1)}\right) &= - 24097 - 0.26\left(T_{M}^{(1)}\right) + 1.69\text{x}10^{-3}\left(T_{M}^{(1)}\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(T_{M}^{(1)}\right)}\\ \Delta H_{r}^{0}\left(T_{M}^{(1)}\right) &= - 24097 - 0.26\left(300\right) + 1.69\text{x}10^{-3}\left(300\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(300\right)}\\ \Delta H_{r}^{0}\left(T_{M}^{(1)}\right) &= -23522.9 \end{array}

  • STEP 2: CALCULATE* :math:`f(T_{L}^{(0)}) \times f(T_{M}^{(1)})`

.. math::

\begin{array}{rl} Check &= f\left(T_{L}^{(0)}\right)\times f\left(T_{M}^{(1)}\right) < 0 ?\\ Check &= \left(173.6000\right)\times\left(-17.9000\right) < 0 ?\\ Check &= \left(-3107.44\right) < 0 ? \;\; \rightarrow \;\; \text{Yes, therefore the root falls between }T_{L}^{(0)}\text{ and }T_{M}^{(1)}. \end{array}

  • :math:`T_{U}^{(1)} = 300`
  • :math:`T_{L}^{(1)} = 200`
  • :math:`f(T_{U}^{(1)}) = -17.9000`
  • :math:`f(T_{L}^{(1)}) = 173.6000`
  • :math:`\Delta H_{r}^{0}(T_{U}^{(1)}) = -23522.9`
  • :math:`\Delta H_{r}^{0}(T_{L}^{(1)}) = -23331.4`
    • ITERATION 2**
  • STEP 1: CALCULATE MID-POINT* :math:`T_{M}^{(2)}`

.. math::

\begin{array}{rl} T_{M}^{(2)} &= \frac{T_{U}^{(1)} + T_{L}^{(1)}}{2}\\ T_{M}^{(2)} &= \frac{\left(300\right) + \left(200\right)}{2}\\ T_{M}^{(2)} &= 250 \end{array}


.. math::

\begin{array}{rl} f\left(T_{M}^{(2)}\right) &= - 24097 - 0.26\left(T_{M}^{(2)}\right) + 1.69\text{x}10^{-3}\left(T_{M}^{(2)}\right)^{2} + \frac{1.5\text{x}10^{5}}{T_{M}^{(2)}} - \left(-23505\right)\\ f\left(T_{M}^{(2)}\right) &= - 24097 - 0.26\left(250\right) + 1.69\text{x}10^{-3}\left(250\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(250\right)} - \left(-23505\right)\\ f\left(T_{M}^{(2)}\right) &= 48.6250\\ & \\ \Delta H_{r}^{0}\left(T_{M}^{(2)}\right) &= - 24097 - 0.26\left(T_{M}^{(2)}\right) + 1.69\text{x}10^{-3}\left(T_{M}^{(2)}\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(T_{M}^{(2)}\right)}\\ \Delta H_{r}^{0}\left(T_{M}^{(2)}\right) &= - 24097 - 0.26\left(250\right) + 1.69\text{x}10^{-3}\left(250\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(250\right)}\\ \Delta H_{r}^{0}\left(T_{M}^{(2)}\right) &= -23456.375 \end{array}


  • STEP 2: CALCULATE* :math:`f(T_{L}^{(1)}) \times f(T_{M}^{(2)})`

.. math::

\begin{array}{rl} Check &= f\left(T_{L}^{(1)}\right)\times f\left(T_{M}^{(2)}\right) < 0 ?\\ Check &= \left(173.6000\right)\times\left(48.6250\right) < 0 ?\\ Check &= \left(8441.3\right) < 0 ? \;\; \rightarrow \;\; \text{No, therefore the root falls between }T_{M}^{(2)}\text{ and }T_{U}^{(1)}. \end{array}

  • :math:`T_{U}^{(2)} = 300`
  • :math:`T_{L}^{(2)} = 250`
  • :math:`f(T_{U}^{(2)}) = -17.9000`
  • :math:`f(T_{L}^{(2)}) = 48.6250`
  • :math:`\Delta H_{r}^{0}(T_{U}^{(2)}) = -23522.9`
  • :math:`\Delta H_{r}^{0}(T_{L}^{(2)}) = -23456.375`
  • STEP 3: CHECK STOPPING CRITERION*

.. math::

\begin{array}{rl} \epsilon_{tol,x}^{(2)} &= \left| \frac{T_{M}^{(2)}-T_{M}^{(1)}}{T_{M}^{(1)}} \right| < \epsilon_{tol} ? \\ \epsilon_{tol,x}^{(2)} &= \left| \frac{\left(250\right)-\left(300\right)}{\left(300\right)} \right| < \left(10^{-6}\right) ? \\ \epsilon_{tol,x}^{(2)} &= 0.1\bar{6} < \left(10^{-6}\right) ? \;\; \rightarrow \;\; \text{No, therefore keep going}. \\ & \\ \epsilon_{tol,f}^{(2)} &= \left| \frac{\Delta H_{r}^{0}\left(T_{M}^{(2)}\right)-\Delta H_{r}^{0}\left(T_{M}^{(1)}\right)}{\Delta H_{r}^{0}\left(T_{M}^{(1)}\right)} \right| < \epsilon_{tol} ? \\ \epsilon_{tol,f}^{(2)} &= \left| \frac{\left(-23456.375\right)-\left(-23522.9\right)}{\left(-23522.9\right)} \right| < \left(10^{-6}\right) ? \\ \epsilon_{tol,f}^{(2)} &= 0.0028 < \left(10^{-6}\right) ? \;\; \rightarrow \;\; \text{No, therefore keep going}. \end{array}

    • ITERATION 3**
  • STEP 1: CALCULATE MID-POINT* :math:`T_{M}^{(3)}`

.. math::

\begin{array}{rl} T_{M}^{(3)} &= \frac{T_{U}^{(2)} + T_{L}^{(2)}}{2}\\ T_{M}^{(3)} &= \frac{\left(300\right) + \left(250\right)}{2}\\ T_{M}^{(3)} &= 275 \end{array}

.. math::

\begin{array}{rl} f\left(T_{M}^{(3)}\right) &= - 24097 - 0.26\left(T_{M}^{(3)}\right) + 1.69\text{x}10^{-3}\left(T_{M}^{(3)}\right)^{2} + \frac{1.5\text{x}10^{5}}{T_{M}^{(3)}} - \left(-23505\right)\\ f\left(T_{M}^{(3)}\right) &= - 24097 - 0.26\left(275\right) + 1.69\text{x}10^{-3}\left(275\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(275\right)} - \left(-23505\right)\\ f\left(T_{M}^{(3)}\right) &= 9.7608\\ & \\ \Delta H_{r}^{0}\left(T_{M}^{(3)}\right) &= - 24097 - 0.26\left(T_{M}^{(3)}\right) + 1.69\text{x}10^{-3}\left(T_{M}^{(3)}\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(T_{M}^{(3)}\right)}\\ \Delta H_{r}^{0}\left(T_{M}^{(3)}\right) &= - 24097 - 0.26\left(275\right) + 1.69\text{x}10^{-3}\left(275\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(275\right)}\\ \Delta H_{r}^{0}\left(T_{M}^{(3)}\right) &= -23495.2392 \end{array}

  • STEP 2: CALCULATE* :math:`f(T_{L}^{(2)})\times f(T_{M}^{(3)})`

.. math::

\begin{array}{rl} Check &= f\left(T_{L}^{(2)}\right)\times f\left(T_{M}^{(3)}\right) < 0 ?\\ Check &= \left(48.6250\right)\times\left(9.7608\right) < 0 ?\\ Check &= \left(474.6189\right) < 0 ? \;\; \rightarrow \;\; \text{No, therefore the root falls between }T_{M}^{(3)}\text{ and }T_{U}^{(2)}. \end{array}

  • :math:`T_{U}^{(3)} = 300`
  • :math:`T_{L}^{(3)} = 275`
  • :math:`f(T_{U}^{(3)}) = -17.9000`
  • :math:`f(T_{L}^{(3)}) = 9.7608`
  • :math:`\Delta H_{r}^{0}(T_{U}^{(3)}) = -23522.9`
  • :math:`\Delta H_{r}^{0}(T_{L}^{(3)}) = -23495.2392`
  • STEP 3: CHECK STOPPING CRITERION*

.. math::

\begin{array}{rl} \epsilon_{tol,x}^{(3)} &= \left| \frac{T_{M}^{(3)}-T_{M}^{(2)}}{T_{M}^{(2)}} \right| < \epsilon_{tol} ? \\ \epsilon_{tol,x}^{(3)} &= \left| \frac{\left(275\right)-\left(250\right)}{\left(250\right)} \right| < \left(10^{-6}\right) ? \\ \epsilon_{tol,x}^{(3)} &= 0.1 < \left(10^{-6}\right) ? \;\; \rightarrow \;\; \text{No, therefore keep going}. \\ & \\ \epsilon_{tol,f}^{(3)} &= \left| \frac{\Delta H_{r}^{0}\left(T_{M}^{(3)}\right)-\Delta H_{r}^{0}\left(T_{M}^{(2)}\right)}{\Delta H_{r}^{0}\left(T_{M}^{(2)}\right)} \right| < \epsilon_{tol} ? \\ \epsilon_{tol,f}^{(3)} &= \left| \frac{\left(-23495.2392\right)-\left(-23456.375\right)}{\left(-23456.375\right)} \right| < \left(10^{-6}\right) ? \\ \epsilon_{tol,f}^{(3)} &= 0.0017 < \left(10^{-6}\right) ? \;\; \rightarrow \;\; \text{No, therefore keep going}. \end{array}

    • ITERATION 4**
  • STEP 1: CALCULATE MID-POINT* :math:`T_{M}^{(4)}`

.. math::

\begin{array}{rl} T_{M}^{(4)} &= \frac{T_{U}^{(3)} + T_{L}^{(3)}}{2}\\ T_{M}^{(4)} &= \frac{\left(300\right) + \left(275\right)}{2}\\ T_{M}^{(4)} &= 287.5 \end{array}

.. math::

\begin{array}{rl} f\left(T_{M}^{(4)}\right) &= - 24097 - 0.26\left(T_{M}^{(4)}\right) + 1.69\text{x}10^{-3}\left(T_{M}^{(4)}\right)^{2} + \frac{1.5\text{x}10^{5}}{T_{M}^{(4)}} - \left(-23505\right)\\ f\left(T_{M}^{(4)}\right) &= - 24097 - 0.26\left(287.5\right) + 1.69\text{x}10^{-3}\left(287.5\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(287.5\right)} - \left(-23505\right)\\ f\left(T_{M}^{(4)}\right) &= -5.3218\\ & \\ \Delta H_{r}^{0}\left(T_{M}^{(4)}\right) &= - 24097 - 0.26\left(T_{M}^{(4)}\right) + 1.69\text{x}10^{-3}\left(T_{M}^{(4)}\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(T_{M}^{(4)}\right)}\\ \Delta H_{r}^{0}\left(T_{M}^{(4)}\right) &= - 24097 - 0.26\left(287.5\right) + 1.69\text{x}10^{-3}\left(287.5\right)^{2} + \frac{1.5\text{x}10^{5}}{\left(287.5\right)}\\ \Delta H_{r}^{0}\left(T_{M}^{(4)}\right) &= -23510.3218 \end{array}

  • STEP 2: CALCULATE* :math:`f(T_{L}^{(3)})\times f(T_{M}^{(4)})`

.. math::

\begin{array}{rl} Check &= f\left(T_{L}^{(3)}\right)\times f\left(T_{M}^{(4)}\right) < 0 ?\\ Check &= \left(9.7608\right)\times\left(-5.3218\right) < 0 ?\\ Check &= \left(-51.945\right) < 0 ? \;\; \rightarrow \;\; \text{Yes, therefore the root falls between }T_{L}^{(3)}\text{ and }T_{M}^{(4)}. \end{array}

  • :math:`T_{U}^{(4)} = 287.5`
  • :math:`T_{L}^{(4)} = 275`
  • :math:`f(T_{U}^{(4)}) = -5.3218`
  • :math:`f(T_{L}^{(4)}) = 9.7608`
  • :math:`\Delta H_{r}^{0}(T_{U}^{(4)}) = -23510.3218`
  • :math:`\Delta H_{r}^{0}(T_{L}^{(4)}) = -23495.2392`
  • STEP 3: CHECK STOPPING CRITERION*

.. math::

\begin{array}{rl} \epsilon_{tol,x}^{(4)} &= \left| \frac{T_{M}^{(4)}-T_{M}^{(3)}}{T_{M}^{(3)}} \right| < \epsilon_{tol} ? \\ \epsilon_{tol,x}^{(4)} &= \left| \frac{\left(287.5\right)-\left(275\right)}{\left(275\right)} \right| < \left(10^{-6}\right) ? \\ \epsilon_{tol,x}^{(4)} &= 0.0\bar{45} < \left(10^{-6}\right) ? \;\; \rightarrow \;\; \text{No, therefore keep going}. \\ & \\ \epsilon_{tol,f}^{(4)} &= \left| \frac{\Delta H_{r}^{0}\left(T_{M}^{(4)}\right)-\Delta H_{r}^{0}\left(T_{M}^{(3)}\right)}{\Delta H_{r}^{0}\left(T_{M}^{(3)}\right)} \right| < \epsilon_{tol} ? \\ \epsilon_{tol,f}^{(4)} &= \left| \frac{\left(-23510.3218\right)-\left(-23495.2392\right)}{\left(-23495.2392\right)} \right| < \left(10^{-6}\right) ? \\ \epsilon_{tol,f}^{(4)} &= 0.00064 < \left(10^{-6}\right) ? \;\; \rightarrow \;\; \text{No, therefore keep going}. \end{array} </rst>

Retrieved from "https://learnche.org/wiki_3E4/index.php?title=Tutorial_5_-_2010_-_Solution/Question_1&oldid=1086"