{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figure 2\n",
    "\n",
    "This notebook recreates Figure 2 in Rein & Tamayo 2016. The figure illustrates the use of second order variational equations together with Newton's method in a simple optimization problem.\n",
    "\n",
    "We start by import the REBOUND, numpy and matplotlib packages. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import rebound\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create an $N$-body simulation of two planets that runs for 10 orbits. The initial semi-major axis of the outer planet is a free parameter. We want to minimize the $x$ position of the outer planet after 10 orbits.\n",
    "\n",
    "We create a function that runs one $N$-body simulation at a time. The function returns the $x$ coordinate of the inner planet as well as the corresponding coordinates of first and second order variational particles at the end of the simulation. The variational particles are initialized with the `vary()` function to correspond to varying the outer planet's semi-major axis. Note that this function is identical to the one for Figure 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def run_sim_var(a):\n",
    "    sim = rebound.Simulation()\n",
    "    sim.add(m=1.)\n",
    "    sim.add(primary=sim.particles[0],m=1e-3, a=1)\n",
    "    sim.add(primary=sim.particles[0],m=1e-3, a=a)\n",
    "    var_da = sim.add_variation()\n",
    "    var_dda = sim.add_variation(order=2,first_order=var_da)\n",
    "    var_da.vary(2,\"a\")\n",
    "    var_dda.vary(2,\"a\")\n",
    "    \n",
    "    sim.integrate(2.*np.pi*10.)\n",
    "    return sim.particles[1].x, var_da.particles[1].x, var_dda.particles[1].x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find the local optimum where $x$ is minimal, we run Newton's method for 14 iterations. We start at $a=1.56$. This is the value indicate by the a red dot in Figure 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N=14\n",
    "a = 1.56\n",
    "x_chain = np.zeros(N)\n",
    "for i in range(N):\n",
    "    x, dxda, ddxdda = run_sim_var(a)\n",
    "    x_chain[i] = x\n",
    "    Delta_a = -dxda/ddxdda\n",
    "    a += Delta_a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To illustrate the fast convergence of Newton's method, we plot the relative error of the inner planet's final position from the *optimum*. We here choose the final value in the iteration as the optimum. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAERCAYAAABVU/GxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYVOWZ/vHvw75EXECEZpG2UdlRhlFUjDTYiriiEgUa\nlygmE4NJDDOamTAw9kxGMZr4c5sAQqIS9xhjECIZaBSME1wQARF/0ICCgLiMUQkiPPNHnSZNU91d\nVV1Vp07V/bmuuug6XXXO3Yj91Dnve97H3B0REZFkNQk7gIiIRJMKiIiIpEQFREREUqICIiIiKVEB\nERGRlKiAiIhISlRAREQkJSogIiKSkmZhB0iGmV0InAscAsx294UhRxIRKVgWxTvRzeww4HZ3nxh2\nFhGRQhXKJSwze8DMtpvZylrbR5rZWjNbZ2Y31bOLHwP3ZjaliIjUJ6wxkDnA2TU3mFkT4J5ge19g\nrJn1Cr43wczuNLMiM7sVeM7dV2Q7tIiI/E0oBcTdlwIf19p8EvCOu29y9z3Ao8CFwesfcvcbgUuA\nEcClZnZdNjOLiMiBcmkQvQvwbo3n7xErKvu5+93A3fXtxMyiN6gjIpID3N2SeX1eTuMdP348GzZs\nwN0j95g6dWroGZQ//BzKH71HlLO7p/a5O5cKyBage43nXYNtSZs7dy5lZWVUVVWlJZiIiBwszAJi\nwaPacqCnmR1tZi2Ay4Hfpbrz9evXM2XKlEZGFBGRuoQ1jffXwEvAcWa22cyudve9wCTgeWA18Ki7\nv9WY42zdurXxYbNs2LBhYUdoFOUPl/KHJ8rZUxXJGwnrU3MQfdy4ccydOzfMOCIikWBmuAbRY5o3\nb87RRx8ddgwRkbyVlwVk/PjxLF68mIceeognnngi7DgiInkpLy9hVf9Mb7zxBmVlZfzmN79h6NCh\nIScTEclduoRVy8CBA3n44Ye59NJLefvtt8OOIyKSVyJVQMysl5ndb2aPmdk1ibznrLPO4ic/+Qnn\nnHMO27dvz3REEZGCEclLWGZmxKb5Xhbnex7vZ5o6dSrz589n8eLFtG3bNhsxRUQiIzKXsBqznLuZ\nnQ/MI7bYYsKmTZtG7969GTt2LHv37k09vIiIANFbzr2zuz/r7qOAq5I5oJkxc+ZMdu3axQ033JDy\n2i8iIhITteXcjzOzu8zsF8DiZI/bokULnnzySV588UXuuOOORv4UIiKFLWrLuS8BljTmIIceeijP\nPfccp5xyCt26deOyyw4aRhERkQTkUgFJm2HDhtGjRw969OjBsGHDDlqjpmvXrvz+97+nrKyMoqIi\nTj/99HCCioiEpLKyksrKSjZu3MjGjRtT2kdos7DM7GjgWXcfEDwfAkxz95HB85sBd/fbktxv3FlY\n8Tz//PNMmDCBJUuW0KtXr+R+ABGRPBKZWViBjC7nnoizzjqLW2+9lVGjRukeERGRJOX1cu6JuPrq\nq7niiis477zz+PzzzzN9OBGRvBHJGwnrk8wlrGruzje/+U127tzJ008/TbNmeTk0JCJSp6hdwsoZ\nZsaMGTPYvXu37hEREUmQCkigefPmPPnkkyxbtozbb7897DgiIjlP12pqaNeuHfPmzePUU0+le/fu\nXH755WFHEhHJWSogtVTfI3LmmWdSVFTE17/+9bAjiYjkJF3CimPAgAHMnTuXMWPG8NZbGZ8IJiIS\nSZErIGbWxsyWm9moTB6nrKyM6dOnM2rUKLZt25bJQ4mIRFIUL2HdBDyWjQNdeeWVbNq0iXPPPZcl\nS5bwta99LRuHFRGJhFDuAzGzB4DzgO3VS5kE20cCPyd2ZvRA7WVMzOxMoD3QCtjp7vPi7Dvp+0Dq\n4+5cc801bNy4kU6dOvH+++/TpUsXKioqKC4uTttxRETClMp9IGEVkKHAZ8CDNdbCagKsA0YAW4kt\nbXK5u681swnAIKAd8L/E+oV84e6j4+w7rQUEYN26dZxwwgns2rVr/7aSkhIWLlyoIiIieSEyNxKm\n2A/kB+5+TdAXZC4wM1t5b7nllgOKB8D69euZMmVKtiKIiOScXBoDabAfSDV3f7C+HTW0nHuytmzZ\nEnf71q1bG7VfEZGwpGM591wqIGlTWVmZ1v116dIl7vaioqK0HkdEJFtqf7g2S+rqFZBb03i3AN1r\nPO8abAtdRUUFJSUlB2wrKSmhoqIipEQiIuEr6H4giSouLmbhwoWMHz+eAQMG0LFjRw2gi0jBC2sW\n1q+BYcSm5G4Hprr7HDM7hwOn8d6awr7TPgurps8++4xOnTqxY8cO2rRpk7HjiIhkU2Sm8WZSpgsI\nwKmnnsott9zCmWeemdHjiIhkS2Sm8Ubd8OHDWbx4cdgxRERCpQKSgtLSUhYtWhR2DBGRUOkSVgp2\n7drFkUceyfvvv88hhxyS0WOJiGSDLmFlSevWrRk8eDBLly4NO4qISGhUQFI0fPhwXcYSkYIWqQJi\nZmeY2Qtmdr+ZhdoqsLS0VAPpIlLQIlVAAAf+ArQktlZWaE466STefvttPv649pqQIiKFIZQCYmYP\nmNl2M1tZa/tIM1trZuvM7Kba73P3F9z9XOBm4JZs5Y2nZcuWnHLKKbzwwgthxhARCU1YZyBzgLNr\nbgj6gdwTbO8LjDWzXsH3JpjZnWbWOXj5J0CLLOaNS9N5RaSQRakfyI3AEDP7L+BXxIpNqHRDoYgU\nslxazr3BfiDu/jTwdEM7Snc/kLr83d/9HZs2beKDDz7gyCOPzMgxREQyQf1A6pDufiB1adasGaef\nfjqVlZWMGTMmK8cUEUkH9QPJAZrOKyKFSv1AGkkD6SJSqMKaxvtr4CXgODPbbGZXu/teYBLwPLAa\neNTd3wojXzIGDhzIjh071B9dRAqOFlNMg9GjRzNmzBjGjRuX1eOKiKSLFlMMidbFEpFCpAKSBhpI\nF5FCpAKSBn379uWzzz5j06ZNYUcREckaFZA0MDOGDRumsxARKSgqIGmi6bwiUmgiNQvLYrdKVgDt\ngOXu/lCc12R9FhbAunXrGDFiBJs3b07pjk4RkTAVwiysC4ndof4lIfcDqe3YY49l3759rF+/Puwo\nIiJZEal+IMDxwDJ3nwx8JythE2Rmms4rIgUlUv1AgK38bRn4vdmLmxhN5xWRQhLaGIiZHQ086+4D\ngudDgKnufk7w/GbA3f22Gu9pDdwNfA6sdff74+w3lDEQgI0bNzJkyBDef/99jYOISKSkMgaSS8u5\nJ9IPZBdwbUM7ylY/kNp69OhB69ateeutt+jTp09Wjikikgr1A6lDtvqBxFM9nVcFRERymfqB5CC1\nuRWRQqF+IGlWWlpKZWUl+/btCzuKiEhGqR9ImnXp0oUOHTqwcuXKhl8sIhJhoYyBuHvcxhnuPh+Y\nn+U4aVc9nfeEE04IO4qISMbk0hhI3tANhSJSCCK1FlYiwrwPpNqOHTs47rjj2LlzJ82a5eVENxHJ\nM4WwFlYkdOzYka5du/L666+HHUVEJGNUQDJEl7FEJN+pgGSI1sUSkXwXqTEQMxsKjCc2e6y3uw+N\n85rQx0AAPvroI3r06MHOnTtp0aJF2HFEROqV92Mg7r7U3f8B+D3wq7Dz1OeII46gZ8+eLF++POwo\nIiIZEbV+INXGAb/ObMrG0ziIiOSzSPUDMbPOZtYN+MTdP8926GRpHERE8lmk+oEE26cBC9z95Tr2\nmxNjIACffvopRUVF7Ny5k1atWoUdR0SkTnnfDwTA3ac1tKOw+oHU1q5dO/r168ef/vQnSktLQ8kg\nIhJPxvuBmFlT4LagB3lkhNkPpLbqy1gqICKSSzLeDyRYIfegqbIZkhf9QGrTQLqI5KsGx0DM7H5i\nl5eeINaLHAB3/02jDmzWg9gYSP/geVPgbWAE8D7wZ2Bssku659IYCMAXX3xBx44d2b59O23btg07\njohIXJm6D6QV8CEwHDg/eJyXfLy/yed+ILW1adOGE088kWXLloUdRUQkrSJ1J3oicu0MBGDq1Kns\n3r2bW2+9NewoIiJxZeQMxMy6mtnTZrYjeDxlZl1Tj1l4dD+IiOSjRMZAFhK76/uhYFM5MN7dyzKc\nLSW5eAby17/+lQ4dOrBlyxYOPfTQsOOIiBwkU2MgR7r7HHf/Knj8EjgypYQFqlWrVpx88sm8+OKL\nYUcREUmbRArIh2ZWbmZNg0c5sUF1SYKm84pIvkmkgHwT+Aawjdj02kuBqzMZKh9pHERE8k29YyDB\nvRk3uPvPshepbmbWBbgb+Ah4p/Y6WcFrcm4MBGDPnj20b9+eqqoq2rdvH3YcEZEDpH0MJLg3Y2yj\nUqXXAOBJd78WOCHsMMlo3rw5p512GkuWLAk7iohIWiRyCWuZmd1jZqeb2aDqR2MO2oh+IMuAb5nZ\nH4EFjckQBl3GEpF8ksg03ni/8dzdh6d80Fhr2s+AB2ss594EWEdsKZOtwHLgcndfa2YTgEHADuBF\nd19qZk+4+5g4+87JS1gAr7zyCldeeSWrV68OO4qIyAFSuYTV0BhIE+BSd3+8seHi7DvpfiBmNgD4\nV+AD4C/u/k9x9puzBWTv3r106NCBtWvXctRRR4UdR0Rkv7T3A3H3fWb2T0DaC0gcDfYDcfeVxGaB\n1StX+oHU1rRpU77+9a+zePFiLr/88rDjiEgBy3g/kMAfzWwy8BgHrsb7UUpHzIJc6gdSW/U4iAqI\niIQpHf1AEikglwV/Xl9jmwPHJH20+uVlP5DaSktLue+++8KOISLSaA0WEHcvztCxLXhUWw70DMZG\n3gcuJ7emEKdF//79+eijj3jvvffo2lVrUopIdCWyGm8bM/uxmc0Inh9rZuoHkqImTZowbNgwTecV\nkchLZBrvY8CrwBXu3s/M2gAvuXtO3siXy7Owqt13330sX76cOXPmhB1FRATI3Gq8Je4+HdgD4O5f\ncOClJ0mSbigUkXyQSAH50sxaExs4x8xKgN0ZTZXnevXqxe7du6mqqgo7iohIyhIpIFOJLRvSzczm\nAv8NHHQDnyTOzCgtLdXy7iISaQ0WEHdfCFwMXAU8Agx298rMxsp/uowlIlHX4CB61ERhEB1g/fr1\nnH766WzZsiWlG3hERNIpU4PoOcPMepvZY2Z2r5ldEnaexjjmmGNo1qwZ69atCzuKiEhKIlVAgHOA\n/+fu1wNXhB2mMcxMbW5FJNISKiBmNtTMrg6+PtLMGnV3eiP6gTwEXG5m04EjGpMhF2gcRESiLJEb\nCacCg4Hj3f04MysCnnD301I+aGr9QE4Ebnf394PXPuXuo+PsOxJjIACbN29m8ODBbNu2jSZNonYy\nKCL5JFNjIKOBCwhW4nX3rcAhycf7G3dfCnxca/NJxPqcb3L3PcCjwIXB6x9y9xuBFmb2C+BXwO2N\nyZALunfvTrt27dRgSkQiKZHVeL90dzez6hsJ22YoSyL9QDYB32poR7naDySe6stY/fv3DzuKiBSQ\nbPUDeTz41H+YmU0EvgnMTOloWZLL/UBqGz58OI899hg33HBD2FFEpICkox9IIjcS/hR4EngKOB74\nV3e/O+kjNawg+oHUNmzYMF544QX27t0bdhQRkaQ0eAZiZjcCjwV3pKdTQfYDqa1z58506tSJN954\ng0GDBoUdR0QkYYkMoh8CPG9mL5rZd83sqMYetJD7gcSjdbFEJIoSXsrEzAYQa297CfCeu5+ZyWCp\nitI03mpPPfUUs2fPZt68eWFHEZEClemlTHYA24APgY7JHETqd8YZZ7B06VL27NkTdhQRkYQl0tL2\nO2ZWSWwZ9/bAxOqb/yQ9OnToQI8ePXj11VfDjiIikrBEpvF2A77v7isyHaaQVa+LNWTIkLCjiIgk\npM4zEDNrF3x5O7DZzI6o+chOvMKhdbFEJGrqHEQ3s9+7+3lmVkWsnW3NwRV392OyETBZURxEB/jk\nk0/o1q0bO3fupGXLlmHHEZECk9ZBdHc/L/iz2N2PCf6sfmS8eJhZsZnNMrPHa2xrY2a/NLNfmNm4\nTGfIpo8//pgWLVpwyimnUF5ern7pIpLzElmN97/dfURD2zLFzB53928EX5cDH7v7PDN71N0vj/P6\nyJ2BVFVVUVZWxvr16/dvKykpYeHChRQXN2rlfBGRhKT1DMTMWgVjHR3M7PAa4x89iC18mGioVHt/\nxNOVvy24mDdrf0yZMuWA4gGxlrdTpkwJKZGISMPqm8b7LeBVoFfwZ/XjGeCeJI4xBzi75oagn8c9\nwfa+wFgz6xV8b4KZ3WlmnatfXuOt7xIrIrW3R9qWLfGX/Nq6dWuWk4iIJK6+MZC73L0YmFxrDGSg\nuydcQBrR+2O3md0PnFDjDOVp4FIzuxd4NuGfMsd16RL/hK6oqCjLSUREEtfgfSDufreZ9QP6AK1q\nbH+wEcdNpPfHR8A/1Nr2BbHl5OsVpX4gABUVFbz88ssHXMY65phjqKioCDGViOSzrPQDCVraDiNW\nQJ4DzgGWAo0pIBkVpX4gAMXFxSxcuJApU6awdetWVq1axeTJkzWALiIZk5V+IMClxPqUb3P3q4GB\nwKFJH+lABdn7oz7FxcU8/PDDLFq0iDvvvJNnnnkm7EgiIvVKpIDscvd9wFfB3ek7iC1vkow6e3+Y\nWQtivT9+l+Q+89Yll1zCK6+8kvJppYhINiRSQF4xs8OItbF9FXgN+FOiB1Dvj+S1bt2acePGMXv2\n7LCjiIjUKeF+IADBPSDt3H1lAy8NTRRvJIznzTff5JxzzmHjxo00a5bImpciIqlL942Eg2o/gCOA\nZsHXkkH9+/enW7duzJ8/P+woIiJx1beYYn1Lw7q7D89MpMbJlzMQgNmzZ/Pb3/6W3/1Ow0Miklmp\nnIEkdQkrCvKpgHz++ed069aNN998s86bDUVE0iEjLW2DFXB/bGYzgufHmtl5qYaUxLVt25bLLrtM\ng+kikpMSmYU1B/gSODV4vgX494wlkgNMnDiRBx54gH379oUdRUTkAIkUkBJ3nw7sgf3LiWR8IcM6\n+oEctC3fDRo0iPbt27Nw4cKwo4iIHCCRAvKlmbUm1pUQMysBdmc0FeDuVe5+bUPbCsF1113HjBkz\nwo4hInKARArIVGAB0M3M5gL/DfxTogdIcz+QgjR27FgWLVrE9u3bw44iIrJfvQXEYqtrrQUuBq4C\nHgEGu3tlEsdIZz8Q6tmWt9q1a8fFF1/ML3/5y7CjiIjsV28BCebDPufuH7r7PHf/vbvvTOYA6ewH\nEnRErN0jpCBcd911zJo1S4PpIpIzElkj4zUz+3t3X57G46baD+SgbfFErR9IIk466SRat25NZWUl\nw4fn5D2cIhIhWekHApwMjDezTcDnxC4fubsPSOmIWRC1fiCJMDMmTpzIzJkzVUBEpNHS0Q8kkQJy\ndsMvSZr6gaSgvLycKVOmsHPnTjp06BB2HBEpcA3OwgrGKQ56JHkc9QNJg8MPP5wLLriABx/M2WaQ\nIlJAEpnG2yjqB5Je1Zex8mW9LxGJLi2mGDHuTp8+fZg5cyZDhw4NO46I5ImMLKYouaXmYLqISJh0\nBhJBO3fupGfPnlRVVXH44YeHHUdE8oDOQApEhw4dGDlyJHPnzg07iogUMBWQiJo4cSIzZszQYLqI\nhEYFJKJKS0v54osv+POf/xx2FBEpUDldQOroCXKhmc0ws0fMrCzMfGFq0qQJ1157rQbTRSQ0kRhE\nN7PH3f0btbYdBtzu7hNrbc/7QfRq27Zto3fv3mzevJlDDjkk7DgiEmE5O4ieoZ4gPwbuTV/K6OnU\nqRPDhw/nkUceCTuKiBSgbF3CSmtPEDO7ldgy8ysynjzHVQ+mi4hkW1YKSJp7gkwCRgCXmtl12cif\ny8rKyvjggw94/fXXw44iIgUmkdV4MyXVniB3A3fXt+N87AdSl6ZNm+4fTL/vvvvCjiMiEZGOfiBZ\nG0Q3s6OBZ6v7iJjZJcDZ7n5d8LwcOMndb2jkcQpmEL3ae++9x4ABA3j33Xdp27Zt2HFEJIJydhC9\nDuoJkiZdu3bltNNO4/HHH2/4xSIiaZLNAqKeIBmkBRZFJNuyNY1XPUEybNSoUWzatIlVq1aFHUVE\nCkQkbiRMRiGOgVSbMmUKn376KXfddVfYUUQkYlIZA1EBySMbN25k8ODBvPfee7Rq1SrsOCISIVEb\nRJc069GjB4MHD+app54KO4qIFAAVkDyjwXQRyRZdwsozX375Jd27d2fJkiUcf/zxYccRkYjQJSyh\nRYsWXHnllcyaNSvsKCKS53L2DMTMioF/AdpVL+UeLLb4PeAI4Hl3fyDO+wr6DATgnXfeYejQoWze\nvJmWLVuGHUdEIiCvzkDcvcrdr621ba27/wOxmw7PCidZ7jv22GPp27cvzzzzTNhRRCSPZbyApLsX\niJmdD8wjtnqv1EGD6SKSadk4A0lrLxB3f9bdRwFXZTp4lI0ePZoVK1awYcOGsKOISJ7KeAFJcy+Q\nM8zsLjP7BbA409mjrFWrVkyYMIEHHjhomEhEJC3C6geSai+QJcCShnZeSP1A6jNx4kRGjBjBtGnT\naN68edhxRCSHpKMfSJgNpTKmsrIy7Ag5oXfv3pSUlDBv3jwuuuiisOOISA6p/eHaLKkJWEB4s7DU\nCyRLNJguIpmSrQKiXiAhufTSS3n55ZfZvHlzRo9TVVVFeXk5paWllJeXU1VVldHjiUj4Mn4jYdAL\nZBjQHtgOTHX3OWZ2DvBzYkXsAXe/NU3HK/gbCWv77ne/S4cOHZg2bVpG9l9VVUVZWRnr16/fv62k\npISFCxdSXFyckWOKSHppOXdUQOJ54403OP/886mqqqJp06Zp3/+4ceN45JFHDto+fvx4Hn744bQf\nT0TSL5UCkpeD6HKggQMH0rlzZ/7whz8watSolPezb98+Nm7cyJtvvsmqVav2P1avXh339Vu3bk35\nWCKS+1RACsRFF13Et771LXr27EmXLl2oqKio8/KSu7Nt2zZWrVp1QLFYs2YN7du3p1+/fvTr149R\no0Zx0003cdttt/HoowcvDFBUVJTpH0tEQqRLWAWgqqqKESNGHDCwXT1Gcdhhh7F69eqDziqaNGlC\n//799xeLfv360bdvXw499NC4+689BtK9e3cqKys1BiISERoDQQUknvLycubOnXvQ9tatW9O0aVP6\n9u17ULHo2LFjUvPCq6qqmDJlClu3buXDDz+kc+fOLFiwIJ0/hohkkMZAJK4tW+LfYjNw4ECWLVtG\nkyaNn81dXFy8f8B8165d9O7dm8rKyoJdBUCkEOTscu5mVmxms8zs8Vrb25jZcjNLfTS4wHTp0iXu\n9pKSkrQUj9pat27N9OnT+f73v8/evXvTvn8RyQ05W0Di9QMJ3AQ8lu08UVZRUUFJSckB20pKSqio\nqMjYMceMGUO7du2YPXt2xo4hIuHKxo2EDwDnAdvdfUCN7SM58EbC2+p4/+M1OhKeSeyGxFbATnef\nF+f1GgOJo+YYRVFRUb2zsNLltdde49xzz2Xt2rVxB99FJHfk5CC6mQ0FPgMerC4gQT+QdcAIYCux\npU0ud/e1ZjYBOBG43d3fN7Mn3H1M8L5/B9oQ6yHyhbuPjnM8FZAccs0119C+fXumT58edhQRqUdO\nFhAAMzsaeLZGARlCbEmTc4LnNwNe8yzEzI4A/gM4E5hV63tXEDsDeS7OsVRAcsi2bdvo168fL7/8\nMj179gw7jojUIUqzsFLqB1Ljew/Wt3P1A8kdnTp1YvLkyUyePJnf/va3YccRkUA6+oGEdQZyCXC2\nu18XPC8HTnL3G9JwLJ2B5Ji//vWv9OnTh5kzZzJixIiw44hIHKmcgagfiGRcq1atuP322/nBD37A\nV199FXYcEUkT9QORrLj44os54ogjmDVrVthRRCRN1A9EsmbFihWMHDmStWvXcthhh4UdR0RqyNlZ\nWNmkApLbJk6cSLt27bjjjjvCjiIiNaiAoAKS67Zv307fvn156aWXOO6448KOIyKBKA2iS4E66qij\nuOmmm5g8eXLYUUSkkVRAJOtuuOEG1qxZw8KFC8OOIiKNoAIiWdeyZUt++tOfalqvSMSpgEgoLrzw\nQo466ihmzJgRdhQRSVHODqKbWTHwL0C7GqvxngFUAKuBR9z9hTjv0yB6RKxcuZKysjLWrl3L4Ycf\nHnYckYKWV4PodfQDceAvQEti62dJhA0YMIDRo0fzb//2b2FHEZEURKofSI1tHYE73b08zut1BhIh\nH3zwAX369OHFF1+kV69eYccRKVi5egYyBzi75oagH8g9wfa+wFgz6xV8b4KZ3WlmnatfHmefnwAt\nMhdZsuXII4/k5ptv5oc//GHYUUQkSRkvIO6+FPi41uaTgHfcfZO77wEeBS4MXv+Qu98I7Daz+4ET\nzOwmADMbbWb/BfyKWAGSPDBp0iTeeecdFixYEHYUEUlCpPqBuPvTwNMN7Vz9QKKlRYsW3HHHHdx4\n442MGDGC5s2bhx1JJO+pH0j8Y2kMJILcnbPOOosLLriASZMmhR1HpODk6hhIPOoHIgcwM372s59R\nUVHBhx9+GHYcEUmA+oFIzujXrx9jxozRtF6RiFA/EMkpO3fupHfv3ixZsoQ+ffqEHUekYGg5d1RA\n8sHPf/5zFixYwPz58zFL6t+ziKQoSmMgInW6/vrr2bhxI/Pnzw87iojUQwVEck7z5s33T+vds2dP\n2HFEpA4qIJKTRo0aRY8ePbjvvvvCjiIiddAYiOSsNWvWMGzYMNasWUOHDh3CjiOS1zSIjgpIvpk0\naRL79u3j3nvvDTuKSF7LqwJSRz8QI9YPpB2w3N0fivM+FZA88uGHH9K7d28WLVpEv379wo4jkrfy\nahZWHf1ALiR21/qX5Gk/kMrKyrAjNEq687dv354pU6Zw4403ko0PBvr7D1eU80c5e6oyXkDM7AEz\n225mK2ttH2lma81sXfVquwk4Hljm7pOB76Q9bA6I+j/CTOT/9re/zYYNGygtLaW0tJTy8nKqqqrS\neoyqqirKy8u56qqrMrp/5a9//1HMH+XsjebuGX0AQ4ETgJU1tjUB/j9wNNAcWAH0Cr43AbgT6Bw8\nf6LG+8YBlwZfP1rH8TzKpk6dGnaERslE/g0bNnjnzp2dWEdKB7ykpMQ3bNiQtv2XlJRo/9p/Tu07\nG/uvKfigF+OgAAAGo0lEQVTdmdTv97BW4x1CbEmTc4LnNwfhb6vxniOA/wDOBGa5+21m1hq4G/gc\nWOvu98c5lmfjZ8qUadOmMW3atLBjpCwT+cvLy5k7d+5B24uKihg0aFCj9//aa6+xdetW7V/7z6l9\n17f/8ePH8/DDDzd6/zXl7CB6tpdzb+w+REQKUbIFJKyGUhmT7F+AiIikRv1AREQkJeoHIiIiKcnG\nNN5fAy8Bx5nZZjO72t33ApOA54HVxGZUvZXpLCIikj4ZLyDuPs7di9y9pbt3d/c5wfb57n68ux/r\n6Wsmlcq9JTnBzLqa2SIzW21mb5pZoycUZJuZNTGz18wskmeTZnaomT1hZm8F/x1ODjtToszsR0Hm\nlWY2Nzizz1nx7g8zs8PN7Hkze9vM/mBmh4aZsT515J8e/NtZYWZPmVm7MDPWp67784Lv/dDM9gUz\nYeuVs3eiJ8vMmgD3AGcDfYGxZtYr3FRJ+Qq40d37AqcA10csP8D3gDVhh2iEu4Dn3L03MBCIxFlx\nMMtxInBiMNOxGbHLwrlsDrH/V2u6Gfijux8PLAJ+lPVUiYuX/3mgr7ufALxD9PJjZl2BMmBTIjvJ\nmwICnAS84+6b3H0P8CixpU8iwd23ufuK4OvPiP3y6hJuqsQF//BGAbPCzpKK4NPi6TXOkL9y909D\njpWoT4kt79PWzJoBbYCDbx7IIe6+FPi41uYLgV8FX/8KuCiroZIQL7+7/9Hd9wVPXyY2OSgn1fH3\nD/Az4B8T3U8+FZAuwLs1nr9HhH4B12RmPYjdvf8/4SZJSvU/vKjeh1MM7DSzOcFluBnBjas5z90/\nBu4ANhObzfiJu/8x3FQp6eju2yH2gQroGHKexvgmEKmWmmZ2AfCuu7+Z6HvyqYDkBTP7GvAk8L3g\nTCTnmdm5wPbgDKr2jLuoaAYMAu5190HAF8QuqeQ8MzsG+AGxpYGKgK+Z2bhwU6VFJD+MmNm/AHvc\n/ddhZ0lU8GHpn4GpNTc39L58KiCRv7ckuPzwJPCQuz8Tdp4knAZcYGYbgEeAUjN7MORMyXqP2Kev\nV4LnTxIrKFEwmNgiox8FMxx/A5wacqZUbDezowDMrBOwI+Q8STOzq4hdyo1aAS8BegBvmFkVsd+f\nr5pZvWeB+VRA8uHektnAGne/K+wgyXD3fw5m2B1D7O99kbtfEXauZASXTt41s+OCTSOIzoSAt4Eh\nZtYq6JkzgmhMAKh9tvo74Krg6yuBXP8QdUB+MxtJ7DLuBe6+O7RUiduf391XuXsndz/G3YuJfaA6\n0d3rLeJ5U0CCT17fJaL3lpjZacB4YLiZvR5chx8Zdq4CcwMw18xWEJuF9ZOQ8yTE3d8AHgReBd4g\n9kthRqihGhDv/jDgVqDMzN4mVgTTMr0/E+rIfzfwNWBh8P/vfaGGrEcd+WtyEriElbMdCUVEJLfl\nzRmIiIhklwqIiIikRAVERERSogIiIiIpUQEREZGUqICIiEhKVEBEajGzpcGfR5vZ2DTv+0e1ni9N\n5/5Fskn3gYjUwcyGAT909/OTeE/T4KbWur7/F3c/JB35RMKmMxCRWszsL8GX/wkMDe4q/l7QMGu6\nmf1P0DRoYvD6M8zsBTN7htgqCJjZ02a2PGgOdm2w7T+B1sH+Hqp1LMzs9uD1b5jZN2rse3GNRlcP\n1Xj9rWa2KsgyPRt/NyI1NQs7gEgOqj4tv5nYGcgFAEHB+MTdTw7WW1tmZs8Hrz2RWDOhzcHzq939\nEzNrBSw3s6fc/Udmdn2w2u8BxzKzS4AB7t4/WMBuuZktCV5zAtAH2BYc81RgLXCRu/cK3p+z3e8k\nf+kMRCRxZwFXmNnrxHq1HAEcG3zvzzWKB8D3gzW1qhsLHUv9TiO2kjHBAnaVwN/X2Pf7HrvevILY\nqqn/C+wys1lmNhrY1cifTSRpKiAiiTNgkrufGDxKajRu+nz/i8zOAIYDJwftTVcArWrsI9FjVau5\nsuteoFkwznISsWXnzwMWJP3TiDSSCojIwap/ef8FqDng/QfgO0HfFszsWDNrE+f9hwIfu/vuoK/9\nkBrf+7L6/bWO9SJwWTDOciRwOvDnOgPGjnuYuy8AbgQGJP7jiaSHxkBEDlY9BrIS2Bdcsvqlu98V\ntBt+Lei7sYP4fbsXAN82s9XEenX8qcb3ZgArzexVd59QfSx3f9rMhhBbjn0f8I/uvsPMeteRrR3w\nTDDGArGOhCJZpWm8IiKSEl3CEhGRlKiAiIhISlRAREQkJSogIiKSEhUQERFJiQqIiIikRAVERERS\n8n9Fht/TbYV2UwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa59562d2b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(6,4))\n",
    "ax = plt.subplot(111)\n",
    "ax.set_xlabel(\"iterations\")\n",
    "ax.set_ylabel(\"relative error\")\n",
    "ax.set_yscale('log')\n",
    "ax.plot(range(N), np.fabs((x_chain-x_chain[-1])/x_chain[-1])+1e-16,\"ro-\", color=\"black\")\n",
    "plt.savefig('paper_test2.pdf',bbox_inches='tight') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Note that the Newton's method will only find a local optimum. We use Newton's method because it is the simplest derivate based optimization method. Better and more reliable derivate based optimization algorithm exist."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}