method_of_bisection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "increasing-round",
   "metadata": {},
   "source": [
    "# Guided Programming 3 - Method of Bisection \n",
    "\n",
    "**Mathematical Methods for Chemical Engineers (MTHS1008)**\n",
    "\n",
    "**Dr Matthew Scase**\n",
    "\n",
    "We have seen in Guided Programming 1 and Guided Programming 2 how to solve two types of equations that we are familiar with, quadratic equations and simultaneous equations.  But we already knew how to do that!\n",
    "\n",
    "In this notebook we will revisit the problem of solving a quadratic equation but we will use a method that means we will be able to **solve *any* equation** of the form\n",
    "\n",
    "$$f(x) = 0.$$\n",
    "\n",
    "This form of equation includes the quadratic equation we looked at previously\n",
    "\n",
    "\\begin{equation*}\n",
    "2x^2 - 7x + 6 = 0,\n",
    "\\label{eq:quadratic} \\tag{1}\n",
    "\\end{equation*}\n",
    "\n",
    "(where $f(x) = 2x^2 - 7x + 6$) but the form $f(x) = 0$ also includes equations that are similar but not exactly quadratic.  We can well imagine that the equation\n",
    "\n",
    "$$2x^{1.9} - 7x + 6 = 0,$$\n",
    "\n",
    "has solutions that are close to the solutions of \\eqref{eq:quadratic}, but we cannot use the quadratic formula to find these solutions anymore.  More generally, even simply expressed mathematical equations, such as \n",
    "\n",
    "$$\\cos x = x,$$\n",
    "\n",
    "do not have simple solutions of the form we are used to being able to find at school.  However, we can define, in this case, $f(x) = \\cos x - x$ and look for a solution of $f(x)=0$ using the methods we shall implement below.\n",
    "\n",
    "In practical engineering applications we are frequently confronted with equations of the form $f(x) = 0$ that we need to solve, but no formula for their solution exists.  In this notebook we will implement a strategy for solving these types of equation known as the [**method of bisection**](https://en.wikipedia.org/wiki/Bisection_method).\n",
    "\n",
    "## The method of bisection\n",
    "\n",
    "We are seeking the special value of $x$ that satisfies $f(x)=0$.  Our strategy is to trap this special value of $x$ between to endpoints, let's call them $x_L$ and $x_R$, and then systematically move our endpoints closer and closer to each other *ensuring that $x$ remains trapped between the endpoints*.  When our endpoints are sufficiently close, say within $10^{-6}$ of each other, then we have a numerical approximation to the solution of our equation that is accurate to 6 decimal places and we can stop.\n",
    "\n",
    "### Trapping the solution\n",
    "\n",
    "We need to provide two endpoints to use our method; $x_L$ and $x_R$, where $x_L$ is less than the true solution $x$, and $x_R$ is greater than the true solution $x$.\n",
    "\n",
    "<div style=\"text-align:center;\">\n",
    "<img src=\"method_of_bisection_figure_1.png\" alt=\"figure\" width = \"460\"/>\n",
    "    <strong><center>Figure 1</center></strong>\n",
    "</div>\n",
    "\n",
    "\n",
    "Figure 1 above shows the curve $y = f(x)$.  The solution of our problem, or desired \"root\", is the special value of $x$ for which $f(x) = 0$, this is indicated by the blue circle on $y = 0$.  The endpoint $x_L$ is less than the root and we can see that $f(x_L)<0$;  $x_L$ is a lower estimate for the desired root.  The endpoint $x_R$ is greater than the root and we can see that $f(x_R)>0$; $x_R$ is an upper estimate for the desired root.  It is *essential* that the root lies inbetween the two endpoints.\n",
    "\n",
    "For the method to work we must ensure that $f(x_L)$ and $f(x_R)$ have opposite signs.  It does not matter whether $f(x_L)<0$ and $f(x_R)>0$, as shown in Figure 1, or $f(x_L)>0$ and $f(x_R)<0$.  \n",
    "\n",
    "Let us assume that we have chosen two suitable values of $x_L$ and $x_R$ so that we can begin our method.\n",
    "\n",
    "<div style=\"text-align:center;\">\n",
    "<img src=\"method_of_bisection_figure_2.png\" alt=\"figure\" width = \"460\"/>\n",
    "    <strong><center>Figure 2</center></strong>\n",
    "</div>\n",
    "\n",
    "- **Step 1**: Calculate the midpoint\n",
    "\n",
    "     We calculate the *midpoint* between $x_L$ and $x_R$ using\n",
    "    $$x_M = \\frac{x_L + x_R}{2}.$$\n",
    "    The midpoint is exactly halfway between $x_L$ and $x_R$, bisecting our endpoints - hence the name of the method.\n",
    "\n",
    "\n",
    "- **Step 2**: Use the midpoint to narrow the distance between the endpoints.\n",
    "\n",
    "    The midpoint of the values of $x_L$ and $x_R$ in Figure 1 has been calculated and plotted in Figure 2.  We can immediately see that the midpoint value $x_M$ is a much better *lower estimate* of the desired root than $x_L$.  So we can redefine our lower limit $x_L$ as being equal to the midpoint value $x_M$.\n",
    "\n",
    "    Observing that $x_M$ is a much better lower estimate than $x_L$ is straightforward from looking at Figure 2 - but how do we get Python to work out that it is the lower estimate that is improved and not the upper estimate?  If we replace the upper limit $x_R$ with the midpoint then we would no-longer be trapping our solution and the method would fail.\n",
    "\n",
    "    We see that it is the lower limit that gets updated in this case since $f(x_M)$ has the *same sign* as $f(x_L)$.  This means that both $x_M$ and $x_L$ are the same side of the desired root.  \n",
    "    \n",
    "    So for the illustrated case we update the lower endpoint, $x_L$, setting it equal to the midpoint $x_M$.\n",
    "    \n",
    "<div style=\"text-align:center;\">\n",
    "<img src=\"method_of_bisection_figure_3.png\" alt=\"figure\" width = \"460\"/>\n",
    "    <strong><center>Figure 3</center></strong>\n",
    "</div>\n",
    "\n",
    "Now we go back to Step 1 and recalculate the midpoint based on our new endpoint.  This is shown in Figure 3.  The original lower estimate is shown in light blue, our new updated lower estimate is at $x_L$ and we have a new midpoint, $x_M$.\n",
    "\n",
    "This time we see that $f(x_M)>0$ meaning that for the case shown $f(x_L)$ and $f(x_M)$ have opposite signs and that therefore the desired root is trapped between $x_L$ and $x_M$.  We have that $x_M$ is an improved upper estimate of the desired root and so in Step 2 we update $x_R$ by setting $x_R = x_M$.\n",
    "\n",
    "<div style=\"text-align:center;\">\n",
    "<img src=\"method_of_bisection_figure_4.png\" alt=\"figure\" width = \"460\"/>\n",
    "    <strong><center>Figure 4</center></strong>\n",
    "</div>\n",
    "\n",
    "As can be seen from Figure 4, when we go back to Step 1 and recalculate the midpoint we now have a very good estimate of our desired root.\n",
    "\n",
    "We can keep iterating around Steps 1 and 2 until our numerical solution is accurate enough for our application.  We know that our root is always trapped between $x_L$ and $x_R$, so if the distance between $x_L$ and $x_R$ is less than a given tolerance that we can specify then we know our solution is accurate.\n",
    "        \n",
    "        \n",
    "## Implementing the method in Python\n",
    "\n",
    "It is **very important** when writing some code to solve a problem that you start on a problem you already know the answer to.  Let's start with \n",
    "\n",
    "\\begin{equation*}\n",
    "    f(x) = 2x^2 - 7x + 6\n",
    "    \\label{eq:f1} \\tag{2}\n",
    "\\end{equation*}\n",
    "\n",
    "and look for solutions $x$ that satisfy $f(x) = 0$.  We know that $x=2$ and $x = \\frac{3}{2}$ satisfy \\eqref{eq:f1} so if we begin our method with $x_L = 1.75$ and $x_R = 3$ we should be trapping the solution $x=2$ between our two endpoints.  Let us set our initial endpoint estimates `XL` and `XR`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "scenic-pointer",
   "metadata": {},
   "outputs": [],
   "source": [
    "XL = 1.75 # ....... This is our initial lower estimate of the actual root\n",
    "XR = 3 # .......... This is our initial upper estimate of the actual root"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "radical-carnival",
   "metadata": {},
   "source": [
    "We now need to iterate over the steps in our bisection algorithm.  We can do this using a `for` loop.  Let us start with a small number of iterations to check that our method is working, we can always increase the number of iterations later if we would like a more accurate solution.\n",
    "\n",
    "        for i in range(5):\n",
    "        \n",
    "- First we will calculate the midpoint, let's call it `XM` between the current lower and upper endpoint estimates `XL` and `XR`\n",
    "\n",
    "            XM = (XL + XR)/2\n",
    "        \n",
    "- Next we will need to calculate $y_L = f(x_L)$ and $y_M = f(x_M)$ so that we can decide whether it is the upper or lower estimate that needs updating.  \n",
    "\n",
    "            YL = 2*XL**2 - 7*XL + 6\n",
    "            YM = 2*XM**2 - 7*XM + 6\n",
    "         \n",
    "- Now we are able to decide whether to update the lower estimate or the upper estimate by calculating the product `YL*YM`.  **If** `YL*YM` is greater than 0 then the sign of $f(x_L)$ is equal to the sign of $f(x_M)$ and we update the lower estimate `XL` to be equal to the midpoint `XM`, **else** we update the upper estimate `XR` to be equal to the midpoint `XM`.\n",
    "\n",
    "            if YL*YM > 0:\n",
    "                XL = XM\n",
    "            else:\n",
    "                XR = XM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "signal-preparation",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(5): # Carry out algorithm 5 times\n",
    "    \n",
    "    XM = (XL + XR)/2 # ......... Calculate the midpoint between XL and XR\n",
    "    YL = 2*XL**2 - 7*XL + 6 # .. Evaluate f(XL)\n",
    "    YM = 2*XM**2 - 7*XM + 6 # .. Evaluate f(XM)\n",
    "    \n",
    "    if YL*YM > 0:\n",
    "        XL = XM # .............. If f(XL) has the same sign as f(XM) update XL\n",
    "    else:\n",
    "        XR = XM # .............. Otherwise update XR"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "grave-basic",
   "metadata": {},
   "source": [
    "Now we can write out our solution to the screen, the current midpoint `XM` is our best estimate of the desired root."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "animal-initial",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = 2.0234375 f(x) = 0.0245361328125\n"
     ]
    }
   ],
   "source": [
    "print('x =', XM, 'f(x) =', YM)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "diagnostic-quarter",
   "metadata": {},
   "source": [
    "The correct answer becomes more and more accurate the greater the number of iterations used.\n",
    "\n",
    "To solve the equation\n",
    "$$\\cos x = x$$\n",
    "we set $f(x) = \\cos x - x$. We need to `import` the cosine function into Python using the Numerical Python package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "hourly-dryer",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "right-fishing",
   "metadata": {},
   "source": [
    "We can now use special functions such as cosine, sine, tan, exp, $\\ldots$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "patient-network",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cos(0) = 1.0\n",
      "exp(1) = 2.718281828459045\n"
     ]
    }
   ],
   "source": [
    "print('cos(0) =', numpy.cos(0))\n",
    "print('exp(1) =', numpy.exp(1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "satisfactory-bleeding",
   "metadata": {},
   "source": [
    "Let us use our bisection method algorithm with $x_L = 0$ and $x_R = 1$ to find the solution of $\\cos x = x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "surprising-criminal",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution: x = 0.7390851080417633  Error: f(x) = 4.213050042167765e-08\n"
     ]
    }
   ],
   "source": [
    "XL = 0 # ....... This is our initial lower estimate of the actual root\n",
    "XR = 1 # ....... This is our initial upper estimate of the actual root\n",
    "\n",
    "for i in range(25): # Carry out algorithm 25 times\n",
    "    \n",
    "    XM = (XL + XR)/2 # ........... Calculate the midpoint between XL and XR\n",
    "    YL = numpy.cos(XL) - XL # .. Evaluate f(XL)\n",
    "    YM = numpy.cos(XM) - XM # .. Evaluate f(XM)\n",
    "    \n",
    "    if YL*YM > 0:\n",
    "        XL = XM # .............. If f(XL) has the same sign as f(XM) update XL\n",
    "    else:\n",
    "        XR = XM # .............. Otherwise update XR\n",
    "\n",
    "print('Solution: x =', XM, ' Error: f(x) =', YM)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adult-order",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "- Try modifying the code above to solve the near quadratic\n",
    "   $$2x^{1.9} - 7x + 6 = 0,$$\n",
    "  using the same initial guess as above.  \n",
    "\n",
    "    [**Hint**: You should find the solution is approximately $x = 2.553905$]\n",
    "\n",
    "\n",
    "- Find an numerical approximation to the square root of 2 by solving $x^2 - 2 = 0$.\n",
    "\n",
    "\n",
    "- Find a numerical solution to $2x^3 - 3x^2 + 4x - 6 = 0$ and hence factorize the cubic.  A cubic equation *always* has at least one real root so if you do not find a solution initially try moving your endpoints further apart.  \n",
    "    [Answer: $(2x-3)(x^2+2)$.]\n",
    "\n",
    "  \n",
    "- Add some code that ensures that we have trapped the root with our initial choices of `XL` and `XR`.\n",
    "\n",
    "\n",
    "- Consider the output of your script when the midpoint is an exact solution of your equation - is this what you would like the code to do?\n",
    "\n",
    "\n",
    "- What happens when you try and solve the equation $\\frac{1}{x} = 0$ with endpoints $x_L < 0$, $x_R > 0$?  The code thinks that $x=0$ is a solution! Disaster!  What has happened?  Write some code that lets the user know that the found solution is not a good solution."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}