hermite_shape_fun.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hermite interpolating polynomials\n",
    "We first construct the Hermite basis functions, expressed as\n",
    "$$\n",
    "H_{n,j}(x) = \\left[1 - 2\\,(x-x_{j})\\,L'_{n,j}(x_{j})\\right]\\,L^{2}_{n,j}(x),\n",
    "$$\n",
    "and \n",
    "$$\n",
    "\\widehat{H}_{n,j}(x) = (x-x_{j})\\,L^{2}_{n,j}(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we re-use the function to generate the Lagrange basis functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "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": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lagrange_sf(x,k,node_x):\n",
    "  '''Evaluates the generic Lagrange shape function\n",
    "\n",
    "  :param x: the locations where to evaluate the shape function\n",
    "  :type x: double or np.array\n",
    "  :param k: the function to be evaluated from 0 to len(node_x)-1. This is 1.0 at x_k and zero at the other points.\n",
    "  :type k: int\n",
    "  :param node_x: the x locations where the values of f are known. \n",
    "  :type node_x: list or np.array containing double x locations\n",
    "  :raises Exception: Generic exception when k has an invalid range\n",
    "  :return: the Lagrange shape functions evaluated at x.\n",
    "  :rtype: double, list or np.array, depending on x and node_x.\n",
    "  '''\n",
    "  # Polynomial order\n",
    "  n = len(node_x)-1\n",
    "  if(k < 0)or(k > n):\n",
    "    raise Exception(\"ERROR: invalid k\")\n",
    "  else:\n",
    "    # Construct Lagrange shape functions    \n",
    "    res = np.ones(len(x))\n",
    "    for j in range(n+1):\n",
    "      if(j != k):\n",
    "        res *= (x-node_x[j])/(node_x[k]-node_x[j])\n",
    "  return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Additionally we need a way to generate the derivative $L'_{n,j}(x_{j})$ as discussed in class. \n",
    "Remember that \n",
    "$$\n",
    "L'_{n,j}(x_{j}) = \\sum_{k=0,k\\ne j}^{n}\\,\\dfrac{1}{x_{j} - x_{k}}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lagrange_prime(j,node_x):\n",
    "  '''Eval the derivative of the j-th Lagrange polynomial at xj\n",
    "\n",
    "  :param j: Basis id\n",
    "  :type j: int\n",
    "  :param node_x: x coordinate at the nodes\n",
    "  :type node_x: list or np.array\n",
    "  :raises Exception: exception for invalid basis id j\n",
    "  :return: the value of the derivative\n",
    "  :rtype: double\n",
    "  '''\n",
    "  n = len(node_x)-1\n",
    "  if(j < 0)or(j > n):\n",
    "      raise Exception(\"ERROR: Invalid shape function number j\")\n",
    "  res = 0.0\n",
    "  for k in range(n+1):\n",
    "    if(k != j):\n",
    "      res += 1.0/(node_x[j]-node_x[k])\n",
    "  # return \n",
    "  return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next function computed the shape function \n",
    "$$\n",
    "H_{n,j}(x) = \\left[1 - 2\\,(x-x_{j})\\,L'_{n,j}(x_{j})\\right]\\,L^{2}_{n,j}(x)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hermite_sf1(x,j,node_x):\n",
    "  '''The Hermite shape functions of first type.\n",
    "  The j-th function should be zero at all nodes except xj.\n",
    "  In addition all first derivatives should be zero at the nodes.\n",
    "\n",
    "  :param x: x coordinates where to evaluate the basis function\n",
    "  :type x: list or np.array\n",
    "  :param j: basis function id\n",
    "  :type j: int\n",
    "  :param node_x: x coordinates at the nodes\n",
    "  :type node_x: list or np.array\n",
    "  :return: the shape function at the locations x\n",
    "  :rtype: list or np.array\n",
    "  '''\n",
    "  l_prime = lagrange_prime(j,node_x)\n",
    "  l_sf = lagrange_sf(x,j,node_x)\n",
    "  return (1.0 - 2.0*(x-node_x[j])*l_prime)*(l_sf**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the shape function $H_{n,j}(x)$. \n",
    "This function should have all derivatives equal to zeros at the nodes and should be equal to zero at all the nodes except node $x_{j}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZwcdZ3/8ddnunvuI9fkzmRCiEgCIYYxgrAaRCAG/UV2WQVd1mV1Izxk1b0UREVXWdhl12O9IGpWViL8UAFZTUBwRfCAMIkJBEJCCLlIyEkymZnM0TOf/aMr0Ex6Mt2Z7unuqffz8ahHV32rqvtT+Wb63XV0l7k7IiISXiX5LkBERPJLQSAiEnIKAhGRkFMQiIiEnIJARCTkovku4ESMGTPGGxsb812GiEhRWbVq1T53r+/bXpRB0NjYSHNzc77LEBEpKma2NVW7Dg2JiIScgkBEJOSyEgRmttTM9pjZun7mzzezQ2a2Jhg+nzRvgZltMLNNZnZtNuoREZH0ZWuP4AfAggGWeczd5wTDPwOYWQT4FvAuYCZwuZnNzFJNIiKShqwEgbs/Chw4gVXnAZvcfbO7dwF3AYuyUZOIiKRnKM8RnG1ma81shZnNCtomAduTltkRtImIyBAZqstHVwNT3b3VzBYC9wEzAEuxbMqfQzWzxcBigIaGhlzVKSISOkOyR+DuLe7eGowvB2JmNobEHsCUpEUnAzv7eY4l7t7k7k319cd8H0JERE7QkASBmY03MwvG5wWvux94EphhZtPMrBS4DLh/KGoSEZGErBwaMrM7gfnAGDPbAdwAxADc/VbgUuBqM4sDR4DLPHFHnLiZXQM8CESApe7+TDZqEhGR9Fgx3qGsqanJ9RMTIiKZMbNV7t7Ut13fLBYRCTkFgYhIyCkIRERCTkEgIhJyCgIRkZBTEIiIhJyCQEQk5BQEIiIhpyAQEQk5BYGISMgpCEREQk5BICIScgoCEZGQUxCIiIScgkBEJOQUBCIiIacgEBEJuawEgZktNbM9Zraun/kfNLOnguH3ZnZG0rwtZva0ma0xM912TERkiGVrj+AHwILjzH8ReLu7zwa+BCzpM/88d5+T6hZqIiKSW1m5eb27P2pmjceZ//ukyceBydl4XRERGbx8nCP4MLAiadqBX5rZKjNb3N9KZrbYzJrNrHnv3r05L1JEJCyyskeQLjM7j0QQnJvUfI677zSzscBDZvacuz/ad113X0JwSKmpqcmHpGARkRAYsj0CM5sNfA9Y5O77j7a7+87gcQ9wLzBvqGoSEZEhCgIzawDuAa5w941J7VVmVnN0HLgQSHnlkYiI5EZWDg2Z2Z3AfGCMme0AbgBiAO5+K/B5YDTwbTMDiAdXCI0D7g3aosCP3P2BbNQkIiLpydZVQ5cPMP8jwEdStG8Gzjh2DRERGSr6ZrGISMgpCEREQk5BICIScgoCEZGQUxCIiIScgkBEJOQUBCIiIacgEBEJOQWBiEjIKQhEREJOQSAiEnIKAhGRkFMQiIiEnIJARCTkFAQiIiGnIBARCbmsBIGZLTWzPWaW8jaTlvCfZrbJzJ4ys7lJ8xaY2YZg3rXZqCeVZcugsRFKShKPy5bl6pUkE+qXwqM+KUw57Rd3H/QAvA2YC6zrZ/5CYAVgwFnAE0F7BHgBOAkoBdYCMwd6vTPPPNMzcccd7pWV7vDaUFmZaJf8Ub8UHvVJYcpWvwDNnuI91RLzBs/MGoGfu/tpKebdBjzi7ncG0xtI3OO4EfiCu18UtF8XhNNNx3utpqYmb25uTru2xkbYuhXKJh8gOqKdeEs5XTtH0jApwpYtaT+NZNnRfulr6lTUL3miPilM2eoXM1vlifvFv05W7lmchknA9qTpHUFbqva3pHoCM1sMLAZoaGjI6MW3bUs8Vs3aQc2cxMt5vITW9RPZvHc6J9VXZ/R8kh1H+yXddsk99UlhynW/DNXJYkvR5sdpP7bRfYm7N7l7U319fUYvfjQ3XvnVLF66bT67f/xmWp+eTPWpu7joa4/y1Yc20tObnT0jSV9/eZ5hzksWqU8KU677ZaiCYAcwJWl6MrDzOO1ZdeONUFkJHo8QP1hFx+axdPz2dK477TwuPn0CX//V83xo6Ur2t3Zm+6XlOI72S7LKykS75If6pDDlvF9SnTg4kYHE8f7+ThZfzOtPFq8M2qPAZmAar50snjXQa2V6stg9cVJl6lR3s8Rj8kmWu1Zu9RnXL/fzbvm17zp4JOPnlhN3vH6R/FCfFKZs9Au5PFlsZneSOPk7BtgN3ADEgqC51cwM+CawAGgHrnT35mDdhcDXSFxBtNTdB8y4TE8Wp+PJLQe48r+eZFRVKXcuPotJIyqy+vwiIvnW38nirF01NJRyEQQAa7Yf5IrvP8G42nJ+etVbqauMZf01RETypb8g0DeLk8yZMoIlVzSxbX87f/PDZjrjPfkuSUQk5xQEfZw9fTS3/PlsVr54gC/9/Nl8lyMiknMKghQWzZnER99+Enc8vo37/vhSvssREckpBUE//unCU5g3bRTX3fM0G3cfznc5IiI5oyDoRzRSwjcvfxNVZVGu+dFqOrp1vkBEhicFwXGMrS3nlj+fzcbdrXz1oY35LkdEJCcUBAM475SxXD5vCkse20zzlgP5LkdEJOsUBGm4/uKZTB5ZwT/8eC3tXfF8lyMiklUKgjRUl0W55dIz2Lq/na8//Hy+yxERySoFQZrOOmk072+awvd++yLrd7XkuxwRkaxREGTguoVvZERFjOvueZpe/Wy1iAwTCoIMjKgs5bPvPpU12w+ybKXu1CEiw4OCIEPvnTOJc04ezb898Bx7DnfkuxwRkUFTEGTIzPjSotPo6O7hlgc25LscEZFBUxCcgJPqq/nrc6bx41U7WLv9YL7LEREZFAXBCbrmHSczprqML/zPMxTjPR1ERI5SEJygmvIYn1pwCn/cdpD71ugXSkWkeGUlCMxsgZltMLNNZnZtivn/ZGZrgmGdmfWY2ahg3hYzezqYl/3bjuXQpXMnM3tyHTeveI62Tn3jWESK06CDwMwiwLeAdwEzgcvNbGbyMu5+i7vPcfc5wHXAb9w9+Yd7zgvmH3MLtUJWUmLc8J5Z7G7p5NuPbMp3OSIiJyQbewTzgE3uvtndu4C7gEXHWf5y4M4svG5BOHPqSC550yS++9iLvHTwSL7LERHJWDaCYBKwPWl6R9B2DDOrBBYAP01qduCXZrbKzBb39yJmttjMms2see/evVkoO3v+8aJTMODfH9TlpCJSfLIRBJairb/LaN4D/K7PYaFz3H0uiUNLHzOzt6Va0d2XuHuTuzfV19cPruIsmzSigr8+dxr3/vEl1r10KN/liIhkJBtBsAOYkjQ9GdjZz7KX0eewkLvvDB73APeSONRUdK6eP51RVaXc+Iv1upxURIpKNoLgSWCGmU0zs1ISb/b3913IzOqAtwM/S2qrMrOao+PAhcC6LNQ05GrLY3zi/Bn8YfN+fr1hT77LERFJ26CDwN3jwDXAg8B64G53f8bMrjKzq5IWvQT4pbu3JbWNA35rZmuBlcAv3P2BwdaULx94SwPTxlRx0/LniPf05rscEZG0WDEexmhqavLm5sL8ysED617mqjtW8S+XnM4H3tKQ73JERF5lZqtSXaavbxZn2UWzxtE0dSRfeWijvmQmIkVBQZBlZsZnLj6Vfa2dLHl0c77LEREZkIIgB+Y2jOTi2RNY8uhmdrfongUiUtgUBDny6YveSLy3l68+tDHfpYiIHJeCIEcaRlfyl2c3cnfzdja8fDjf5YiI9EtBkEN/+46TqS6LctOK9fkuRUSkXwqCHBpRWcrHz5/BIxv28tjzhfX7SCIiRykIcuyKs6cyZVQFN/5iPT29xfedDREZ/hQEOVYWjfCpi97Icy8f5p7VO/JdjojIMRQEQ+DdsycwZ8oI/v2XGzjS1ZPvckREXkdBMATMjOsvPpXdLZ187zF9yUxECouCYIi8uXEUC2aN5zu/eYE9h/UlMxEpHAqCIfTpd72RrngvX3v4+XyXIiLyKgXBEJo2poq/OGsqd63cxvO79SUzESkMCoIh9vHzZ1BVGuWmFc/luxQREUBBMORGVZXysXeczP8+t4ffbdqX73JERLITBGa2wMw2mNkmM7s2xfz5ZnbIzNYEw+fTXXc4+qu3NjJpRAX/snw9vfqSmYjk2aCDwMwiwLeAdwEzgcvNbGaKRR9z9znB8M8ZrjuslMcifGrBKTyzs4X71ryU73JEJOSysUcwD9jk7pvdvQu4C1g0BOsWtffMnsjsyXX82wMbaO/SncxEJH+yEQSTgO1J0zuCtr7ONrO1ZrbCzGZluC5mttjMms2see/e4v8Bt5IS43PvnsnLLR1855EX8l2OiIRYNoLAUrT1PfC9Gpjq7mcA3wDuy2DdRKP7Endvcvem+vr6Ey62kLy5cRSL5kzktkc3s21/e77LEZGQykYQ7ACmJE1PBnYmL+DuLe7eGowvB2JmNiaddYe76951KtES48blz+a7FBEJqWwEwZPADDObZmalwGXA/ckLmNl4M7NgfF7wuvvTWXe4G19XzsfOO5kHn9mtexaISF4MOgjcPQ5cAzwIrAfudvdnzOwqM7sqWOxSYJ2ZrQX+E7jME1KuO9iais2Hz53G1NGVfPF/nqW7pzff5YhIyJh78V3H3tTU5M3NzfkuI6sefnY3H/nvZj737pl8+Nxp+S5HRIYhM1vl7k192/XN4gJx/qljedsb6vnawxvZe7gz3+WISIgoCAqEmXHDe2bS2d3Ll3+hE8ciMnQUBAVken01V8+fzs/W7OTRjTpxLCJDQ0FQYK6eP52TxlTx2fvW0dGt21qKSO4pCApMeSzCly85jW0H2vnG/+oGNiKSewqCAvTW6WP4s7mTue03m9moG9iISI4pCArU9RefSk15lM/c87R+qlpEckpBUKBGVZXymYWn0rz1FX60clu+yxGRYUxBUMAuPXMy5548hpuWr2f7Af0onYjkhoKggJkZ/3rpbMyMT/3kKR0iEpGcUBAUuEkjKvjcu0/lD5v388PHt+a7HBEZhhQEReB9TVOYf0o9N694ji372vJdjogMMwqCImBm3Pyns4lGjH/88Vp6dIhIRLJIQVAkxteV88X/N4vmra/w7V9vync5IjKMKAiKyCVvmsR750zkqw9vZOWLB/JdjogMEwqCImJmfPmS02kYVckn7vojB9u78l2SiAwDWQkCM1tgZhvMbJOZXZti/gfN7Klg+L2ZnZE0b4uZPW1ma8xseN1tJgeqy6J84/K57Gvt5J9+8hTFeGMhESksgw4CM4sA3wLeBcwELjezmX0WexF4u7vPBr4ELOkz/zx3n5PqzjlyrNMn1/HpBW/koWd3c/vvt+S7HBEpctnYI5gHbHL3ze7eBdwFLEpewN1/7+6vBJOPA5Oz8Lqh9uFzp3H+G8dy4/L1NG/R+QIROXHZCIJJwPak6R1BW38+DKxImnbgl2a2yswW97eSmS02s2Yza967VzdtMTO+8v45TBpRwVV3rOblQx35LklEilQ2gsBStKU8cG1m55EIgk8nNZ/j7nNJHFr6mJm9LdW67r7E3Zvcvam+vn6wNQ8LdRUxlvxlE0e64nz0jlW6kY2InJBsBMEOYErS9GRgZ9+FzGw28D1gkbvvP9ru7juDxz3AvSQONUma3jCuhv943xms3X6Qz963TiePRSRj2QiCJ4EZZjbNzEqBy4D7kxcwswbgHuAKd9+Y1F5lZjVHx4ELgXVZqClUFpw2gU+cP4OfrNrBN/5XXzYTkcxEB/sE7h43s2uAB4EIsNTdnzGzq4L5twKfB0YD3zYzgHhwhdA44N6gLQr8yN0fGGxNYfTJd85g+4F2vvLQRiaOqODSM3U+XkTSY8V4KKGpqcmbm/WVg7664r1c+YOVPLH5AD+4ch7nzhiT75JEpICY2apUl+nrm8XDSGm0hO/8xZmcPLaaj/6wmdXbXhl4JREJPQXBMFNbHuP2v55HfU0ZH1q6kqd3HMp3SSJS4BQEw9C42nJ+9DdnUVcR44qlT7B+V0u+SxKRAqYgGKYmjqjgzr85i4pYhA9893Ge2nEw3yWJSIFSEAxjU0ZVctfis6gqi3L5ksf5/aZ9+S5JRAqQgmCYmzq6ip9e/VYmjazgr/7rSR5YtyvfJYlIgVEQhMC42nLu/ujZzJpUy9XLVnPrb17QN5BF5FUKgpAYUVnKso+8hYWnTeDmFc/xd/9/jX6bSEQABUGoVJZG+eYH3sQ/XvgG7luzk/fd9ge27m/Ld1kikmcKgpAxM655xwyWXHEmW/a1sfDrj3HP6h35LktE8khBEFIXzhrPik++jVkT6/j7u9fy8Tv/yIE23QNZJIwUBCE2aUQFdy4+i7+/4A0sf3oX5//HI9zdvF0nkkVCRkEQcpES4+Pnz+AXH/8TptdX86mfPMX7b3uctdv1BTSRsFAQCACnjK/h7o+ezb/+2em8sLeVRd/6HVf9cBWb9hzOd2kikmODvh+BDB8lJcb739zAxbMn8v3HXuS7j23mwWdf5oJTx/GRPzmJNzeOJLh3hEhRcXdaO+McOtLN4Y44nfFeunt66Y730hk89jrEIkY0UkK0xIiWGLFoCTVlUWorYtSWxyiPlQzLvwHdj0D6daCti+//djPLntjGwfZuTp9Ux/uaJvPu2RMZWVWa7/Ik5Dq6e9h7uJM9hzvY09LJ7pYO9hzuZHdLou1gezeHjnTT0tFNy5FuerPwVheLGHUVMepryhlfW8b4unLG1ZYzvracySMraRxTycS6CkpKCjMs+rsfgYJABnSkq4efrt7BD/+wlQ27DxOLGPNPGcuFM8fx9lPqGVtTnu8ST4i7c6S7h9aOOK2dcdo6ezjc2U1bZw9tnXHauuJ0dvcS7+2lu8fp7ukl3uN09yYee92JlhglJUbEjEhJMJgRiRgVsQiVpREqSqNJ44nHqtIoteUxasqjBfumkS/uzuHOOC8f6mDnwSO8fKiDXYc62HXoCLsOdbC7pYPdLZ0cOtJ9zLrREqO+poyxNWWMqiqltiJGXfBpvq4iRm1FlJrgk30sUkJppIRYNPFoBvEeJx70b7zX6erppbUjHoRJ4vFgexe7Wzp5Oahlf5+r7cqiJUwdXUnj6CqmjUkM08dWM72+mlF5/gCV0yAwswXA10ncqvJ77n5zn/kWzF8ItAN/5e6r01k3lRMKgmXL4PrrYds2aGiAG2+ED34ws+cIOXfn2V0t3Lv6Jf7nqZ3sbukE4LRJtcxrHM2ZU0cyd+oIxteWp7/7nMV+6el1Xmnv4pW2Lva3pXhs7+JAWxf7WxPj+9u66Ir3ZvQaZhCLlBArMUrM6PHEG0Zvr9PjTqZ/TmZQXRalLnjD6jvUBkPf9hFBeyQXIZLDv5WueC/72zrZd7iLfa2dvNwSvMkfPPK68bau13/r3Qzqq8uYEHwCH1dbztiaMsbVllNfW8a4mnLG1pYxqrJ0yIO1K97L7pYOtr/SzpZ97by4r5UX97WzZX8b2/a309Xz2v+xkZUxptdXc1J9FdPrE+EwfWw1U0ZWEI0McMo2C/2SsyAwswiwEbgA2EHiZvaXu/uzScssBP6WRBC8Bfi6u78lnXVTyTgIli2DxYuhvf21tspKWLJEYXCC3J31uw7z6w17+M3GvazdfpDO4E21tjzK9LHVnFxfzbT6KsbWJP5o62vKGFEZoyKW+GRcetdd2EeP7ZfOW5fQ9qd/TltnnMMdiU/mrZ1xWo50s7818WZ+oL2LA8H4/rZODrR1cfBId79vxDVlUUZVlzKyspTRVaWMCoaRVaVUl0WpLotSVRalqizy6nh1WZSyaMmrx4xjkZIB33iPBkJPr3Okq4f27h6OdMVp7+qhvauHI109ib2QYHtajnTT0pE4dp08tASPnQME1dHj18eESOVrQZL4RJzYKyl/dSihIhahLBhPfCK24/6t9Fz+AbrivXTFe+ns6aEr3ktHdy+HOxLH3VuOPh55bXp/Wxf7Dneyr7WTfa1dKT/Fm8HYmjLG11Uwsa6c8XXlTKgrZ0JdBROC6XG15cQGeqMsQD29zs6DR9i0t5UX9rTywt42Xtjbyua9bexr7Xx1uVjEmDq6iulBQJxUnwiHo/9HR973EyIp/lYyfQ/LZRCcDXzB3S8Kpq8DcPebkpa5DXjE3e8MpjcA84HGgdZNJeMgaGyErVuPbZ86FbZsSf95pF9d8V7W72phzfaDPL/nMJuC//R7D3f2u06kt4dYTze9VoKb0WMl9JZEBnytEoORla+9mY/u+wZfXcboqqAtmFcaLb43EUgcB2/pExKphlTLdHRnvrdT0tODuWM45mDei+F0R2L0pNE3R5UY1JTHGF1VypjqMsbUBI/VZYyufm18fF3iQ0IxvskP1qH2bl7YlwiIzfvagqBoZev+duIpTmiUxrso7Ynznfv+hT/ZsibRmOF7WH9BkI2rhiYB25Omd5D41D/QMpPSXBcAM1sMLAZoaGjIrMJt2zJrl4yVRks4Y8oIzpgy4nXt7V3x4IReJ3taOjnc0Z34ZNzdw5Eb/pmuSBRzJ+K9lARDaU8P1bfc9Oqn8qpgqKuIMqqqjLpcHRIpQEc/wY+tzfw8TGe8Jykk4nR09wRDb+IxnthD6Yz30tndgwP+5RvpNcMx3AgejVhvD6Wf+yyl0RLKoiWUBkNZtITa8teOvdeUJ859VJZGhuXVNdlUVxljbsNI5jaMfF17d08v2w60s+tgBweCQ50HPn09nZEY3ZEYE1qS7iuSpfewbARBqt7uG2f9LZPOuolG9yXAEkjsEWRSIA0NqfcIMg0UyVhlaZSpo6NMHV117MwrH+1/T+2cabkvbpgri0YYWxPJ7GT+4uP0yTtvz15x0q9YpOTV8wev2v67nL6HZWN/bAcwJWl6MrAzzWXSWXfwbrwxcTwtWWVlol3yR/1SeNQnhSnX/eLugxpI7FVsBqYBpcBaYFafZS4GVpDYAzgLWJnuuqmGM8880zN2xx3uU6e6myUe77gj8+eQ7FO/FB71SWHKQr8AzZ7iPTVbl48uBL5G4hLQpe5+o5ldFQTNrcHlo98EFpC4fPRKd2/ub92BXk/fIxARyZy+UCYiEnL9BUH4rtkSEZHXURCIiIScgkBEJOQUBCIiIacgEBEJOQWBiEjIKQhEREJOQSAiEnIKAhGRkFMQiIiEnIJARCTkFAQiIiGnIBARCTkFgYhIyCkIRERCTkEgIhJygwoCMxtlZg+Z2fPB48gUy0wxs1+b2Xoze8bMPpE07wtm9pKZrQmGhYOpR0REMjfYPYJrgV+5+wzgV8F0X3HgH9z9VBL3K/6Ymc1Mmv9Vd58TDMsHWY+IiGRosEGwCLg9GL8deG/fBdx9l7uvDsYPA+uBSYN8XRERyZLBBsE4d98FiTd8YOzxFjazRuBNwBNJzdeY2VNmtjTVoaWkdRebWbOZNe/du3eQZYuIyFEDBoGZPWxm61IMizJ5ITOrBn4KfNLdW4Lm7wDTgTnALuA/+lvf3Ze4e5O7N9XX12fy0iIichzRgRZw93f2N8/MdpvZBHffZWYTgD39LBcjEQLL3P2epOfenbTMd4GfZ1K8iIgM3mAPDd0PfCgY/xDws74LmJkB3wfWu/tX+sybkDR5CbBukPWIiEiGBhsENwMXmNnzwAXBNGY20cyOXgF0DnAF8I4Ul4n+m5k9bWZPAecBfzfIekREJEMDHho6HnffD5yfon0nsDAY/y1g/ax/xWBeX0REBk/fLBYRCTkFgYhIyCkIRERCTkEgIhJyCgIRkZBTEIiIhJyCQEQk5BQEIiIhpyAQEQk5BYGISMgpCEREQk5BICIScgoCEZGQUxCIiIScgkBEJOQUBCIiITeoIDCzUWb2kJk9HzyO7Ge5LcGdyNaYWXOm64uISO4Mdo/gWuBX7j4D+FUw3Z/z3H2Ouzed4PoiIpIDgw2CRcDtwfjtwHuHeH0RERmkwQbBOHffBRA8ju1nOQd+aWarzGzxCawvIiI5MuDN683sYWB8ilnXZ/A657j7TjMbCzxkZs+5+6MZrE8QIIsBGhoaMllVRESOY8AgcPd39jfPzHab2QR332VmE4A9/TzHzuBxj5ndC8wDHgXSWj9YdwmwBKCpqckHqltERNIz2END9wMfCsY/BPys7wJmVmVmNUfHgQuBdemuLyIiuTXYILgZuMDMngcuCKYxs4lmtjxYZhzwWzNbC6wEfuHuDxxvfRERGToDHho6HnffD5yfon0nsDAY3wyckcn6IiIydPTNYhGRkFMQiIiEnLkX3wU4ZrYX2HqCq48B9mWxnHzSthSe4bIdoG0pVIPZlqnuXt+3sSiDYDDMrLnPz1wULW1L4Rku2wHalkKVi23RoSERkZBTEIiIhFwYg2BJvgvIIm1L4Rku2wHalkKV9W0J3TkCERF5vTDuEYiISBIFgYhIyA3bIDCzBWa2wcw2mdkxdz6zhP8M5j9lZnPzUWc60tiW+WZ2KLgV6Boz+3w+6hyImS01sz1mtq6f+UXRJ2lsR1H0B4CZTTGzX5vZejN7xsw+kWKZYumXdLal4PvGzMrNbKWZrQ2244splslun7j7sBuACPACcBJQCqwFZvZZZiGwAjDgLOCJfNc9iG2ZD/w837WmsS1vA+YC6/qZXyx9MtB2FEV/BLVOAOYG4zXAxiL+W0lnWwq+b4J/5+pgPAY8AZyVyz4ZrnsE84BN7r7Z3buAu0jcFjPZIuC/PeFxYERwT4RCk862FAVP3IzowHEWKYo+SWM7ioa773L31cH4YWA9MKnPYsXSL+lsS8EL/p1bg8lYMPS9qierfTJcg2ASsD1pegfH/odIZ5lCkG6dZwe7kivMbNbQlJZ1xdIn6Si6/jCzRuBNJD6BJiu6fjnOtkAR9I2ZRcxsDYmbdT3k7jntk0H9DHUBsxRtfRM1nWUKQTp1ribxGyKtZrYQuA+YkfPKsq9Y+mQgRdcfZlYN/BT4pLu39J2dYpWC7ZcBtqUo+sbde4A5ZjYCuNfMTnP35HNSWe2T4bpHsAOYkjQ9Gdh5AssUggHrdPeWo7uS7r4ciJnZmKErMWuKpU+Oq9j6w8xiJN44l7n7PSkWKZp+GWhbitIOruUAAAEYSURBVK1v3P0g8AiwoM+srPbJcA2CJ4EZZjbNzEqBy0jcFjPZ/cBfBmffzwIOufuuoS40DQNui5mNNzMLxueR6Nf9Q17p4BVLnxxXMfVHUOf3gfXu/pV+FiuKfklnW4qhb8ysPtgTwMwqgHcCz/VZLKt9MiwPDbl73MyuAR4kcdXNUnd/xsyuCubfCiwnceZ9E9AOXJmveo8nzW25FLjazOLAEeAyDy4tKCRmdieJqzbGmNkO4AYSJ8KKqk/S2I6i6I/AOcAVwNPBMWmAzwANUFz9QnrbUgx9MwG43cwiJILqbnf/eS7fv/QTEyIiITdcDw2JiEiaFAQiIiGnIBARCTkFgYhIyCkIRERCTkEgIhJyCgIRkZD7PwrNvBWrjvptAAAAAElFTkSuQmCC",
      "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,2,3])\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,500)\n",
    "sf_id = 0\n",
    "hermite_y = hermite_sf1(x,sf_id,node_x)\n",
    "\n",
    "plt.axis('equal')\n",
    "plt.plot(node_x,np.zeros(len(node_x)),'ro')\n",
    "plt.plot(node_x,np.ones(len(node_x)),'bo')\n",
    "plt.plot(x,hermite_y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also compute the second type of shape function\n",
    "$$\n",
    "\\widehat{H}_{n,j}(x) = (x-x_{j})\\,L^{2}_{n,j}(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hermite_sf2(x,j,node_x):\n",
    "  '''The Hermite shape functions of second type.\n",
    "  The j-th function should be zero at all nodes.\n",
    "  In addition, all first derivatives should be zero at the nodes except 1 at xj.\n",
    "\n",
    "  :param x: x locations for shape function evaluation\n",
    "  :type x: list or np.array\n",
    "  :param j: shape function id\n",
    "  :type j: int\n",
    "  :param node_x: x coordinates at the nodes\n",
    "  :type node_x: list or np.arrray\n",
    "  :return: the shape functions at the locations x\n",
    "  :rtype: list or np.array\n",
    "  '''\n",
    "  l_sf = lagrange_sf(x,j,node_x)\n",
    "  return (x-node_x[j])*(l_sf**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAZiElEQVR4nO3de5Bc5X3m8e/TPTOSRgJd0EjoLtYI1oADgYmMi9QGxciWlWRle+1aWBchVKpUpGBjwHGVdr1lO1VLLRWvc10WotjKKmUCiddgtASMgYrNel22NcICSYCQrCA0GlkzEqALus70b//oM6NWay496p5L630+VV3nnPd9T/d75pXO0+f06dOKCMzMLF25se6AmZmNLQeBmVniHARmZolzEJiZJc5BYGaWuIax7sD5mDlzZixevHisu2FmVlc2bdp0ICJaysvrMggWL15MW1vbWHfDzKyuSNrdX7lPDZmZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniahIEktZJ6pS0dYB6SfpLSTslvSrp+pK6FZK2Z3VratEfMzOrXK2OCP4XsGKQ+k8AS7LHauBhAEl54KGs/irgNklX1ahPZmZWgZoEQUS8BLwzSJNVwN9F0U+AaZLmAEuBnRGxKyJOAY9nbc3MbJSM1mcE84A9JcvtWdlA5eeQtFpSm6S2rq6uEeuomVlqRisI1E9ZDFJ+bmHE2ohojYjWlpZz7plkZmbnabRuOtcOLChZng90AE0DlJuZ2SgZrSOCDcDvZlcP3Qgcioh9wEZgiaTLJDUBt2ZtzcxslNTkiEDSY8DNwExJ7cBXgEaAiHgEeAZYCewEjgF3ZnXdku4BngPywLqI2FaLPpmZWWVqEgQRcdsQ9QHcPUDdMxSDwszMxoC/WWxmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSWuJkEgaYWk7ZJ2SlrTT/0XJW3OHlsl9UiakdW9JWlLVtdWi/6YmVnlqv6FMkl54CFgOcUfqd8oaUNEvNbbJiK+Bnwta/87wH0R8U7J0yyLiAPV9sXMzIavFkcES4GdEbErIk4BjwOrBml/G/BYDV7XzMxqoBZBMA/YU7LcnpWdQ1IzsAL4TklxAN+XtEnS6hr0x8zMhqEWP16vfspigLa/A/y/stNCN0VEh6RZwPOS3oiIl855kWJIrAZYuHBhtX02M7NMLY4I2oEFJcvzgY4B2t5K2WmhiOjIpp3AkxRPNZ0jItZGRGtEtLa0tFTdaTMzK6pFEGwElki6TFITxZ39hvJGkqYCvwE8VVI2WdJFvfPAx4CtNeiTmZlVqOpTQxHRLeke4DkgD6yLiG2S7srqH8mafgr4fkS8X7L6bOBJSb19+fuI+F61fTIzs8opYqDT+eNXa2trtLX5KwdmZsMhaVNEtJaX+5vFZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJq0kQSFohabuknZLW9FN/s6RDkjZnjy9Xuq6ZmY2sqn+zWFIeeAhYDrQDGyVtiIjXypr+34j47fNc18zMRkgtjgiWAjsjYldEnAIeB1aNwrpmZlYDtQiCecCekuX2rKzcRyS9IulZSVcPc10krZbUJqmtq6urBt02MzOoTRCon7IoW34ZWBQR1wJ/BXx3GOsWCyPWRkRrRLS2tLScd2fNzOxstQiCdmBByfJ8oKO0QUQcjoij2fwzQKOkmZWsa2ZmI6sWQbARWCLpMklNwK3AhtIGki6VpGx+afa6BytZ18zMRlbVVw1FRLeke4DngDywLiK2Sborq38E+AzwB5K6gePArRERQL/rVtsnMzOrnIr74/rS2toabW1tY90NM7O6ImlTRLSWl/ubxWZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeJqEgSSVkjaLmmnpDX91H9O0qvZ48eSri2pe0vSFkmbJfnXZszMRlnVP1UpKQ88BCyn+GP0GyVtiIjXSpr9C/AbEfGupE8Aa4EPl9Qvi4gD1fbFzMyGrxZHBEuBnRGxKyJOAY8Dq0obRMSPI+LdbPEnwPwavK6ZmdVALYJgHrCnZLk9KxvI7wPPliwH8H1JmyStHmglSasltUlq6+rqqqrDZmZ2RtWnhgD1Uxb9NpSWUQyCXy8pvikiOiTNAp6X9EZEvHTOE0aspXhKidbW1n6f38zMhq8WRwTtwIKS5flAR3kjSb8CfANYFREHe8sjoiObdgJPUjzVZGZmo6QWQbARWCLpMklNwK3AhtIGkhYCTwC3R8SbJeWTJV3UOw98DNhagz6ZmVmFqj41FBHdku4BngPywLqI2Cbprqz+EeDLwCXA/5QE0B0RrcBs4MmsrAH4+4j4XrV9MjOzyimi/k63t7a2Rlubv3JgZjYckjZlb8LP4m8Wm5klzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklrhZ3H03SsVPdvPP+Kd47dprDx09TKPmCdvOEPBdPbODiiY3MmNxEQ955a2bjl4OgAu++f4of/+IgG996hx2dR9jZeZT9h09WtG5OcOnFE5kzbRJzp01i/vRJLJk1hStmX8Tls6YwsTE/wr03Mxucg2AAB46e5ImX29nwSgfbOg4TAc1NeZbMmsJNl8/kAy1TaJkyganNjVw8sZGGfPFnGSLg/VPdHD5+msMnuuk8fIK97x2n473jvNr+Hs9u2Ud3dvggwcIZzSyZdRFXXjqFD82byjXzpjJv2iSyG/GZmY04B0GZn7/9Ln/9w1288Pp+ugvBdQumcd8tV3DT5TP5lflTaazyNM/pngK7D77Pm/uP8ub+I+zIpj/Y3tkXEDMmN/GheVOLj/nF6ZypEx0OZjYifPfRzM7OI3ztue08t20/05sb+cwN8/n3v7aAy2ddVNPXGciJ0z1s/+URXt17iC3t77Fl72He3H+EniwcZk45Ew7XZA+Hg5kNx0B3H03+iOBkdw///bntfPNH/0JzUwP3L7+C3//1y5g8YXT/NBMb81y7YBrXLpgGLAKK4fDavsNsaT/Elr2H2NJ+iB++2dX3wfSMyU1cPfdirp47lWvmFaeLZjSTyzkczKxySQfBG788zL2Pb+aNXx7htqUL+aOPXcElUyaMdbf6TGzMc/3C6Vy/cHpf2bFT3by+7wjbOg6xde8htnUc5ps/2sXpnmI6TJnQwFVzL+aqORdzxeyLWDJ7Cpe3TGH65Kax2gwzG+eSDYJHf7qbP97wGhdPauRvf+/XWPavZ411lyrS3NTADYumc8OiM+FwsruHHfuP9gXD1o5D/MPGPRw/3dPXZuaUJi6fNYUls4rhcNnMySyY3szcaZNoavDlrWYpq0kQSFoB/AXFn6r8RkQ8WFavrH4lcAz4vYh4uZJ1a+bRR+FLX6Lw9h7+66p7WXflb3LzlS18/bPXjqujgPMxoSHf97lBr0Ih2PvecXZ2HWXn/qPs6DzCjs6jfPfnezlysruvXU4wZ2rxstaFM5qZN30Ssy6ayOyLJ/RNL5kygfxInW7KxoW334aFC+GBB+BznxuZ17LKXIBj0lMITvcUONVT4HR3oe8Iuj/lH7sJaMznaGzI0ZgXjbnc2Jx+HcFxqToIJOWBh4DlQDuwUdKGiHitpNkngCXZ48PAw8CHK1y3eo8+CqtXc/rESb74W/fx3SuXcefmf+K/XHMT+SlLa/pS40UuJxbMaGbBjGaWXXnmaCci2H/4JLsPvs/b7xxjz7vH2fPOMfa8c4wfvNlF15Fzvx+RE8yYPIGpkxqYOqmRac1NTJ3UeNZj8oQ8ExvzTGjIM6kpz8SGHBMb89kjR04ilxM5QV5CErknvkPu/vvJvf8+ucaJFH55gJ4//AI9p0Xh05+mpxD0FIJClE45q6y70DstnFnuCXrOqi/Q3VPePivP2ndnz9lX31NSX7Je78UVkhAU9xKAEFJxUeVlWUMJGnIin1M2zZ21nOsrz6b5XFn7wdfNZ3/fnHrnVTJf/DeRz8pyOfpv87+/TXz+Xjh+nMKEKUTnu8QffoHoFvHpf0dEEBQvkw6CQhT/TUUUr4jrzna43T3F6eme4tj0LpfWdxeK9b3Lp3oK2TrF8lPdZct9O/GS5d5Hd8mOPqs/3X1muVDja2LyORVDIZ/LHsX5poYcExqK/+YnZP8HJpxVNvB0Qsl0Yvn06Q1M+OIfMeHIIZoR+d27YfXqYmdqEAZVXzUk6SPAVyPi49nyfwKIiP9W0uavgR9ExGPZ8nbgZmDxUOv2Z9hXDS1eDLt3s2bFf+Txaz/OF3+4nrt/8m1YtAjeeqvy50nAqe4CB46epPPISfYfPkHn4RN0HjnJgaMnOXT89FmP946d5siJ7qGftE7kBA25XMmOuHTHnCOXK+5ke3eCvf91zkyLO8nest42vWWFwpmAKg2Xnlrvperc2TvYHE150VC2sz1r59vbNnvH3lSybmND2XJeNDXkaMjlznnnP9CusBBBd1n49AbbqdIgKhQ42V3g5OkCJ7t7+qYnepe7C5w4fWZ6vsP+t9/+Cst2bSouDHMfNpJXDc0D9pQst1N81z9Um3kVrguApNXAaoCFCxcOr4dvvw3AHZv+Dx/at4PPvfK9s8rtjKaGHHOzb0FXoqcQHDlxmvdP9XDidO+jcM58IXvX2BPFd96FgLj7bnqUoyARyqEokC8UyEeQe/ih4rvX7F1sPlc6X9whN+SHeledO+vddLH9ueX53nfKY3S1VZQd3fTEmaOUc452CuVBcuZop9D79y07iioEJfPFo6q+Nn3tszb33w8R5CIQgSKKRzkEeuihviMeUTyC6J2XiqdPGvp24sW/c+9yQ+7Mzr0hO73SkFfffOnOPYVLoiM7+iwPh/Lw6AuRO+7kZL6Rkw1NXH6gZJdZo31YLYKgv1Erz7qB2lSybrEwYi2wFopHBMPpIAsXwu7dfLDrLT7Y9dbZ5VaVfE5Ma25iWvN5rHxwK+zefW75okXw4UVV961eKAu1hvFwt5HOnw88JjemMyYjTTpzamlKJZeqH3qz/3Gp0T6sFpeLtAMLSpbnAx0Vtqlk3eo98AA0l+2pmpuL5TZ2PC7jj8dkfBrpcSl+0HP+D4pHFbuAy4Am4BXg6rI2vwU8S/EI4EbgZ5Wu29/jhhtuiGH71rciFi2KkIrTb31r+M9htedxGX88JuNTDcYFaIt+9qk1ucWEpJXAn1O8BHRdRDwg6a4saB7JLh/9H8AKipeP3hkRbQOtO9TrjcQtJszMLnQDfVjsew2ZmSVioCDwV0rNzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEVRUEkmZIel7Sjmw6vZ82CyT9s6TXJW2T9PmSuq9K2itpc/ZYWU1/zMxs+Ko9IlgDvBgRS4AXs+Vy3cAXIuKDFH+v+G5JV5XU/1lEXJc9nqmyP2ZmNkzVBsEqYH02vx74ZHmDiNgXES9n80eA14F5Vb6umZnVSLVBMDsi9kFxhw/MGqyxpMXArwI/LSm+R9Krktb1d2qpZN3VktoktXV1dVXZbTMz6zVkEEh6QdLWfh6rhvNCkqYA3wHujYjDWfHDwAeA64B9wNcHWj8i1kZEa0S0trS0DOelzcxsEA1DNYiIWwaqk7Rf0pyI2CdpDtA5QLtGiiHwaEQ8UfLc+0va/A3w9HA6b2Zm1av21NAG4I5s/g7gqfIGkgR8E3g9Iv60rG5OyeKngK1V9sfMzIap2iB4EFguaQewPFtG0lxJvVcA3QTcDvxmP5eJ/omkLZJeBZYB91XZHzMzG6YhTw0NJiIOAh/tp7wDWJnN/wjQAOvfXs3rm5lZ9fzNYjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHFVBYGkGZKel7Qjm04foN1b2S+RbZbUNtz1zcxs5FR7RLAGeDEilgAvZssDWRYR10VE63mub2ZmI6DaIFgFrM/m1wOfHOX1zcysStUGweyI2AeQTWcN0C6A70vaJGn1eayPpNWS2iS1dXV1VdltMzPrNeSP10t6Abi0n6ovDeN1boqIDkmzgOclvRERLw1jfSJiLbAWoLW1NYazrpmZDWzIIIiIWwaqk7Rf0pyI2CdpDtA5wHN0ZNNOSU8CS4GXgIrWNzOzkVPtqaENwB3Z/B3AU+UNJE2WdFHvPPAxYGul65uZ2ciqNggeBJZL2gEsz5aRNFfSM1mb2cCPJL0C/Az4p4j43mDrm5nZ6Bny1NBgIuIg8NF+yjuAldn8LuDa4axvZmajx98sNjNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLXFVBIGmGpOcl7cim0/tpc6WkzSWPw5Luzeq+KmlvSd3KavpjZmbDV+0RwRrgxYhYAryYLZ8lIrZHxHURcR1wA3AMeLKkyZ/11kfEM+Xrm5nZyKo2CFYB67P59cAnh2j/UeAXEbG7ytc1M7MaqTYIZkfEPoBsOmuI9rcCj5WV3SPpVUnr+ju11EvSakltktq6urqq67WZmfUZMggkvSBpaz+PVcN5IUlNwL8Fvl1S/DDwAeA6YB/w9YHWj4i1EdEaEa0tLS3DeWkzMxtEw1ANIuKWgeok7Zc0JyL2SZoDdA7yVJ8AXo6I/SXP3Tcv6W+ApyvrtpmZ1Uq1p4Y2AHdk83cATw3S9jbKTgtl4dHrU8DWKvtjZmbDVG0QPAgsl7QDWJ4tI2mupL4rgCQ1Z/VPlK3/J5K2SHoVWAbcV2V/zMxsmIY8NTSYiDhI8Uqg8vIOYGXJ8jHgkn7a3V7N65uZWfX8zWIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxVQWBpM9K2iapIKl1kHYrJG2XtFPSmpLyGZKel7Qjm06vpj9mZjZ81R4RbAU+Dbw0UANJeeAhij9efxVwm6Srsuo1wIsRsQR4MVs2M7NRVFUQRMTrEbF9iGZLgZ0RsSsiTgGPA6uyulXA+mx+PfDJavpjZmbDNxqfEcwD9pQst2dlALMjYh9ANp010JNIWi2pTVJbV1fXiHXWzCw1Q/54vaQXgEv7qfpSRDxVwWuon7KoYL2zV4hYC6wFaG1tHfb6ZmbWvyGDICJuqfI12oEFJcvzgY5sfr+kORGxT9IcoLPK1zIzs2EaMghqYCOwRNJlwF7gVuA/ZHUbgDuAB7NpJUcYbNq06YCk3efZn5nAgfNcd7zxtow/F8p2gLdlvKpmWxb1V6iI8z/LIulTwF8BLcB7wOaI+LikucA3ImJl1m4l8OdAHlgXEQ9k5ZcA/wgsBN4GPhsR75x3hyrrc1tEDHipaz3xtow/F8p2gLdlvBqJbanqiCAingSe7Ke8A1hZsvwM8Ew/7Q4CH62mD2ZmVh1/s9jMLHEpBsHase5ADXlbxp8LZTvA2zJe1XxbqvqMwMzM6l+KRwRmZlbCQWBmlrgLNggGuuNpSb0k/WVW/6qk68ein5WoYFtulnRI0ubs8eWx6OdQJK2T1Clp6wD1dTEmFWxHXYwHgKQFkv5Z0uvZnYQ/30+behmXSrZl3I+NpImSfibplWw7/rifNrUdk4i44B4Uv6/wC+BfAU3AK8BVZW1WAs9SvAXGjcBPx7rfVWzLzcDTY93XCrbl3wDXA1sHqK+XMRlqO+piPLK+zgGuz+YvAt6s4/8rlWzLuB+b7O88JZtvBH4K3DiSY3KhHhEMdsfTXquAv4uinwDTsttcjDeVbEtdiIiXgMG+MFgXY1LBdtSNiNgXES9n80eA1zlzU8he9TIulWzLuJf9nY9mi43Zo/yqnpqOyYUaBIPd8XQ4bcaDSvv5kexQ8llJV49O12quXsakEnU3HpIWA79K8R1oqbobl0G2BepgbCTlJW2meP+15yNiRMdkNO41NBYqueNpTe6KOgoq6efLwKKIOJrdzuO7wJIR71nt1cuYDKXuxkPSFOA7wL0Rcbi8up9Vxu24DLEtdTE2EdEDXCdpGvCkpGsiovQzqZqOyYV6RDDYHU+H02Y8GLKfEXG491AyirfzaJQ0c/S6WDP1MiaDqrfxkNRIccf5aEQ80U+TuhmXobal3sYmIt4DfgCsKKuq6ZhcqEHQd8dTSU0U73i6oazNBuB3s0/fbwQORfYjOePMkNsi6VJJyuaXUhzXg6Pe0+rVy5gMqp7GI+vnN4HXI+JPB2hWF+NSybbUw9hIasmOBJA0CbgFeKOsWU3H5II8NRQR3ZLuAZ7jzB1Pt0m6K6t/hOJN8FYCO4FjwJ1j1d/BVLgtnwH+QFI3cBy4NbJLC8YTSY9RvGpjpqR24CsUPwirqzGpYDvqYjwyNwG3A1uyc9IA/5niHYHralyobFvqYWzmAOtV/L33HPCPEfH0SO6/fIsJM7PEXainhszMrEIOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS9/8Bloq5Jzb+1oAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.axis('equal')\n",
    "plt.plot(node_x,np.zeros(len(node_x)),'ro')\n",
    "plt.plot(x,hermite_sf2(x,0,node_x))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assemble the Hermite interpolatory polynomial from the two shape functions defined above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hermite_int(x,node_x,node_y,node_y_prime):\n",
    "  '''Generates the Hermite osculating polynomial\n",
    "\n",
    "  :param x: locations where the Hermite polynomial need to be evaluated\n",
    "  :type x: list or np.array\n",
    "  :param node_x: x coordinates at the nodes\n",
    "  :type node_x: list or np.array\n",
    "  :param node_y: y values at the nodes\n",
    "  :type node_y: list ot np.array\n",
    "  :param node_y_prime: function derivatives at the nodes\n",
    "  :type node_y_prime: list ot np.array\n",
    "  :return: the Hermite osculating polynomial at the locations x\n",
    "  :rtype: list or np.array\n",
    "  '''\n",
    "  n = len(node_x)-1\n",
    "  res = np.zeros(len(x))\n",
    "  for j in range(n+1):\n",
    "    res += hermite_sf1(x,j,node_x)*node_y[j] + hermite_sf2(x,j,node_x)*node_y_prime[j]\n",
    "  return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare Hermite interpolating polynomial and original function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "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_int(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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}