3.2 - Finding Roots, Built-in Methods.ipynb

{
 "metadata": {
  "name": "3.2 - Finding Roots, Built-in Methods"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### 3.2 - Solving Roots, Built-in Methods\n",
      "\n",
      "Use the built-in functions for finding the roots."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Functional forms\n",
      "\n",
      "The function of interest is $f: \\mathbb{R}^2 \\to \\mathbb{R}^2$, where:\n",
      "\n",
      "$$f(x,y) = \\begin{bmatrix}\n",
      "x^3(y+3)^2 - sin(x) ~\\\\\\\n",
      "x^2y^2 + x + y - 1\n",
      "\\end{bmatrix}$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import sympy as sp\n",
      "import timeit\n",
      "from sympy.utilities.lambdify import lambdify\n",
      "from scipy import optimize\n",
      "\n",
      "x,y = sp.symbols('x y');\n",
      "# The function of interest\n",
      "f = sp.Matrix(\n",
      "    [  x**3*(y + 3)**2 - sp.sin(x),\n",
      "         x**2*y**2 + x + y - 1      ]\n",
      ");\n",
      "# The matrix of first partials\n",
      "Df = f.jacobian([x,y]);\n",
      "\n",
      "# Evaluatable (by numpy) versions of f(.), Df(.)\n",
      "_f = lambdify([x,y], f, modules='numpy')\n",
      "__f = lambda vec: _f(vec[0],vec[1]).flatten().tolist()[0]\n",
      "_j = lambdify([x,y], Df, modules='numpy')\n",
      "__j = lambda vec: _j(vec[0],vec[1])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 71
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### scipy.optimize.root"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "time1 = timeit.timeit('optimize.root(__f, [0.5,1])', 'from __main__ import optimize, __f', number=1)\n",
      "time2 = timeit.timeit('optimize.root(__f, [0.5,1], jac=__j)', 'from __main__ import optimize, __f, __j', number=1)\n",
      "\n",
      "sol = optimize.root(__f, [0.5,1], jac=__j)\n",
      "print \"scipy.optimize.root took %s seconds (Jacobian not provided)\" % time1\n",
      "print \"scipy.optimize.root took %s seconds (Jacobian provided)\" % time2\n",
      "print sol.x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "scipy.optimize.root took 0.00413799285889 seconds (Jacobian not provided)\n",
        "scipy.optimize.root took 0.00321006774902 seconds (Jacobian provided)\n",
        "[ 0.2689308   0.69603126]\n"
       ]
      }
     ],
     "prompt_number": 72
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### scipy.optimize.fsolve"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "time1 = timeit.timeit('optimize.fsolve(__f, [0.5,1])', 'from __main__ import optimize, __f', number=1)\n",
      "time2 = timeit.timeit('optimize.fsolve(__f, [0.5,1], fprime=__j)', 'from __main__ import optimize, __f, __j', number=1)\n",
      "\n",
      "x,y = optimize.fsolve(__f, [0.5,1], fprime=__j)\n",
      "print \"scipy.optimize.fsolve took %s seconds (Jacobian not provided)\" % time1\n",
      "print \"scipy.optimize.fsolve took %s seconds (Jacobian provided)\" % time2\n",
      "print [x,y]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "scipy.optimize.fsolve took 0.00342202186584 seconds (Jacobian not provided)\n",
        "scipy.optimize.fsolve took 0.00328707695007 seconds (Jacobian provided)\n",
        "[0.26893079771661266, 0.69603126166710938]\n"
       ]
      }
     ],
     "prompt_number": 73
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### scipy.optimize.newton_krylov"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "time = timeit.timeit('optimize.newton_krylov(__f, [0.5,1])', 'from __main__ import optimize, __f', number=1)\n",
      "\n",
      "sol = optimize.newton_krylov(__f, [0.5,1])\n",
      "print \"scipy.optimize.newton_krylov took %s seconds\" % time\n",
      "print sol"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "scipy.optimize.newton_krylov took 0.020437002182 seconds\n",
        "[ 0.26893111  0.69603089]\n"
       ]
      }
     ],
     "prompt_number": 74
    }
   ],
   "metadata": {}
  }
 ]
}