{
 "metadata": {
  "name": "",
  "signature": "sha256:1c1ede3b022ac427065ea22b63e82351d3ce2b6d01ce93081fcbd3cb880814bc"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Newton method for system"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For $x = (x_0,x_1,x_2) \\in R^3$, define $f : R^3 \\to R^3$ by\n",
      "$$\n",
      "f(x) = \\begin{bmatrix}\n",
      "x_0 x_1 - x_2^2 - 1 \\\\\n",
      "x_0 x_1 x_2 - x_0^2 + x_1^2 - 2 \\\\\n",
      "\\exp(x_0) - \\exp(x_1) + x_2 - 3\n",
      "\\end{bmatrix}\n",
      "$$\n",
      "The gradient of $f$ is given by\n",
      "$$\n",
      "f'(x) = \\begin{bmatrix}\n",
      "x_1 & x_0 & -2 x_2 \\\\\n",
      "x_1 x_2 - 2 x_0 & x_0 x_2 + 2 x_1 & x_0 x_1 \\\\\n",
      "\\exp(x_0) & -\\exp(x_1) & 1\n",
      "\\end{bmatrix}\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy import array,zeros,exp\n",
      "from numpy.linalg import norm,solve"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def f(x):\n",
      "    y = zeros(3)\n",
      "    y[0] = x[0]*x[1] - x[2]**2 - 1.0\n",
      "    y[1] = x[0]*x[1]*x[2] - x[0]**2 + x[1]**2 - 2.0\n",
      "    y[2] = exp(x[0]) - exp(x[1]) + x[2] - 3.0\n",
      "    return y"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def df(x):\n",
      "    y = array([[x[1],               x[0],               -2.0*x[2]], \\\n",
      "               [x[1]*x[2]-2.0*x[0], x[0]*x[2]+2.0*x[1], x[0]*x[1]], \\\n",
      "               [exp(x[0]),          -exp(x[1]),        1.0]])\n",
      "    return y"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 35
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def newton(fun,dfun,x0,M=100,eps=1.0e-14,debug=False):\n",
      "    x = x0\n",
      "    for i in range(M):\n",
      "        g = fun(x)\n",
      "        J = dfun(x)\n",
      "        h = solve(J,-g)\n",
      "        x = x + h\n",
      "        if debug:\n",
      "            print i,x,norm(g)\n",
      "        if norm(h) < eps * norm(x):\n",
      "            return x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x0 = array([1.0,1.0,1.0])\n",
      "x = newton(f,df,x0,debug=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0 [ 2.1893261   1.59847516  1.39390063] 2.44948974278\n",
        "1 [ 1.85058965  1.44425142  1.278224  ] 2.52438353632\n",
        "2 [ 1.7801612   1.42443598  1.23929244] 0.412336169632\n",
        "3 [ 1.77767471  1.42396093  1.23747382] 0.0148129839519\n",
        "4 [ 1.77767192  1.4239606   1.23747112] 1.83106597611e-05\n",
        "5 [ 1.77767192  1.4239606   1.23747112] 2.43794280754e-11\n",
        "6 [ 1.77767192  1.4239606   1.23747112] 6.66133814775e-16\n"
       ]
      }
     ],
     "prompt_number": 37
    }
   ],
   "metadata": {}
  }
 ]
}