20230324 Rational vs Float inverse Laplace works.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "76f2437d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9820cf17",
   "metadata": {},
   "outputs": [],
   "source": [
    "s=symbols(\"s\", positive=True)\n",
    "t=symbols(\"t\", real=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ee7dbf8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "H=1/s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b3cbe250",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\theta\\left(t\\right)$"
      ],
      "text/plain": [
       "Heaviside(t)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h=inverse_laplace_transform(H,s,t)\n",
    "h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "de023273",
   "metadata": {},
   "outputs": [],
   "source": [
    "G=1/s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f24b0893",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle t \\theta\\left(t\\right)$"
      ],
      "text/plain": [
       "t*Heaviside(t)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F=G*H\n",
    "f=inverse_laplace_transform(F,s,t)\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9b88c811",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\theta\\left(t\\right) - e^{- t} \\theta\\left(t\\right)$"
      ],
      "text/plain": [
       "Heaviside(t) - exp(-t)*Heaviside(t)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H=1/(s+1)\n",
    "F=G*H\n",
    "f=inverse_laplace_transform(F,s,t).simplify().expand()\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "027efeac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1.0 \\theta\\left(t\\right) - 1.0 e^{- t} \\theta\\left(t\\right)$"
      ],
      "text/plain": [
       "1.0*Heaviside(t) - 1.0*exp(-t)*Heaviside(t)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H=1.0/(s+1.0)\n",
    "F=G*H\n",
    "f=inverse_laplace_transform(F,s,t).simplify().expand()\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "028b9ba0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\theta\\left(t\\right) - \\frac{10 e^{- 10 t} \\theta\\left(t\\right)}{9} + \\frac{e^{- 100 t} \\theta\\left(t\\right)}{9}$"
      ],
      "text/plain": [
       "Heaviside(t) - 10*exp(-10*t)*Heaviside(t)/9 + exp(-100*t)*Heaviside(t)/9"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H=1000/(s+10)/(s+100)\n",
    "F=G*H\n",
    "f=inverse_laplace_transform(F,s,t).simplify().expand()\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "fcadd4d4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\theta\\left(t\\right) - \\frac{2 \\sqrt{280173} e^{- \\frac{t}{2}} \\sin{\\left(\\frac{\\sqrt{280173} t}{366} \\right)} \\theta\\left(t\\right)}{1531}$"
      ],
      "text/plain": [
       "Heaviside(t) - 2*sqrt(280173)*exp(-t/2)*sin(sqrt(280173)*t/366)*Heaviside(t)/1531"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sympy as sympy\n",
    "#wo^2=2.34153145345 chosen to break the inverse laplace transform.\n",
    "H=(s**2+2.34153145345)/(s**2+s+2.34153145345)\n",
    "F=G*H\n",
    "#Use Apart first.  Apart can turn rationals into floats\n",
    "F=apart(G*H,s)\n",
    "# We ahve to turn the floats into rationals or the inverse laplace might not work.\n",
    "for i in F.atoms(sympy.Float):\n",
    "        r = sympy.Rational(i).limit_denominator(1000)\n",
    "        F = F.subs(i, r)\n",
    "f=inverse_laplace_transform(F,s,t).simplify().expand()\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "65744f39",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\theta\\left(t\\right) - 0.691461406407671 e^{- \\frac{t}{2}} \\sin{\\left(\\frac{\\sqrt{280173} t}{366} \\right)} \\theta\\left(t\\right)$"
      ],
      "text/plain": [
       "Heaviside(t) - 0.691461406407671*exp(-t/2)*sin(sqrt(280173)*t/366)*Heaviside(t)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f=sympy.N(f)\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "96b23506",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVQAAAD1CAYAAAAPpQ1KAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAyi0lEQVR4nO2deXwV5fX/34dAANkFDEhYxQ1QQMCtLMFdakVt1aLiShHXassPa7WVqt/Woha1LriA4oLWKq5FVJSg1AUQARHEImtkk00MAiHh/P545pLJ9d5kkszdyHm/XvO6M/PMzPO5c2fOPc92HlFVDMMwjOpTK9UCDMMw9hXMoBqGYYSEGVTDMIyQMINqGIYREmZQDcMwQsIMqmEYRkiYQTUqRESuEpH1IlIoIs2TkN9TInJniNdbISInhXW9RCEiKiKd00DHpSIyM9U6qoOI9BORJb7tpDwDGWVQRaSviHwkIt+LyGYR+a+I9PHSEvoQiEi+iOz0jMpGEZksIq0TlV+6ICJ1gH8Ap6hqQ1XdFJWeJyIFMc7LF5FhydK5r+N7/tr69p0kIisCnh/qn1S6ICIdvD+iQt8yX1U/VNVD45wzWkSeTYSejDGoItIYeBP4J7A/0Ab4C7AriTKuVdWGQGegIXBPEvNOFTlAPeDLVAupCBGpnWoNCWY78KdUi0hTmnp/+A1VtXsiMyrvOcsYgwocAqCqz6tqiaruUNV3VHWBiBwOjAOO8/6htgKISF0RuUdEVnlF1nEiUt9LyxORAhH5o+dxrhCRC4MIUdWtwKtAj8g+ETlMRN71POclInKeL+0pL+93ReQHEZkhIu196ceLyGzP854tIsf70vJF5A7PG/9BRN4RkRZeWj0ReVZENonIVu/cHC+tiYiMF5G1IvKtiNwpIlmxvo93n+4TkTXecp+37xAgUmzaKiLvB7k/Ma4/WkReFJGnve/wpYj09qX3FJG5Xtq/cAbcf/4ZIjLP+44ficiRvrQVInKTiCwAtotIbREZKiIrvftyS9S1jhaRj71rrRWRB0Uk25euInK1iPzP03OHiBzknbPN+x7Z3rHlPkPlPX9e+v/zNKwRkcsD3MoHgCESp1pARA73npet3j0+09s/HLgQGOW9H294+9uKK2l9592rB6Oud4+IbBGR5SJyejxR3ve+WUQWecc/KSL1vLRmIvKml8cWbz3Xd+6lIrLMu9fLI/dPRDqLe0++9+7tvwLcH7+meCWn04A/Aud792K+tz/u++Jp/K+IjBWRzcDouBmrakYsQGNgEzAROB1oFpV+KTAzat99wOs4j7YR8AbwNy8tDyjGFWfrAgNwHsChcfLPB4Z5682BacBr3nYDYDVwGVAbOArYCHT10p8CfgD6e3ndH9HqadsCDPXOHeJtN/fl+w3uD6W+t32Xl3al9532A7KAXkBjL+1V4FFP2wHALODKON/tduAT77iWwEfAHV5aB0CB2nHOzQMKKrhfo4GdwCBP59+AT7y0bGAlcCNQB/gVsBu400s/CtgAHOOdewmwAqjrpa8A5gFtvfvTBSj03et/eL/zSd7xvYBjvXvdAVgM3ODTrbhnpjHQFVcCeg/oBDQBFgGXBHmGKP/5Ow1YD3TzfqNJXt6dy3v+vLye9fadBKzw1usAS3HGIhs4AffMRbQ8Fbmn3nYWMB8Y6+VfD+jre5d2A7/xjrsKWANIHG0rgIXeb7A/8F/f79cc+CXuGW0E/Bt41ffebPNpbE3pO/M8cAvO6durLUbeHYjxfBL1XHoaI8/A6Mg99KW/Spz3xbsfxcB1uOemflw7lWpDWZkFONx7MAq8L/g6kBPLoAKCe7gP8u07Dlge9TI08KW/CPypnAf6R+B77wecB7Tz0s4HPow6/lHgNt/D/IIvrSFQ4j2AQ4FZUed+DFzqy/dWX9rVwFRv/XKc8Tsy6vwcnCGo79s3BJge57t9AwzybZ9K6Ysa84GN9+BG3S+/QZ3mS+sC7PDW+xP1snrfKfJCPoJn3H3pS4ABvhflcl/an6PudQOgCO9liqHzBuAV37YCP/Ntfwbc5Nu+F7ivomeIip+/CXh/jN72IQQzqC1xz2BXyhrUfsA6oJbvnOeB0b5n8M4oLd/F+l1x79JS3/Z+nrZWcbStAEb4tgcB38Q5tgewxffbbMUZ3PpRxz0NPAbkxrqO77gOnratvmUklTCoVPC+ePdjVXk6IksmFflR1cWqeqmq5uL+2Q/EeQGxaIl7ED7zikBbgane/ghbVHW7b3uld814XK+qTYAjgWZApOjSHjgmko+X14VAK9+5q33foxDY7OV1oJevn5W4OuII63zrP+IMMsAzwNvAC16xcYy4RqT2OI9lrU/Po7h/3lhEa6joPvgp9vKKpg7Oy4n3HeqJq4s6EPhWvSfXl3+E9sDvo+5t2yh9q33rB1L2Xm/HlWwAEJFDvGLnOhHZBvwVaBGlfb1vfUeM7Ya+7XjPUEXPXxmd/PQZiImqfgc8iCtV+DkQWK2qe6Ku2YbYtAVWqmpxnPS9v5eq/uitNoxzLPz0uxwIICL7icijXhXMNuADoKmIZHn37XxgBO5Z/Y+IHOZdYxTuT2mWV31RUZVIC1Vt6i2VbdsI8r6sjnViNBllUP2o6le4f91ukV1Rh2zEPfxdfTe6ibpGpQjNRKSBb7sdzluqKO8vgDuBh0REcDd7hi+fpuoqx6/yneZvnW2IKxqt8Zb2lKUd8G0AHbtV9S+q2gU4HjgDuNjTs4uyD1ljVe0a51LRGgLdB49VQAvvOwHg3ZP2BDMSa4E23jn+/COsBv4v6t7up6rP+47x//ZrKXuv98MVOyM8AnwFHKyqjXFFZH/elSXeM1TR81dGJ2W/c0XcDQzEVV9EWAO0FRH/O+1/jqLfj9VAOwmvIS/6u0Sen98DhwLHePe7v7dfAFT1bVU9GVfc/wp43Nu/TlV/o6oH4qq2Ho5Xd1wFYt2Lit6X6HNikjEGVVyjz+8jFdriuo8MwdX9gfMiciMNBt4/9ePAWBE5wDunjYicGnXpv4hItoj0wxmkfweUNBH3D3YmrvfBIeIaQ+p4Sx9xjWURBonr9pUN3AF8qqqrgSneuReIa1A5H1ckfjPAPRkoIkd4lefbcB5hiaquBd4B7hWRxiJSS1zDyoA4l3oeuFVEWopr8PozEKhbiaquAj4F/i4iDUWkLvD/cJ7rJ+We7PjYO/Z67/ufAxztS38cGCEix4ijgYj8XEQaxbneS8AZvnt9O2Wf80a4e1XoeUNXxbhGZfnJMxTg+XsRuFREunhG/7agmalrFL0X58VF+BRXxTDKe/7ygF8AL3jp63H1wBFm4Yz6Xd49rSciP6vMl47iGhHJFZH9cX9SkUakRrg/lq1e2t7vKSI5InKm94e0C1f3XeKlnSuljVdbcAatpBr6/KwHOkT+fKrwvsQlYwwqroL9GOBTEdmOe1kX4v4BAd7Hde1ZJyIbvX034SrqP/GKG9Nw/5YR1uF+rDXAc7h6oK+CiFHVIlyr659U9QfgFODX3rXWAX/HNVREmIR7mDbjPIsLvetswr2Ev8cVTUcBZ6jqRiqmFc6AbMM1rsyg1BBejGucWOR9x5dwXkAs7gTmAAuAL4C53r6gnI/7c1mK84hOxNXJ7qzoRO8+noOrp9riXWuyL30OrnHkQS99qXdsvOt9CVyDu99rvXP8rb0jgQtwz9PjlL74VaW8Zyju86eqb+Gqq973jqlsD4r78RkY7z6eiWuw3Qg8DFzs0zIe6OIVaV9V1RKcwe2MK2UU4O59VZmEM0rLvCXy/NyHayzciHtnp/rOqYV77tfg3osBuDYCgD64d70Q11byW1VdXg19fiJO0yYRmeutV+Z9iYuUrbqqOXj/4M969bGJzuspXAX5rYnOy0geyXyG0hlxgwuGqeq0VGtJNZnkoRqGYaQ1STGoIjJBRDaIyMI46ReKyAJv+UhEuidDl2EYRpgkpcgvIv1xFc5Pq2q3GOnHA4tVdYu4ERmjVfWYhAszDMMIkaSMfVbVD0SkQznpH/k2P6G0f6dhGEbGkI7BJK4A3oqXKG5c8nCA+vXr92rbtm28Q5POnj17qFUrM6ulM1W76U4+maq9Krq//vrrjarasuIjPYIMpwpjwQ0RW1jBMQNx3X+aB7lmr169NJ2YPn16qiVUmUzVbrqTT6Zqr4puYI5Wws6ljYcqLoLQE8DpGhVz0zAMIxNIC79dRNrhOnMPVdWvU63HMAyjKiTFQxWR53HRX1qIi1F4G15ADVUdhxvq2Bw3XhegWFV7x76aYRhGepKsVv4hFaQPw4Umqza7d++moKCAnTsrHPUYOk2aNGHx4sVJzzcM/Nrr1atHbm4uderECiJlGEY80qYONSwKCgpo1KgRHTp0oGwAo8Tzww8/0KhRvJgd6U1Eu6qyadMmCgoK6NixY6plGUZGsc8Z1J07d6bEmO4riAjNmzfnu+++q9Z1duyAhQth+XLYvh0aN4a2baF7d6hbt+LzDSMT2ecMKmDGtJpU9f6pwhtvwPjxMHUqFBX99Ji6daFvX7jwQvjVryBDHXrDiMk+aVCN5PPJJ3D11fD5525bBLp2hcMOgwYNYNs2+PprWLQI3nvPLb//PdxwA1x/PTRtmkr1hhEOadFtyshcSkrgllvgZz9zxrR1a7jnHvj2W1fkf+klmDgRXnkFvvwSNm2CJ56A446DLVvgttvg4IPhqaech2sYmYwZVKPKFBbCWWfBX//qtm+6Cb75xnmereOE5t1/f7jiCvjoI8jPh379YONGuOwyyMuDFSuSo90wEoEZ1CSzY8cOBgwYQEmJC7ZeUFDAv/7lgsYXFRXRv39/iovLzpumnus2evToMtsA06ZNY+jQoaHlFZTt2+HnP4c333RGcto0uOsuqF+/4nMjDBgAM2bAM8/AAQfABx+4RqsXXqj4XMNIR8ygJpkJEyZwzjnnkJWVBcB7773H3LluFobs7GxOPPHEvUYvwnPPPceYMWPYuXMnY8aM4bnnntubNn/+fHr27BlaXkHYtQsGD3YGsE0b+PhjGDiw0pcBXF3rRRfB4sVw9tmurnXIEBg+3OVjGJmEGdQEsXjxYvr378+RRx7J3XffTefObsLG5557jsGDBwMwc+ZMfve73/HSSy/Ro0cPli9fzllnnVXGYAJcdNFFtG3bljFjxtCuXTsuuuiivWnz589n3bp19OvXj1atWjFtWuksFFXJqyJUYcQI16jUqhW8/z4cckiVblEZ9t8fXn4Zxo2DevXg8cfhhBNg/fqKzzWMdGGfNqgiiVkqori4mAsvvJD777+fBQsWsGzZMrp160ZRURHLli2jQ4cOAPTt25c+ffrw2muvMW/ePDp27Ei3bt2YPXt2metNmjSJ1atXM2rUKFatWsWkSZP2ps2fP58WLVrw4Ycf8vDDD+81kFXNqyIeeMA1IO23H7z1VjjGNIIIXHklzJwJubmunrV3b/jss/DyMIxEYt2mEsDkyZPp3r373qJ4ly5dOOCAA9i4cSNNo/oHLVmyhEMPLZ2INSsri+zs7DKjroYMGYKIMHr0aEaNGrW3DnX37t1s3ryZkSNHAs6QR65f2byCMHu2a3ACZ1R79Ah0WqXp1QvmzIFzznFGtV8/V6965pmJyc8wwmKf9lBVE7NUxIIFC+jhszYLFy6kR48e1K9fv0yMgU2bNtGkSZOfjJnftWsX9erV27sd6WgfaZSKbC9atIju3bvvDZq7YMECunVzM8xUNa94/Pijq+ssKYEbb4Rzz63wlGqRk+OqEy67zI26OvtseOSRxOZpGNVlnzaoqaJ58+Z8/bWLQjhv3jyeffZZunfvTrNmzSgpKdlr6JYvX86BBx5Y5txNmzbRsmXLQIFJ5s+fT/fupfMZLliwgCOPPBIg9Lxuusl1zO/atbSbVKKpW9eNuho9GvbscQMHbr7ZrRtGOmIGNQEMHTqUOXPm0KdPHyZMmECHDh3o1KkTAKeccgozZ84E4LDDDmPjxo1069aNjz5y02pNnz6dQYMGBcpn/vz5ew0oOE844qGGmdd//wsPPgh16rguTgEc2tAQcZ3/x4+HrCzXNevii2MPazWMlFOZ8P7ptsSaAmXRokWVnOQgPLZt26aqqj/88MPefWPGjNFbbrll7/bcuXP1oosuinuNs88+W7/66qtQ9FQmr4j2CJH7uHu3avfurrLj1ltDkVVlpk5VbdjQaRk4UHXLlpo1HUe6kKnakzEFinmoCWDs2LF07dqVHj16sGLFCv70pz/tTevZsycDBw7c29neT1FREWeddVaZhqPqEEZe48bB/PnQvr0rbqeSU091AwFatYLp012QlXXrLHSVkUZUxvqm25KuHmomEstD/e471SZNnEf4yispkRWT5ctVu3Rxupo126WzZqVaUeXJVC9PNXO1m4dqpJS77oLvv4dTTnEjo9KFDh1cve6JJ8KWLdkMGOCCrxhGqjGDasSkuNg1RIEzrOkWYrZpUzewYNCgtezYAb/8Jdx7r0WsMlLLPmlQ1d6qaqGqbN3qxtKffz7ECRWQcurUgZEjl/DXvzpDOnKkxQAIgqr7w7TXJHySNevpBOAMYIOqdouRLsD9wCDgR+BSVZ1blbzq1avHpk2baN68uUXurwKqytq1m5g3rx5ZWXDHHalWVD4irrGsUye45BIXa3X+fBeHtV27VKtLDps2wbJlbrqZFStg5UrYsMHtjyyFha6r2e7dzphGyM4uXerXd55/s2aln82aQfPmriGwVSsXlnHDhroUFblzjLIka+jpU8CDwNNx0k8HDvaWY4BHvM9Kk5ubS0FBQbXnRKoKO3fuDDTqKB3xa//qq3rcdlsul1/ugj9nAuefD507u2lVZs+Go46CSZNc/e++gir8738waxYsWOD+OBYsgHXrKn+tWrXcAImiorJ9er/9NsjZxwEuoE3r1qXGtlUrN8ItertFC5dfTSBZ00h/ICIdyjlkMPC016r2iYg0FZHWqrq2snnVqVMnJbN1rlsHL774Oddff3jS8w6D/Px8evbsyapVcN557uVNdTepytKrlwukctFFrn711FPdMNn/+7/KxWlNF4qL3feZOdM1ws2cCbH8hIYN3Z9Jx46uwa5DB2fM9t/feZfNm7u5uyKeaO3azrP3G9SiIje8eOtWt2zZ4patW12e69fD2rXuOV+1ahdbttRl82bYvNnNxFAeWVku3q3f4LZs6bzgJk3cp3+9SRP3e9Wt6waRRPRWh5ISKCqqRUmJ05Mo0iU4ShtgtW+7wNv3E4MqIsOB4QA5OTnk5+cnQ1+F/OEPR/Dppz1p1GgWHTv+mGo5laawsJD8/HweeKAzxcW5nHjielauXMzKlalWVj4R3X5GjoTWrdvz1FMdGDtWmDx5Ozff/BWHHhosCEwyiKUb4Lvvspk9e39mzdqfzz5rRmFh2WHBzZoV0a3b9xx0UKG3bCcnZ2dcD/D7791SWZo0cUs87fXrN2Tbtjps3pzNpk3ZbN6czZYt7tO/vmVLNtu21WHtWmeQq0KtWkqdOnvIzt5DnTp7KvR2S0qE4uLIUoviYkFVgP78859z6dZtW9WEBCBdDGqs/5+YVeaq+hjwGEDv3r01Ly8vgbKCUzq+/GjSRFKlyM/Pp2vXPKZOddtjx+ZwxBE5qRUVgPz8fGI9AyecAFddBUOHwldfNeDqq3tx9dVw++2uXjDVRHQXFTnv86233EyxX3xR9rjOnd3UMH37uqVTp2xEWgItUyEbiH/P41FU5Op0161znu66dc7rjRj7iFe8dWvp9s6drnFxxw7Ys0fYtSuLXbuq7lqKQO3ae+je/Sj69avyZSokXQxqAdDWt50LrEmRlioRaTGtqPiTzjzwgHuAzzgDjjgi1WqqT+/eMHeum0Tw/vtdN7Dnn3cNbZdf7oqUqWDDBpg6tRUPPQRvvw3+6IkNGrg/g9NOc1UWBx2UGo1hkp3t4tvm5lbt/OJiZ1x37nRLeb0TVF3vD/+Sne2K+fn5H9CvX17VRAQkXQzq68C1IvICrjHq+6rUn6YDCxemWkHV2LGj1t5+p5lWd1oe9evDP/4Bl17qpqueMcNFrfrrX10Ercsuc0YskajCvHlu/q0333SNZqqH7U3v0gUGDXJGtG/f1Bn6dKV2bbck+ncKg2R1m3oeyANaiEgBcBtQB0BVxwFTcF2mluK6TV2WDF2JIFM91HffzWHrVje98/HHp1pN+Bx5pBv///LL8Je/uD++665z3usFFziD26dPOK3RqrB0qcvv/ffd54YNpel160L37pu4+OLm/PznrhHJ2DdIViv/kArSFbgmGVoSRaQYUlDg6oCiguWnNarwyiuuPHbddSkWk0BEXLeqc86BV1+Fe+5xEwyOG+eWnBznKfbv73oMHH6484zKY/duWLXKTTL42WeuimHOHFgTVWF14IFultgzznBDZmfP/qJS9ZBGZpAuRf59ikWLMsvLy8+HFSsa0KqVG8K5r1OrljOq55zjPNXx453nuno1PPmkW8AZ0zZt3NK4sfMsa9VydZ7ff+8aWAoKYge8btHCNSadcIKbEfbQQ9Nv+K4RPmZQQ8JfUb5wYWYZ1H/+032OGFHzRr906wZjx7p61i++cI1Es2c7L3P5cjfqqLyuYyKusaVzZzeYoFcvtxx8cM3pzG6UYgY1AWRSPerKlfDaa65LyZVX1lwLIOLqWX0TILBzpyu6FxSUDt0sKXGd5Bs3dh3m27WzRiSjFDOoIRHtoWYKjzziiqwDB35Hq1bp3+80mdSr52IEeLPXGEaF1FyXJIFkioe6c6cLJgJwzjmBBnEbhlEOZlBDwu+hrl8fe8x1uvHqqy4SUc+ecPjhiRuOZxg1BTOoIRMJNpUJXur48e5z2DBrgTaMMDCDGjKHe8Gm0r0edcUKmDbN/QFccEGq1RjGvoEZ1JCIFPkjY+AXLEidliBE+lr+8peZNQjBMNIZM6ghE5kuZN68lMool5KSUoN6xRWp1WIY+xJmUEMi4qFGDOoXX5SdaiKdmDbNjQrq1AkGDEi1GsPYdzCDGjJNm7rI6Tt3wldfpVpNbCKNUZdfbqN5DCNMAr1OItJFRHK89YYi8hcR+bOI7JdYeZmDv9tUxEv9/PPUaCmPjRtdd6latdykdoZhhEdQ/2QS0NRbvwfoj5up69EEaMpoRNLboD77rIuQdOqpVQ/4axhGbIIOPe2gqku86Z7PBroCO4DlCVOWYWSCh6paWty3xijDCJ+gBnWXiDQCugCrVXWjiNQGMnPO5ATjb+lXTZ9O87Nnu/6xLVvCL36RajWGse9RmSL/+8BE4Clv31GYh7qXiIcq4uYqP+AAF2h6xYpUqipLxDsdOrTmhekzjGQQyKCq6o3ALcBVqurNPMQe4MZECctk0rEedft2N0EdWHHfMBJF0Fb+B1T1HVWdHtmnqnOAMxOmLEOJFO/TzaC+9JKLNH/ssW5SOMMwwidokf/SOPuHhqQj44me2rZ3b/c5a1bytcTCGqMMI/GU2yglIpdHjvOtR+gEbAyakYicBtwPZAFPqOpdUelNgGeBdp6ue1T1yaDXTxciHuoxx7jPWbNcAOdUdqD/+mv48EM3De/556dOh2Hs61TUyh/xQLMp640qsB4I1DVcRLKAh4CTgQJgtoi8rqqLfIddAyxS1V+ISEtgiYg8p6pFQfJINdEeam6um9zt22+dQTvssNjnJYPIuP3zznPTdxiGkRjKNaiqOhBARO5U1Vurkc/RwFJVXeZd7wVgMOA3qAo08vq6NgQ2A2k6Gj4+/i5SxxwDkyfDp5+mzqAWF8PEiW7divuGkVjiGlQREdW9ftefRSRmoVVVY0yi+xPaAKt92wXAMVHHPAi8DqwBGgHnx7q2iAwHhgPk5OSQn58fIPvEs317H6ABs2bNYsOGHwFo0aItcBCTJ39L+/b/S4mujz5qztq1R9C27Y8UFc0i3u0qLCxMm3tZGUx38slU7UnRraoxF2Cbb30PUBK17AFK4p0fda1zcfWmke2hwD+jjvkVMBYQoDOuj2vj8q7bq1cvTRcOO0wVVBctKt2Xn+/29eyZOl2DBzsNY8aUf9z06dOTISd0THfyyVTtVdENzNEANi6ylFfk7+pb71hNu10AtPVt5+I8UT+XAXd5X2KpiCwHDgPSpJ28fKLrUMG19Neq5YJN//gj7JfkUDJr18Kbb0Lt2q4zv2EYiSVu27Oqrvatr4y3BMxnNnCwiHQUkWzg17jivZ9VwIkAXmSrQ4Fllfky6UaDBi6Cf0kJfPZZ8vOfONHl/YtfQKtWyc/fMGoagcbyi8j+wEigB67BaC+q2r+i81W1WESuBd7GdZuaoKpfisgIL30ccAfwlIh8gSv236SqgbtlpQvR4/aPOw7mz4eZM6Ffv+TpUC2dInrYsOTlaxg1maDBUSYBdYEXgR+rkpGqTgGmRO0b51tfA5xSlWunA7GK/OAi4o8bBzNmwM03J0/PjBnwzTeu69appyYvX8OoyQQ1qMcDLVV1VyLF7AtEe6iRKUZmznRxSOvUSY6OiHd6+eWQlZWcPA2jphN0/M4CXEOSEYd4Hmrr1nDIIS44ydy5ydGyZQu8/LIz7pdHj28zDCNhBPVQ3wemisiTwDp/gqpOCF1VBhMr9umAAW601IwZpUNSE8mkSW5Oq5NPhg4dEp+fYRiOoB5qP1zXp5NxfUgjy0UJ0pVxxPNQobTYn4y+0Krw+ONu3UZGGUZyCeShqjcE1aiYeB4quHrU4mLXLzRRfPaZ61Ww//5w1lmJy8cwjJ8SNB5qrXhLogVmCuV5qLm50Lmzi0c6e3ZidTz0kPu85BKoWzexeRmGUZagBrEY2B1nMXzEmz/qtNPc51tvJS7vTZtKo/JffXXi8jEMIzZBDWpHXPzTyPIz4A28ICVGxZx+uvtMpEEdPx527XJ5de6cuHwMw4hN0DrU6CGmK0XkEtyQ0vGhq8pAyivyA+TluSL4nDmwfj3k5ISbf0kJPPywW7/mmnCvbRhGMKpTB9oYaBmWkH2FeEX+/faDgV7T3ttvh5/vf/4DK1dCp06l1QuGYSSXoGP5n8EFgI6wH9AfN2WJQcUeKrii+NSp8MYbcPHF4eb/oDcX7dVX28gow0gVQT3UpcA3vuUT4AJVvS5RwjKVeB4qwODB7nPKFDdyKiw+/xzefdd5wZddFt51DcOoHEHrUP+SaCGZThAPtX17N43zJ5+4Ivp554WT9913u8/f/Mb1PzUMIzVYP9KQKc9DhVIj+uKL4eS3fDn8619usMDvfhfONQ3DqBpmUEMiiIcKcO657vM//3Ed/avLvfe6aaqHDIF27ap/PcMwqo4Z1JCpyEPNzYW+fV3wkup6qatXl4bpGzWqetcyDKP6BB162jzRQjKdoB4qlEbQf+yx6uV5++2uI/9550G3btW7lmEY1Seoh7paRF4TkV95c0IZ1eDcc6FJE5g1C+bNq9o1vv4annzSdZG6445Q5RmGUUWCGtT2wHvATcA6EXlMRPomTlbmUlGRH1z3pkg/1Pvuq1o+N9/sRkddeqkLYG0YRuoJZFBV9TtVfUBV+wDHARuAZ0RkmYjcLiLtE6oyA6hMkR/ghhucd/nss7BiReXOnToVJk92s6qOHl25cw3DSBxVaZRq5S2NcZ382wCfi8gfyjtJRE4TkSUisjTesSKSJyLzRORLEZlRBW0pJ4iHCm6I6JAhzsv829+CX3/HDrj2Wrc+erRr5DIMIz0I2ijVVUT+JiKrgEeA/wFHqurJqnoFcBTwx3LOzwIeAk4HugBDRKRL1DFNgYeBM1W1K3BuFb5Pyqishwrwxz86L/WJJ4LXpY4a5WYz7doVfvvbyudpGEbiCOqhfgA0An6lql1U9e+q+m0kUVVXAPeVc/7RwFJVXaaqRcALwOCoYy4AJqvqKu+aGwJqSyuCeqgAhx/uvM09e9wY/JKS8o9/9VU3Zr9OHXj66eTNoGoYRjAqHHrqeZePArer6s54x6nqn8u5TBtgtW+7AIieru4QoI6I5OOM9/2q+nQMPcPx4rDm5OSQn4yJmgKwc+dxQF0+/vhjvvkm+GzbJ52UxTPPHM3HH9flkktWMmzY8pjHLV7ciBtv7AFkMWzYUrZtKwh1jqrCwsK0uZeVwXQnn0zVnhTdqlrhAmwEagU5Ns755wJP+LaHAv+MOuZBXNCVBkALXLXCIeVdt1evXpoutGmjCqqrV1f+3PffV61Vy51/zz2qe/aUTX/nHdXGjV365Zf/ND0Mpk+fHv5Fk4DpTj6Zqr0quoE5WglbF7TI/zQwohp2uwBo69vOBdbEOGaqqm5X1Y24aobu1cgzqVSlDjXCwIGlwaFHjoQzzoCXXoLXXoOhQ+GUU2DbNtd/ddy4ylUrGIaRPILOv3k0cJ2IjMIV3feaD1XtH+D82cDBItIR+Bb4Na7O1M9rwIMiUhvIxlUJjA2oL22oqrG78kpo3BhGjHDh/aZMKU2rXRtuvRX+9CeoZYOFDSNtCWpQH/eWKqGqxSJyLfA2kAVMUNUvRWSElz5OVReLyFRgAbAHV0WwsKp5ZiJDhripUsaPdyH+SkrgqKPgiitcNyvDMNKboPFQJ1Y3I1WdAkyJ2jcuavtu4O7q5pUKqlPk99O6tfNGDcPIPIJ6qIhIDq7o3wLYW7BV1QkJ0JWxWP2mYdRcgs4pdRZu/qj/AV2BL4FuwEzADCrheaiGYWQuQZs47gQuU9WewHbvczjwWcKUZSjmoRpGzSWoQW2nqv+O2jcRCHnuzszFPFTDMIIa1A1eHSrAChE5DjgI12Jv+DAP1TBqLkEN6uNAJP7pWGA6MB8XzMTAPFTDMIK38t+tqnsAVPVpb7x9A1VdnDBlGYp5qIZRcwkaHKVQRJqq6i4A9SJCGYZhGKVUWORX1RLga8Am6iuHSJHfPFTDqLkELfI/B7wpIvfjgpj4x/K/nwhhhmEYmUZQg3qV9zk6ar8CNsoca5QyDCP4WP6OiRayr2BFfsOouQSdU+q1OPsnhysnczEP1TCMoP1QB8bZnxeSjn0G81ANo+ZSbpFfRG73VrN96xE6ASsToioDMQ/VMIyK6lAj05bUouwUJoqL3D86AZoyGvNQDaPmUq5BVdXLAETkI1WtcsR+wzCMmkBcgyoiHVR1hbf5nojE7B6lqssSISzTsI79hmGU56F+ATTy1pfiivnR5kKxiFOGYRhAOQZVVRv51m2uzQowD9UwjAoNpYhkicg3IlK3OhmJyGkiskRElorIH8o5ro+IlIjIr6qTn2EYRrIJGhylBKhf1Uy8iFUPAacDXYAhItIlznF/x003nVFYtynDMIIW5e8D/iUiA0TkIBHpFFkCnn80sFRVl6lqEfACMDjGcdcBLwMbAl437bAiv2HUXIIGR3nQ+zw5an/QRqk2uH6rEQqAY/wHiEgb4GzgBKBPvAuJyHDcBIHk5OSQn58fIPvEU1zcF6jNzJkzadiwONVyKk1hYWHa3MvKYLqTT6ZqT4buoMFRqtsoFctviy4k3wfcpKolUo6bp6qPAY8B9O7dW/Py8qopLRyyvL+Vfv360qRJarVUhfz8fNLlXlYG0518MlV7MnQH9VCBvV7kgcC3qrqmEqcWUHakVS4QfX5v4AXPmLYABolIsaq+WhmNqcLqUA3DCBptqp2IfIgbu/8fYJWIzBSR9gHzmQ0cLCIdRSQb+DXwuv8AVe2oqh1UtQPwEnB1phhTP1aHahg1l6BF+YnAZ0ATVT0AaIozkhODnKyqxcC1uNb7xcCLqvqliIwQkRGVVm0YhpGGBC3y9wJOUdXdAKpaKCI3AZuCZqSqU4ApUfvGxTn20qDXTResY79hGEE91E9wXZ/89AY+DleOYRhG5lJecBR//NNvgCki8h9c96e2wCBgUmLlZQ7moRqGUV6Rv23UdmS6kwOAXcArQL1EiDIMw8hEyguOEomFWgs31cl/VXVXknRlHNZtyjCMIGP59wCvmTENhhX5DaPmErRR6gMROTahSjIc81ANwwjabWol8JY3nfRqfMNGVfXPiRCWqZiHahg1l6AGtT7wqreemxgphmEYmU3Q4CiXxdrvNVgZWLcpwzCC16GWQUSOEJG7cUFPDMMwDCphUEWkpYj8VkTmAvNwI6d+myhhmYZ5qIZhlFvkF5E6wJnApcCpuNlPnwfaA+eqasZG1jcMwwibijzU9cCjwBLgWFXtoqp3AEUJV5ZhmIdqGEZFBnUBLlTfMUAfEWmWcEWGYRgZSrkGVVXzgIOAd4CRwDoReQNoANRJuLoMwjr2G4YRZOjpSlW9Q1UPBk4E1gJ7gPkiMibRAjMNK/IbRs2lUt2mVHWmqg4HWuGmfD4iIaoMwzAykCr1Q1XVnar6vKqeHragTMUapQzDsJFOhmEYIWEGNSTMQzUMI2kGVUROE5ElIrJURP4QI/1CEVngLR+JSPdkaTMMwwiDpBhUEckCHgJOB7oAQ0SkS9Rhy4EBqnokcAfwWDK0hYV5qIZhJMtDPRpYqqrLVLUIeAEY7D9AVT9S1S3e5idYmEDDMDKMoPFQq0sbXGDqCAW40VfxuAJ4K1aCiAwHhgPk5OSQn58fksTqoToAEGbMyKdWBtZMFxYWps29rAymO/lkqvZk6E6WQY1VEI45tkhEBuIMat9Y6ar6GF51QO/evTUvLy8kieGQl5eXkQY1Pz+fdLuXQTDdySdTtSdDd7IMagFlp6XOBdZEHyQiRwJPAKer6qYkaQsFG3pqGEayfKnZwMEi0lFEsoFfA6/7DxCRdsBkYKiqfp0kXaFjjVKGUXNJioeqqsUici3wNpAFTFDVL0VkhJc+Dvgz0Bx4WJxVKlbV3snQZxiGEQbJKvKjqlOAKVH7xvnWhwHDkqUnUZiHahg1lwxsPjEMw0hPzKCGgDVIGYYBZlANwzBCwwxqCJQOOzVX1TBqMmZQDcMwQsIMaghYYBTDMMAMqmEYRmiYQTUMwwgJM6ghYN2mDMMAM6iGYRihYQY1BKzblGEYYAbVMAwjNMyghoB1mzIMA8ygGoZhhIYZ1BAwD9UwDDCDahiGERpmUA3DMELCDGoIWMd+wzDADKphGEZoJM2gishpIrJERJaKyB9ipIuIPOClLxCRo5KlrbpYx37DMCBJBlVEsoCHgNOBLsAQEekSddjpwMHeMhx4JBnaDMMwwiJZs54eDSxV1WUAIvICMBhY5DtmMPC0qirwiYg0FZHWqro2DAG7d8MDD4RxpdjXBus2ZRg1nWQZ1DbAat92AXBMgGPaAKEY1OJiGDkyjCvFp27dPUBWYjMxDCNtSZZBjeW7RVc4BjkGERmOqxIgJyeH/Pz8QAKKi4Vzz+0U6Niq0rXrGvLzdyQ0j0RRWFgY+F6mE6Y7+WSq9mToTpZBLQDa+rZzgTVVOAZVfQx4DKB3796al5cXWMRJJwU+tErk539DZfSkE/n5+Rmp3XQnn0zVngzdyWrlnw0cLCIdRSQb+DXwetQxrwMXe639xwLfh1V/ahiGkQyS4qGqarGIXAu8jatknKCqX4rICC99HDAFGAQsBX4ELkuGNsMwjLBIVpEfVZ2CM5r+feN86wpckyw9hmEYYWMjpQzDMELCDKphGEZImEE1DMMICdEMDpUkIt8BK1Otw0cLYGOqRVSRTNVuupNPpmqviu72qtoy6MEZbVDTDRGZo6q9U62jKmSqdtOdfDJVezJ0W5HfMAwjJMygGoZhhIQZ1HB5LNUCqkGmajfdySdTtSdct9WhGoZhhIR5qIZhGCFhBtUwDCMkzKCGREVzZqUjItJWRKaLyGIR+VJEfptqTZVBRLJE5HMReTPVWiqDNxvFSyLylXfvj0u1piCIyI3ec7JQRJ4XkXqp1hQPEZkgIhtEZKFv3/4i8q6I/M/7bBZ2vmZQQyDgnFnpSDHwe1U9HDgWuCZDdEf4LbA41SKqwP3AVFU9DOhOBnwHEWkDXA/0VtVuuKhxv06tqnJ5Cjgtat8fgPdU9WDgPW87VMyghsPeObNUtQiIzJmV1qjqWlWd663/gHux26RWVTBEJBf4OfBEqrVUBhFpDPQHxgOoapGqbk2pqODUBuqLSG1gP2IEgE8XVPUDYHPU7sHARG99InBW2PmaQQ2HePNhZQwi0gHoCXyaYilBuQ8YBexJsY7K0gn4DnjSq654QkQapFpURajqt8A9wCrcPG/fq+o7qVVVaXIiQeu9zwPCzsAMajgEmg8rXRGRhsDLwA2qui3VeipCRM4ANqjqZ6nWUgVqA0cBj6hqT2A7CSh6ho1X3zgY6AgcCDQQkYtSqyr9MIMaDoHmw0pHRKQOzpg+p6qTU60nID8DzhSRFbjqlRNE5NnUSgpMAVCgqpGSwEs4A5vunAQsV9XvVHU3MBk4PsWaKst6EWkN4H1uCDsDM6jhEGTOrLRDRARXl7dYVf+Raj1BUdWbVTVXVTvg7vX7qpoR3pKqrgNWi8ih3q4TgUUplBSUVcCxIrKf99ycSAY0pkXxOnCJt34J8FrYGSRtCpR9mXhzZqVYVhB+BgwFvhCRed6+P3rT1RiJ4zrgOe/PdxkZMH+aqn4qIi8Bc3G9Qz4njYegisjzQB7QQkQKgNuAu4AXReQK3B/EuaHna0NPDcMwwsGK/IZhGCFhBtUwDCMkzKAahmGEhBlUwzCMkDCDahiGERJmUI2MwYt0lJekvLqIyJwAx10vInclQ5OR/li3KSNtEJFC3+Z+wC6gxNu+UlWfS6KWl4F/q+oLFRxXD1gKHKWqoY+8MTILM6hGWuINKx2mqtNSkHdr4EvgQFXdGeD4x4ElqnpPwsUZaY0V+Y2MQURWiMhJ3vpoEfm3iDwrIj+IyBcicoiI3OwFFl4tIqf4zm0iIuNFZK2IfCsid3pxbGNxMjDXb0xF5CbvvB+8QOIn+o7Px4USNGo4ZlCNTOYXwDNAM9xQyLdxz3Qb4HbgUd+xE3FDJjvjwhSeAgyLc90jgCWRDW/c/bVAH1VtBJwKrPAdvxgXKNqo4ZhBNTKZD1X1bVUtBv4NtATu8qIhvQB08KYbycHNpnCDqm736jrHEj/ifFPgB992CVAX6CIidVR1hap+40v/AWgS6jczMhILjmJkMut96zuAjapa4tsGaIiL31kHWOsCJQHOmfAHBfezBWgU2VDVpSJyAzAa6CoibwO/U9VIiMZGwPfV+ibGPoF5qEZNYDWux0ALVW3qLY1VtWuc4xcAh/h3qOokVe0LtMcFD/+7L/lwYH4CdBsZhhlUY5/Hm+7iHeBeEWksIrVE5CARGRDnlHeBoyKzeorIoSJygojUBXbivN8S3/EDgLcS+BWMDMEMqlFTuBjIxgVz3oKLlN861oGquh54n9KJFuviYmluBNbh5iL6I+zthzqI0snfjBqM9UM1jBh402lPBI7Wcl4SEbkOaKuqo5ImzkhbzKAahmGEhBX5DcMwQsIMqmEYRkiYQTUMwwgJM6iGYRghYQbVMAwjJMygGoZhhIQZVMMwjJD4/2GFTcFNhikTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "N=1000\n",
    "time=np.linspace(-1,10,N)\n",
    "flam = lambdify(t, f, \"numpy\")\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "fig.subplots_adjust(left=0.25, bottom=0.25)\n",
    "legend='$g(t)*h(t) $'\n",
    "[line] = ax.plot(time, flam(time), linewidth=2, color='blue',label=legend)\n",
    "ax.set_xlabel('Time (s)',fontsize=12)\n",
    "ax.set_ylabel('Arbritrary units',fontsize=12) \n",
    "ax.legend()\n",
    "ax.set_title('Step Reponse of Underdamped Notch pass Filter')\n",
    "ax.grid()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a77c1fec",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}