Root-finding algorithms.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Root-Finding Algorithms\n",
    "_By Dhruv Jain_\n",
    "\n",
    "### **Objective: Find a value x* for a function f(x), such that f(x*) = 0**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Key libraries: Numpy(for mathematical procedures) and matplotlib(to create plots)\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bisection Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bisection Method Function\n",
    "def bisection(func, ll, ul, tol=1e-10, divtol=10, Nmax=50):\n",
    "    \"\"\"\"Dhruv Jain, 9 Sept, 2021\n",
    "    Obj: Apply Bisection Method to find a root of a function in a given interval for a certain tolerance and/or a certain number of iterations\n",
    "    Args:\n",
    "        func: function\n",
    "            Evaluation function on which the bisection method is applied\n",
    "        ll: float\n",
    "            Lower limit of the interval\n",
    "        ul: float\n",
    "            Upper limit of the interval\n",
    "        tol: float\n",
    "            Convergence tolerance value (used as stopping condition)\n",
    "        divtol: float \n",
    "            Divergence tolerance value (used as stopping condition)\n",
    "        Nmax: int\n",
    "            # Allowed iterations (used as stopping condition)\n",
    "            \n",
    "    Output:\n",
    "        midpt: numpy ndarray, float \n",
    "            Stores midpoint value calculated after each iteration, last element is {x*: f(x*)=0} basedd on # iteration and/or tolerance\n",
    "        count: int \n",
    "            # Iterations used to compute the midpoint\n",
    "    \"\"\"\n",
    "    \n",
    "    ful = func(ul) # f(x) of lower limit of the interval\n",
    "    fll = func(ll) # f(x) of upper limit of the interval\n",
    "    \n",
    "    if ful*fll > 0: # check if zeros exist by checking if the sign of ful and fll is opposite, if none exist then exit from the function\n",
    "        print('No zero of the function exists in the given interval')\n",
    "        return None, None\n",
    "    \n",
    "    midpt = np.zeros(Nmax)\n",
    "    midpt.fill(np.nan)\n",
    "    count = 0\n",
    "    \n",
    "    # Make sure valid Nmax\n",
    "    if count >= Nmax:\n",
    "        print('Enter Nmax more than 0')\n",
    "        return None, None\n",
    "    \n",
    "    # Stopping Condition: #iterations<Nmax, Interal below convergence tolerance or above divergence tolerance\n",
    "    while count < Nmax and np.abs(ul-ll) > tol and np.abs(ul-ll) < divtol:\n",
    "        \n",
    "        midpt[count] = ll + (ul-ll)/2 # mid point of the interval\n",
    "        fmidpt = func(midpt[count]) # f(x) @ mid point\n",
    "        \n",
    "        #Update limit of interval\n",
    "        if fmidpt*fll < 0:\n",
    "            ul = midpt[count]\n",
    "        else:\n",
    "            ll = midpt[count]\n",
    "        count = count + 1\n",
    "        \n",
    "        if count == Nmax:\n",
    "            print('The maximum number of iterations :',Nmax,'has been met')\n",
    "\n",
    "    midptcalc = midpt[:count]\n",
    "    \n",
    "    # Everything 0 when exceeded divergence tolerance\n",
    "    if np.abs(ul-ll) >= divtol:\n",
    "        midptcalc = 0\n",
    "        count = 0\n",
    "    \n",
    "    return midptcalc, count"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Newton's Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Newton's Method\n",
    "def newton(func, dfunc, x0, tol=1e-12,divtol=1e2, Nmax=50): \n",
    "    \"\"\"Dhruv Jain, 7th Sept 2021   \n",
    "    Obj: Apply Newton's Method to obtain the zero of the funcion using an inital guess for a certain tolerance and/or # iterations \n",
    "    Args:\n",
    "        func: function\n",
    "            Evaluation function on which the Newton's method is applied\n",
    "        dfunc: function \n",
    "            Derivative of the func\n",
    "        x0: float\n",
    "            Initial guess\n",
    "        tol: float\n",
    "            Convergence tolerance value (used as stopping condition)\n",
    "        divtol: float \n",
    "            Divergence tolerance value (used as stopping condition)\n",
    "        Nmax: int\n",
    "            # Allowed iterations (used as stopping condition)\n",
    "            \n",
    "    Output:\n",
    "        root: numpy ndarray, float \n",
    "            Store intermediate x* approximation after every iteration, last element is {x*: f(x*)=0} basedd on # iteration and/or tolerance\n",
    "        count: int \n",
    "            # Iterations used to compute the zero\n",
    "    \"\"\"\n",
    "    \n",
    "    # Update Step\n",
    "    def g(val):\n",
    "        return val - func(val)/dfunc(val)\n",
    "    \n",
    "    xk = np.zeros(Nmax+1)\n",
    "    count = 0\n",
    "    \n",
    "    xk[0] = x0 # First element is the initial guess\n",
    "    \n",
    "    # Make sure valid Nmax\n",
    "    if count >= Nmax:\n",
    "        print('Enter Nmax more than 0')\n",
    "        return None, None\n",
    "\n",
    "    # Stopping Condition: #iterations<Nmax, Interal below convergence tolerance or above divergence tolerance\n",
    "    while count < Nmax and np.abs(xk[count] - xk[count-1]) > tol and np.abs(xk[count] - xk[count-1]) < divtol:\n",
    "        xk[count+1] = g(xk[count])\n",
    "        count = count + 1\n",
    "        \n",
    "    roots = xk[:count+1]\n",
    "    \n",
    "    # Everything 0 when exceeded divergence tolerance\n",
    "    if np.abs(xk[count] - xk[count-1]) >= divtol:\n",
    "        roots = 0 \n",
    "        count = 0\n",
    "        \n",
    "    return roots, count  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Secant's Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Secant's Method\n",
    "def secant(func, x0,x1, tol=1e-12,divtol=1e2, Nmax=50): \n",
    "    \"\"\"Dhruv Jain, 7th Sept 2021\n",
    "    Obj: Apply Secant's Method to obtain the zero of the funcion using a set of inital guesses for a certain tolerance and/or # iterations \n",
    "    Args:\n",
    "        func: function\n",
    "            Evaluation function on which the Newton's method is applied\n",
    "        x0 float\n",
    "            First initial guess\n",
    "        x1: float\n",
    "            Second initial guess\n",
    "        tol: float\n",
    "            Convergence tolerance value (used as stopping condition)\n",
    "        divtol: float \n",
    "            Divergence tolerance value (used as stopping condition)\n",
    "        Nmax: int\n",
    "            # Allowed iterations (used as stopping condition)\n",
    "            \n",
    "    Output:\n",
    "        root: numpy ndarray, float \n",
    "            Store intermediate x* approximation after every iteration, last element is {x*: f(x*)=0} basedd on # iteration and/or tolerance\n",
    "        count: int \n",
    "            # Iterations used to compute the zero  \n",
    "    \"\"\"\n",
    "    \n",
    "    # Update Step\n",
    "    def secant_update(val0, val1):\n",
    "        return val1 - func(val1)*(val1-val0)/(func(val1)-func(val0))\n",
    "    \n",
    "    xk = np.zeros(Nmax+1)\n",
    "    count = 0\n",
    "    \n",
    "    # Make sure valid Nmax\n",
    "    if count >= Nmax:\n",
    "        print('Enter Nmax more than 0')\n",
    "        return None, None\n",
    "    \n",
    "    xk[0] = x0 # First element is the first initial guess\n",
    "    count = count + 1\n",
    "    xk[count] = x1 # Second element is the second initial guess\n",
    "    \n",
    "    # Stopping Condition: #iterations<Nmax, Interal below convergence tolerance or above divergence tolerance\n",
    "    while count < Nmax and np.abs(xk[count] - xk[count-1]) > tol and np.abs(xk[count] - xk[count-1]) < divtol:\n",
    "        xk[count+1] = secant_update(xk[count-1],xk[count])\n",
    "        count = count + 1\n",
    "        \n",
    "    roots = xk[:count+1]\n",
    "    \n",
    "    # Everything 0 when exceeded divergence tolerance\n",
    "    if np.abs(xk[count] - xk[count-1]) >= divtol:\n",
    "        roots = 0 \n",
    "        count = 0\n",
    "    \n",
    "    return roots, count  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fixed Point Iteration Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fixed_pt_iter(func, x0, tol=1e-12,divtol=1e3, Nmax=50): \n",
    "    \"\"\"Dhruv Jain, 7 Nov 2021        \n",
    "    Obj: Apply Fixed Point Iteration Method to obtain the zero of the funcion using an inital guesses for a certain tolerance and/or # iterations \n",
    "    Args:\n",
    "        func: function\n",
    "            Evaluation function on which the Newton's method is applied\n",
    "        x0: float\n",
    "            Initial guess\n",
    "        tol: float\n",
    "            Convergence tolerance value (used as stopping condition)\n",
    "        divtol: float \n",
    "            Divergence tolerance value (used as stopping condition)\n",
    "        Nmax: int\n",
    "            # Allowed iterations (used as stopping condition)\n",
    "            \n",
    "    Output:\n",
    "        xcalc: numpy ndarray, float \n",
    "            Store intermediate x* approximation after every iteration, last element is {x*: f(x*)=0} basedd on # iteration and/or tolerance\n",
    "        count: int \n",
    "            # Iterations used to compute the zero  \n",
    "    \"\"\"\n",
    "    xcalc=np.zeros(Nmax)\n",
    "    xcalc[0]= x0\n",
    "\n",
    "    count  = 0\n",
    "    \n",
    "    for i in range(1,Nmax):\n",
    "        \n",
    "        xcalc[i] = func(xcalc[i-1])\n",
    "        if (abs(xcalc[i]-xcalc[i-1]) < tol):\n",
    "            print('The convergence tolerance:',tol,'has been met')\n",
    "            xfinal = xcalc[:count+1]\n",
    "            break\n",
    "        \n",
    "        count = count + 1\n",
    "            \n",
    "        if ( abs(xcalc[i]-xcalc[i-1]) > divtol):    \n",
    "            print (\"No solution found\")\n",
    "            xfinal = np.zeros(count)\n",
    "            xfinal.fill(np.nan)\n",
    "            break\n",
    " \n",
    "    return xcalc, count "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: $f(x) = 1-ln(x)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Recall that the last element in the array is the converged solution and the other numbers show the estimated solution after every iteration\n",
      "The maximum number of iterations : 10 has been met\n",
      "\n",
      "Bisection Method: Root Estimate after 10 iteraions is:  [2.5        2.75       2.625      2.6875     2.71875    2.703125\n",
      " 2.7109375  2.71484375 2.71679688 2.71777344]\n",
      "\n",
      " Newton's Method: Root Estimate after 5 iteraions is:  [2.5        2.70927317 2.71826688 2.71828183 2.71828183]\n",
      "\n",
      " Secant's Method: Root Estimate after 6 iteraions is:  [2.5        2.6        2.71343127 2.71817548 2.71828173 2.71828183]\n",
      "\n",
      " Fixed Point Iteration Method: Root Estimate after 10 iteraions is:  [ 2.5         0.08370927  3.48040558 -0.24714883         nan         nan\n",
      "         nan         nan         nan         nan]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-22-e70c2ba7497a>:3: RuntimeWarning: invalid value encountered in log\n",
      "  return 1-np.log(x)\n"
     ]
    }
   ],
   "source": [
    "# Example function \n",
    "def f1(x):\n",
    "    return 1-np.log(x)\n",
    "#Derivative of example function\n",
    "def df1(x):\n",
    "    return -1/x\n",
    "\n",
    "# Input paramters\n",
    "ul = 3 # Upper limit of interval\n",
    "ll = 2 # Lowerr limit of interval\n",
    "x0 = 2.5 # First Initial guess\n",
    "x1 = 2.6 # Second Initial guess\n",
    "tol = 1e-5 # Convergence tolerance\n",
    "Nmax = 10 # Maximum number of allowed iterations\n",
    "\n",
    "print('Recall that the last element in the array is the converged solution and the other numbers show the estimated solution after every iteration')\n",
    "\n",
    "#Bisection Method\n",
    "bisection_res, bisec_count = bisection(f1,ll,ul,tol=tol,Nmax=Nmax)\n",
    "\n",
    "print('\\nBisection Method: Root Estimate after',len(bisection_res),'iteraions is: ',bisection_res)\n",
    "\n",
    "# Newton's Method\n",
    "new_res, new_count = newton(f1,df1,x0,tol=tol,Nmax=Nmax)\n",
    "print('\\n Newton\\'s Method: Root Estimate after',len(new_res),'iteraions is: ',new_res)\n",
    "\n",
    "#Secant's Method\n",
    "sec_res, sec_count = secant(f1,x0,x1,tol=tol,Nmax=Nmax)\n",
    "print('\\n Secant\\'s Method: Root Estimate after',len(sec_res),'iteraions is: ',sec_res)\n",
    "\n",
    "#Fixed Point Iteration Method\n",
    "fix_res, fix_count = fixed_pt_iter(f1,x0,tol=tol,Nmax=Nmax)\n",
    "print('\\n Fixed Point Iteration Method: Root Estimate after',len(fix_res),'iteraions is: ',fix_res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Method | Benefit | Disadvantage|\n",
    "|----|--------|-----|\n",
    "|Bisection Method | Does not require an intigal guess or derivative of a function | Linear convergence |\n",
    "|Newton's Method | Quadratic Convergence | Require first derivative of the function|\n",
    "|Secant's Method | Nearly quadratic convegence and does not require derivative of a function| Requires 2 inital guesses |\n",
    "|Fixed point iteration Method | Convergence can be faster than other orders | Does not gurantee convergence |"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}