hermite_div_diff.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hermite interpolation with divided differences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "  return np.sin(x)\n",
    "def f_prime(x):\n",
    "  return np.cos(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hermite_div_diff(x,node_x,node_y,node_y_prime):\n",
    "  '''Build the Hermite osculating polynomial using the divided difference approach\n",
    "\n",
    "  :param x: locations where to evaluate the resulting Hermite polynomial\n",
    "  :type x: numpy array of doubles\n",
    "  :param node_x: x coordinate of the locations to interpolate\n",
    "  :type node_x: list or numpy array of doubles\n",
    "  :param node_y: values f(x) to interpolate\n",
    "  :type node_y: list or numpy array of doubles\n",
    "  :param node_y_prime: values of the first derivative f'(x) to construct the osculating polynomial\n",
    "  :type node_y_prime: list or numpy array of doubles\n",
    "  :return: the osculating Hermite polynomial at the locations x\n",
    "  :rtype: numpy array of doubles\n",
    "  '''\n",
    "  n = len(node_x)-1\n",
    "  # Initialize storage matrix\n",
    "  z = np.zeros(2*n+2)\n",
    "  Q = np.zeros((2*n+2,2*n+2))\n",
    "\n",
    "  # Fill the first two columns of the Q matrix with the \n",
    "  # zero-th and first order divided differences\n",
    "  for i in range(n+1):\n",
    "    # Set the z sequecne\n",
    "    z[2*i] = node_x[i]\n",
    "    z[2*i+1] = node_x[i]\n",
    "    # The first column of Q has the function values\n",
    "    Q[2*i,0] = node_y[i]\n",
    "    Q[2*i+1,0] = node_y[i]\n",
    "    # The second column has the derivatives every other entry\n",
    "    Q[2*i+1,1] = node_y_prime[i]\n",
    "    if(i != 0):\n",
    "      Q[2*i,1] = (Q[2*i,0]-Q[2*i-1,0])/(z[2*i]-z[2*i-1])\n",
    "  \n",
    "  # Complete the matrix with the other divided differences\n",
    "  for i in range(2,2*n+2):\n",
    "    for j in range(2,i+1):\n",
    "      Q[i,j] = (Q[i,j-1]-Q[i-1,j-1])/(z[i]-z[i-j])\n",
    "\n",
    "  # Generate the interpolating polynomial\n",
    "  poly = np.zeros(len(x))\n",
    "  for i in range(2*n+2):\n",
    "    # Compute the product of the x\n",
    "    prod_x = 1.0\n",
    "    for j in range(i):\n",
    "      prod_x *= (x-z[j])\n",
    "    poly += Q[i,i]*prod_x\n",
    "  # Return the polynomial\n",
    "  return poly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd1yV5f/H8dd12AguQEVRcKCJ4SRH5kpzZbYsd64yrUzNr2Va2bLxdWZpRpqjnGmZX1fmSM3KxD1w4EARFMTFEBDO9fsD84dGctAD9zmHz/Px8CHn3Dfnft/c+vb2Pve5LqW1RgghhP0zGR1ACCGEdUihCyGEg5BCF0IIByGFLoQQDkIKXQghHISzURv29fXVQUFBRm1eCCHs0s6dOy9orf1yW2ZYoQcFBREREWHU5oUQwi4ppaL/bZlcchFCCAchhS6EEA5CCl0IIRyEYdfQhRC27/r168TExJCWlmZ0lCLH3d2dgIAAXFxcLP4eKXQhxL+KiYnB29uboKAglFJGxykytNYkJiYSExND5cqVLf6+PC+5KKW+UUrFK6UO/MtypZSaqpSKUkrtU0rVz0duIYQNS0tLw8fHR8q8kCml8PHxyff/jCy5hj4HaH+H5R2A4Bu/BgJf5iuBEMKmSZkb425+7nlectFab1FKBd1hlceBeTp7HN4/lVIllVL+Wuu4fKcRwgZkZmYRHRvP6bh4Eq9cwdndCx//imSaNb9vWoc2a9zcXHF3c8PL04MqQZWoXjWI4u4ueLo6SQEKw1jjGnoF4EyOxzE3nvtHoSulBpJ9Fk+lSpWssGkh7k3chUscO5tAplsJzl66xsLwzzh/LhZzlvnmOv7BodRp/SQA635aiTkr65bXCAxtRM2m7TBnZfLbgqmULVuGCv7+hFSvQlhoTerdVwVnZ6dC3S9H4uXlRXJy8s3Hc+bMISIigi+++MLq23rnnXdo3rw5bdq0YcqUKQwcOBBPT0+rb6egWKPQczsdyXXWDK11OBAOEBYWJjNriEKXmpbOum072RaxlwOHIjkXF0fJcgE0ebI/vsVcCa5enQZ17qdSBX8CyvrhW7ok/n4++JYugZNJ8UHn+ZhMipS0dFKupXPlajLpWmFy8+Jc4hUywuoTExvHrj17+G3bNsKBkAfb0OaRtgSWdMXHOZ3Q4CBMJrljuLBkZmbi7GxZ1b3//vs3v54yZQq9evUqcoUeA1TM8TgAiLXC6wphFdfSr3M0PoUDsVdYMHsmsScO4+ziTNWq1WjxUBOa1AulUZ2auLs4QdsaFr2mt6cH3p4elPMpefO5+yuUoE3tIQCYzWaOnznHX/sPk+Xpw4WUDLZu38nutYvx9fOjccMwOrduRoOQqlLu9yAhIYFBgwZx+vRpILuEmzZtyrvvvktsbCynTp3C19eXtm3bsnz5crKysjhw4AAjRowgIyODb7/9Fjc3N1avXk3p0qXp27cvnTp1IjY2ltjYWFq1aoWvry+bNm1i3bp1jB07lvT0dKpWrcrs2bPx8vIy+CdwK2sU+grgFaXUIqARcEWunwujmc1mftt1kKWrN/BXxC6aPDOIMj6l6NShHUGlOtOqUR3c3VwLbPsmk4ngwPIEB5a/+Vx0qA9ryrvx6+9/sWrNWlauWkNZf39eGzacpiEVKeZm23cRr9wXS9xl696P7l/SnU61y99xnWvXrlG3bt2bjy9evEjnzp0BGDp0KMOHD+ehhx7i9OnTtGvXjsjISAB27tzJb7/9hoeHB3PmzOHAgQPs3r2btLQ0qlWrxqeffsru3bsZPnw48+bNY9iwYTe38eqrrzJp0iQ2bdqEr68vFy5c4MMPP2T9+vUUK1aMTz/9lEmTJvHOO+9Y9edxr/L8E6SUWgi0BHyVUjHAWMAFQGs9A1gNdASigFSgX0GFFSIvqWnpLFq9iR/+9zNnz8bg6upKWFgDejWsQP3qlQx9wzLQ35dB3R5jULfHOJ94mR9/2cqWHfvYeCKJrdGH8bxyikfCahBWK9iwjLbIw8ODPXv23Hz89zV0gPXr13Po0KGby65evUpSUhIAnTt3xsPD4+ayVq1a4e3tjbe3NyVKlOCxxx4DIDQ0lH379t0xw59//smhQ4do2rQpABkZGTRp0sQ6O2hFltzl0j2P5Rp42WqJhLgLadez+PNEIr/sjmL17G8oW64cgwb0pXunVhQvZnvXQMv6lLxZ7nFXrrH9RCLjP/yRJfOvUr1GDfp1fYK2D9a3qcsxeZ1JG8FsNvPHH3/cUtx/K1as2C2P3dzcbn5tMpluPjaZTGRmZt5xO1prHnnkERYuXGiF1AXHdv60CHEXLielMO7Lb+k24mN+PnieahX9+WL8OFbPnsKLXR+1yTK/nX8JD56oF8CKWVPo0bUL586d4833P6XLS2+yZVck2edMIjdt27a95W6XnGfy98rb2/vm2X7jxo3Ztm0bUVFRAKSmpnL06FGrbctapNCFXUrPuM4X3y2nU59XWPrjCkq4agY+FEi/ppVpWjvYps5sLeVT0puRA7ry83fTeb5PLxIvXmRxRAyzfjvJ2UupRsezSVOnTiUiIoLatWsTEhLCjBkzrPbaAwcOpEOHDrRq1Qo/Pz/mzJlD9+7dqV27No0bN+bw4cNW25a1KKP+9Q8LC9MywYW4G7/uimTcxM+5kJBAteBghr/Qmwfr1jQ6ltWlX89kZ/RlNhyOZ/u6HwjyKcbIgb0I9M91spoCERkZSc2ajveztRe5/fyVUju11mG5rW/bb6sLkUNicjqr9sex71QSHl7Fef+FPjza/AG7PBu3hJuLMw9W86VuxZKMO1qOn39Zz7MRO+n17NMM7t5ZPqwk/sEx/yYIh5KZmcX4WYt59uVRHD+fROcGlflh+kc81rKRw5Z5Tp5uzox77XnmT59A1apV+ebbBTz54kh2Rh43OpqwMY7/t0HYtd2Rx3lq0OssWLyUSmV9GdwskObV/XB2Knp/dGtUDuC7SWN5fehLpKSmsWDnOX4+eI7MHMMUiKJNLrkIm5SZmcWkOUtZvOwHPDw8GTXsFZ5p36xInJHficlkovujrXjikWas2n+OTYfjmT9/IUO6P8oDodWNjicMVrT/dgibdCX1Ol9vieJ/6zdTp04dvg+fRNeOLYp8mefk4epMlwYBPH5fcU4e2c+g199m8pylmM1ytl6UyRm6sCnLN/7BrssemJxd+OTd0TSpXl6K/A6a1Arih5mTGfnx58xbsJgde/Yz8a1h+PuWMjqaMID8TRE2IeN6JqMnfs17n0wiPnI7rzxcjab3BUiZW8CvVAm++WQ0z/fpxdEjR+gz4j1OzloAQUFgMmX/Pn++0THv2rhx46hVqxa1a9embt26bN++HYDnn3/+lo/95yU9PZ02bdpQt25dFi9ebLV8y5cvvyXHO++8w/r16632+vkhZ+jCcNFxCbz2/kROHD/Ow61aMu61fgU6cJYjMplMvNzzcRrXC2HZd78w84ftdDQXp4nWmKKjYeDA7BV79jQ2aD798ccfrFy5kl27duHm5saFCxfIyMgAYObMmfl6rd27d3P9+nWrfpoUsgu9U6dOhISEALcOwVvY5PRHGOrPfUfoPWQUMWfOMGLIICa++bKU+T1oEBLMO4snc19cFBN9KjA4IIRk5QSpqTBmjNHx8i0uLg5fX9+b4674+vpSvnz2mDItW7a8OUiXl5cXY8aMoU6dOjRu3Jjz58/f8jrx8fH06tWLPXv2ULduXY4fP05QUBAXLlwAICIigpYtWwLw7rvv0r9/f1q2bEmVKlWYOnXqzdeZN28etWvXpk6dOvTu3Zvff/+dFStWMHLkyJuv27dvX5YuXQrAhg0bqFevHqGhofTv35/09HQAgoKCGDt2LPXr1yc0NNRqnzqVM3RhmN2nL/H9ngSKeZfgsw9GU69mVaMjOQT3k8fppaM4W64aS5zd6BFQk2nnjlPxxpjh92LixIn/eK5Bgwa0bNmSjIwMPv/8838sb9KkCQ8++CDJycl89dVXtywbMWLEHbfXtm1b3n//fapXr06bNm3o2rUrLVq0+Md6KSkpNG7cmHHjxvH666/z9ddf89Zbb91cXqZMGWbOnMmECRNYuXJlnvt5+PBhNm3aRFJSEjVq1GDw4MEcPXqUcePGsW3bNnx9fbl48SKlS5emc+fOdOrUiS5dutzyGmlpafTt25cNGzZQvXp1nnvuOb788subw/T6+vqya9cupk+fzoQJE/L9P47cyBm6KHRms5mpi9ayeMcZ7gsqz/dffiplbk2VKqGAUeeieCvlEue1ple5quyoep/RyfLNy8uLnTt3Eh4ejp+fH127dmXOnDn/WM/V1ZVOnToB2f/AnDp16p62++ijj+Lm5oavry9lypTh/PnzbNy4kS5duuDr6wtA6dKl7/gaR44coXLlylSvnn07aZ8+fdiyZcvN5U899ZTV8v5NztBFoUpLz+DV9yezY0cEvfs/T/+H2uJkkkmVrWrcuOxr5qmpPHUplsCMVEb4VmJIuSC+PHGeelXK3vVL3+mM2tXV9Y7Lvby88jwjz42TkxMtW7akZcuWhIaGMnfuXPr27XvLOi4uLjfHundycspzOFwAZ2fnm7d5pqXdOnFHzqF2/349rXW+xtPPa5ysv7dhaV5LyBm6KDQJl67Qa/hYduyIoEfXLgx79hEp84LQsyeEh0NgIChFA98SzH/zFVr3G8z3e+P56+RFoxNa7MiRIxw7duzm4z179hAYGGiV1w4KCmLnzp0ALFu2LM/1W7duzZIlS0hMTASyZ06CW4fZzem+++7j1KlTN4fc/fbbb3O9XGRNUuiiUETHJdDr1bc4dfIk/3l1MCMHdJVbEgtSz55w6hSYzXDqFBVefIF3nmtH9TJeTFuylg+mfWsXH0JKTk6mT58+hISEULt2bQ4dOsS7775rldceO3YsQ4cOpVmzZjg55T3QWa1atRgzZgwtWrSgTp06vPbaawB069aN8ePHU69ePY4f///xddzd3Zk9ezbPPPMMoaGhmEwmBg0aZJXs/0aGzxUFLiEpnf8u3shvKxfz4aihtG5cz+hIRVaWWTN83Ods3bqVVq1aMOGNl+74D6sMn2us/A6fK6dIokBFnb1A+JbjlCgXyHczJkuZG8zJpJgy5hU6tm/Hpk2befndyWRmZhkdS1iJFLooML/viaT3S69x7ngkL7aoStVy8nF0W2AymfhgWH+effpJ/vzzT14Y/QkZ163zppwwlhS6KBC/7TzI8LfH4enpybAnHsTP2y3vbxKFxmQy8eaLPejfuwdXTd7M336a6/8yDK/MaWqMu/m5S6ELq9u4fS8jxn5E8eLF+WbC+1St5G90JPEvhvR+kjde7MXR+BSmr44gNS39luXu7u4kJiZKqRcyrTWJiYm4u7vn6/vkPnRhVXujzvDmB/+lRKlSzJrwLhXL+hodSeShYeXSpKdf4z9vjGb7hiBmfjIaT/fs/1EFBAQQExNDQkKCwSmLHnd3dwICAvL1PXKXi7h38+fDmDGcuZzGrFa9SChfnNFjX6ZCGR+jk4l8mLV0DdO+nk1ISAgzPxktY+rYKLnLRRSc+fNh4ED+upDE5PsfplhiPB8t+JwKv6w1OpnIpwFdOjCwX28OHjzIS2MnyN0vdkgKXdybMWPYr10YVro8B2IiGfDXj5S4fMEuR/YTMKjbY/Tp0Y3du3YzdsYiuXZuZ+Qaurgnx88n8opfEE7A5+ejKJ2enL3ACiP7CWMM6/s07l7enDD7smJvLJ3rlM/XGCbCOHKGLu7a2fhEBlW8jwxgysUYaqQl///CSpUMyyXu3YtPP0Krmv5sjTzLF4tWGx1HWEgKXdyVtOtZvDX1Wy77+vLptcs0SLn8/ws9PbNH/BN2SylF+/vL4XzuIN/MnsP0BSuMjiQsIIUu8i0zy8z87acpV7cF/506keZTPrk5sh+Bgdkj/dnZVGfin5RSfDjkOUJrhzJz7ncsXr3Z6EgiD3LbosgXs9nMmM+/I8O3Bj0eCqZB4J0H+Rf2LzUtnV7D3uF0dDQT3htNy4a1jY5UpN3zbYtKqfZKqSNKqSil1KhclpdQSv1PKbVXKXVQKdXvXkML2zT28zmsXfU/KmSelTIvIjzd3fjq49GULu3DB5/PJP5qWt7fJAyRZ6ErpZyAaUAHIATorpQKuW21l4FDWus6QEtgolJKPpXgYL5avIqVq9bQskVzhvV+wug4ohD5lSrBtI/G0KRTD+b+cYqktOtGRxK5sOQMvSEQpbU+obXOABYBj9+2jga8Vfa9TV7ARUCGb3MgP2/byVez51IzJITxeYyhLRxTcGB5XnykNldTMxgzbcE/xn0RxrPkb2UF4EyOxzE3nsvpC6AmEAvsB4Zqrf8xdJtSaqBSKkIpFSFjQ9iP81eu8d8Z8yhTthzTP3gdZ+e8Z3cRjqliaU8a+2ay9eeVDPtgil3MelSUWFLouX2i4PZ3UtsBe4DyQF3gC6VU8X98k9bhWuswrXWYn59fvsOKwpeSnsm8P6Np2rkX0z8cTUnvYkZHEgZ7tFl9uj/7NDt2RPDxV/ONjiNysKTQY4CKOR4HkH0mnlM/4AedLQo4CdxnnYjCKBnXM3lr2kIuJ6fxQutaVKlYzuhIwkb8p/+zNG7cmKU/rmDhqk1GxxE3WFLoO4BgpVTlG290dgNu/5TBaaA1gFKqLFADOGHNoKLwvTnhK35du4K6xVOoWNrT6DjChphMJiaPGULlKlX4LHw2J85dMjqSwIJC11pnAq8APwORwBKt9UGl1CCl1N9TWH8APKiU2g9sAN7QWl8oqNCi4M1Y9D82bvqVju3b8cwjTY2OI2yQu5srn7//Oi2e7MP3e86TnC73QRhNPlgk/uHXv/bx2jsfElIzhDnj35Y3QcUdxVxKJXzLCdTF04zp1RZXFxnzryDJeOjCYpdS0nlv0nR8fHz5bOwIKXORp4BSnjTyM7Ni4WxGTww3Ok6RJoUubsrMMrNoRwz123Vl4tiR+JT0NjqSsBOPPlib9u3asmHjJub99IvRcYosKXRx08xVvxOdmELfR+pTu3plo+MIO/PB0P5UCw7m8/Bv2HnomNFxiiQpdAHA3OXr+OqLzyifGUdoQAmj4wg75OzsxJSxI/D0LMYbH0ziaoqM+VLYpNAFe4+c4IuvZ1MtOJiXnm5jdBxhxyqU8eHDUcOoVL8Fy/edkynsCpkUehF3OSmFke9PxN3dg8nyJqiwgmYN7uf5p9oQGZfEur3RRscpUqTQizCtNa9/8gUXLiTw7sghBJTxMTqScBBNqvhQMv0cb40Zw5aI/UbHKTKk0Iuw7ScvkulTlZ7dn6V143pGxxEORCnFwEebUKJ4cd7+9DPOXZBPkhYGKfQiKuZiMqv2xdGsUQNe6/O00XGEAypV3IuPx7xGSnIyr30oIzMWBin0IuhqSioDhr3J+WO7eSasItnD2AthfQ/cH8yAPj2JPHSIyXOWGh3H4UmhF0FvfPolCefP07XZ/Xi5yce0RcF68dlHqVe/Pn8eO8eZi6lGx3FoUuhFzNzl6/jzzz95svNjct1cFAqTycQX740krEVbFu04Tdr1LKMjOSwp9CLkyMkYps+cS9Vq1XjzxR5GxxFFiKebM90eqMjxqGOMmhAu19MLiBR6EZFl1oSv3IaTszPjRw+V+81FoQv0KUagcxKbN23k2xXrjY7jkKTQi4gNkedxrRDC5PEfUzlAZh4Sxhg5oCtVq1Vj+qx5HIu+feIzca+k0IuAX3fsZ9nGv2gQWIpGweWNjiOKMGdnJz4d9SpKKd74+DMyM+V6ujVJoTu4S1eTeW/C5xz/Yy2PhpY1Oo4QVK3kz5CB/Th54gRffr/W6DgORe5Zc3BjJoZz5fJlpn40Fg9XF6PjCAFAz04PE5dsJs61DGcupsqctVYiZ+gObMnazfzxxx888VgnHmpQy+g4QtziladbUdzDlXm/HiQp5ZrRcRyCFLqDOptwkSkzvqFipUqMklsUhQ1yd3Gi432l+Gn257z3+Wyj4zgEKXQHpLVm3ZHLVAtrycdvDpVJe4XNql25LK2aP8iGjZtYs3WH0XHsnhS6A/rrxAWOxqcwqNtj1Kpayeg4QtzR2CH9KFO2HJ9M/YrEy0lGx7FrUugO5vjpOEa8/ibuSWdpUlXGNxe2z9PdjfdGvkxS0lXenhRudBy7JoXuQMxmM6PHT+N6+jV6taojoygKu9G49n080akTpxJTOHj2stFx7JYUugOZtmAFR48coX/vHlSt5G90HCHy5c1BPWn/VHd+2htHakam0XHskhS6gzgWHcu3C5dQMySEF57pYHQcIfLNxdmJLg0COBcXy3tffGt0HLskhe4AtNZMX/oLymTi45EvYzLJYRX2qXxJD8pej2Pd2tUsXbfF6Dh2R/7mO4C/Tl7EvXJ9PvrwAwIrlDE6jhD3ZOSAblSoEMBnX83hwqWrRsexK1Lodu7EmXMs3ryXqn7FeLh2kNFxhLhnri7OjB3xEikpyYz9bKbRceyKFLodM5vNjJkwnT+Wz6VTLT+5q0U4jAfuD6bzox35/fc/WPPbTqPj2A2LCl0p1V4pdUQpFaWUGvUv67RUSu1RSh1USm22bkyRm7nLf+FwZCQ9n32KcqW9jY4jhFW9MbAHjR7uyK5LrjJtnYXyLHSllBMwDegAhADdlVIht61TEpgOdNZa1wKeKYCsIoez8YmEz51P5SpVeKn7Y0bHEcLqPNxcGf1CV1IyYfW+s0bHsQuWnKE3BKK01ie01hnAIuDx29bpAfygtT4NoLWOt25Mcbuxk7/mekYG7/9nsNzVIhxWJR9PqntdZ+LHH7AlYr/RcWyeJU1QATiT43HMjedyqg6UUkr9qpTaqZR6LrcXUkoNVEpFKKUiEhIS7i6x4ODZK1wye/Ds009yf7Ugo+MIUaCebHIfrs4mxn32FWnpGUbHsWmWFHpu77Tp2x47Aw2AR4F2wNtKqer/+Catw7XWYVrrMD8/v3yHFZB2PYsV+2Jp/khHXusnV7aE4yvh5cmIwQOIP3+e8TMXGx3HpllS6DFAxRyPA4DbZ3eNAdZqrVO01heALUAd60QUOX0cvpATRw/zVP0KODvJpRZRNHRq0ZBGjRqyfOVK9h45YXQcm2VJI+wAgpVSlZVSrkA3YMVt6/wENFNKOSulPIFGQKR1o4rfdh7kfytWUCLtnEzZJYqcd4e9gJubO18sXInWt18kEGBBoWutM4FXgJ/JLuklWuuDSqlBSqlBN9aJBNYC+4C/gJla6wMFF7voybieybipX1G8RAneeinXtyiEcGjlfEry3tuj8QltyY5Tl4yOY5MsmspGa70aWH3bczNuezweGG+9aCKnKXOXci4ujrdHDqOkdzGj4whhiDb1q3Mq5SQ//HmEAK8alPcrbXQkmyIXYe3A8TPn+P7Hn6hbty5PPdLU6DhCGEYpRcdafmxe/BXvTvna6Dg2Rwrdxmmt2Rydyv0PdeDdoS8YHUcIw1Xw8aZTu9bs2BHBqs1/GR3Hpkih27gDZ69wLD6F55/pKCMpCnHDiP5d8fMrw6QZ35Calm50HJshhW7DLielMHzUW2TFn6BJFZkfVIi/ubu5MuKl/lxMTGSC3Jt+kxS6Dfvoy3lcTIjnicbVMZlkJEUhcmrXtAGNGjVkZ9RZ4q+mGR3HJkih26gdB46xfsMmWrZozoN1axodRwib9Ombr1K/dWdW7I2Ve9ORQrdJZrOZcVPD8SxWjDEv9TE6jhA2q4SnG21DyrH78HF+3Ljd6DiGk0K3Qcs2bif61CkG9OqGT0kZ51yIO2kYVIqzO9YxZXo4l5NSjI5jKCl0G5OakcnhtBJ06f8SfZ54xOg4Qtg8JycTr7/Un+TkZD79ar7RcQwlhW5jftx+jLTrWQzo0ETGORfCQg/WrUnz5s34ef169hw+bnQcw0hj2JDf90Qy6aP38M+Kp1wJd6PjCGFXxrz0HO7uHnw4dSZms9noOIaQQrcRZrOZT6bNolixYvRq19joOELYHb9SJRjwXA8o4c+u6ItGxzGERYNziYI3a9lazkRHM/zlgTL4lhB3qf+TbUkrVZV1kQncH1AKdxcnoyMVKjlDtwEJl64we/5iqlStRq/HWhsdRwi7pZSic53ynDhyiInfLDE6TqGTQrcB363bTnpGBqOH9Jc3QoW4RxVLe+Jx9QzLfvyR/UdPGh2nUEl7GCzuyjUSXP15/e33aRASbHQcIRzC6Jf74ObmzrgvvilSb5BKoRvIbDYze+12PFyc6NygitFxhHAY/r6l6PHMUxw5fJil67YaHafQSKEbaPGazXw/axrBLpfxcC1ab94IUdAGd3uMcv7lmTFnQZEZYlcK3SBXU1L5au4CKlQIoEvrRkbHEcLhODs78Z+X+lGh9kP8cbJozEEqty0aZPLs77ly+TLvjXwVZ2c5OxeiILRuVJdzlGLz0QuEBfpQwtPF6EgFSs7QDXD8dBz/W7WGRo0a0uKBUKPjCOHQOt7vT/TBCMZMmG50lAInhW6ApdsO4OZZjDcH9zU6ihAOr1QxV6qWcGLbb7+xafteo+MUKCn0QnYiIZnLbuV474MPCSzvZ3QcIYqE4f26UKJkSSbOmENmZpbRcQqMFHohyszMYvr36yju7kSLGmWNjiNEkeHt6cELz3Xn7NkY5i5fZ3ScAiOFXohmLVvLhuULqe5yCRcn+dELUZi6d2xJUOUqzFn4PUmpjjkHqbRKIbl0NZlvFy+larVqPNm6idFxhChyTCYTo14ZQEirJ9l24rLRcQqE3LZYSCbMXERqSgqvD+4n47UIYZBGodWJTvdgW9QFHggqhY+Xm9GRrEqapRAcORnD2l9+4cEHm9AwtLrRcYQo0trVKsexvzby+rjJRkexOin0QrB6zylK+pbjjUG9jY4iRJFXwsOFOoE+7Nq5kw3b9xgdx6qk0AvYiYRk4nUJxoweTcWyvkbHEUIArz6XfRvjlPB5DjUaoxR6AcrMzGLSnGV4OplpVl3uORfCVhTzcGNAr67EnDnDvJ9+MTqO1VhU6Eqp9kqpI0qpKKXUqDus94BSKksp1cV6Ee3XvBW/sO2XlVRxuii3KQphY3p2epiKgYHMWbSMaxnXjY5jFXm2jFLKCZgGdFzRGIgAABDnSURBVABCgO5KqZB/We9T4Gdrh7RHV1NSmbPgewKDgujesYXRcYQQtzGZTIx+9UVqt+/B78cdY1JpS04bGwJRWusTWusMYBHweC7rDQGWAfFWzGe3Pp+3jKSrV3lt4HNym6IQNqpxaDCNQ6qw5WgCF5Pt/8NGljRNBeBMjscxN567SSlVAXgSmHGnF1JKDVRKRSilIhISEvKb1W5ExyXw08o11G9Qn+ZhMpqiELasXa2y7Fi7hLETvzQ6yj2zpNBVLs/p2x5PAd7QWt9x1ButdbjWOkxrHebn57hvEm44GEup8oG88eJzRkcRQuTBx8uNOtUqsu33P4g4eMzoOPfEkkKPASrmeBwAxN62ThiwSCl1CugCTFdKPWGVhHYm7so1opJMDH31VaoHVcj7G4QQhvvPgG64e3gwfsZcu76N0ZJC3wEEK6UqK6VcgW7AipwraK0ra62DtNZBwFLgJa31cquntQPjZy5BpyXRqkYZo6MIISzkU9Kbrk89ztEjR1i1ZYfRce5anoWutc4EXiH77pVIYInW+qBSapBSalBBB7Qnqzb/xaafV1Lu+jmZ9FkIO/Ni1074+vkRPn8ZZvPtV5Xtg9LamOBhYWE6IiLCkG0XhMzMLB4bMJyMjAxWzZmKu5ur0ZGEEPm0fsch1kYl82zjajSsXNroOLlSSu3UWofltkzup7OSucvXcS4ujoG9u0uZC2GnWofVpFp5H9YdjCMlLd3oOPkmhW4FV1NSmbdoGUGVK/NM+2ZGxxFC3CWlFA8Hl2LdghlMmb3E6Dj5JoVuBduOxlO8QlWGv9BbPkQkhJ2rXr401Sr5s2LVWs7GJxodJ1+kfe5Rcnom28+k0LV7D/kQkRAOYuSLvbmeeZ0JXy8wOkq+SKHfoynzfiD+7Gna1SpndBQhhJXUrFyRh1u1YvOWrew/dtLoOBaTQr8Hh06c5odlS3G/fBw/b8eaykqIom7EgG64uLoyfcH/jI5iMSn0ezAh/DtcnF0YMaC70VGEEFbm71uK4cNfo3hIc6ITU4yOYxEp9Lv0286D7N61m86PtqdCGR+j4wghCsCTzetQspgbP+2MtoshAaTQ74LZbGbyzG/x8vZmyHNPGx1HCFFA3JydqOenmD/tvyxdt9XoOHmSQr8LkXFXcfKpRJ/uz1K8mKfRcYQQBahNvWBKFvfm628Xk3E90+g4dySFnk9ms+aXyHgaNW9Nv6faGR1HCFHAnJ2dGNy3OxcSEvhm2Rqj49yRFHo+LVyzmT27dtE2pCxOptyGihdCOJrOrRpTpWpV5n//I1dTUo2O86+k0PMhNS2d8DnfcfFYBLXKFzc6jhCikJhMJoYO6ElyUhJzV2wyOs6/kkLPh2nzl3P1yhWGDugpH/EXoohpHhbKc4OHccEzkJR027yWLq1kocTLSfywYhW1atWideN6RscRQhigZ+sGXDdr1u0/k/fKBpBCt9CUOUtIT0tj+As9jY4ihDBImeLulEo9y8fvvc3x03FGx/kHKXQLXLl2nXNZXrRu05YGIcFGxxFCGOjZVg3Q5iwmzLS9gbuk0C2w8fB5ylYN4c2XnjM6ihDCYIHl/WjXpg3bt29nd+Rxo+PcQgo9DweiTrH4x5U0CChO6WIyE5EQAob1ewY3d3cmz5pvdJRbSKHnYdLMBZzctZXGQXKbohAim2/J4jzxWEf27z9ARKTtDK/rbHQAW/b7nkh279rNU493xt+3lNFxhBA25OUeT3DZzZ+IeDMN7tMoZfwHDeUM/Q4+mzUfz2LFGPLcU0ZHEULYGC9Pd55oej/RiansO20bU9VJof+Ltb9FcPTIEZ59sjMlvYsZHUcIYYPCgkoTt+dX3hz7AVlZxg+vK4WeC601O86kUKVmKC927WR0HCGEjXIyKdo2vJ+42Fjmr9xgdBwp9NwcjL1KqltpRr82BHc3ubNFCPHvunVsQYUKAcxdtIy09AxDs0ih3yYzM4vJsxbibcqgXsWSRscRQtg4k8nEoL7duZiYyKylxg6vK4V+m7nL17Hrtw0Ee6RgkuFxhRAW6NgsjKrVqrF05VrSDJwEQwo9h9S0dL5d8gMVAwN5+pGHjI4jhLATJpOJUUNeoF6nPvxx4qJxOQzbsg36eskqrly+zMt9u8vwuEKIfAmrWYW6Vcqx5UgCSanphmSQ1rrhclIKS35cQfUaNWjXtIHRcYQQdqhVdR82L53J+K+NGRLAokJXSrVXSh1RSkUppUblsrynUmrfjV+/K6XqWD9qwdp65Bwlyldh6AAZHlcIcXcq+ngTXMmftet+4Wx84X/YKM9CV0o5AdOADkAI0F0pFXLbaieBFlrr2sAHQLi1gxak1IxMdsam0b13Hx6sW9PoOEIIOzZsQHcys7KYPHtJoW/bkjP0hkCU1vqE1joDWAQ8nnMFrfXvWutLNx7+CQRYN2bBmrFkDRfOx9I2pKzRUYQQdq5m5Yq0aNaMTb/+yrHo2ELdtiWFXgHIOd9SzI3n/s0AINebMZVSA5VSEUqpiISEBMtTFqDouAQWLlxA1pl9lC3ubnQcIYQDGNa/GyaTienzfyzU7VpS6LndjK1zXVGpVmQX+hu5Lddah2utw7TWYX5+fpanLECfzVlMVlYWrw3obnQUIYSDCPT3Y+DgV3ANbkL81bRC264lhR4DVMzxOAD4x/8jlFK1gZnA41pr2xh6LA/HomP5dfMWWrZoTnBgeaPjCCEcSI9HGuHh5sbPBwtv7lFLCn0HEKyUqqyUcgW6AStyrqCUqgT8APTWWh+1fsyCMWnWQpycnBjat6vRUYQQDqaYmzPVPdP4csI4tu87UijbzLPQtdaZwCvAz0AksERrfVApNUgpNejGau8APsB0pdQepVREgSW2knNXrnEh05WO7dsT6G8bl3+EEI6lQ8MakHWdqd8UzoTSFs1YpLVeDay+7bkZOb5+HnjeutEK1vrIeGo1fpiR7WoYHUUI4aBKFffi6c6P8t2i79m4fS8PNyrYj+gUyU+K7jhwjE2/R9C0qg/F3GQWPiFEwXmxW2e8ixdn+txFmM0FOwlGkSz0yTO/I3LLCh6o5G10FCGEg/PydKfrk505HhXF+u17C3RbRe70dNP2vUQeOkTPbs/I1HJCiEIxoEsHjl3RHL3mxSO64CaULlJn6GazmWlzF+Hl7c2gbp2NjiOEKCLc3Vzp/Whzzkbs50Ddh8BkgqAgmG/dQbyKVKGv3RrB8agouj71OF6e8qlQIUThqbd1DVcXfc3oKylkag3R0TBwoFVLvcgUutaarcfO418piOe7dDA6jhCiiDG9NYaHTh8gJiuTVX5B2U+mpsKYMVbbRpG5hh4ZlwSlg3j/raYy8bMQovCdPk1XrYmpG0q5C6dved5aisQZutls5uulP1PKw0T9SqWMjiOEKIoqVcIE/GfPGhrFHLzleWspEoW+cPWv/LpyCZX0BZn4WQhhjHHjwNPz1tEOPT2zn7cShy/0zMws5ixcSjn/8jzbvpnRcYQQRVXPnhAeDoGBoFT27+Hh2c9bicNfQ5+34hcuJCTw1n9elYmfhRDG6tnTqgV+O4duuLT0DL5b8iMVAwN5sk1To+MIIUSBcuhC33wgmiwnNwb17ipn50IIh+ewl1yuZ5nZFZ9J9xdepf1DVYyOI4QQBc5hT1uXbtzBhUtXaVurnJydCyGKBIdsuqspqUz7cgbnd/5MFT8vo+MIIUShcMhCn77gJ1KSkxncq4vRUYQQotA4XKEnXk7ip5VrCK0dStN6IUbHEUKIQuNwhT7tux9Ju3aNIX27Gx1FCCEKlUMVekp6JjuPnqZe/Xo8cH+w0XGEEKJQOdRti1uPJVCr1eMMaSm3KQohih6HOUOPTbjIuogj1A0oiX8pmVpOCFH0OEyhT5m9mM2LZ9CwoqfRUYQQwhAOUejRZ+PZuGkzDzZuTFA5H6PjCCGEIRyi0KfMWYxGM6xfV6OjCCGEYey+0I9Fx7J561ZaNm9O1Ur+RscRQgjD2H2hL9u8C1dXN17t+6zRUYQQwlB2XegJSelc9Q5k1DsfEOjvZ3QcIYQwlF0X+tKt+3AxKVrfX8HoKEIIYTi7LfS9R07w9Wfj8bp0DG93F6PjCCGE4ey20D/7ZiEurq70ffxho6MIIYRNsKjQlVLtlVJHlFJRSqlRuSxXSqmpN5bvU0rVt35UYP58CApih7cPu79fSqeA8viVKlEgmxJCCHuTZ6ErpZyAaUAHIATorpS6fVzaDkDwjV8DgS+tnDO7zAcOhOhoppb0xz0zk1e++yr7eSGEEBadoTcEorTWJ7TWGcAi4PHb1nkcmKez/QmUVEpZ96bwMWMgNZXDxf046uzC49fTKJV8Nft5IYQQFhV6BeBMjscxN57L7zoopQYqpSKUUhEJCQn5S3r6NACeysSAUuV5Of7kLc8LIURRZ0mhq1ye03exDlrrcK11mNY6zM8vn/eNV6qU/duV8wzcvRpvc9YtzwshRFFnSaHHABVzPA4AYu9inXszbhx43jaSoqdn9vNCCCEsKvQdQLBSqrJSyhXoBqy4bZ0VwHM37nZpDFzRWsdZNWnPnhAeDoGBoFT27+Hh2c8LIYTIe8YirXWmUuoV4GfACfhGa31QKTXoxvIZwGqgIxAFpAL9CiRtz55S4EII8S8smoJOa72a7NLO+dyMHF9r4GXrRhNCCJEfdvtJUSGEELeSQhdCCAchhS6EEA5CCl0IIRyEyn4/04ANK5UARN/lt/sCF6wYx0iyL7bJUfbFUfYDZF/+Fqi1zvWTmYYV+r1QSkVorcOMzmENsi+2yVH2xVH2A2RfLCGXXIQQwkFIoQshhIOw10IPNzqAFcm+2CZH2RdH2Q+QfcmTXV5DF0II8U/2eoYuhBDiNlLoQgjhIGy60G1mcmorsGBfWiqlriil9tz49Y4ROfOilPpGKRWvlDrwL8vt6ZjktS/2ckwqKqU2KaUilVIHlVJDc1nHLo6LhftiL8fFXSn1l1Jq7419eS+Xdax7XLTWNvmL7KF6jwNVAFdgLxBy2zodgTVkz5jUGNhudO572JeWwEqjs1qwL82B+sCBf1luF8fEwn2xl2PiD9S/8bU3cNSO/65Ysi/2clwU4HXjaxdgO9C4II+LLZ+h28bk1NZhyb7YBa31FuDiHVaxl2Niyb7YBa11nNZ6142vk4BI/jmnr10cFwv3xS7c+Fkn33jocuPX7XehWPW42HKhW21yahtgac4mN/57tkYpVatwolmdvRwTS9nVMVFKBQH1yD4bzMnujssd9gXs5LgopZyUUnuAeOAXrXWBHheLJrgwiNUmp7YBluTcRfYYDclKqY7AciC4wJNZn70cE0vY1TFRSnkBy4BhWuurty/O5Vts9rjksS92c1y01llAXaVUSeBHpdT9Wuuc79lY9bjY8hm6bUxObR155tRaX/37v2c6e4YoF6WUb+FFtBp7OSZ5sqdjopRyIbsA52utf8hlFbs5Lnntiz0dl79prS8DvwLtb1tk1eNiy4VuG5NTW0ee+6KUKqeUUje+bkj2sUks9KT3zl6OSZ7s5ZjcyDgLiNRaT/qX1eziuFiyL3Z0XPxunJmjlPIA2gCHb1vNqsfFZi+5aFuanPoeWbgvXYDBSqlM4BrQTd94G9yWKKUWkn2Xga9SKgYYS/abPXZ1TMCifbGLYwI0BXoD+29crwUYDVQCuzsuluyLvRwXf2CuUsqJ7H90lmitVxZkh8lH/4UQwkHY8iUXIYQQ+SCFLoQQDkIKXQghHIQUuhBCOAgpdCGEcBBS6EII4SCk0IUQwkH8H2ruSPTSdHORAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Get values and derivatives from function f\n",
    "node_x = np.array([0,1.0,2.0,3.0])\n",
    "node_y = f(node_x)\n",
    "node_y_prime = f_prime(node_x)\n",
    "\n",
    "# Set the x values\n",
    "x = np.linspace(0,3,100)\n",
    "sf_id = 2\n",
    "\n",
    "plt.plot(node_x,f(node_x),'ro')\n",
    "plt.plot(x,hermite_div_diff(x,node_x,node_y,node_y_prime),label='Hermite',alpha=0.6)\n",
    "plt.plot(x,f(x),'k--',label='Sin function',alpha=0.6)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
  },
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "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.8.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}