Newton-Raphson

NR.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I want to solve this system using the NR method\n",
    "$$\\left\\{ \\begin{array}{l} y + \\sin x + 1 = 0\\\\ y^2 + x = 0 \\end{array}\\right.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def F(x,y):\n",
    "    return np.array([ y + np.sin(x) + 1, y**2 + x ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First I calculate the Jacobian\n",
    "\n",
    "$$\n",
    "    J = \\left(\\begin{array}{cc} \n",
    "        \\frac{\\partial f_1}{\\partial x}& \\frac{\\partial f_2}{\\partial x} \\\\ \n",
    "        \\frac{\\partial f_1}{\\partial y}& \\frac{\\partial f_2}{\\partial y}  \n",
    "    \\end{array}\\right) = \n",
    "    \\left(\\begin{array}{cc} \n",
    "        \\cos x & 1 \\\\ \n",
    "        1 & 2 y  \n",
    "    \\end{array}\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def J(x,y):\n",
    "    return np.array([[ np.cos(x) , 1   ], \n",
    "                     [         1 , 2*y ]  ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Newton-Raphson method is an iterative method expressed with this formula:\n",
    "$$ \\vec{r}_{n+1} = \\vec{r}_n - J^{-1}\\cdot F(\\vec{r}_n) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def step(r):\n",
    "    x = r[0]\n",
    "    y = r[1]\n",
    "    \n",
    "    j = J(x,y)\n",
    "    f = F(x,y)\n",
    "    \n",
    "    jinv = np.linalg.inv(j)\n",
    "    rnew = r - jinv.dot(f)\n",
    "    \n",
    "    return rnew"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One starts with some initial guess ([x=0,y=0] in this case) and applies this formula itertively until it converges to some value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0., -1.])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "step([0,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.33333333, -0.66666667])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "step([ 0., -1.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.3861205 , -0.62292371])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "step([-0.33333333, -0.66666667])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.38728575, -0.62232316])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "step([-0.3861205 , -0.62292371])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.38728607, -0.62232312])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "step([-0.38728575, -0.62232316])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.38728607, -0.62232312])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "step([-0.38728607, -0.62232312])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets check wheter these values solve our equations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  3.62356367e-10,  -4.31346558e-09])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(-0.38728607, -0.62232312)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}