bisection.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Root finding algorithms\n",
    "\n",
    "## Bisection method"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are finding the roots of a continuous function $f: [a,b] \\rightarrow F$.\n",
    "\n",
    "The bisection method works by repeatedly dividing the interval (in which the root is known to lie) at midpoints $x_1 = a + \\frac{b-a}{2}$. \n",
    "\n",
    "At each step one determine the subinterval in which the function changes sign, here $[a,x_1]$ or $[x_1, b]$. This process is repeated until the root is found to within a certain level of accuracy (or tolerance).\n",
    "\n",
    "The bisection method is guaranteed to converge to the root of the function, but it is relatively slow and may not be the most efficient method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import *\n",
    "\n",
    "def bissection(f, a, b, tol):\n",
    "    niter = 0\n",
    "    zero = 0.0\n",
    "    res = 0.0\n",
    "    err = 0.0\n",
    "    inc = []\n",
    "    y = (a+b)/2\n",
    "    \n",
    "    if f(a)*f(b)<0:\n",
    "        while abs(b-a)>tol:\n",
    "            x=(a+b)/2\n",
    "            inc.append(abs(x-y))\n",
    "            y = x\n",
    "            niter += 1\n",
    "            if f(a)*f(x)<=0:\n",
    "                b = x\n",
    "            else:\n",
    "                a = x        \n",
    "        zero = (a+b)/2\n",
    "        res = f(zero)\n",
    "        err = abs(a-b)\n",
    "        return zero, res, niter, inc, err\n",
    "    elif f(a)*f(b)>0:\n",
    "        print(\"The bisection method cannot be applied to the interval [a,b]\")\n",
    "    else:\n",
    "        if f(a)==0:\n",
    "            print(f\"One root of f on [a,b] is {a=}\")\n",
    "        else:\n",
    "            print(f\"One root of f on [a,b] is {b=}\")\n",
    "    "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First tests"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us look for the root of the polynomial of degrees two $f(x) = x^2 - 10$.\n",
    "Knowing the solution, it is reasonable to choose as search interval $[a,b] = [3,4]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "a = 3\n",
    "b = 4\n",
    "tol = 10**(-6)\n",
    "def f(x):\n",
    "    return x**2-10\n",
    "\n",
    "zero, res, niter, inc, err = bissection(f, a, b, tol)\n",
    "print(f\"An estimate of a zero of the function x**2-10 : {zero  = }\")\n",
    "print(f\"The value of the function at the point obtained : {res = }\")\n",
    "print(f\"The number of iterations : {niter = }\")\n",
    "print(f\"The vector of residuals at each iteration : {inc = }\")\n",
    "print(f\"The length of the final bisection interval : {err = }\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now solve the equation \n",
    "\n",
    "$$\\frac{x}{2}-\\sin(x)+\\frac{\\pi}{6}-\\frac{\\sqrt{3}}{2}=0, \\quad x \\in [-\\frac{\\pi}{2}, \\pi ] $$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us set\n",
    "$$f(x)=\\frac{x}{2}-\\sin(x)+\\frac{\\pi}{6}-\\frac{\\sqrt{3}}{2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from matplotlib.pylab import *\n",
    "\n",
    "def F(x):\n",
    "    return x/2-sin(x)+pi/6-((3**(1/2))/2)\n",
    "\n",
    "x = linspace(-pi/2, pi, 100)\n",
    "y = F(x)\n",
    "plot(x,y)\n",
    "axhline(0, color='black')\n",
    "xlabel('x')\n",
    "ylabel('y')\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Nombre de zeros de f:** f possede 2 recines sur $[-\\frac{\\pi}{2}, \\pi ]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Zeros pour lesquelles on peut appliquer la dichotomie**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous observons qu'il y a un zero de f dans $[-\\frac{\\pi}{2}, 0 ]$ et un zero dans $[\\frac{\\pi}{2}, \\pi ]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "i. $\\forall x \\in [-\\frac{\\pi}{2}, 0 ], \\quad f(x) \\leq 0$. Par suite, $\\forall x,y \\in [-\\frac{\\pi}{2}, 0 ], \\quad f(x)f(y) \\geq 0$. Par consequent on ne peut appliquer la dichotomie pour ce zero"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ii. $f$ est continue sur $[\\frac{\\pi}{2}, \\pi ]$. De plus $f(\\frac{\\pi}{2})f(\\pi) < 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"f(pi/2)= \" + str(F(pi/2)))\n",
    "print(\"f(pi)= \" + str(F(pi)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Par suite, d'après le théorème des valeurs intermédiaires, f a au moins un zéro dans l’intervalle $[\\frac{\\pi}{2}, \\pi ]$ et on peut appliquer la methode de dichotomie"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Determinons le nombre minimal d'iterations possibles pour $tol=10^{-10}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "i. D'apres la question precedente, il n'est pas possible d'appliquer la methode de dichotomie sur $I_1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ii. $I_2 = [\\frac{\\pi}{2}, \\pi ], \\quad tol=10^{-10}.$ Soit $n$ le nombre d'itérations nécessaires pour satisfaire cette tolérance. Apres $n$ iterations, l'intervalle mesure $\\frac{b-a}{2^n}$. Par suite, $n$ est telle que $\\frac{b-a}{2^n} < tol$. La fonction $\\ln$ etant croissante:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "\\frac{b-a}{2^n} < tol &\\Rightarrow & \\ln (b-a) - \\ln (2^n) < \\ln (tol) \\\\\n",
    "&\\Rightarrow & n \\ln 2 > \\ln (b-a)-\\ln (tol) \\\\\n",
    "&\\Rightarrow & n > \\frac{\\ln (b-a)-\\ln (tol)}{\\ln 2} \\\\\n",
    "&\\Rightarrow & n > \\frac{\\ln (\\frac{\\pi}{2})+10\\ln (tol)}{\\ln 2} \\\\\n",
    "&\\Rightarrow & n > 33.870777078 \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**On conclut que le nombre minimal d'itérations nécessaires pour satisfaire la tolérance $tol=10^{-10}$ est $n=34$**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Verifions en utilisant la fonction bissection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import *\n",
    "\n",
    "a = pi/2\n",
    "b = pi\n",
    "tol = 10**(-10)\n",
    "\n",
    "zero, res, niter, inc, err = bissection(F, a, b, tol)\n",
    "\n",
    "print(\"Une estimation d’un zero de la function : zero  = \" + str(zero))\n",
    "print(\"Le reste de la fonction au point zero : f(zero) = \" + str(res))\n",
    "print(\"Le nombre d’itérations : niter = \" + str(niter))\n",
    "print(\"Un vecteur contenant les residus a chaque iteration : inc = \" + str(inc))\n",
    "print(\"La longueur de l'intervalle final : err = \" + str(err))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Nous avons bien $niter=34$**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Verifions que l'erreur commise satisfait la tolerance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = 2.246005589297974 #bonne approximation\n",
    "err = abs(x0-zero)\n",
    "print(err<tol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $x_0$ une bonne approximation du zero.\n",
    "On a $err=|x_0-zero|<tol$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**On conclut donc que l'erreur commise avec la solution numerique obtenue satisfait la tolerance**"
   ]
  }
 ],
 "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.10.4"
  },
  "vscode": {
   "interpreter": {
    "hash": "3ad933181bd8a04b432d3370b9dc3b0662ad032c4dfaa4e4f1596c548f763858"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}