{ "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": {} } ] }