{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module root_finding\n",
    "\n",
    "Module generated from eponymous notebook.\n",
    "\n",
    "**Author: Brenton LeMesurier**, brenton.lemesurier@unco.edu and lemesurierb@cofc.edu\n",
    "\n",
    "**Last Revised: February 25, 2021**\n",
    "\n",
    "A collection of functions for rootfinding, as used for example in MATH375 Assignments 2 and 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bisection(f, a, b, errorTolerance, maxIterations=100, demoMode=False):\n",
    "    \"\"\"Solve f(x) = 0 to within absolute error errorTolerance by the Bisection Method,\n",
    "    but avoiding infinite loops by giving up after maxIterations iterations.\n",
    "    Note: There is an optional input parameter \"demoMode\" which controls whether to\n",
    "    - print intermediate results (for \"study\" purposes), or to\n",
    "    - work silently (for \"production\" use).\n",
    "    \"\"\"\n",
    "    functionEvaluations = 0\n",
    "    # Whenever a  (resp. c) is updated, also update fa (resp. fc).\n",
    "    fa = f(a)\n",
    "    functionEvaluations += 1\n",
    "    for iteration in range(1, maxIterations):\n",
    "        c = (a + b)/2\n",
    "        fc = f(c)\n",
    "        functionEvaluations += 1\n",
    "        if fa * fc < 0:\n",
    "            b = c\n",
    "        else:\n",
    "            a = c\n",
    "            fa = fc  # Copy rather than redundantly evaluating f(a) for the new a, which is the old c\n",
    "        errorBound = (b-a)/2\n",
    "        if demoMode: print(f\"{iteration=}, approximate root {c}, {errorBound=:0.3}\")\n",
    "        if errorBound <= errorTolerance: break\n",
    "    root = c             \n",
    "    return (root, errorBound, iteration, functionEvaluations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def newton(f, Df, x0, errorTolerance, maxIterations=100, demoMode=False):\n",
    "    \"\"\"\n",
    "    Solve f(x) = 0 to within absolute error errorTolerance by Newton's Method\n",
    "    but avoiding infinite loops by giving up after maxIterations iterations.\n",
    "    Note: There is an optional input parameter \"demoMode\" which controls whether to\n",
    "    - print intermediate results (for \"study\" purposes), or to\n",
    "    - work silently (for \"production\" use).\n",
    "    The default is silence.\n",
    "    \"\"\"\n",
    "    functionEvaluations = 0\n",
    "    x = x0\n",
    "    for iteration in range(1, maxIterations + 1):\n",
    "        fx = f(x)\n",
    "        Dfx = Df(x)\n",
    "        functionEvaluations += 2\n",
    "        # Note: a more careful, robust code would check for the possibility of division by zero here.\n",
    "        dx = fx/Dfx\n",
    "        x -= dx  # Aside: this is shorthand for \"x = x - dx\"\n",
    "        errorEstimate = abs(dx)\n",
    "        if demoMode: print(f\"{iteration=}, {x=}, {errorEstimate=:0.3}\")\n",
    "        if errorEstimate <= errorTolerance: break\n",
    "    return (x, errorEstimate, iteration, functionEvaluations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fpi(g, x_0, errorTolerance=1e-6, maxIterations=20, demoMode=False):\n",
    "    \"\"\"Fixed point iteration for approximately solving x = f(x),\n",
    "    x_0: the initial value\n",
    "    \"\"\"\n",
    "    x = x_0\n",
    "    for iteration in range(maxIterations):\n",
    "        x_new = g(x)\n",
    "        errorEstimate = np.abs(x_new - x)\n",
    "        x = x_new\n",
    "        if demoMode: print(f\"x_{iteration} = {x}, {errorEstimate=:0.3}\")\n",
    "        if errorEstimate <= errorTolerance: break\n",
    "    return (x, errorEstimate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def false_position(f, a, b, errorTolerance=1e-15, maxIterations=15, demoMode=False):\n",
    "    \"\"\"Solve f(x)=0 in the interval [a, b] by the Method of False Position.\n",
    "    This code also illustrates a few ideas that I encourage, such as:\n",
    "    - Avoiding infinite loops, by using for loops sand break\n",
    "    - Avoiding repeated evaluation of the same quantity\n",
    "    - Use of descriptive variable names\n",
    "    - Use of \"camelCase\" to turn descriptive phrases into valid Python variable names\n",
    "    - An optional \"demonstration mode\" to display intermediate results.\n",
    "    \"\"\"\n",
    "    fa = f(a)\n",
    "    fb = f(b)\n",
    "    functionEvaluations = 2\n",
    "    for iteration in range(1, maxIterations+1):\n",
    "        if demoMode: print(f\"\\nIteration {iteration}:\")\n",
    "        c = (a * fb - fa * b)/(fb - fa)\n",
    "        fc = f(c)\n",
    "        functionEvaluations += 1\n",
    "        if fa * fc < 0:\n",
    "            b = c\n",
    "            fb = fc  # N.B. When b is updated, so must be fb = f(b)\n",
    "        else:\n",
    "            a = c\n",
    "            fa = fc\n",
    "        errorBound = b - a\n",
    "        if demoMode:\n",
    "            print(f\"The root is in interval [{a}, {b}]\")\n",
    "            print(f\"The new approximation is {c}, with error bound {errorBound:0.3}, backward error {abs(fc):0.3}\")\n",
    "        if errorBound < errorTolerance: break\n",
    "    # Whether we got here due to accuracy of running out of iterations,\n",
    "    # return the information we have, including an error bound:\n",
    "    root = c  # the newest value is probably the most accurate\n",
    "    iterations = iteration\n",
    "    return (root, errorBound, iterations, functionEvaluations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def secant(f, a, b, errorTolerance=1e-15, maxIterations=15, demoMode=False):\n",
    "    \"\"\"Solve f(x)=0 in the interval [a, b] by the Secant Method.\"\"\"\n",
    "    # Some more descriptive names\n",
    "    x_older = a\n",
    "    x_more_recent = b\n",
    "    f_x_older = f(x_older)\n",
    "    f_x_more_recent = f(x_more_recent)\n",
    "    functionEvaluations = 2\n",
    "    for iteration in range(1, maxIterations+1):\n",
    "        x_new = (x_older * f_x_more_recent - f_x_older * x_more_recent)/(f_x_more_recent - f_x_older)\n",
    "        f_x_new = f(x_new)\n",
    "        functionEvaluations += 1\n",
    "        (x_older, x_more_recent) = (x_more_recent, x_new)\n",
    "        (f_x_older, f_x_more_recent) = (f_x_more_recent, f_x_new)\n",
    "        errorEstimate = abs(x_older - x_more_recent)\n",
    "        if demoMode: print(f\"{iteration=}, x={x_more_recent}, {errorEstimate=:0.3}\")\n",
    "        if errorEstimate < errorTolerance: break\n",
    "    # Whether we got here due to accuracy of running out of iterations,\n",
    "    # return the information we have, including an error estimate:\n",
    "    root = x_new\n",
    "    iterations = iteration\n",
    "    return (root, errorEstimate, iterations, functionEvaluations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "This work is licensed under [Creative Commons Attribution-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-sa/4.0/)"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}