{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Root-finding solver based on Brent's method\n",
    "\n",
    "Now we start with our [Brent root-finding algorithm](https://en.wikipedia.org/wiki/Brent%27s_method)\n",
    "example code.\n",
    "\n",
    "We start with untyped, script-like code, and modularize it bit by bit.\n",
    "\n",
    "This document is a Jupyter notebook! Thus it shows how to use \n",
    "[MyST-NB](https://myst-nb.readthedocs.io/en/latest/index.html) to include notebooks in a Sphinx\n",
    "site."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 3.0, -0.9899924966004454, 3.0, 1.9899924966004454, None)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from math import cos\n",
    "\n",
    "x_rel_tol = 1e-12  # X relative tolerance\n",
    "x_abs_tol = 1e-12  # X absolute tolerance\n",
    "y_tol = 1e-12  # Y tolerance\n",
    "verbose = True  # Whether to display progress\n",
    "\n",
    "f = cos\n",
    "a = 0.0  # first point that brackets the root\n",
    "b = 3.0  # second point that brackets the root\n",
    "fa = f(a)  # f(a)\n",
    "fb = f(b)  # f(b)\n",
    "\n",
    "if abs(fa) < abs(fb):\n",
    "    # force abs(fa) >= abs(fb), make sure that b is the best root approximation known so far\n",
    "    # and a is the contrapoint\n",
    "    b, a = a, b\n",
    "    fb, fa = fa, fb\n",
    "\n",
    "c = a  # Previous iterate\n",
    "fc = fa  # f(c)\n",
    "d = a  # Iterate before the previous iterate\n",
    "fd = fa  # f(d)\n",
    "last_step = None\n",
    "step = None\n",
    "iter = 1  # Current iteration number (1-based)\n",
    "\n",
    "# display the initial state\n",
    "(iter, b, fb, abs(a-b), abs(fa-fb), last_step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\t1.500e+00\t7.074e-02\t1.500e+00\t1.061e+00\tNone\n",
      "3\t1.600e+00\t-2.923e-02\t1.000e-01\t9.997e-02\tbisection\n",
      "4\t1.571e+00\t1.434e-05\t2.925e-02\t2.924e-02\tsecant\n",
      "5\t1.571e+00\t-2.042e-09\t1.434e-05\t1.434e-05\tsecant\n",
      "Met y_tol criterion\n",
      "The approximate root is 1.5707963267948966\n"
     ]
    }
   ],
   "source": [
    "while True: # loop forever\n",
    "    iter = iter + 1\n",
    "    last_step = step\n",
    "    if fa != fc and fb != fc:\n",
    "        # perform a quadratic interpolation step\n",
    "        # https://en.wikipedia.org/wiki/Inverse_quadratic_interpolation\n",
    "        L0 = (a * fb * fc) / ((fa - fb) * (fa - fc))\n",
    "        L1 = (b * fa * fc) / ((fb - fa) * (fb - fc))\n",
    "        L2 = (c * fb * fa) / ((fc - fa) * (fc - fb))\n",
    "        s = L0 + L1 + L2\n",
    "        step = \"quadratic\"\n",
    "    else:\n",
    "        # perform a secant step\n",
    "        s = b - fb * (b - a) / (fb - fa)\n",
    "        step = \"secant\"\n",
    "    perform_bisection = False\n",
    "    if a <= b and not ((3 * a + b) / 4 <= s <= b):\n",
    "        perform_bisection = True\n",
    "    elif b <= a and not (b <= s <= (3 * a + b) / 4):\n",
    "        perform_bisection = True\n",
    "    elif last_step == \"bisection\" and abs(s - b) >= abs(b - c) / 2:\n",
    "        perform_bisection = True\n",
    "    elif last_step != \"bisection\" and abs(a - b) >= abs(c - d) / 2:\n",
    "        perform_bisection = True\n",
    "    elif last_step == \"bisection\" and abs(b - c) < x_abs_tol:\n",
    "        perform_bisection = True\n",
    "    elif last_step != \"bisection\" and abs(c - d) < x_abs_tol:\n",
    "        perform_bisection = True\n",
    "    if perform_bisection:\n",
    "        # perform a bisection step\n",
    "        s = min(a, b) + abs(b - a) / 2\n",
    "        step = \"bisection\"\n",
    "    fs = f(s)\n",
    "    d = c\n",
    "    c = b\n",
    "    fc = fb\n",
    "    # check which point to replace to maintain (a,b) have different signs\n",
    "    if f(a) * f(s) < 0:\n",
    "        b = s\n",
    "        fb = fs\n",
    "    else:\n",
    "        a = s\n",
    "        fa = fs\n",
    "    # keep b as the best guess\n",
    "    if abs(fa) < abs(fb):\n",
    "        b, a = a, b\n",
    "        fb, fa = fa, fb\n",
    "\n",
    "    # Checks convergence\n",
    "    if fb == 0:\n",
    "        print(\"Exact root found\")\n",
    "        break\n",
    "    x_delta = abs(a - b)\n",
    "    if x_delta <= x_abs_tol:\n",
    "        print(\"Met x_abs_tol criterion\")\n",
    "        break\n",
    "    if x_delta / max(a, b) <= x_rel_tol:\n",
    "        print(\"Met x_rel_tol criterion\")\n",
    "        break\n",
    "    y_delta = abs(fb)\n",
    "    if y_delta <= y_tol:\n",
    "        print(\"Met y_tol criterion\")\n",
    "        break\n",
    "\n",
    "    dx = abs(a - b)\n",
    "    dy = abs(fa - fb)\n",
    "    print(f\"{iter}\\t{b:.3e}\\t{fb:.3e}\\t{dx:.3e}\\t{dy:.3e}\\t{last_step}\")\n",
    "\n",
    "print(f\"The approximate root is {b}\")"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "56a39a9d3534467ed83694461484cafc606d79ce7db33d76d260741ac1a4cb13"
  },
  "kernelspec": {
   "display_name": "Python 3.9.12 ('.venv': poetry)",
   "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.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}