{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Physics 503: Scientific Computing\n",
    "### Lecture 20: Root Finding\n",
    "#### J.R. Gladden, Univ. of Mississippi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finding the root of an equation means determining the value for the independent variable for which the function is 0. If this happens multiple times, then there are multiple roots.  Mathematically, the roots $x_i$ of a function $F(x) are such that $F(x_i)=0$. This is useful and necessary in all kinds of analysis such as solving optimization problems, determining critical behavior in a system, ...  More specific examples are given in the slides.\n",
    "\n",
    "### Simplest Methods\n",
    "The most obvious methods are to (a) simply plot the function and zoom in on the visual $x=0$ crossings or (b) a \"brute force\" method where you march along $x$ values and evaluate $F(x)$ and look for the case where $F(x)$ switches sign. These can work, but are not very efficient or easy to automate.  Here are two better options..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bisection Method\n",
    "Based on the Binary Search logic we discussed earlier, for a range between $x_{low}$ to $x_{high}$, test if $F(x_{low})\\cdot F(x_{high}) <0$.\n",
    "- If **no** , then there are no roots (or an even number!)\n",
    "- If **yes**, then test if $F(x_{low})\\cdot F(x_{mid} <0$\n",
    "     - if yes, then set $x_{high} = x_{mid}$ and repeat\n",
    "     - if no, then set $x_{low} = x_{mid}$ and repeat\n",
    "\n",
    "Keep repeating until $|x_{high} - x_{low}|<\\epsilon$ (within some tolerance value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib wx\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bisect(f,xlow,xhigh,switch=0,tol=1.0e-8):\n",
    "    '''\n",
    "    Supplied functions:\n",
    "1. bisect()\n",
    "Usage: root=bisect(f,xlow,xhigh,switch=0,tol=1.0e-8)\n",
    "Returns: the root in the bracket (float)\n",
    "Parameters:\n",
    "f - function f(x) to find the root of\n",
    "xlow,xhigh - lower and upper bounds of root\n",
    "switch - whether to check whether |f(x)| is getting smaller with each step \n",
    "tol - accuracy of root\n",
    "    '''\n",
    "    \n",
    "    from math import ceil,log\n",
    "    #for single numbers, these are more efficient than those supplied by numpy\n",
    "    \n",
    "    f1=f(xlow)\n",
    "    if f1 == 0. : return xlow\n",
    "    f2=f(xhigh)\n",
    "    if f2 == 0. : return xhigh\n",
    "    if f1*f2 > 0. :\n",
    "        print \"An even number of roots, or no roots, lie inside the bracket.\"\n",
    "        return None\n",
    "    n= int(ceil(log(abs(xhigh - xlow)/tol)/log(2.0)))\n",
    "    print \"Number of iteration: \", n\n",
    "    for i in range(n):\n",
    "        xmid = (xhigh + xlow)/2.0\n",
    "        f3 = f(xmid)\n",
    "        if (switch == 1) and (abs(f3) > abs(f1)) and (abs(f3) > abs(f2)):\n",
    "            return None\n",
    "        if f3 == 0. : return xmid\n",
    "        if f2*f3 < 0. :\n",
    "            xlow = xmid\n",
    "            f1 = f3\n",
    "        else:\n",
    "            xhigh = xmid\n",
    "            f2 = f3\n",
    "        print \"Intermediate root is: \", xmid\n",
    "    return (xhigh + xlow)/2.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of iteration:  28\n",
      "Intermediate root is:  1.0\n",
      "Intermediate root is:  1.5\n",
      "Intermediate root is:  1.25\n",
      "Intermediate root is:  1.125\n",
      "Intermediate root is:  1.0625\n",
      "Intermediate root is:  1.03125\n",
      "Intermediate root is:  1.046875\n",
      "Intermediate root is:  1.0546875\n",
      "Intermediate root is:  1.05078125\n",
      "Intermediate root is:  1.048828125\n",
      "Intermediate root is:  1.0478515625\n",
      "Intermediate root is:  1.04833984375\n",
      "Intermediate root is:  1.04858398438\n",
      "Intermediate root is:  1.04870605469\n",
      "Intermediate root is:  1.04876708984\n",
      "Intermediate root is:  1.04879760742\n",
      "Intermediate root is:  1.04878234863\n",
      "Intermediate root is:  1.04878997803\n",
      "Intermediate root is:  1.04879379272\n",
      "Intermediate root is:  1.04879188538\n",
      "Intermediate root is:  1.04879283905\n",
      "Intermediate root is:  1.04879331589\n",
      "Intermediate root is:  1.04879355431\n",
      "Intermediate root is:  1.0487934351\n",
      "Intermediate root is:  1.0487934947\n",
      "Intermediate root is:  1.0487934649\n",
      "Intermediate root is:  1.0487934798\n",
      "Intermediate root is:  1.04879348725\n"
     ]
    }
   ],
   "source": [
    "def f(x): return x**3*np.sin(x)-1 #-(x-2.0)**2+0.5\n",
    "\n",
    "xlow=0\n",
    "xhigh=2\n",
    "root = bisect(f,xlow,xhigh)\n",
    "\n",
    "x=np.linspace(xlow,xhigh,100)\n",
    "plt.close('all')\n",
    "plt.plot(x,f(x))\n",
    "plt.grid()\n",
    "if root:\n",
    "    plt.annotate(\"Root\",xytext=(root,1),xy=(root,0),arrowprops={'arrowstyle':'->'})\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Newton-Raphson\n",
    "Newton-Raphson is more sophisticated and can usually jump to the root in many fewer iterations.  The theory is laid out in the slides, but the fundamental concept is to us the derivative of the function to point \"downhill\" to find a new $x$ value and repeat until a tolerance is obtained. Problems can arrise if you happen to land on an $x$ value at (or near) an extremum so the derviative is very close to zero. The code below also adds interesting visualization which shows you that path to the root from the starting point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def newtonraphson(f,xlow,xhigh,xStart,maxIter=10,tol=1.0e-8,plotit=True):\n",
    "    \n",
    "    x0 = xStart\n",
    "    def dfdx(x0,f,dx=tol/2.0):\n",
    "        return (f(x0+dx)-f(x0))/dx  #Use euler to get derivative\n",
    "\n",
    "    if plotit:\n",
    "        line=[]\n",
    "    #Loop until x0 changes less than tol\n",
    "    for Step in range(maxIter+1):\n",
    "        xnew=x0-f(x0)/dfdx(x0,f,dx=abs(xhigh-xlow)/1000.)\n",
    "        if abs(xnew - x0) < tol:\n",
    "            print 'Found root in ', Step, ' iterations.'\n",
    "            break\n",
    "        m=-f(x0)/(xnew-x0)\n",
    "        if plotit: line.append((x0,f(x0),xnew,0))\n",
    "        x0=xnew\n",
    "        print \"For iteration \", Step, \", the root is: %2.6f\" % x0\n",
    "    \n",
    "    ## Print out final root and number of required iters.\n",
    "    print '-'*20\n",
    "    if  Step == maxIter:\n",
    "            print \"No roots were found after\", maxIter,\" iterations.\"\n",
    "    else:\n",
    "            print 'From initial guess ', xStart, ', the root found was: ', x0\n",
    " \n",
    "\n",
    "    if plotit:\n",
    "        x=np.linspace(xlow,xhigh,200)\n",
    "        plt.figure(1)\n",
    "        plt.close('all')\n",
    "        ##Plot function with initial guess (black diamond) and root found (red circle)\n",
    "        plt.plot(x,f(x),'b-',[xStart],[0],'kD',[x0],[0],'ro',linewidth=2)\n",
    "        plt.plot([xStart,xStart],[0,f(xStart)],'--k')\n",
    "        ## Also plot the tangent lines showing the path to the root\n",
    "        for l in range(Step):\n",
    "            x0=line[l][0]\n",
    "            y0=line[l][1]\n",
    "            x1=line[l][2]\n",
    "            y1=line[l][3]\n",
    "            dx=x1-x0\n",
    "            dy=y1-y0\n",
    "            #arrow(x0,y0,dx,dy,facecolor='grey',alpha=0.5)\n",
    "            plt.annotate('', (x1,y1), xytext=(x0,y0),arrowprops=dict(facecolor='grey',alpha=0.5,width=4))\n",
    "            plt.plot([x1,x1],[0,f(x1)],'--k')\n",
    "        plt.grid()\n",
    "        plt.xlabel('X')\n",
    "        plt.ylabel('f(x)')\n",
    "        plt.title('Newton-Raphson Method')\n",
    "        plt.show()\n",
    "    return x0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For iteration  0 , the root is: 1.701176\n",
      "For iteration  1 , the root is: 1.214515\n",
      "For iteration  2 , the root is: 1.072848\n",
      "For iteration  3 , the root is: 1.049545\n",
      "For iteration  4 , the root is: 1.048798\n",
      "For iteration  5 , the root is: 1.048794\n",
      "For iteration  6 , the root is: 1.048793\n",
      "Found root in  7  iterations.\n",
      "--------------------\n",
      "From initial guess  0.6 , the root found was:  1.04879348388\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.0487935109404454"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def f(x): return x**3*np.sin(x)-1 #-(x-2.0)**2+0.5\n",
    "\n",
    "xlow=0\n",
    "xhigh=5\n",
    "newtonraphson(f,xlow,xhigh,0.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}