root-finding.ipynb

{
 "metadata": {
  "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.5-final"
  },
  "orig_nbformat": 2,
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Root finding problem\n",
    "\n",
    "Finding root or solution for an equation of the form `f(x) = 0`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding simple roots\n",
    "\n",
    "We will first look at finding \"simple\" roots, i.e. roots with multiplicuty 1,\n",
    "of the function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bisection method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bisection(f, a, b, tol):\n",
    "    \"\"\"\n",
    "    Arguments\n",
    "    ---------\n",
    "    f: function\n",
    "        function 'f' to find root of\n",
    "    a: float\n",
    "    b: float\n",
    "        a and b such that f(a)*f(b) < 0\n",
    "    tol: float\n",
    "        tolerance\n",
    "    \"\"\"\n",
    "    \n",
    "    p = (a + b)/2.0\n",
    "    while b - p > tol:\n",
    "        if f(p) == 0:\n",
    "            break\n",
    "        elif f(a)*f(p) < 0:\n",
    "            b = p\n",
    "        else:\n",
    "            a = p\n",
    "        p = (a + b)/2.0\n",
    "    \n",
    "    return p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Using Bisection method\n",
    "\n",
    "from math import cos, sin, inf\n",
    "\n",
    "f = lambda x: sin(x) + 3*cos(x) - 2\n",
    "a = 1\n",
    "b = 2\n",
    "tol = 1e-6\n",
    "\n",
    "bisection(f, a, b, tol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fixed Point method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fixed_point(g, x, tol, n):\n",
    "    \"\"\"\n",
    "    Arguments\n",
    "    ---------\n",
    "    g: function\n",
    "        function 'g' such that x = g(x)\n",
    "    x: float\n",
    "        Initial guess of root\n",
    "    tol: float\n",
    "        tolerance\n",
    "    n: int\n",
    "        maximum number of iterations\n",
    "    \"\"\"\n",
    "    next_sequence_item = lambda x: g(x)\n",
    "\n",
    "    x_old = x\n",
    "    x = next_sequence_item(x)\n",
    "    n = n - 1\n",
    "    while abs(x - x_old) > tol and n > 0:\n",
    "        x_old = x\n",
    "        x = next_sequence_item(x)\n",
    "        n = n - 1\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "## Using Fixed Point method\n",
    "\n",
    "from math import cos, sin, acos, inf\n",
    "\n",
    "f = lambda x: sin(x) + 3*cos(x) - 2\n",
    "g = lambda x: acos((2 - sin(x))/3)\n",
    "x = 1.2\n",
    "tol = 1e-6\n",
    "n = inf\n",
    "\n",
    "fixed_point(g, x, tol, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Newton's method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def newton(f, f_der, x, tol, n):\n",
    "    \"\"\"\n",
    "    Arguments\n",
    "    ---------\n",
    "    f: function\n",
    "        function `f` to find root of\n",
    "    f_der: function\n",
    "        Derivative of `f`\n",
    "    x: float\n",
    "        Initial guess of root\n",
    "    tol: float \n",
    "        tolerance\n",
    "    n: int\n",
    "        maximum number of iterations\n",
    "    \"\"\"\n",
    "    next_sequence_item = lambda x: x - f(x)/f_der(x)\n",
    "    \n",
    "    x_old = x\n",
    "    x = next_sequence_item(x)\n",
    "    n = n - 1\n",
    "    print(x, f(x))\n",
    "    while abs(x - x_old) > tol and n > 0:\n",
    "        x_old = x\n",
    "        x = next_sequence_item(x)\n",
    "        n = n - 1\n",
    "        print(x, f(x))\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Using newton's method\n",
    "\n",
    "from math import cos, sin, inf\n",
    "\n",
    "f = lambda x: sin(x) + 3*cos(x) - 2\n",
    "f_der = lambda x: cos(x) - 3*sin(x)\n",
    "x = 1\n",
    "tol = 1e-6\n",
    "n = inf\n",
    "\n",
    "newton(f, f_der, x, tol, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Secant Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def secant(f, x0, x1, tol, n):\n",
    "    \"\"\"\n",
    "    Arguments\n",
    "    ---------\n",
    "    f: function\n",
    "        function `f` to find root of\n",
    "    x0, x1: float\n",
    "        Initial two guess of root\n",
    "    tol: float \n",
    "        tolerance\n",
    "    n: int\n",
    "        maximum number of iterations\n",
    "    \"\"\"\n",
    "    next_sequence_item = lambda x0, x1: x1 - f(x1)*(x1-x0)/(f(x1)-f(x0))\n",
    "\n",
    "    x2 = next_sequence_item(x0, x1)\n",
    "    n = n - 1\n",
    "    while abs(x2 - x1) > tol and n > 0:\n",
    "        x0 = x1\n",
    "        x1 = x2\n",
    "        x2 = next_sequence_item(x0, x1)\n",
    "        n = n - 1\n",
    "        # print(x)\n",
    "\n",
    "    return x1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Using secant method\n",
    "\n",
    "from math import sin, cos, inf\n",
    "\n",
    "f = lambda x: sin(x) + 3*cos(x) - 2\n",
    "x0 = 0\n",
    "x1 = 1.5\n",
    "tol = 1e-6\n",
    "n = inf\n",
    "\n",
    "secant(f, x0, x1, tol, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## False Position Method "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import inf\n",
    "\n",
    "def false_position(f, x0, x1, tol, n):\n",
    "    \"\"\"\n",
    "    Arguments\n",
    "    ---------\n",
    "    f: function\n",
    "        function `f` to find root of\n",
    "    x0, x1: float\n",
    "        Initial two guess of root\n",
    "    tol: float \n",
    "        tolerance\n",
    "    n: int\n",
    "        maximum number of iterations\n",
    "    \"\"\"\n",
    "    next_sequence_item = lambda x0, x1: x1 - f(x1)*(x1-x0)/(f(x1)-f(x0))\n",
    "    x2 = next_sequence_item(x0, x1)\n",
    "    n = n - 1\n",
    "    # print(x2, f(x2))\n",
    "\n",
    "    while abs(x2 - x1) > tol and n > 0:\n",
    "        if f(x1)*f(x2) < 0:\n",
    "            x3 = next_sequence_item(x1, x2)\n",
    "            x0 = x1\n",
    "            x1 = x2\n",
    "            x2 = x3\n",
    "        else:\n",
    "            x3 = next_sequence_item(x0, x2)\n",
    "            x1 = x2\n",
    "            x2 = x3\n",
    "        n = n - 1\n",
    "        # print(x2, f(x2))\n",
    "    \n",
    "    return x2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "## Using false position\n",
    "\n",
    "from math import sin, cos, inf\n",
    "\n",
    "f = lambda x: sin(x) + 3*cos(x) - 2\n",
    "x0 = 0\n",
    "x1 = 1.5\n",
    "tol = 1e-6\n",
    "n = inf\n",
    "\n",
    "false_position(f, x0, x1, tol, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding non-simple roots\n",
    "\n",
    "Finding roots with multiplicity `m > 1`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modified Netwon's method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def modified_newton(f, f_der, f_der2, x, tol, n):\n",
    "    \"\"\"\n",
    "    Arguments\n",
    "    ---------\n",
    "    f: function\n",
    "        function `f` to find root of\n",
    "    f_der: function\n",
    "        Derivative of `f`\n",
    "    f_der2: function\n",
    "        Derivative of `f_der`\n",
    "    x: float\n",
    "        Initial guess of root\n",
    "    tol: float \n",
    "        tolerance\n",
    "    n: int\n",
    "        maximum number of iterations\n",
    "    \"\"\"\n",
    "    numerator = lambda x: f(x)*f_der(x)\n",
    "    denominator = lambda x: f_der(x)**2 - f(x)*f_der2(x)\n",
    "    next_sequence_item = lambda x: x - numerator(x)/denominator(x)\n",
    "    \n",
    "    x_old = x\n",
    "    x = next_sequence_item(x)\n",
    "    n = n - 1\n",
    "    # print(x, f(x))\n",
    "    while abs(x - x_old) > tol and n > 0:\n",
    "        x_old = x\n",
    "        x = next_sequence_item(x)\n",
    "        n = n - 1\n",
    "        # print(x, f(x))\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "## Using modifed newton's method\n",
    "\n",
    "from math import inf\n",
    "\n",
    "f = lambda x: x**3 - 7*x**2 + 11*x - 5\n",
    "f_der = lambda x: 3*x**2 - 14*x + 11\n",
    "f_der2 = lambda x: 6*x - 14\n",
    "x = 0\n",
    "tol = 1e-6\n",
    "n = inf\n",
    "\n",
    "modified_newton(f, f_der, f_der2, x, tol, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Muller's method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cmath import sqrt\n",
    "\n",
    "def muller(f, x0, x1, x2, tol, n):\n",
    "    \"\"\"\n",
    "    Arguments\n",
    "    ---------\n",
    "    f: function\n",
    "        function `f` to find root of\n",
    "    x0, x1, x2: complex/float\n",
    "        Initial 3 guesses of root\n",
    "    tol: float \n",
    "        tolerance\n",
    "    n: int\n",
    "        maximum number of iterations\n",
    "    \"\"\"\n",
    "    # print(x0, x1, x2)\n",
    "    while abs(x2 - x1) > tol and n > 0:\n",
    "        h1 = x1 - x0\n",
    "        h2 = x2 - x1\n",
    "        δ1 = (f(x1)-f(x0))/h1\n",
    "        δ2 = (f(x2)-f(x1))/h2\n",
    "\n",
    "        # Coefficients of parabola (as quadratic eqn)\n",
    "        a = (δ2 - δ1)/(h2 + h1)\n",
    "        b = δ2 + h2*a\n",
    "        c = f(x2)\n",
    "\n",
    "        D = sqrt(b**2 - 4*a*c)\n",
    "        E = b+D if abs(b-D) < abs(b+D) else b-D\n",
    "        \n",
    "        # Root approximation\n",
    "        x = x2 - 2*c/E\n",
    "\n",
    "        # Update sequence variables\n",
    "        x0 = x1\n",
    "        x1 = x2\n",
    "        x2 = x\n",
    "        n = n - 1\n",
    "        # print(x0, x1, x2)\n",
    "    \n",
    "    return x2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "## Using muller's method\n",
    "\n",
    "from math import inf\n",
    "\n",
    "f = lambda x: x**3 - 5*x - 1\n",
    "x0 = 0\n",
    "x1 = 0.5\n",
    "x2 = 1\n",
    "tol = 1e-6\n",
    "n = inf\n",
    "\n",
    "muller(f, x0, x1, x2, tol, n)"
   ]
  }
 ]
}