{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Secant method"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The aim of the lab is to code a function which finds the zeroes of a function using the *secant* method.\n",
      "\n",
      "#### Secant method\n",
      "\n",
      "The idea for this method is simple. We are **given** a function $f$ whose roots we are supposed to find.\n",
      "\n",
      "+ Start with two distinct points $x_1$ and $x_2$.\n",
      "+ Find the line joining $(x_1, f(x_1))$ with $(x_2, f(x_2))$.\n",
      "+ Find the point $x_3$ where this line meets the $x$-axis.\n",
      "+ Repeat the procedure with $x_1 = x_2$ and $x_2 = x_3$.\n",
      "\n",
      "Starting with $x$ and $y$, the line joining $(x, f(x))$ and $(y, f(y))$ is\n",
      "\n",
      "$\\frac{t - f(x)}{s - x} = \\frac{f(y) - f(x)}{y - x}$.\n",
      "\n",
      "It is easy to see that $t = 0$ is attained for\n",
      "\n",
      "$s = x - \\frac{y - x}{f(y) - f(x)} f(x) = \\frac{x (f(y) - f(x)) - f(x)(y - x)}{f(y) - f(x)}$.\n",
      "\n",
      "#### The algorithm\n",
      "\n",
      "1. Start with some arbitrary x, y (preferably such that f(x) and f(y) have\n",
      "   opposite signs.\n",
      "1. Set `count_iterations` to zero.\n",
      "1. Repeat until f(y) is negigible :\n",
      "    1. z = x - ((y - x)/f(y) - f(x)) * f(x)\n",
      "    1. x = y\n",
      "    1. y = z\n",
      "    1. Increase `count_iterations` by 1.\n",
      "1. If f(y) is negligible,\n",
      "    1. return the result.\n",
      "    1. Otherwise, return None."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Problem 1\n",
      "\n",
      "Code the above algorithm. You can use a random generator to decide on the starting values. Note that, _it is not necessary for f(x) and f(y) to have opposite signs_."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Problem 2\n",
      "\n",
      "Modify the program so that if, in case, f(x) and f(y) have opposite signs, after compute z repeat the procedure on (x, z) if f(x) and f(z) have opposite signs; otherwise repeate the procedure on (z, y). This is called the *method of false position*."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Solutions (this is only for the T.A.s)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Secant method\n",
      "\n",
      "def next_point(f, x, y, e=1E-7) :\n",
      "    \"\"\"Given x and y as above finds z.\"\"\"\n",
      "    if abs(y - x) < e :\n",
      "        print \"x = %g and y = %g are too close.\" % (x, y)\n",
      "        retval = None\n",
      "    else :\n",
      "        slope = (f(y) - f(x))/(y - x)\n",
      "        if abs(slope) < e :\n",
      "            print \"Slope is too small to continue.\"\n",
      "            retval = None\n",
      "        else :\n",
      "            z = x - f(x)/slope\n",
      "            retval = z\n",
      "    return retval"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "if __name__ == '__main__' :\n",
      "    print next_point(lambda x : x*x - 100, 0, 1)\n",
      "    print next_point(lambda x : x*x*x*x - 16, -1, 1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100\n",
        "Slope is too small to continue.\n",
        "None\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def solve_by_secant(f, x = -10, y = 15, e = 1E-7, max_iter = 10000) :\n",
      "    \"\"\"Solves for f(x) = 0 using the secant method.\n",
      "    \n",
      "    f : function\n",
      "    x, y : initial starting point.\n",
      "    e : no. below which everything is considered to be 0.\n",
      "    max_iter : Maximum no of iterations to try.\n",
      "    \"\"\"\n",
      "    \n",
      "    x = float(x)\n",
      "    y = float(y)\n",
      "    z = next_point(f, x, y, e)\n",
      "    if z == None :\n",
      "        print \"Cannot continue with present values of x and y.\"\n",
      "        retval = None\n",
      "    else :\n",
      "        if abs(f(z)) < e :\n",
      "            # Yaay! We found a solution\n",
      "            retval = z\n",
      "        else :\n",
      "            cont_loop = True\n",
      "            no_iter = 0\n",
      "            while cont_loop :\n",
      "                no_iter += 1\n",
      "                x = y\n",
      "                y = z\n",
      "                z = next_point(f, x, y, e)\n",
      "                if z == None :\n",
      "                    print \"Cannot continue with present values of x and y.\"\n",
      "                    retval = None\n",
      "                    cont_loop = False\n",
      "                elif abs(f(z)) < e :\n",
      "                    # Eureka!\n",
      "                    retval = z\n",
      "                    cont_loop = False\n",
      "                elif no_iter > max_iter :\n",
      "                    print \"Maximum no of iterations reached.\"\n",
      "                    retval = None\n",
      "                    cont_loop = False\n",
      "    return retval"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "if __name__ == '__main__' :\n",
      "    print solve_by_secant(lambda x : x*x - 4, 999.00001, 1000, 1E-7, 1000)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2.00000000001\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# To solve the second question regarding method of false position.\n",
      "# We'll also give an alternate code for secant method\n",
      "\n",
      "def next_pair(f, x, y, e, method=\"false_position\") :\n",
      "    z = next_point(f, x, y, e)\n",
      "    if z == None :\n",
      "        print \"Cannot find next pair.\"\n",
      "        retval = None\n",
      "    elif abs(f(z)) < e :\n",
      "        retval = (z, \"root\")\n",
      "    else :\n",
      "        if method == \"false_position\" :\n",
      "            if f(x) * f(y) > 0 :\n",
      "                print \"For method of false position, f(x) and f(y)\",\n",
      "                print \"should have opposite signs.\"\n",
      "                retval = None\n",
      "            else :\n",
      "                if f(x) * f(z) < 0 :\n",
      "                    retval = ((x, z), \"pair\")\n",
      "                else :\n",
      "                    retval = ((z, y), \"pair\")\n",
      "        elif method == \"secant\" :\n",
      "            retval = ((y, z), \"pair\")\n",
      "        else :\n",
      "            print \"Unknown method.\"\n",
      "            retval = None\n",
      "    return retval"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def find_root(f, x=-10, y=15, e=1E-5, max_iter=10000, method=\"false_position\", show_steps=False) :\n",
      "    x = float(x)\n",
      "    y = float(y)\n",
      "    cont_loop = True\n",
      "    no_iter = 0\n",
      "    while cont_loop :\n",
      "        no_iter += 1\n",
      "        z = next_pair(f, x, y, e, method)\n",
      "        if z == None :\n",
      "            print \"Cannot continue.\"\n",
      "            retval = None\n",
      "            cont_loop = False\n",
      "        elif z[1] == \"root\" :\n",
      "            retval = z[0]\n",
      "            cont_loop = False\n",
      "        elif z[1] != \"pair\" :\n",
      "            print \"If you see this message, there is a bug\",\n",
      "            print \"in the code.\"\n",
      "            retval = None\n",
      "            cont_loop = False\n",
      "        elif no_iter > max_iter :\n",
      "            print \"Maximum of iterations reached.\"\n",
      "            retval = None\n",
      "            cont_loop = False\n",
      "        else :\n",
      "            (x, y) = z[0]\n",
      "            if show_steps :\n",
      "                print \"%7d : %9g, %9g\" %(no_iter, x, y)\n",
      "            \n",
      "    return retval"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "if __name__ == '__main__' :\n",
      "    import math\n",
      "    print find_root(lambda t : math.sin(t * math.pi/180.0), 31.234, 1000.876876, show_steps=True, method=\"secant\")\n",
      "    print find_root(lambda t : math.sin(t * math.pi/180.0), 31.234, 1000.876876, show_steps=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "      1 :   1000.88,   366.302\n",
        "      2 :   366.302,   430.105\n",
        "      3 :   430.105,   357.869\n",
        "      4 :   357.869,   360.617\n",
        "359.999898995\n",
        "      1 :   366.302,   1000.88\n",
        "      2 :   430.105,   1000.88\n",
        "      3 :   430.105,   709.297\n",
        "      4 :   430.105,   663.252\n",
        "      5 :   430.105,   553.505\n",
        "      6 :   528.955,   553.505\n",
        "      7 :   528.955,   540.019\n",
        "539.999883922\n"
       ]
      }
     ],
     "prompt_number": 7
    }
   ],
   "metadata": {}
  }
 ]
}