Lec09-accurate-deriv.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b37b10e3-cb22-4fc5-a428-e70da4e313ee",
   "metadata": {},
   "source": [
    "### More accurate derivatives\n",
    "\n",
    "( see https://en.wikipedia.org/wiki/Five-point_stencil )\n",
    "\n",
    "Taylor series \n",
    "\n",
    "$\n",
    "f(x+h)=f(x) + f'(x)h + \\frac{f''(x)h^{2}}{2!} + \\frac{f'''(x)h^3}{3!}+ \\ldots\n",
    "$\n",
    "\n",
    "lets discretize $x$ with a step size $h$ i.e.\n",
    "\n",
    "$\n",
    "x_j=x_0 + h(j-1)\n",
    "$\n",
    "\n",
    "we can the write\n",
    "\n",
    "$ f(x_{j-2}) \\approx f(x_j) - f'(x_j)2h + \\frac{1}{2}f''(x_j)4h^{2} - \\frac{1}{6}f'''(x_j)8h^3 + \\frac{1}{24}f''''(x_j)16h^4 - \\frac{1}{120}f'''''(x_j)32h^5 + \\ldots $\n",
    "\n",
    "$ f(x_{j-1}) \\approx f(x_j) - f'(x_j)h  + \\frac{1}{2}f''(x_j)h^{2} - \\frac{1}{6}f'''(x_j)h^3 + \\frac{1}{24}f''''(x_j)h^4 - \\frac{1}{120}f'''''(x_j)h^5 + \\ldots $\n",
    "\n",
    "$ f(x_{j-1}) \\approx f(x_j) + f'(x_j)h  + \\frac{1}{2}f''(x_j)h^{2} + \\frac{1}{6}f'''(x_j)h^3 + \\frac{1}{24}f''''(x_j)h^4 + \\frac{1}{120}f'''''(x_j)h^5 + \\ldots $\n",
    "\n",
    "$ f(x_{j+2}) \\approx f(x_j) + f'(x_j)2h + \\frac{1}{2}f''(x_j)4h^{2} + \\frac{1}{6}f'''(x_j)8h^3 + \\frac{1}{24}f''''(x_j)16h^4 + \\frac{1}{120}f'''''(x_j)32h^5 + \\ldots $\n",
    "\n",
    "and then rearrange for $f'(x_j)$\n",
    "\n",
    "$\n",
    "f'(x_j) \\approx \\frac{f(x_{j-2})-8f(x_{j-1})+8f(x_{j+1})+f(x_{j+2})}{12h}\n",
    " + \\ldots \n",
    "$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eeffe989-6676-4bbd-b35f-08a83ab89c38",
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Lets test these formulas\n",
    "#\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "h=0.005;\n",
    "L=2*np.pi;N=np.round(L/h)+1;h=L/(N-1)\n",
    "x=np.linspace(0,L,int(N) )\n",
    "\n",
    "def ff(x,L):\n",
    "    return np.cos(x/L*np.pi*2)\n",
    "\n",
    "def fop(xvec, L, f):\n",
    "    ans=np.zeros(np.size(xvec))\n",
    "    i=0\n",
    "    for x in xvec:\n",
    "        ans[i]=f(x, L)\n",
    "        i+=1\n",
    "    return ans\n",
    "\n",
    "fx=fop(x, L, ff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b795fcb-be1c-4bc9-b4fd-d5482d7b7adb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def deriv1(fx,h):\n",
    "    deriv=fx[1:]-fx[0:-1]\n",
    "    deriv=deriv/h\n",
    "    return deriv\n",
    "\n",
    "def deriv2(fx,h):\n",
    "    deriv=fx[2:]-fx[0:-2]\n",
    "    deriv=deriv/(2*h)\n",
    "    return deriv\n",
    "\n",
    "def deriv5(fx,h):\n",
    "    fxm2=fx[0:-4]\n",
    "    fxm1=fx[1:-3]\n",
    "    fx_0=fx[2:-2]\n",
    "    fxp1=fx[3:-1]\n",
    "    fxp2=fx[4:]\n",
    "    deriv=fxm2-8*fxm1+8*fxp1-fxp2\n",
    "    deriv=deriv/(12*h)\n",
    "    return deriv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0030a694-c80c-4f23-93fb-fbc88a39f398",
   "metadata": {},
   "outputs": [],
   "source": [
    "d1=deriv1(fx,h)\n",
    "d2=deriv2(fx,h)\n",
    "d5=deriv5(fx,h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd85f28d-232b-4779-9518-c45f848f3b8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(x,fx)\n",
    "plt.plot(x[1:],d1)\n",
    "plt.plot(x[1:],-np.sin(x[1:]/L*np.pi*2))\n",
    "plt.show()\n",
    "plt.plot(x[1:],d1+np.sin(x[1:]/L*np.pi*2))\n",
    "plt.show()\n",
    "e1=np.linalg.norm( d1+np.sin(x[1:]/L*np.pi*2) )\n",
    "print(e1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cfaaabf-fe50-4470-9865-2242a137419b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(x[1:-1],d2)\n",
    "plt.show()\n",
    "plt.plot(x[1:-1],d2+np.sin(x[1:-1]/L*np.pi*2))\n",
    "plt.show()\n",
    "e2=np.linalg.norm( d2+np.sin(x[1:-1]/L*np.pi*2) )\n",
    "print(e2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c501492-8a94-4dbb-84f9-a3e6ea44dc0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(x[2:-2],d5)\n",
    "plt.show()\n",
    "plt.plot(x[2:-2],d5+np.sin(x[2:-2]/L*np.pi*2))\n",
    "plt.show()\n",
    "e5=np.linalg.norm( d5+np.sin(x[2:-2]/L*np.pi*2) )\n",
    "print(e5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f8473f1-e1b9-4298-9d54-36b764c69d8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now lets look at convergence\n",
    "def get_x(h):\n",
    "    L=2*np.pi;N=np.round(L/h)+1;h=L/(N-1)\n",
    "    x=np.linspace(0,L,int(N) )\n",
    "    return x, h\n",
    "\n",
    "h=0.05\n",
    "nh=5\n",
    "e1vec=np.zeros(nh)\n",
    "e2vec=np.zeros(nh)\n",
    "e5vec=np.zeros(nh)\n",
    "hvec=np.zeros(nh)\n",
    "\n",
    "for i in [*range(0,nh)]:\n",
    "    x,h=get_x(h)\n",
    "    fx=fop(x, L, ff)\n",
    "    d1=deriv1(fx,h)\n",
    "    e1vec[i]=np.linalg.norm( d1+np.sin(x[1:]/L*np.pi*2) )\n",
    "    d2=deriv2(fx,h)\n",
    "    e2vec[i]=np.linalg.norm( d2+np.sin(x[1:-1]/L*np.pi*2) )\n",
    "    d5=deriv5(fx,h)\n",
    "    e5vec[i]=np.linalg.norm( d5+np.sin(x[2:-2]/L*np.pi*2) )\n",
    "    hvec[i]=h\n",
    "    h=h/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbb227ec-9ad2-4056-99c2-225ac5de6207",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(x[1:], d1+np.sin(x[1:]/L*np.pi*2) )\n",
    "plt.show()\n",
    "plt.plot(x[1:-1],d2+np.sin(x[1:-1]/L*np.pi*2) )\n",
    "plt.show()\n",
    "plt.plot(x[2:-2],d5 + np.sin(x[2:-2]/L*np.pi*2) )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4444e4e-92de-4686-9218-64b62828cab1",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.semilogy(e1vec/e1vec[0],'.')\n",
    "plt.show()\n",
    "plt.semilogy(e2vec/e2vec[0],'.')\n",
    "plt.show()\n",
    "plt.semilogy(e5vec/e5vec[0],'.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42101544-350c-493d-9396-43eb24c8064f",
   "metadata": {},
   "outputs": [],
   "source": [
    "hvec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c3220dc-dd6f-493a-8b7a-014a4c5a8f28",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}