01_00_numerical_differentiation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A primer on numerical differentiation\n",
    "========================\n",
    "\n",
    "In order to numerically evaluate a derivative $y'(x)=dy/dx$ at point $x_0$, we approximate is by using finite differences:\n",
    "Therefore we find: $$\\begin{eqnarray}\n",
    "&& dx \\approx \\Delta x &=&x_1-x_0, \\\\\n",
    "&& dy \\approx \\Delta y &=&y_1-y_0 = y(x_1)-y(x_0) = y(x_0+\\Delta_x)-y(x_0),\\end{eqnarray}$$\n",
    "\n",
    "Then we re-write the derivative in terms of discrete differences as:\n",
    "$$\\frac{dy}{dx} \\approx \\frac{\\Delta y}{\\Delta x}$$\n",
    "\n",
    "#### Example\n",
    "\n",
    "Let's look at the accuracy of this approximation in terms of the interval $\\Delta x$. In our first example we will evaluate the derivative of $y=x^2$ at $x=1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 1e+00   3.0000000000000000   1.0000000000000000\n",
      " 1e-01   2.1000000000000019   0.1000000000000019\n",
      " 1e-02   2.0100000000000007   0.0100000000000007\n",
      " 1e-03   2.0009999999996975   0.0009999999996975\n",
      " 1e-04   2.0000999999991720   0.0000999999991720\n",
      " 1e-05   2.0000100000139298   0.0000100000139298\n",
      " 1e-06   2.0000009999243669   0.0000009999243669\n",
      " 1e-07   2.0000001010878061   0.0000001010878061\n",
      " 1e-08   1.9999999878450576  -0.0000000121549424\n",
      " 1e-09   2.0000001654807416   0.0000001654807416\n",
      " 1e-10   2.0000001654807416   0.0000001654807416\n"
     ]
    }
   ],
   "source": [
    "dx = 1.\n",
    "x = 1.\n",
    "while(dx > 1.e-10):\n",
    "    dy = (x+dx)*(x+dx)-x*x\n",
    "    d = dy / dx\n",
    "    print(\"%6.0e %20.16f %20.16f\" % (dx, d, d-2.))\n",
    "    dx = dx / 10.\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why is it that the sequence does not converge? This is due to the round-off errors in the representation of the floating point numbers. To see this, we can simply type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0002000099999999172"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((1.+0.0001)*(1+0.0001)-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try using powers of 1/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.00000e+00   3.0000000000000000   1.0000000000000000\n",
      "5.00000e-01   2.5000000000000000   0.5000000000000000\n",
      "2.50000e-01   2.2500000000000000   0.2500000000000000\n",
      "1.25000e-01   2.1250000000000000   0.1250000000000000\n",
      "6.25000e-02   2.0625000000000000   0.0625000000000000\n",
      "3.12500e-02   2.0312500000000000   0.0312500000000000\n",
      "1.56250e-02   2.0156250000000000   0.0156250000000000\n",
      "7.81250e-03   2.0078125000000000   0.0078125000000000\n",
      "3.90625e-03   2.0039062500000000   0.0039062500000000\n",
      "1.95312e-03   2.0019531250000000   0.0019531250000000\n",
      "9.76562e-04   2.0009765625000000   0.0009765625000000\n",
      "4.88281e-04   2.0004882812500000   0.0004882812500000\n",
      "2.44141e-04   2.0002441406250000   0.0002441406250000\n",
      "1.22070e-04   2.0001220703125000   0.0001220703125000\n",
      "6.10352e-05   2.0000610351562500   0.0000610351562500\n",
      "3.05176e-05   2.0000305175781250   0.0000305175781250\n",
      "1.52588e-05   2.0000152587890625   0.0000152587890625\n",
      "7.62939e-06   2.0000076293945312   0.0000076293945312\n",
      "3.81470e-06   2.0000038146972656   0.0000038146972656\n",
      "1.90735e-06   2.0000019073486328   0.0000019073486328\n",
      "9.53674e-07   2.0000009536743164   0.0000009536743164\n",
      "4.76837e-07   2.0000004768371582   0.0000004768371582\n",
      "2.38419e-07   2.0000002384185791   0.0000002384185791\n",
      "1.19209e-07   2.0000001192092896   0.0000001192092896\n",
      "5.96046e-08   2.0000000596046448   0.0000000596046448\n",
      "2.98023e-08   2.0000000298023224   0.0000000298023224\n",
      "1.49012e-08   2.0000000149011612   0.0000000149011612\n",
      "7.45058e-09   2.0000000000000000   0.0000000000000000\n",
      "3.72529e-09   2.0000000000000000   0.0000000000000000\n",
      "1.86265e-09   2.0000000000000000   0.0000000000000000\n",
      "9.31323e-10   2.0000000000000000   0.0000000000000000\n",
      "4.65661e-10   2.0000000000000000   0.0000000000000000\n",
      "2.32831e-10   2.0000000000000000   0.0000000000000000\n",
      "1.16415e-10   2.0000000000000000   0.0000000000000000\n"
     ]
    }
   ],
   "source": [
    "dx = 1.\n",
    "x = 1.\n",
    "while(dx > 1.e-10):\n",
    "    dy = (x+dx)*(x+dx)-x*x\n",
    "    d = dy / dx\n",
    "    print(\"%8.5e %20.16f %20.16f\" % (dx, d, d-2.))\n",
    "    dx = dx / 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition, one could consider the midpoint difference, defined as:\n",
    "$$ dy \\approx \\Delta y = y(x_0+\\frac{\\Delta_x}{2})-y(x_0-\\frac{\\Delta_x}{2}).$$\n",
    "\n",
    "For a more complex function we need to import it from the math module. For instance, let's calculate the derivative of $sin(x)$ at $x=\\pi/4$, including both the forward and midpoint differences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.00000e+00   0.2699544827129282  -0.4371522984736194   0.6780100988420897  -0.0290966823444578\n",
      "5.00000e-01   0.5048856975964859  -0.2022210835900616   0.6997640691250939  -0.0073427120614536\n",
      "2.50000e-01   0.6118351194488110  -0.0952716617377366   0.7052667953545546  -0.0018399858319930\n",
      "1.25000e-01   0.6611301360648314  -0.0459766451217162   0.7066465151141266  -0.0004602660724210\n",
      "6.25000e-02   0.6845566203276618  -0.0225501608588857   0.7069916978116613  -0.0001150833748863\n",
      "3.12500e-02   0.6959440534591224  -0.0111627277274252   0.7070780092891873  -0.0000287718973603\n",
      "1.56250e-02   0.7015538499518499  -0.0055529312346977   0.7070995881463560  -0.0000071930401916\n",
      "7.81250e-03   0.7043374663312676  -0.0027693148552800   0.7071049829223881  -0.0000017982641595\n",
      "3.90625e-03   0.7057239167465070  -0.0013828644400405   0.7071063316202526  -0.0000004495662950\n",
      "1.95312e-03   0.7064157978737740  -0.0006909833127736   0.7071066687949497  -0.0000001123915979\n",
      "9.76562e-04   0.7067614018394579  -0.0003453793470897   0.7071067530886239  -0.0000000280979237\n",
      "4.88281e-04   0.7069341196006462  -0.0001726615859013   0.7071067741617298  -0.0000000070248177\n",
      "2.44141e-04   0.7070204574165473  -0.0000863237700003   0.7071067794299779  -0.0000000017565697\n",
      "1.22070e-04   0.7070636210573866  -0.0000431601291609   0.7071067807464715  -0.0000000004400761\n",
      "6.10352e-05   0.7070852015604032  -0.0000215796261444   0.7071067810775276  -0.0000000001090200\n",
      "3.05176e-05   0.7070959914854029  -0.0000107897011447   0.7071067811557441  -0.0000000000308035\n",
      "1.52588e-05   0.7071013863605913  -0.0000053948259563   0.7071067811848479  -0.0000000000016996\n",
      "7.62939e-06   0.7071040837909095  -0.0000026973956381   0.7071067811775720  -0.0000000000089756\n",
      "3.81470e-06   0.7071054324915167  -0.0000013486950309   0.7071067811921239   0.0000000000055763\n",
      "1.90735e-06   0.7071061068563722  -0.0000006743301754   0.7071067811921239   0.0000000000055763\n",
      "9.53674e-07   0.7071064440533519  -0.0000003371331957   0.7071067811921239   0.0000000000055763\n",
      "4.76837e-07   0.7071066126227379  -0.0000001685638097   0.7071067811921239   0.0000000000055763\n",
      "2.38419e-07   0.7071066969074309  -0.0000000842791167   0.7071067811921239   0.0000000000055763\n",
      "1.19209e-07   0.7071067392826080  -0.0000000419039395   0.7071067811921239   0.0000000000055763\n",
      "5.96046e-08   0.7071067616343498  -0.0000000195521977   0.7071067802608013  -0.0000000009257463\n",
      "2.98023e-08   0.7071067690849304  -0.0000000121016172   0.7071067802608013  -0.0000000009257463\n",
      "1.49012e-08   0.7071067765355110  -0.0000000046510366   0.7071067690849304  -0.0000000121016172\n",
      "7.45058e-09   0.7071067690849304  -0.0000000121016172   0.7071067690849304  -0.0000000121016172\n",
      "3.72529e-09   0.7071067690849304  -0.0000000121016172   0.7071067988872528   0.0000000177007052\n",
      "1.86265e-09   0.7071068286895752   0.0000000475030276   0.7071067690849304  -0.0000000121016172\n",
      "9.31323e-10   0.7071068286895752   0.0000000475030276   0.7071067094802856  -0.0000000717062619\n",
      "4.65661e-10   0.7071068286895752   0.0000000475030276   0.7071065902709961  -0.0000001909155515\n",
      "2.32831e-10   0.7071065902709961  -0.0000001909155515   0.7071070671081543   0.0000002859216067\n",
      "1.16415e-10   0.7071075439453125   0.0000007627587649   0.7071065902709961  -0.0000001909155515\n"
     ]
    }
   ],
   "source": [
    "from math import sin, sqrt, pi\n",
    "dx = 1.\n",
    "while(dx > 1.e-10):\n",
    "    x = pi/4.\n",
    "    d1 = sin(x+dx) - sin(x); #forward\n",
    "    d2 = sin(x+dx*0.5) - sin(x-dx*0.5); # midpoint\n",
    "    d1 = d1 / dx;\n",
    "    d2 = d2 / dx;\n",
    "    print(\"%8.5e %20.16f %20.16f %20.16f %20.16f\" % (dx, d1, d1-sqrt(2.)/2., d2, d2-sqrt(2.)/2.) )\n",
    "    dx = dx / 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A more in-depth discussion about round-off errors in numerical differentiation can be found <a href=\"http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h10/kompendiet/kap11.pdf\">here</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Special functions in **numpy**\n",
    "\n",
    "numpy provides a simple method **diff()** to calculate the numerical derivatives of a dataset stored in an array by forward differences. The function **gradient()** will calculate the derivatives by midpoint (or central) difference, that provides a more accurate result. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.  2.  4.  6.  8. 10. 12. 14. 16. 17.]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl41OX19/H3EWRHLTVqBSG4YRGLaFCUusStWKlaF0QBraK4Ii4VFdzXStVq+1NsRFQwohEXqLWCCipWxQQICqLiAgFBCVIRDFvI/fxxyBOWQGYmM5nMzOd1XVzJTCbfuc0lJzfne+5zLISAiIikvu2SvQAREYkPBXQRkTShgC4ikiYU0EVE0oQCuohImlBAFxFJEwroIiJpQgFdRCRNKKCLiKSJhnX5ZjvvvHPIzs6uy7cUEUl506ZNWxpCyKrpdXUa0LOzsykqKqrLtxQRSXlmNj+S1ynlIiKSJhTQRUTShAK6iEiaUEAXEUkTCugiImlCAV1E0s6wYTB58qbPTZ7sz9e5/HzIzobttvOP+fkJeysFdBFJO127Qq9eVUF98mR/3LVrHS8kPx8GDID58yEE/zhgQMKCutXlCLqcnJygOnQRqQtvvAG9e8Pll8Pw4VBQALm5dbyI7GwP4ptr1w7mzYv4MmY2LYSQU9Pr6vRgkYhIXZk92zfFd94JN99cx8F88WJ49tnqgzlASUlC3lYpFxFJK8uW+ceOHWH9eg/mw4dvmVOPu59/9lRKjx7Qpg38+c/QpEn1r23bNiFLUEAXkbTw88/Qvz906QLjx0OfPvDKK3DHHZ5u2TinHjcVFf7PAPDfHH37wmefwZAh/nHECGjWbNPvadYM7r47zgtxCugikvJmzoScHHjySY+ps2dvmjPPzfXHhYVxesNPP4Ubb/Qc+ZQp/txll8E778DXX3uep0MH/62Sl+c5czP/mJfnzyeAcugikrJCgEce8exGq1Z+I/TYY6t/bW5uLfPoq1Z5MB49GqZNgwYN4He/g8aN/et77+1/NtenT8IC+Oa0QxeRlBUCjBvnQXzmzK0H85itWuXbffAAftdd/qZ/+xt8+y38+99w6KFxftPYaYcuIilnyhTYc09o3RpeegmaN/dzO3FRUQHvvec78YIC3/p/9RU0auTBfZdd4vRG8acduoikjPXr4fbb4eij4aab/LmWLeMYzAsKYK+94KijYMwYOPVUePzxqq/X42AO2qGLSIpYuNBT0e++6zc+//73OFz0hx/g+ee91HDPPb0CZd99/abmH//oW/8UUuPvNTMbaWZLzGzWZs8PNLPPzGy2mSWjQ4KIZIipU6FzZ78XOWqUZ0NatozxYmvWwMsve8D+1a/8KOm//uVf69kTJkzw3xgpFswhspTLU0CPjZ8ws1zgFKBzCGF/4P74L01ExHXo4FmQGTOgX78Iv6m6plhr10L79nDaafDhhzBwIBQXw6BBCVx93akx5RJCeNfMsjd7+lLgLyGENRtesyT+SxORTPbZZ3DPPZ7C3mknv/kZscqmWGVl/riyKRbA9df7b4jjjoOG6ZV1jvW/Zl/gCDO7G1gN/DmEEK+SfRHJYCH4AaGBAz2l/cUXcMABUV5kyJCqYF6prAyGDo2qKVaqifXecEOgFdANuA4oMDOr7oVmNsDMisysqLS0NMa3E5FMsHw5nHOOH+E/9FCvLY86mP/ww9abXyWoKVZ9EWtAXwi8FNxHQAWwc3UvDCHkhRByQgg5WVlZsa5TRDLAuefCCy/4+Z033oDdd4/hIo0awfbbV/+1BDXFqi9iDeivALkAZrYv0AhYGq9FiUjmqKjwA5kA997rZYlDh/rBzIitW+enN1et8vKXJ56o06ZY9UUkZYtjgA+ADma20Mz6AyOBPTeUMj4HnBfqclKGiKSF776DE0+sul/ZsSMcfniUF/nmGzjySLjmGi9HBC+FqcOmWPVFJFUuZ2/lS33jvBYRySATJ3rc/ekneOghvxla/Z24bSgogIsu8s+few7OOqvqa3XYFKu+0NF/EalTa9fC4MHeqDAry1vaXnxxDMH8vvs8gP/6115LvnEwz1DpVYQpIvXekiU+9+Hii+HBB7dMdUfsj3+EFSvg1lu3fhM0w2iHLiJ1YtIkvwHapo3Ph3jssSiDeQg+S+5Pf/LP993Xy2EUzP8/BXQRSajK0XDHHutzkwF22y3KiyxbBqef7lOBvv8eVq+O+zrTgQK6iMTNsGGbzu2cOdNT3CNHeili794xXPS99+DAA+HVV+GBB3yoRNOmcVtzOlEOXUTipmtXH8ZcUOCHMi+6yHuY338/XHttDBdctQrOOANatID33/fBobJVCugiEjeVw5h79fL6cjMYO9bvX0bl+++9BKZpU9+Z77sv7LBDQtacTpRyEZG4efddb3F76aXes/z662MI5v/6F+y/v2/rwXflCuYR0Q5dRGqtvNwLTu680+d8rloFN9/sRSm5uf6nRmvWeIH63//uOfNTTkn4utONdugiUisLFsAxx/isz+OO8y61BQVwxx1V6ZeNb5RW6/PPoVs3D+aDBvnwiQ4d6mT96UQ7dBGJWVmZt7ldscJTLIsWwQ03VO3IK3PqhYU17NJLS2HxYk+39OxZJ2tPR1aXPbVycnJCUVFRnb2fiCRGeXnVsJ/nn4eDDoJ99onyIj/9BK+9VlXLWFZWi2Oj6c3MpoUQaizxUcpFRKIyZw4cfLBXr4C3UIk6mBcW+m+Bfv28WyIomMeBArqIRCQEbzOek+OplRYtYrhIRYVXrxx+uHfpmjzZhzZLXCiHLiI1Wr7cm2k9/7zfAB09OoZpQiH4IaGXX/ZaxhEjoFWrhKw3Uymgi0iN3njDUyx33+215VFNE6pk5qeNjj8eLrkkhn65UpNIJhaNNLMlG6YTbf61a80smFm180RFJHVVVHibcfCN9Zw5MGRIBME8Px+ys2G77XxS0MknV3XluugiP3WkYJ4QkeTQnwJ6bP6kme0BnACk9xhtkQz03XfQowccdhjMm+fPRXTjMz/f58nNn+8plpISL0UcPTqRy5UNagzoIYR3gWXVfOlvwGBAs0RF0sjEidC5M0yZAg8/7JvsiA0d6uWHm5szJ27rk62LqcrFzE4Bvg0hzIzgtQPMrMjMikpLS2N5OxGpAyF4frxyNFxRkW+2o8qOlGzlH+xbe17iKuqAbmbNgCHALZG8PoSQF0LICSHkZGVlRft2IlJHzLzV7SWXeJn4/vtHeYHly6Ft2+q/trXnJa5i2aHvBbQHZprZPKANMN3Mop1BIiL1wJgx3moc4K9/9YZaUc2PCAEefdRzMwMGbHlAqFkzL4+RhIs6oIcQPgkh7BJCyA4hZAMLgYNCCN/FfXUikjArV8IFF8A558A//uHPRV18Ujka7vLL/bDQhRdCXp4HdzP/mJcHffrEff2ypUjKFscAHwAdzGyhmfVP/LJEJJGKi/3E51NPwU03xViEMmXKpqPhXn0VdtnFg/e8eV73OG+egnkdqvFgUQjh7Bq+nh231YhIwk2b5pvpnXeGt96KsFd5dV57DRo31mi4ekS9XEQyRGVj1S5dvMVtcXEMwXzhQi9/AW94Pn26gnk9ooAukgHeecebGy5a5Ac4b7/dSxOjMn68F6j37evlMNtvDy1bJmS9EhsFdJE0Vl4Ot97qDbV+/hn+978YLrJ6NVx5pY+Ea9cOxo2LsZmLJJqac4mkqQUL/H7klClw3nleyRL1hnrpUm+mVVzso+Huu8/z5lIvaYcukgaGDdtybufFF8PUqV7B8tRTMWZHWrWCTp28H8tDDymY13MK6CJpoGtXH8Y8YQJ8+60H948+8oEUfftGebGffvLfBiUlnnAfPVpzPlOEUi4iaSA31095nnSSl4KvWwcvvBBDFUthoc/4nD8fjjgiht8GkkzaoYukuBB8+M9ll3lGZPFibzkeVTCvqPDfCIcf7ndS33lHwTwFKaCLpLAVK+Dss31uxH77QZMmcPPN3o9l85z6Nj3wAAwe7JUsxcXQvXvC1iyJo5SLSApr2BA+/xz69/dqwrFjfWeem+s59YKCGnbqa9b4tv7ii2HXXaFfP00TSmHaoYukmIoKeOQR3503beo3P/fdd9PgnZvrjwsLt3KRdev8uGi3bl5nvsMOcO65CuYpTjt0kRTy3Xe+iX7zTX98+eV+YHPw4C1fW7lT38I333ieZupUb3cbNHQsXSigi6SI11/3TfTKld6R9sILY7jI889XjSEqKIAzz4z7OiV5FNBFUsA//+mThDp1grffho4dY7hIebmf9Nx/f3j2WcjOjvMqJdkU0EVSQI8ecPXVPvgnqmlCAJ984iPgdtzRW97+8peep5G0E8mAi5FmtsTMZm303F/N7DMz+9jMXjaznRK7TJHM8+yznuquqPCeWA8+GMNouEce8WOkN97oz+22m4J5GoukyuUpoMdmz70BdAoh/Ab4ArgxzusSyVgrV8L553tjrQUL/CR+RPLzPY2y3Xawxx4eyK+4wlst3nZbAlcs9UWNAT2E8C6wbLPnJoYQyjc8/BAfFC0itVRcDAcfDE8/7QeE3n4bdork37/5+X6zc/5835kvXOijic45p2o0nKS9eOTQLwCej8N1RDJaeTmcdpqf9Zk0CY4+OopvHjoUysq2fP6///Udu2SEWgV0MxsKlAP523jNAGAAQNu2bWvzdiJpadkyP9fTsKGf9Gzb1ud9RqWkJLrnJS3F/KvbzP4E9AT6hLD1kwkhhLwQQk4IIScr6plXIunt7bfhN7/xkXDgY+KiDubjxm39hKc2URklpoBuZj2AwcDJIYRq/p0nIttSXg633OL3K5s181RL1FavhoED4dRT/SZokyabfr1ZM69zlIwRSdniGOADoIOZLTSz/sD/AS2BN8ys2MweS/A6RdLGggV+JP/OO/3k5/Tp0KVLDBfq1Qv+7//gqqu8Q9eIEV7faOYf8/K8VEYyhm0jWxJ3OTk5oaioqM7eT6Q+mj4dTjgBHn44hngbghemN2gA77/vU59POikh65T6w8ymhRByanqdToqK1IFVq+CVV/yg0EEHeXVh8+ZRXuSnn/z8f5s2PkT08MMTslZJXapnEkmwTz+FQw/1kvCPP/bnog7mH33keZmCgggL0yUTKaCLJEjlaLicHG97+9prXtESlcrRcN27+53Ud9+FIUMSsl5JfQroIgly0UX+p3t3mDkTTjwxhot8/bUfGa0cDac0i2yDcugitTBsmLdM2XiQxOTJPikoNxf23tuHT0R9WHPWLO+Vu/fefhf117/WNCGpkXboIrXQtatXD06e7NmRiy6Ck0/25/v08SlvUQXztWvh+uvhgAP8Lip483MFc4mAdugitVA5u/OMM7zd+DffeEniNgczb83XX3sZzEcf+dDmE06I+3olvSmgi9TSmjV+aHPZMujZE8aPj+EiY8dC//6+E3/hBf8NIRIlpVxEauGLL/xm55o1XiL+4YfenyVq69f7aLjiYgVziZkCukgMVqzwj99+650SX30Vhg/39EtlTr1GM2f6WCKAs86CKVM051NqRQFdJEr5+d4q5f33vZrllVd85idU5dQLC7dxgRC8B8shh3gf8zVr/PkGDRK+dklvyqGLRGjlSp/o9vTTXlvepo2XJG4uN3cbN0V/+MFz5ePGwe9/D089BY0bJ3LZkkEU0EUiMGMG9O4Nc+f6OZ9bbvGBFFFZudIbuSxe7BOfr7pK5YgSVwroIhF4/XWPx1GPhgNPsZhBixZwzTXw29/64FCROFMOXWQrli71qhXw1Monn8QQzBcsgGOPhXfe8ceDBimYS8Johy5Sjbff9pOeZvDVV57mbtUqyouMGwcXXOA3PUtLE7FMkU1EMrFopJktMbNZGz3XyszeMLO5Gz7+IrHLFKkbG4+Ga9EC/vWvGO5ZbjwaLjvbE/CqLZc6EEnK5Smgx2bP3QC8FULYB3hrw2ORlLZiRdVouPPOg2nTIhwNl5/vgXu77fzjwIFelnj11V7buM8+CV65iKsxoIcQ3gWWbfb0KcDTGz5/Gjg1zusSqXMtWnjsfeYZePJJf1yj/HwYMMBHEIXgH599Fm691StZVJIodSjWm6K7hhAWb/j8O2DXOK1HpE6tWuUb6blzPV8+cmSUcz6HDoWysk2fKyvz+nKROlbrm6IhhGBmW500bWYDgAEAbdu2re3bicTNp596bfknn3jb8ZgyIyUl0T0vkkCx7tC/N7NfAWz4uGRrLwwh5IUQckIIOVlZWTG+nUj8hACPP+6j4b7/Hv7zH7j88hguMmyYf6yONi+SBLEG9PHAeRs+Pw8YF5/liCTeE0942rtyNFyPzW/5R6q42CdZNG266fPNmsHdd9d6nSLRiqRscQzwAdDBzBaaWX/gL8DxZjYXOG7DY5F6be1a/9inDzz2GEyYALvtFuVFJkyAOXM84f7kkzB1qm/327Xz59q1g7y8KBPxIvFhYWv/ZEyAnJycUFRUVGfvJwI+Gm7YMBg92k9+tmwZw0XWroWbboK//tVb3T73XNzXKbI1ZjYthJBT0+t0UlTS2uLF0K8fvPWW9ymPaf/y1Vc+Gq6w0KdYPPhg3NcpEg8K6JK2/vMfPyC0ciWMGOGn8KNublhY6L1YGjTwMXGnn56QtYrEg5pzSVoKwe9L7rYbFBVVjeuM2gEHeG1jcbGCudR7CuiSVubO9T5YZvDii37PsmPHKC9SXOzDJ5YvhyZN/CZnu3YJWa9IPCmgS0oaNmzLuZ1DhkCnTj43AmDXXbesKNymytFwhx7qQf2bb+K2XpG6oIAuKalr16phzCtXwu9+B/feCx06wF9iKaL94QfvjjhwIBx/vBeoH3hg3Nctkki6KSopqXIY8+mne3pl2TKvZhk5MobRcODDQl9/HR56CK68UqPhJCVphy4pKzcXzj+/KpiPGhVlMC8v9zw5eH35Bx/4RCEFc0lRCuiScpYuhRtvhDff9CB+001eorh5Tn2bFizwKRaVxelt2vgAZ5EUppSLpJTJk/1UfWmpH99/6SXfqVfG5oICf7xNlaPh1q6F4cO1I5e0oR26pITycrj5Zj/js8MO3h2xMphDVU69sHAbF1m92nPlp54K7dvD9OnQt2+drF+kLmiHLimhf39Pr5x/PvzjH9C8+Zavyc2tYXdeVgbjx8M113hJTKNGCVuvSDIooEu9VlHhozoHDfI2t2efHeUFQvCt/B/+AK1awaxZvsUXSUNKuUi9tGoVXHpp1eCJgw6KIZgvX+7fdMYZ3uoWFMwlrSmgS70zezYccojf9GzZMsYOiVOnQpcu3lDrnnvgwgvjvk6R+kYBXeqNELxtSteusGSJn/MZNiyGIpQnn4Tf/tbzNVOmeI1jgwYJWbNIfVKrgG5mV5vZbDObZWZjzKxJvBYmmWfRIrj2Wo/FM2f6cf5tys+H7GxPsmdn+2OAgw/2IRTFxXDYYQletUj9EfNNUTNrDVwJdAwhrDKzAqA38FSc1iYZYs4c2G8/aN3aD2t27Ogxepvy830waFmZP54/H/70J/+8Tx945plELlmkXqptyqUh0NTMGgLNgEW1X5JkivXrPb19wAE+Hg68W2KNwRxg6NCqYF6pvNzTKyIZKuYdegjhWzO7HygBVgETQwgTN3+dmQ0ABgC0bds21reTNLNokfdfmTTJsyOnnBLlBUpKqn9+4cJar00kVcW8QzezXwCnAO2B3YHmZrbFsbsQQl4IISeEkJOVlRX7SiVtTJwInTt7emXECBgzBnbcMcqL7LFH9c9r0yAZrDYpl+OAb0IIpSGEdcBLwOHxWZaks3XrPF8+bVoMo+FWrfJymHvugcaNN/1as2Y+d04kQ9UmoJcA3cysmZkZcCwwJz7LknQzdy48/bR/ftJJHsx//esoL1Jc7LXljz7qNz6feMJHw5n5x7w8f14kQ8Uc0EMIU4GxwHTgkw3XyovTuiSNjB7tJz0HD4YVK/y5qMrCQ/AGLoce6heoHBLapw/Mm+f15vPmKZhLxqtVlUsI4dYQwn4hhE4hhH4hhDXxWpikvhUr4Nxz/U+XLlBU5Cc/o1I5Gu7KK+GEE7xAvcb+uCKZSc25JCHWrvUN9eefw223eZVhTKPhiothwgSNhhOJgAK6xFUIHnMbNfLW4506wZFHRnmR8nJ47z04+mhvgD5vHuy2WwJWK5Je1MtF4qa0FE4+2cfBAVx2WQzBvKTEUyrHHguffebPKZiLREQBXeJi8mSvLZ84ERYvjvEiL78MBx7oaZZRo7wfgIhETAFdojZsWNVA5vJyH9J8zDFeX/7RRz6uM2pXXQWnnQZ77QUzZqhiRSQGyqFL1Lp2rRrIvGyZn+Vp0sTrzDt3jvGie+zhrRbvuUej4URipIAuUcvNhUce8aB+ySV+bP/ll6OsJgzBDwbtuquPh7v22oStVyRTKOUiUVm1yoN4//7QuzfcdZdXE0YVzJcv92++6CK1uRWJIwV0idjs2Z5u+ec//fj+c8/BzTfD8OFVOfUaffih3/h88UVPr4wZk9A1i2QSBXSJSF4e5OR4aeJ998Fbb3kO/Y47/GOvXhEE9VmzfBxRCFWj4SJqfi4ikdDfJonItGlwxBF+8h48iFemWXJz/XFh4Va+ubzcP+6/v5/41Gg4kYSwENNI9djk5OSEoqKiOns/qZ333/eOtAceCGvWwPbbx7Chfv11P2H02muqKxeJkZlNCyHk1PQ67dBlC+vXeynikUfCDTf4c40bRxnM166F666DE0+E5s09zSIiCaWyRdnEokXQt6/nw3v3hscei+EiX30FZ5/tOZhLL4UHHoCmTeO+VhHZlAK6/H9z5viuvKzMS8TPPz/G5obDh/tEixdf9NOfIlInapVyMbOdzGysmX1mZnPMTHe6Utg++3jr8aIiP74fVTBfudJ75YIXp3/8sYK5SB2rbQ79YeD1EMJ+QGc0gi7lzJ3rBzWXLvV+5Y8/HuFouPx8yM72xPruu/tvg549vaFLkyZbH+IsIgkTc8rFzHYEjgT+BBBCWAusjc+ypC6MHu0FKI0a+eZ6550j/Mb8fBgwwHMz4O0Vzfz46PbbJ2y9IrJttdmhtwdKgSfNbIaZjTCz5nFalyTQxqPhDjrIa8u7d4/iAkOHVgXzSiHoGL9IktUmoDcEDgKGhxC6AD8DN2z+IjMbYGZFZlZUWlpai7eTeBk82DfZt90GkyZBmzZRXqCkJLrnRaRO1CagLwQWhhCmbng8Fg/wmwgh5IUQckIIOVlZWbV4O6mNEODHH/3zO+6At9+GW2+FBg2iuEh5Odx7L7RuXf3X27at7TJFpBZiDughhO+ABWbWYcNTxwKfxmVVElelpX6/8qSTPCZnZfkx/qiUlPiMzyFD4Kij/Ajpxpo189NIIpI0ta1yGQjkm9nHwIHAPbVfksTTpEk+dOKtt+Ccc6LckVd66SW/yMcfe578mWe8W1e7dn4ztF07f6wpQyJJpV4uaWrdOs+R33svdOjgrW5jmiY0YoT3Lc/J8YvstVe8lyoiNYi0l4tOiqap1au9A+IFF8DDD3s7laiE4Lvv007zfgA33KDRcCL1nJpzpZl//9unCrVs6a1URoyIMpiH4N907LG+zW/VCm65RcFcJAUooKeJsjK4+GK/+fn3v/tzO+0U5UV+/LFqNFyDBl6wLiIpQymXNDBrlsfh2bO9xvyaa2K4yIcfeofEBQs88T54sKYJiaQY/Y1NIcOGbTnm7fbboUsX78UyYYKPh4v69H1FhR/lB3jvPc+XK5iLpBz9rU0hXbtuOrtz8mS/4XnYYX58/4QTorzg4sXw888evF98EWbMgG7d4r5uEakbCugppHJ252mnedzt1cvj8Lvvwq67Rnmx//zH6xj//Gd/vM8+MSTdRaQ+UUBPIevX+5zP5cth6lRvrlU5qDlia9fCtdfC738Pv/oVXHllQtYqInVPAT1FLFoExx8PN93kFYTXXQejRm2ZU9+mr7/2tooPPgiXX+6/FSJqfi4iqUBVLimgogKOO87jcYsWMG4cHHOMz1/u1cvTMBHt1M387unLL/toIhFJK9qh12Nr1niaZbvt4NFH4YorYPx4D+ZQlVMvLNzGRVauhH/8ww8MtW8PX3yhYC6SptTLpZ764guvLT/zTLjxxhgvMmOGX+TLL+GDD+CQQ+K6RhGpG5H2ctEOvR4aNconCc2fD/vvH8MFQvB6xm7dvCxx0iQFc5EMoIBej6xYAf36wXnnwcEHe235ySfHcKELL4SrroLf/Q6Ki71/uYikPd0UrUdmz/ac+O23+9jOmHqXA5x1lteYDxzoN0JFJCMooCdZRQVMmeKb6G7dvJJlaxPetqq83H8LNG7sdY0nnBDDsVERSXW1TrmYWQMzm2Fmr8ZjQZlkyRLvjnj00VWVKlEH88rRcHfdBfPmef5cRDJSPHLog4A5cbhORnnrLc+KTJoEjzziA4Eikp8P2dley5iVBfvt56Ph8vO9j7lSLCIZq1YB3czaACcBI+KznMxwxx1+6vMXv4CPPoLLLoswDufne1fE+fN9J750qRer3367DwwVkYxW2x36Q8BgoCIOa8kYu+7qo+EKC+E3v4niG4cO9UkWG6uo8BJFEcl4MQd0M+sJLAkhTKvhdQPMrMjMikpLS2N9u5Q3dqxXsIBPFoppNNz8+dV/raSk1usTkdRXmx16d+BkM5sHPAccY2bPbP6iEEJeCCEnhJCTlZVVi7dLTZWj4c48Ex5/PMZ7lj/+6KWIW9O2bczrE5H0EXNADyHcGEJoE0LIBnoDk0IIfeO2sjQwa5YPpcjLg+uvh9dei+Ge5QcfwIEHekOts86CZs02/XqzZnD33XFbs4ikLp0UTZB58zyY//ADTJwIf/lLDKPhAObM8YqW996D557z3w7t2vlvhnbt/HGfPvFevoikIDXnirP166tOeD7yCJxxRgzThBYv9nP/PXp4jqasLMqEu4ikEzXnSoL//tfnRcyY4Y8vv7wWo+HOPdcDuZmCuYhERAE9Dtav94OaRx3lVYQx/aNnzRq45pqq0XDvvLNlvlxEZBvUy6WWFi2Cvn19FNw558Dw4bDDDlFeZNUqOOIImDbNt/X33w9NmiRkvSKSvhTQa2nECB/N+eST3vY2ppP3TZvBkzMSAAAJAUlEQVR6vvymmzRNSERipoAegWHDvGKlcm7nmjVecPL99z5NqE8f2GuvKC+6ciUMGgSXXOIXv+uuuK9bRDKLcugR6NrVhzFPngyff+5ThC64AA44wEsRow7m06f7SKKnnvLtvYhIHGiHHoHcXHj+eTjlFFi92tuP33knnHhilBeqHA13/fXeKXHSJE0TEpG40Q49AqtXw8iRPiJu3Tq48krvkxW1/Hy4+mofDTdzpoK5iMSVAnoEGjXySULNmnkgz8/39EvEVqzwj717e/J93Dj45S8TslYRyVwK6FtR2ZW2pMRLwufOhVdf9XuXBQVVOfVtKi/3ypUOHfwOasOG3o9FQyhEJAGUQ6/GkiVegvj667BsmR/ULCioqnLJzfXHhYVVz21h/nwvTH//fb+D2qJFna1fRDKTAvpm3nwT+vWD//3Pe7Fcemn1G+rc3G0E8xdfhAsv9COkzz4LZ5+d0DWLiIAC+ibGjvVUyn77wYQJUU4TqhQCjB4N++4LY8bAnnvGfZ0iItVRQMdjsJkXn9x4IwwZEkM/rNmz/a5p+/YwapQf3W/UKCHrFRGpTsbfFC0ogCOP9HYqLVv6rIioR8Pl5fnpo0GD/LkddlAwF5E6l7EBvawMBgzwopPycli+PIaL/Pij52guvhh++1sP7CIiSVKbIdF7mNlkM/vUzGab2aB4LiyRPvnEN9QjRniK5d13YbfdorzIZ5/5aLhXXoH77vOSmKgvIiISP7XZoZcD14YQOgLdgMvNrGN8lpU4IXg/rGXLfDTcPfdEOBouPx+ys30cXHa2lyPuv7+Phhs82J8XEUmimG+KhhAWA4s3fL7CzOYArYFP47S2uFq2zEfD7bijF6E0bx7FNKH8fM/PlJX54/nzYeBAT7EcemjC1iwiEo24bCvNLBvoAtTL1oHvvefZkcsu88d77hlFMF+9Gq66qiqYVyori7Ghi4hIYtQ6oJtZC+BF4KoQwk/VfH2AmRWZWVFpaWlt3y4q69d7V8SjjvKik6uvjvICS5f6OLilS6v/eklJrdcoIhIvtQroZrY9HszzQwgvVfeaEEJeCCEnhJCTlZVVm7eLyqJFcNxxcMstflBz+nTIqWlm9pdfwq23wrXX+uOdd/ZSxF12qf71bdvGdc0iIrURcw7dzAx4ApgTQngwfkuKjxBg3jyfIXHuudvoh7VsmTc7Hz0aPvjAX9izZ9Vpo9tug3322TSHDn6I6O67E/8fIiISodrs0LsD/YBjzKx4w5/fx2ldMVmzxvuvVFRA69Y+XajaOZ9r13rxOcCDD3py/aeffNbcggUwfvym39Snj98AbdfOn2/Xzh/36VNn/20iIjWxEEKdvVlOTk4oKipKyLU//9zbjRcXezni8cdv9oIQfNzb6NHek3zUKDjpJPj2W2+veOCBamsrIvWSmU0LIdSUNE79Xi4hwNNPwxVXePuU8eM3C+arVsEDD3gAnzvXX3TqqVWHgFq39j8iIiku5U/DDB4M55/vNzxnzoQ//AE/kv/hh/6Cxo39SGjr1vDEEz5oYswYOPjgpK5bRCTeUn6Hfuqp3gtryHXraPDmBLhmtI9423FHWLjQj4HOnh1D+0QRkdSScgG9osLvY/7wA9x7L3TvDt0XPAdtr4TSUp/VedFFXtrScMN/noK5iGSAeh3Qhw2DrsvfJDf/QigpYUnrLvRsNIHCr3fm9E6fsb53OQ06d4Ldd/ceuOeeCz16qHWtiGSkeh3Quy5/k173dKaA9qxnb3otfJ7/0YqreJAHZ12LTf4bdO7kwfzII5O9XBGRpKrXAT03/0IKaM+ZvMBydiCwHSPoT/8dX4JpX8JeeyV7iSIi9Ub9rnIpKSGXt7mMRymnEdfxV/rzpB8CUjAXEdlE/Q7obdsymaMZzqXczB2M4EImc7R6qIiIVKNep1wm9xmxIYfei1zeJpfJ9KKAgj4zyU324kRE6pl6vUMv3PE4CobMJLfdN2BGbrtvKBgyk8Idj0v20kRE6p206eUiIpKuIu3lUq936CIiEjkFdBGRNKGALiKSJhTQRUTShAK6iEiaqNMqFzMrBebH+O07A0vjuJxUp59HFf0sNqWfx6bS4efRLoSQVdOL6jSg14aZFUVStpMp9POoop/FpvTz2FQm/TyUchERSRMK6CIiaSKVAnpeshdQz+jnUUU/i03p57GpjPl5pEwOXUREti2VdugiIrINKRHQzayHmX1uZl+a2Q3JXk+ymNkeZjbZzD41s9lmNijZa6oPzKyBmc0ws1eTvZZkM7OdzGysmX1mZnPM7LBkrylZzOzqDX9PZpnZGDNrkuw1JVq9D+hm1gB4BDgR6AicbWYdk7uqpCkHrg0hdAS6AZdn8M9iY4OAOcleRD3xMPB6CGE/oDMZ+nMxs9bAlUBOCKET0ADondxVJV69D+jAIcCXIYSvQwhrgeeAU5K8pqQIISwOIUzf8PkK/C9r6+SuKrnMrA1wEjAi2WtJNjPbETgSeAIghLA2hPBjcleVVA2BpmbWEGgGLEryehIuFQJ6a2DBRo8XkuFBDMDMsoEuwNTkriTpHgIGAxXJXkg90B4oBZ7ckIIaYWbNk72oZAghfAvcD5QAi4HlIYSJyV1V4qVCQJfNmFkL4EXgqhDCT8leT7KYWU9gSQhhWrLXUk80BA4ChocQugA/Axl5z8nMfoH/S749sDvQ3Mz6JndViZcKAf1bYI+NHrfZ8FxGMrPt8WCeH0J4KdnrSbLuwMlmNg9PxR1jZs8kd0lJtRBYGEKo/FfbWDzAZ6LjgG9CCKUhhHXAS8DhSV5TwqVCQC8E9jGz9mbWCL+xMT7Ja0oKMzM8PzonhPBgsteTbCGEG0MIbUII2fj/F5NCCGm/C9uaEMJ3wAIz67DhqWOBT5O4pGQqAbqZWbMNf2+OJQNuEDdM9gJqEkIoN7MrgAn4neqRIYTZSV5WsnQH+gGfmFnxhueGhBBeS+KapH4ZCORv2Px8DZyf5PUkRQhhqpmNBabj1WEzyIATozopKiKSJlIh5SIiIhFQQBcRSRMK6CIiaUIBXUQkTSigi4ikCQV0EZE0oYAuIpImFNBFRNLE/wMnwsv8RU3TjgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot\n",
    "\n",
    "y = lambda x: x*x\n",
    "\n",
    "x1 = np.arange(0,10,1)\n",
    "x2 = np.arange(0,10,0.1)\n",
    "\n",
    "y1 = np.gradient(y(x1), 1.)\n",
    "print (y1)\n",
    "\n",
    "pyplot.plot(x1,np.gradient(y(x1),1.),'r--o');\n",
    "pyplot.plot(x1[:x1.size-1],np.diff(y(x1))/np.diff(x1),'b--x');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice above that **gradient()** uses forward and backward differences at the two ends."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFgFJREFUeJzt3X+MZWd93/H3Z9dO2jVpMPHUMf6xi1rLlYNqQCMDhSIIxDEuimlFqV3jblJHKypoIaKKCCtBRbQS/RHSpI6gW3Bx2IkDAlys1gVWDpKDBC5jx4B/QO3S3fUui3fAxIbQKvX62z/uGXF3PD/unfv7nvdLGt17n/Pcc54j29/5+JlzzpOqQpLUHjsmPQBJ0nhZ+CWpZSz8ktQyFn5JahkLvyS1jIVfklrGwi9JLWPhl6SWsfBLUsucNekBrOe8886rPXv2THoYkjQz7r333u9V1UIvfaey8O/Zs4fl5eVJD0OSZkaSo732dapHklrGwi9JLbNl4U9ycZIvJnkoyYNJ3tG0Py/J4SSPNK/nbvD9vU2fR5LsHfYJSJL600vifxp4V1VdDrwMeFuSy4F3A3dV1aXAXc3nMyR5HvA+4KXAlcD7NvoFIUkajy0Lf1WdrKr7mvc/BB4GLgSuBW5tut0KvHGdr/8ycLiqnqiqHwCHgauHMXBJ0vb0NcefZA/wYuAe4PyqOtls+i5w/jpfuRB4rOvz8aZtvX3vS7KcZHllZaWfYUnSTFtagj17YMeOzuvS0miP13PhT/Ic4NPAO6vqqe5t1VnGa6ClvKrqYFUtVtXiwkJPl6JK0sxaLfYJ3HgjHD0KVZ3XfftGW/x7KvxJzqZT9Jeq6jNN8+NJLmi2XwCcWuerJ4CLuz5f1LRJUmstLXWK+9Hmyvu1K+D++Mewf//ojt/LVT0BPgo8XFUf7Np0B7B6lc5e4LPrfP3zwFVJzm3+qHtV0yZJrbOa8t/ylk5x38yxY6MbRy937r4CuBH4RpL7m7b3AB8APpnkJuAo8GaAJIvAW6vq16vqiSS/DXy1+d77q+qJoZ6BJM2A1ZS/VcFfdckloxvLloW/qr4EZIPNr12n/zLw612fbwFu2e4AJWmWLS11pm2O9vxABdi1Cw4cGN2YvHNXkkZk7Vz+ZtLE69274eBBuOGG0Y3Lwi9JQ9bPXD50iv3HP975I++RI6Mt+jClT+eUpFnVz1z+rl2jT/frMfFL0hBsJ+VPouiDiV+SBjYLKb+biV+StmmWUn43E78kbcOspfxuJn5J6sOspvxuJn5J6tEsp/xuJn5J2sI8pPxuJn5J2sS8pPxuJn5JWqN7YZS9e+cj5Xcz8UtSl7UJ//TpzfvPSsrvZuKXJPqfx4fZSvndTPySWq/fZ+XPYsrvZuKX1Fr9pPydOzuPTp7VlN9ty8Sf5BbgDcCpqnph0/YJ4LKmy3OBP6+qF63z3SPAD4HTwNNVtTikcUvSQObxap1e9TLV8zHgZuAPVxuq6h+tvk/yO8CTm3z/NVX1ve0OUJKGqd8VsXbv7qyGNS9FH3pbevHuJHvW29YsxP5m4BeHOyxJGr42p/xug87x/13g8ap6ZIPtBXwhyb1J9m22oyT7kiwnWV5ZWRlwWJL0E/N25+2gBr2q53rgtk22v7KqTiT568DhJN+sqrvX61hVB4GDAIuLizXguCQJMOWvZ9uJP8lZwD8APrFRn6o60byeAm4Hrtzu8SSpH6b8jQ2S+F8HfLOqjq+3Mck5wI6q+mHz/irg/QMcT5J6Ysrf3JaJP8ltwJeBy5IcT3JTs+k61kzzJHl+kjubj+cDX0ryNeB/AP+tqj43vKFL0plM+b3p5aqe6zdo/9V12r4DXNO8/zZwxYDjk6SemPJ75527kubC/v2m/F5Z+CXNtNXpna1uyNq1Cw4dgiNH2l30wcIvaYatTu9sVfRN+Wfy6ZySZk6vj11o+1z+Rkz8kmaKKX9wJn5JM6Gfh6vt3t2Zy9f6LPySpl6/l2oeODD6Mc0yp3okTS1vyBoNE7+kqeQNWaNj4pc0VUz5o2filzQ1TPnjYeKXNHGm/PEy8UuaKFP++Jn4JU2EKX9yTPySxs6UP1kmfkljsZrwd+yAvXtN+ZPUywpctyQ5leSBrrZ/leREkvubn2s2+O7VSb6V5NEk7x7mwCXNju7n61TB6dOb9/cRyqPVS+L/GHD1Ou2/W1Uvan7uXLsxyU7gD4DXA5cD1ye5fJDBSpot/c7jgyl/HHpZevHuJHu2se8rgUebJRhJ8sfAtcBD29iXpBnTzzw+OJc/ToPM8b89ydebqaBz19l+IfBY1+fjTZukOdZPyt+5ExJT/rhtt/B/CPgbwIuAk8DvDDqQJPuSLCdZXllZGXR3kiag12flQyfh33orPPOMc/njtq3CX1WPV9XpqnoG+E90pnXWOgFc3PX5oqZto30erKrFqlpcWFjYzrAkTYjX5M+WbV3Hn+SCqjrZfPz7wAPrdPsqcGmSF9Ap+NcB/3hbo5Q0tbwmf/b0cjnnbcCXgcuSHE9yE/BvknwjydeB1wC/0fR9fpI7AarqaeDtwOeBh4FPVtWDIzoPSWNmyp9dqapJj+FZFhcXa3l5edLDkLQBU/70SXJvVS320tc7dyX1zJQ/H3xWj6SemPLnh4lf0qZM+fPHxC9pQ6b8+WTil/Qspvz5ZuKXdAZT/vwz8UsCTPltYuKXZMpvGRO/1GKm/HYy8UstZcpvLxO/1DKmfJn4pRYx5QtM/FIrmPLVzcQvzTlTvtYy8UtzypSvjZj4pTlkytdmelmB65Ykp5I80NX2b5N8M8nXk9ye5LkbfPdIs1LX/UlcWUUaodWEv2MH7N1rytfGepnq+Rhw9Zq2w8ALq+pvA/8T+K1Nvv+aqnpRryvDSOrfasI/ehSq4PTpzfvv2gWHDsGRIxb9Ntqy8FfV3cATa9q+0KypC/AV4KIRjE3SFvqdxwdTvoYzx/9PgU9ssK2ALyQp4D9W1cEhHE8S/c3jg3P5+omBrupJsh94GljaoMsrq+olwOuBtyV51Sb72pdkOcnyysrKIMOS5lo/KX/nTkhM+TrTtgt/kl8F3gDcUFW1Xp+qOtG8ngJuB67caH9VdbCqFqtqcWFhYbvDkuZa91z+VnbtgltvhWeecS5fZ9pW4U9yNfCbwK9U1bqZI8k5SX5m9T1wFfDAen0lbc5r8jVMW87xJ7kNeDVwXpLjwPvoXMXz08DhJABfqaq3Jnk+8JGqugY4H7i92X4W8EdV9bmRnIU0x7wmX8O2ZeGvquvXaf7oBn2/A1zTvP82cMVAo5NabGkJ9u/vbVoHOin/wAGLvrbmnbvSFDLla5R8Vo80RZzL1ziY+KUpYcrXuJj4pQkz5WvcTPzSBJnyNQkmfmkCTPmaJBO/NGamfE2aiV8aE1O+poWJXxoDU76miYlfGiFTvqaRiV8aEVO+ppWJXxoyU76mnYlfGiJTvmaBiV8aAlO+ZomJXxqQKV+zxsQvbZMpX7Oqp8Kf5JYkp5I80NX2vCSHkzzSvJ67wXf3Nn0eSbJ3WAOXJqnftW8PHXLdW02PXhP/x4Cr17S9G7irqi4F7mo+nyHJ8+gs1fhSOgutv2+jXxDSLDDlax70VPir6m7giTXN1wK3Nu9vBd64zld/GThcVU9U1Q+Awzz7F4g0E0z5mheDzPGfX1Unm/ffpbO4+loXAo91fT7etEkzYTXh79gBe/ea8jUfhnJVT1VVkhpkH0n2AfsALrnkkmEMSxrI2qt1Tp/evL9X7GhWDJL4H09yAUDzemqdPieAi7s+X9S0PUtVHayqxapaXFhYGGBY0mD6nccHU75myyCF/w5g9SqdvcBn1+nzeeCqJOc2f9S9qmmTplI/8/jgXL5mU6+Xc94GfBm4LMnxJDcBHwB+KckjwOuazyRZTPIRgKp6Avht4KvNz/ubNmmq9JPyd+6ExJSv2ZWqgabmR2JxcbGWl5cnPQy1hHfeah4kubeqFnvp6527ai2vyVdb+awetZIpX21m4lermPIlE79axJQvdZj4NfdM+dKZTPyaa6Z86dlM/JpLpnxpYyZ+zR1TvrQ5E7/mhilf6o2JX3PBlC/1zsSvmWbKl/pn4tfMMuVL22Pi18wx5UuDMfFrppjypcGZ+DUTTPnS8Jj4NfVM+dJwmfg1tUz50mhsu/AnuSzJ/V0/TyV555o+r07yZFef9w4+ZLVBP2vfuu6t1J9tT/VU1beAFwEk2QmcAG5fp+ufVtUbtnsctcvSEuzf3/ti57t3w4EDFnypH8Oa438t8L+qqsf/XKVncy5fGo9hzfFfB9y2wbaXJ/lakv+e5Bc22kGSfUmWkyyvrKwMaViadqvz+Dt2wN69zuVL45CqGmwHyU8B3wF+oaoeX7PtrwHPVNWPklwD/F5VXbrVPhcXF2t5eXmgcWn69ZPwwZQvbSbJvVW12EvfYST+1wP3rS36AFX1VFX9qHl/J3B2kvOGcEzNsH6v1gFTvjRMw5jjv54NpnmS/DzweFVVkivp/KL5/hCOqRllypcmb6DEn+Qc4JeAz3S1vTXJW5uPbwIeSPI14PeB62rQuSXNpH5S/s6dkJjypVEZKPFX1V8AP7em7cNd728Gbh7kGJp9Xq0jTRfv3NXIeOetNJ18Vo9GwpQvTS8Tv4bKlC9NPxO/hsaUL80GE78GZsqXZouJXwMx5Uuzx8SvbTHlS7PLxK++mfKl2WbiV89M+dJ8MPGrJ6Z8aX6Y+LUpU740f0z82pApX5pPJn49iylfmm8mfp3BlC/NPxO/AFO+1CYmfpnypZYZOPEnOZLkG0nuT/KsFdLT8ftJHk3y9SQvGfSYGg5TvtROw0r8r6mq722w7fXApc3PS4EPNa+aIFO+1F7jmOO/FvjD6vgK8NwkF4zhuNrE/v2mfKmthlH4C/hCknuT7Ftn+4XAY12fjzdtZ0iyL8lykuWVlZUhDEvrWZ3eOXp08367dsGhQ3DkiEVfmjfDKPyvrKqX0JnSeVuSV21nJ1V1sKoWq2pxYWFhCMPSWqvTO1sVfVO+NN8GLvxVdaJ5PQXcDly5pssJ4OKuzxc1bRqTXv+Ia8qX2mGgwp/knCQ/s/oeuAp4YE23O4B/0lzd8zLgyao6Ochx1TtTvqS1Br2q53zg9iSr+/qjqvpckrcCVNWHgTuBa4BHgR8DvzbgMbWFpaXOH2+PHYMdO+D06c37797dSfmS2mGgwl9V3wauWKf9w13vC3jbIMdR79ZeprlV0d+1Cw4cGP24JE0PH9kwJ/q9GQuc3pHaykc2zIF+bsYCb8iS2s7EP8P6Sfk7d0Jiypdk4p9ZPnJB0naZ+GeMD1aTNCgT/wwx5UsaBhP/DDDlSxomE/+UM+VLGjYT/5Qy5UsaFRP/FDLlSxolE/8UMeVLGgcT/5Qw5UsaFxP/hJnyJY2biX+CTPmSJsHEPwGmfEmTZOIfM1O+pEnbduJPcnGSLyZ5KMmDSd6xTp9XJ3kyyf3Nz3sHG+7sMuVLmhaDJP6ngXdV1X3Nurv3JjlcVQ+t6fenVfWGAY4z80z5kqbJthN/VZ2sqvua9z8EHgYuHNbA5oEpX9I0Gsocf5I9wIuBe9bZ/PIkXwO+A/zLqnpwGMecdqZ8SdNq4Kt6kjwH+DTwzqp6as3m+4DdVXUF8B+A/7LJfvYlWU6yvLKyMuiwJsaUL2naDVT4k5xNp+gvVdVn1m6vqqeq6kfN+zuBs5Oct96+qupgVS1W1eLCwsIgw5qY1ZR/9OjWfXftgkOH4MgRi76k8Rrkqp4AHwUerqoPbtDn55t+JLmyOd73t3vMaWXKlzRLBpnjfwVwI/CNJPc3be8BLgGoqg8DbwL+WZKngf8DXFdVNcAxp45z+ZJmzbYLf1V9CcgWfW4Gbt7uMabV0hLs3w/HjsGOHXD69Nbf2b0bDhyw6EuaPO/c7dPahL9V0TflS5o2PqunR/3O44Nz+ZKmk4m/B/3M44MpX9J0M/Fvop+Uv3MnJKZ8SdPPxL8Br9aRNK9M/Gt4Tb6keWfi72LKl9QGJn5M+ZLapfWJ35QvqW1am/hN+ZLaqpWJ35Qvqc1alfhN+ZLUosRvypekjrlP/KZ8STrTXCd+U74kPdtcJn5TviRtbNA1d69O8q0kjyZ59zrbfzrJJ5rt9yTZM8jxNrNa7BO48UbXvZWkjQyy5u5O4A+A1wOXA9cnuXxNt5uAH1TV3wR+F/jX2z3eZtYuct7L4o6mfEltNUjivxJ4tKq+XVV/CfwxcO2aPtcCtzbvPwW8dnXx9WHav7+/Z+Wb8iW12SCF/0Lgsa7Px5u2dftU1dPAk8DPDXDMdR071ls/U74kTdEfd5PsS7KcZHllZaWv715yyebbTfmS9BODFP4TwMVdny9q2tbtk+Qs4GeB76+3s6o6WFWLVbW4sLDQ10AOHOgU926rE0qmfEk60yCF/6vApUlekOSngOuAO9b0uQPY27x/E/AnVb386bU/N9zQKe67d/9k+cOPf7zzR15TviSdads3cFXV00neDnwe2AncUlUPJnk/sFxVdwAfBT6e5FHgCTq/HEbihhss8JLUi4Hu3K2qO4E717S9t+v9/wX+4SDHkCQN19T8cVeSNB4WfklqGQu/JLWMhV+SWiYjuLpyYElWgB4es7au84DvDXE4s6CN5wztPO82njO087z7PefdVdXTTVBTWfgHkWS5qhYnPY5xauM5QzvPu43nDO0871Ges1M9ktQyFn5Japl5LPwHJz2ACWjjOUM7z7uN5wztPO+RnfPczfFLkjY3j4lfkrSJuSn8W63/O4+SXJzki0keSvJgkndMekzjkmRnkj9L8l8nPZZxSfLcJJ9K8s0kDyd5+aTHNGpJfqP5d/uBJLcl+SuTHtMoJLklyakkD3S1PS/J4SSPNK/nDut4c1H4e1z/dx49Dbyrqi4HXga8rSXnDfAO4OFJD2LMfg/4XFX9LeAK5vz8k1wI/AtgsapeSOcpwCN7wu+EfQy4ek3bu4G7qupS4K7m81DMReGnt/V/505Vnayq+5r3P6RTCNYufzl3klwE/D3gI5Mey7gk+VngVXQedU5V/WVV/flkRzUWZwF/tVnIaRfwnQmPZySq6m46j67v1r1m+a3AG4d1vHkp/L2s/zvXkuwBXgzcM9mRjMW/B34TeGbSAxmjFwArwH9uprg+kuScSQ9qlKrqBPDvgGPASeDJqvrCZEc1VudX1cnm/XeB84e143kp/K2W5DnAp4F3VtVTkx7PKCV5A3Cqqu6d9FjG7CzgJcCHqurFwF8wxP/1n0bNnPa1dH7pPR84J8lbJjuqyWhWLhzaJZjzUvh7Wf93LiU5m07RX6qqz0x6PGPwCuBXkhyhM6X3i0kOTXZIY3EcOF5Vq/9H9yk6vwjm2euA/11VK1X1/4DPAH9nwmMap8eTXADQvJ4a1o7npfD3sv7v3EkSOnO+D1fVByc9nnGoqt+qqouqag+df85/UlVznwKr6rvAY0kua5peCzw0wSGNwzHgZUl2Nf+uv5Y5/4P2Gt1rlu8FPjusHQ+09OK02Gj93wkPaxxeAdwIfCPJ/U3be5olMTV//jmw1ISbbwO/NuHxjFRV3ZPkU8B9dK5g+zPm9A7eJLcBrwbOS3IceB/wAeCTSW6i87TiNw/teN65K0ntMi9TPZKkHln4JallLPyS1DIWfklqGQu/JLWMhV+SWsbCL0ktY+GXpJb5/yEj2gNwya29AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pyplot.plot(x2,np.gradient(y(x2),0.1),'b--o');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More discussion about numerical differenciation, including higher order methods with error extrapolation can be found <a href=\"http://young.physics.ucsc.edu/115/diff.pdf\">here</a>. \n",
    "\n",
    "The module **scipy** also includes methods to accurately calculate derivatives:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 1e+00   2.0000000000000000   0.0000000000000000\n",
      " 1e-01   2.0000000000000004   0.0000000000000004\n",
      " 1e-02   2.0000000000000018   0.0000000000000018\n",
      " 1e-03   1.9999999999998352  -0.0000000000001648\n",
      " 1e-04   1.9999999999992246  -0.0000000000007754\n",
      " 1e-05   2.0000000000020002   0.0000000000020002\n",
      " 1e-06   2.0000000000019997   0.0000000000019997\n",
      " 1e-07   2.0000000000575109   0.0000000000575109\n",
      " 1e-08   1.9999999933961727  -0.0000000066038273\n",
      " 1e-09   2.0000000544584391   0.0000000544584391\n",
      " 1e-10   2.0000001654807416   0.0000001654807416\n"
     ]
    }
   ],
   "source": [
    "from scipy.misc import derivative\n",
    "\n",
    "y = lambda x: x**2\n",
    "\n",
    "dx = 1.\n",
    "x = 1.\n",
    "\n",
    "while(dx > 1.e-10):\n",
    "    d = derivative(y, x, dx, n=1, order=3)\n",
    "    print(\"%6.0e %20.16f %20.16f\" % (dx, d, d-2.))\n",
    "    dx = dx / 10."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to improve the roundoff errors is by simply using the **decimal** package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 1e+00   3.0000000000000000   1.0000000000000000\n",
      " 1e-01   2.1000000000000001   0.1000000000000000\n",
      " 1e-02   2.0099999999999998   0.0100000000000000\n",
      " 1e-03   2.0009999999999999   0.0010000000000000\n",
      " 1e-04   2.0001000000000002   0.0001000000000000\n",
      " 1e-05   2.0000100000000001   0.0000100000000000\n",
      " 1e-06   2.0000010000000001   0.0000010000000000\n",
      " 1e-07   2.0000000999999998   0.0000001000000000\n",
      " 1e-08   2.0000000099999999   0.0000000100000000\n",
      " 1e-09   2.0000000010000001   0.0000000010000000\n",
      " 1e-10   2.0000000001000000   0.0000000001000000\n"
     ]
    }
   ],
   "source": [
    "from decimal import Decimal\n",
    "\n",
    "dx = Decimal(\"1.\")\n",
    "while(dx >= Decimal(\"1.e-10\")):\n",
    "    x = Decimal(\"1.\")\n",
    "    dy = (x+dx)*(x+dx)-x*x\n",
    "    d = dy / dx\n",
    "    print(\"%6.0e %20.16f %20.16f\" % (dx, d, d-Decimal(\"2.\")))\n",
    "    dx = dx / Decimal(\"10.\")"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}