Bisection Method for 1D functions-Lesson 2 (Edited Newton's Method).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bisection Method and (Fixed) Newton's method\n",
    "\n",
    "## by Dion Ho\n",
    "\n",
    "### Last updated: January 20, 2019\n",
    "\n",
    "\n",
    "# Import Python Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Defining the Bisection Routine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bisection_method(a,b,f,tol,N):\n",
    "    FA,FB = f(a), f(b)\n",
    "    iteration = 1\n",
    "    if FA*FB<0:\n",
    "        print (\"iter      a           b            root:(x,_)           root:(_,y)\")\n",
    "        print (\"------------------------------------------------------------------\")\n",
    "        p = (a + b)/2\n",
    "        FP = f(p)\n",
    "        while (np.abs((b-a)/2)>tol) & (iteration <= N):\n",
    "            if FA*FP<0:\n",
    "                b=p\n",
    "                FB=FP\n",
    "            else:\n",
    "                a=p\n",
    "                FA=FP\n",
    "            p = (a+b)/2\n",
    "            FP = f(p)\n",
    "            print('{:>3d}  {:> 10.5f}  {:> 10.5f}  {:> 22.16f}  {:>13.6e}'.format(iteration, a, b, p, abs(FP)))\n",
    "            iteration = iteration + 1\n",
    "        if np.abs((b-a)/2)>tol:\n",
    "            print(\"Method failed to converge\")\n",
    "        else:\n",
    "            return p\n",
    "    else:\n",
    "             print(\"Cannot ensure existence of root\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Testing the Bisection Method\n",
    "Here, let's try to solve the equation $3x^2 + 10x + 8 = 0$ on the interval $\\left[-1.5,0\\right]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda x: 3*x**2 + 10*x + 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter      a           b            root:(x,_)           root:(_,y)\n",
      "------------------------------------------------------------------\n",
      "  1    -1.50000    -0.75000     -1.1250000000000000   5.468750e-01\n",
      "  2    -1.50000    -1.12500     -1.3125000000000000   4.296875e-02\n",
      "  3    -1.50000    -1.31250     -1.4062500000000000   1.298828e-01\n",
      "  4    -1.40625    -1.31250     -1.3593750000000000   5.004883e-02\n",
      "  5    -1.35938    -1.31250     -1.3359375000000000   5.187988e-03\n",
      "  6    -1.33594    -1.31250     -1.3242187500000000   1.847839e-02\n",
      "  7    -1.33594    -1.32422     -1.3300781250000000   6.542206e-03\n",
      "  8    -1.33594    -1.33008     -1.3330078125000000   6.513596e-04\n",
      "  9    -1.33594    -1.33301     -1.3344726562500000   2.274752e-03\n",
      " 10    -1.33447    -1.33301     -1.3337402343750000   8.133054e-04\n",
      " 11    -1.33374    -1.33301     -1.3333740234375000   8.137524e-05\n",
      " 12    -1.33337    -1.33301     -1.3331909179687500   2.848916e-04\n",
      " 13    -1.33337    -1.33319     -1.3332824707031250   1.017330e-04\n",
      " 14    -1.33337    -1.33328     -1.3333282470703125   1.017260e-05\n",
      " 15    -1.33337    -1.33333     -1.3333511352539062   3.560289e-05\n",
      " 16    -1.33335    -1.33333     -1.3333396911621094   1.271554e-05\n",
      " 17    -1.33334    -1.33333     -1.3333339691162109   1.271565e-06\n",
      " 18    -1.33333    -1.33333     -1.3333311080932617   4.450495e-06\n",
      " 19    -1.33333    -1.33333     -1.3333325386047363   1.589459e-06\n",
      " 20    -1.33333    -1.33333     -1.3333332538604736   1.589457e-07\n",
      " 21    -1.33333    -1.33333     -1.3333336114883423   5.563098e-07\n",
      " 22    -1.33333    -1.33333     -1.3333334326744080   1.986821e-07\n",
      " 23    -1.33333    -1.33333     -1.3333333432674408   1.986821e-08\n",
      " 24    -1.33333    -1.33333     -1.3333332985639572   6.953876e-08\n",
      " 25    -1.33333    -1.33333     -1.3333333209156990   2.483527e-08\n",
      " 26    -1.33333    -1.33333     -1.3333333320915699   2.483527e-09\n",
      " 27    -1.33333    -1.33333     -1.3333333376795053   8.692343e-09\n",
      " 28    -1.33333    -1.33333     -1.3333333348855376   3.104407e-09\n",
      " 29    -1.33333    -1.33333     -1.3333333334885538   3.104397e-10\n",
      " 30    -1.33333    -1.33333     -1.3333333327900618   1.086543e-09\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-1.3333333327900618"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bisection_method(-1.5,0,f,1e-9,50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Newton's Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def newton_method(x0,f,fp,tol,N):\n",
    "    F = f(x0)\n",
    "    Fp = fp(x0)\n",
    "    iteration = 1\n",
    "    print (\"iter      grad         root:(x,_)           root:(_,y)\")\n",
    "    print (\"------------------------------------------------------\")\n",
    "    while (iteration<=N) & (np.abs(F)>tol):\n",
    "        x = (Fp*x0 - F)/Fp\n",
    "        x0 = x\n",
    "        F = f(x0)\n",
    "        Fp = fp(x0)\n",
    "        print('{:>3d}  {:> 10.5f}  {:> 22.16f}  {:>13.6e}'.format(iteration, Fp.real, x.real, abs(F.real)))\n",
    "        iteration = iteration + 1\n",
    "    if np.abs(F)<=tol:\n",
    "        return x0\n",
    "    else:\n",
    "        print(\"Method failed to converge. Try harder!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter      grad         root:(x,_)           root:(_,y)\n",
      "------------------------------------------------------\n",
      "  1   302.00331     48.6672185430463600   7.600167e+03\n",
      "  2   151.00828     23.5013796788473357   1.899958e+03\n",
      "  3    75.51738     10.9195638905965779   4.749063e+02\n",
      "  4    37.78518      4.6308626067112613   1.186433e+02\n",
      "  5    18.94552      1.4909197718979010   2.957772e+01\n",
      "  6     9.57833     -0.0702791370981974   7.312026e+00\n",
      "  7     4.99797     -0.8336721066328142   1.748306e+00\n",
      "  8     2.89915     -1.1834756070969692   3.670875e-01\n",
      "  9     2.13943     -1.3100947637572193   4.809723e-02\n",
      " 10     2.00454     -1.3325760791700947   1.516229e-03\n",
      " 11     2.00001     -1.3333324751321609   1.716405e-06\n",
      " 12     2.00000     -1.3333333333322286   2.208900e-12\n",
      " 13     2.00000     -1.3333333333333330   8.881784e-16\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-1.333333333333333"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = lambda x: 3*x**2 + 10*x + 8\n",
    "fp = lambda x: 6*x + 10\n",
    "newton_method(99,f,fp,1e-15,20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter      grad         root:(x,_)           root:(_,y)\n",
      "------------------------------------------------------\n",
      "  1    -0.52640      0.1733333333333333   1.101658e+00\n",
      "  2     1.40661     -0.9385562279973102   6.666576e-01\n",
      "  3    -0.48970     -0.5322071222832284   4.378873e-01\n",
      "  4    -1.75342     -0.4569866275367542   7.839953e-03\n",
      "  5    -1.49342     -0.4986483767110790   7.330294e-03\n",
      "  6    -1.50001     -0.5000036235897229   1.799275e-05\n",
      "  7    -1.50000     -0.5000000000252197   1.063500e-10\n",
      "  8    -1.50000     -0.4999999999999999   0.000000e+00\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-0.49999999999999994+0.8660254037844387j)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = lambda x: x**3 - 1\n",
    "fp = lambda x: 3*(x**2)\n",
    "newton_method(0.5-1j,f,fp,1e-10,20) #Newton's method converges quadratically fast.\n",
    "#bisection_method(0,-1-1j,f,1e-10,50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Secant Method (Edited Newton's Method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def secant_method(x0,x1,f,tol,N):\n",
    "    F0 = f(x0)\n",
    "    F1 = f(x1)\n",
    "    iteration = 2\n",
    "    print (\"iter       root:(x,_)           root:(_,y)\")\n",
    "    print (\"------------------------------------------\")\n",
    "    while (iteration<=N) & (np.abs(F1)>tol):\n",
    "        x = x1 - F1*(x1-x0)/(F1-F0)\n",
    "        x0 = x1\n",
    "        x1 = x\n",
    "        F0 = f(x0)\n",
    "        F1 = f(x1)\n",
    "        print('{:>3d}  {:> 22.16f}  {:>13.6e}'.format(iteration, x.real, abs(F1.real)))\n",
    "        iteration = iteration + 1\n",
    "    if np.abs(F1)<=tol:\n",
    "        return x1\n",
    "    else:\n",
    "        print(\"Method failed to converge. Try harder!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter       root:(x,_)           root:(_,y)\n",
      "------------------------------------------\n",
      "  2      0.2500000000000000   9.843750e-01\n",
      "  3      0.4657534246575342   8.989659e-01\n",
      "  4      2.7366456314994219   1.949537e+01\n",
      "  5      0.5658525355572288   8.188202e-01\n",
      "  6      0.6533524304561957   7.211038e-01\n",
      "  7      1.2990633142633885   1.192254e+00\n",
      "  8      0.8967070634963772   2.789726e-01\n",
      "  9      0.9730014529021270   7.882856e-02\n",
      " 10      1.0030506947105988   9.180033e-03\n",
      " 11      0.9999163071920949   2.510574e-04\n",
      " 12      0.9999997451830924   7.644505e-07\n",
      " 13      1.0000000000213276   6.398282e-11\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.0000000000213276"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = lambda x: x**3 - 1\n",
    "secant_method(0,2,f,1e-10,50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## False Position + Secant Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def false_position_secant_method(x0,x1,f,tol,N):\n",
    "    F0 = f(x0)\n",
    "    F1 = f(x1)\n",
    "    iteration = 2\n",
    "    print (\"iter       root:(x,_)           root:(_,y)\")\n",
    "    print (\"------------------------------------------\")\n",
    "    while (iteration<=N) & (np.abs(F1)>tol):\n",
    "        x2 = x1 - F1*(x1-x0)/(F1-F0)\n",
    "        if f(x2)*f(x1) < 0:\n",
    "            x0 = x1\n",
    "            x1 = x2\n",
    "        else:\n",
    "            x0 = x0\n",
    "            x1 = x2\n",
    "        F0 = f(x0)\n",
    "        F1 = f(x1)\n",
    "        print('{:>3d}  {:> 22.16f}  {:>13.6e}'.format(iteration, x2.real, abs(F1.real)))\n",
    "        iteration = iteration + 1\n",
    "    if np.abs(F1)<=tol:\n",
    "        return x2\n",
    "    else:\n",
    "        print(\"Method failed to converge. Try harder!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter       root:(x,_)           root:(_,y)\n",
      "------------------------------------------\n",
      "  2      0.2500000000000000   9.843750e-01\n",
      "  3      0.4657534246575342   8.989659e-01\n",
      "  4      0.6403630266802741   7.374097e-01\n",
      "  5      0.7699424960309564   5.435693e-01\n",
      "  6      0.8585771251953758   3.670959e-01\n",
      "  7      0.9154532149617496   2.328002e-01\n",
      "  8      0.9503612361259737   1.416466e-01\n",
      "  9      0.9711796483435732   8.399316e-02\n",
      " 10      0.9833781169037836   4.904138e-02\n",
      " 11      0.9904509283841171   2.837453e-02\n",
      " 12      0.9945266192248776   1.633043e-02\n",
      " 13      0.9968668475591338   9.370038e-03\n",
      " 14      0.9982078234034287   5.366900e-03\n",
      " 15      0.9989753090132928   3.070924e-03\n",
      " 16      0.9994142694126268   1.756163e-03\n",
      " 17      0.9996652337877771   1.003962e-03\n",
      " 18      0.9998086844367171   5.738369e-04\n",
      " 19      0.9998906700980669   3.279538e-04\n",
      " 20      0.9999375235748394   1.874176e-04\n",
      " 21      0.9999642984686835   1.071008e-04\n",
      " 22      0.9999795988908509   6.120208e-05\n",
      " 23      0.9999883421468972   3.497315e-05\n",
      " 24      0.9999933383446933   1.998483e-05\n",
      " 25      0.9999961933316738   1.141996e-05\n",
      " 26      0.9999978247582949   6.525711e-06\n",
      " 27      0.9999987570038709   3.728984e-06\n",
      " 28      0.9999992897162139   2.130850e-06\n",
      " 29      0.9999995941234582   1.217629e-06\n",
      " 30      0.9999997680705173   6.957883e-07\n",
      " 31      0.9999998674688572   3.975934e-07\n",
      " 32      0.9999999242679152   2.271962e-07\n",
      " 33      0.9999999567245219   1.298264e-07\n",
      " 34      0.9999999752711550   7.418653e-08\n",
      " 35      0.9999999858692313   4.239231e-08\n",
      " 36      0.9999999919252750   2.422417e-08\n",
      " 37      0.9999999953858714   1.384239e-08\n",
      " 38      0.9999999973633551   7.909935e-09\n",
      " 39      0.9999999984933458   4.519963e-09\n",
      " 40      0.9999999991390547   2.582836e-09\n",
      " 41      0.9999999995080312   1.475906e-09\n",
      " 42      0.9999999997188750   8.433750e-10\n",
      " 43      0.9999999998393572   4.819285e-10\n",
      " 44      0.9999999999082041   2.753877e-10\n",
      " 45      0.9999999999475452   1.573645e-10\n",
      " 46      0.9999999999700259   8.992240e-11\n",
      "iter      a           b            root:(x,_)           root:(_,y)\n",
      "------------------------------------------------------------------\n",
      "  1     1.00000     2.00000      1.5000000000000000   2.375000e+00\n",
      "  2     1.50000     2.00000      1.7500000000000000   4.359375e+00\n",
      "  3     1.75000     2.00000      1.8750000000000000   5.591797e+00\n",
      "  4     1.87500     2.00000      1.9375000000000000   6.273193e+00\n",
      "  5     1.93750     2.00000      1.9687500000000000   6.630829e+00\n",
      "  6     1.96875     2.00000      1.9843750000000000   6.813961e+00\n",
      "  7     1.98438     2.00000      1.9921875000000000   6.906616e+00\n",
      "  8     1.99219     2.00000      1.9960937500000000   6.953216e+00\n",
      "  9     1.99609     2.00000      1.9980468750000000   6.976585e+00\n",
      " 10     1.99805     2.00000      1.9990234375000000   6.988287e+00\n",
      " 11     1.99902     2.00000      1.9995117187500000   6.994142e+00\n",
      " 12     1.99951     2.00000      1.9997558593750000   6.997071e+00\n",
      " 13     1.99976     2.00000      1.9998779296875000   6.998535e+00\n",
      " 14     1.99988     2.00000      1.9999389648437500   6.999268e+00\n",
      " 15     1.99994     2.00000      1.9999694824218750   6.999634e+00\n",
      " 16     1.99997     2.00000      1.9999847412109375   6.999817e+00\n",
      " 17     1.99998     2.00000      1.9999923706054688   6.999908e+00\n",
      " 18     1.99999     2.00000      1.9999961853027344   6.999954e+00\n",
      " 19     2.00000     2.00000      1.9999980926513672   6.999977e+00\n",
      " 20     2.00000     2.00000      1.9999990463256836   6.999989e+00\n",
      " 21     2.00000     2.00000      1.9999995231628418   6.999994e+00\n",
      " 22     2.00000     2.00000      1.9999997615814209   6.999997e+00\n",
      " 23     2.00000     2.00000      1.9999998807907104   6.999999e+00\n",
      " 24     2.00000     2.00000      1.9999999403953552   6.999999e+00\n",
      " 25     2.00000     2.00000      1.9999999701976776   7.000000e+00\n",
      " 26     2.00000     2.00000      1.9999999850988388   7.000000e+00\n",
      " 27     2.00000     2.00000      1.9999999925494194   7.000000e+00\n",
      " 28     2.00000     2.00000      1.9999999962747097   7.000000e+00\n",
      " 29     2.00000     2.00000      1.9999999981373549   7.000000e+00\n",
      " 30     2.00000     2.00000      1.9999999990686774   7.000000e+00\n",
      " 31     2.00000     2.00000      1.9999999995343387   7.000000e+00\n",
      " 32     2.00000     2.00000      1.9999999997671694   7.000000e+00\n",
      " 33     2.00000     2.00000      1.9999999998835847   7.000000e+00\n",
      " 34     2.00000     2.00000      1.9999999999417923   7.000000e+00\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.9999999999417923"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = lambda x: x**3 - 1\n",
    "false_position_secant_method(0,2,f,1e-10,50)\n",
    "bisection_method(0,2,f,1e-10,50)"
   ]
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}