{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fe293564",
   "metadata": {},
   "source": [
    "# Equations non-linéaires\n",
    "\n",
    "**Objectif :** trouver les zéros (ou racines) d'une fonction $f:[a,b] \\rightarrow \\mathbf R$ : \n",
    "$$\\alpha \\in [a,b] \\,:\\, f(\\alpha) = 0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61d4ed28",
   "metadata": {},
   "outputs": [],
   "source": [
    "# importing libraries used in this book\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "349ed263",
   "metadata": {},
   "outputs": [],
   "source": [
    "# defining the fonction\n",
    "def f(x):\n",
    "    return x*np.sin(x*2.*np.pi) + 0.5*x - 0.25\n",
    "\n",
    "[a,b] = [-2,2]\n",
    "\n",
    "# points used to plot the graph \n",
    "z = np.linspace(a, b, 100)\n",
    "\n",
    "plt.plot(z, f(z),'-')\n",
    "\n",
    "# labels, title, legend\n",
    "plt.xlabel('x'); plt.ylabel('$f(x)$'); #plt.title('data')\n",
    "plt.legend(['$f(t)$'])\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "563afd65",
   "metadata": {},
   "source": [
    "## Méthode de dichotomie ou bissection\n",
    "\n",
    "Si $f$ est continue et elle change de signe dans $[a,b]$, alors il existe au moins un $\\alpha$ tel que $f(\\alpha) = 0$.\n",
    "\n",
    "On peut alors définir l'algorithme suivant :\n",
    "\n",
    "$a^{(0)}=a$, $b^{(0)}=b$. Pour $k=0,1,...$\n",
    "\n",
    "1. $x^{(k)}=\\frac{a^{(k)}+b^{(k)}}{2}$\n",
    "2. si $f(x^{(k)})=0$, alors $x^{(k)}$ est le zéro cherché.\n",
    "\n",
    "    Autrement:\n",
    "\n",
    "    1. soit $f(x^{(k)})f(a^{(k)})<0$, alors le zéro\n",
    "    $\\alpha\\in[a^{(k)},x^{(k)}]$. \n",
    "\n",
    "    On pose $a^{(k+1)}=a^{(k)}$ et $b^{(k+1)}=x^{(k)}$\n",
    "\n",
    "    2. soit $f(x^{(k)}f(b^{(k)})<0$, alors le zéro\n",
    "    $\\alpha\\in[x^{(k)},b^{(k)}]$. \n",
    "\n",
    "    On pose $a^{(k+1)}=x^{(k)}$ et $b^{(k+1)}=b^{(k)}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6926310",
   "metadata": {},
   "source": [
    "### Exercice\n",
    "Comprenez et completez la fonction suivante qui effectue l'algorithme de dichotomie \n",
    "\n",
    "Ensuite testez-la pour trouver la racine de la fonction $f(x) = x\\sin(2\\pi x) + \\frac12 x - \\frac14$ dans l'intervalle $[-1.5,1]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef2d6b71",
   "metadata": {},
   "outputs": [],
   "source": [
    "def bisection(a,b,fun,tol,maxIterations) :\n",
    "    # [a,b] interval of interest\n",
    "    # fun function\n",
    "    # tolerance desired accuracy\n",
    "    # maxIterations : maximum number of iteration\n",
    "    # returns:\n",
    "    # zero, residual, number of iterations, x_sequence\n",
    "    # x_sequence : the sequence computed by the fixed point iterations\n",
    "    \n",
    "    if (a >= b) :\n",
    "        print(' b must be greater than a (b > a)')\n",
    "        return 0,0,0,[0]\n",
    "    \n",
    "    # what we consider as \"zero\"\n",
    "    eps = 1e-12\n",
    "    \n",
    "    # evaluate f at the endpoints\n",
    "    fa = fun(a)\n",
    "    fb = fun(b)\n",
    "\n",
    "    if abs(fa) <  eps : # a is the solution\n",
    "        zero = a\n",
    "        esterr = fa\n",
    "        k = 0\n",
    "        return zero, esterr, k, [zero]\n",
    "    \n",
    "    if abs(fb) <  eps : # b is the solution\n",
    "        zero = b\n",
    "        esterr = fb\n",
    "        k = 0\n",
    "        return zero, esterr, k, [zero]\n",
    "\n",
    "    if fa*fb > 0 :\n",
    "        print(' The sign of FUN at the extrema of the interval must be different')\n",
    "        return 0,0,0,[0]\n",
    "\n",
    "\n",
    "\n",
    "    # We want the final error to be smaller than tol, \n",
    "    # i.e. k > log( (b-a)/tol ) / log(2) - 1\n",
    "\n",
    "    nmax = int(np.ceil(np.log( (b-a)/tol ) / np.log(2))) - 1\n",
    "\n",
    "    \n",
    "    # but nmax shall be smaller the the nmaximum iterations asked by the user\n",
    "    if ( maxIterations < nmax ) :\n",
    "        nmax = int(round(maxIterations))\n",
    "        print('Warning: maxIterations is smaller than the minimum number of iterations necessary to reach the tolerance wished');\n",
    "\n",
    "    # vector of intermadiate approximations etc\n",
    "    x = np.zeros(nmax+1)\n",
    "\n",
    "    # initial error is the length of the interval.\n",
    "    esterr  = (b - a)\n",
    "\n",
    "    # do not need to store all the a^k and b^k, so I call them with a new variable name:\n",
    "    ak = a\n",
    "    bk = b\n",
    "    # the values of f at those points are fa and fk\n",
    "\n",
    "    for k in range(nmax+1) :\n",
    "\n",
    "        # approximate solution is midpoint of current interval\n",
    "        # COMPLETE the code below\n",
    "        #> x[k] =  \n",
    "        #> fx = \n",
    "        \n",
    "        # error estimator is the half of the previous error\n",
    "        #> esterr  =\n",
    "\n",
    "        # COMPLETE the code below\n",
    "\n",
    "        # if we found the solution, stop the algorithm\n",
    "        if np.abs(fx) <  eps :\n",
    "            # error is zero\n",
    "            #> zero = \n",
    "            #> esterr = \n",
    "            #> return  \n",
    "\n",
    "        if fx*fa < 0  :  # alpha is in (a,x)\n",
    "            #> bk = \n",
    "            #> fb = \n",
    "        elif fx*fb < 0 : # alpha is in (x,b)\n",
    "            #> ak =  \n",
    "            #> fa = \n",
    "        else  :\n",
    "            error('Algorithm not operating correctly')\n",
    "\n",
    "\n",
    "\n",
    "    zero = x[k];\n",
    "\n",
    "    if esterr > tol :\n",
    "        print('Warning: bisection stopped without converging to the desired tolerance because the maximum number of iterations was reached');\n",
    "\n",
    "    return zero, esterr, k, x \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23fd5adb",
   "metadata": {},
   "outputs": [],
   "source": [
    "[a,b] = [-1.5,1]\n",
    "tol = 1e-10\n",
    "maxIter = 4\n",
    "zero, esterr, k, x = bisection(a,b,f,tol,maxIter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d60e8132",
   "metadata": {},
   "outputs": [],
   "source": [
    "from NonLinearEquationsLib import plotBisectionIterations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32ba3753",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'font.size': 16})\n",
    "plt.figure(figsize=(8, 5))\n",
    "\n",
    "plotBisectionIterations(a,b,f,x)\n",
    "\n",
    "# plt.savefig('Bisection-iterations.png', dpi=600)\n",
    "\n",
    "print(f'The estimated root is {zero:2.12f}, the estimated error {esterr:1.3e} and the residual is {f(zero):1.3e}, after {k} iterations')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  },
  "unianalytics_cell_mapping": [
   [
    "fe293564",
    "fe293564"
   ],
   [
    "61d4ed28",
    "61d4ed28"
   ],
   [
    "349ed263",
    "349ed263"
   ],
   [
    "563afd65",
    "563afd65"
   ],
   [
    "b6926310",
    "b6926310"
   ],
   [
    "ef2d6b71",
    "ef2d6b71"
   ],
   [
    "23fd5adb",
    "23fd5adb"
   ],
   [
    "d60e8132",
    "d60e8132"
   ],
   [
    "32ba3753",
    "32ba3753"
   ]
  ],
  "unianalytics_notebook_id": "349ecc8e-3236-4d72-8858-c6f20ced3308"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}