{
 "cells": [
  {
   "cell_type": "markdown",
   "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",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import array,zeros,exp\n",
    "from numpy.linalg import norm,solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 [2.1893261  1.59847516 1.39390063] 2.449489742783178\n",
      "1 [1.85058965 1.44425142 1.278224  ] 2.5243835363236373\n",
      "2 [1.7801612  1.42443598 1.23929244] 0.4123361696320246\n",
      "3 [1.77767471 1.42396093 1.23747382] 0.014812983951911496\n",
      "4 [1.77767192 1.4239606  1.23747112] 1.8310659761129026e-05\n",
      "5 [1.77767192 1.4239606  1.23747112] 2.4379428075383574e-11\n",
      "6 [1.77767192 1.4239606  1.23747112] 6.661338147750939e-16\n",
      "Root =  [1.77767192 1.4239606  1.23747112]\n"
     ]
    }
   ],
   "source": [
    "x0 = array([1.0,1.0,1.0])\n",
    "x = newton(f,df,x0,debug=True)\n",
    "print('Root = ',x)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}