stationary_phase_approx.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ec5b9e42-79cc-4f89-a60b-aa00f07ecc4f",
   "metadata": {},
   "source": [
    "# **Illustration of the stationary phase approximation to an integral.**\n",
    "<i class=\"fa fa-home fa-2x\"></i><a href=\"./index.ipynb\" style=\"font-size: 20px\"> Go back to index</a>\n",
    "\n",
    "**Source code:** https://github.com/Taylor-96/SP_tutorial/tree/osscar\n",
    "\n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18b95e21-a64f-45f5-a94e-c40055686318",
   "metadata": {},
   "source": [
    "# **Goals**\n",
    "* Appreciate the difficulty in numerically computing the integral of a highly oscillatory function.\n",
    "* Observe how the cancellation of the contributions to an integral from the regions in which the integrand has rapid oscillations leads to the idea of the stationary phase approximation. \n",
    "* Note how the stationary phase approximation becomes better and better as the strength of the integrand's oscillations increases and as the size of the integration range enclosing the stationary point is extended."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6132b69b",
   "metadata": {},
   "source": [
    "\n",
    "## **Background Theory**\n",
    "<p style=\"text-align: justify;font-size:15px;width: 90%\">\n",
    "\n",
    "For a brief review of the background theory underlying the stationary phase approximation along with a sketch of its derivation see <a href=\"./stationary_phase_approx_intro.ipynb\" style=\"font-size: 20px\"> stationary phase approximation background theory.</a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc44a9cf-2888-4693-b2ed-bb07d406be5c",
   "metadata": {},
   "source": [
    "# **Tasks and exercises** \n",
    "<ol style=\"text-align: justify;font-size:15px\">\n",
    "    <li> Investigate the difficulty associated with trying to approximate numerically the integral $\\int  \\limits_{-L}^{L}dx e^{i\\alpha(x^2-a)}$. To do so, vary the slider which controls the number of grid points in the numerical quadrature scheme (i.e how fine the mesh is), and  the slider which \n",
    "    <details style=\"color: red\">\n",
    "    <summary>Solution</summary>\n",
    "        It should be clear that as we increase the strength of the oscillations in the integrand (controlled by the parameter $\\alpha$) we are forced to introduce more and more grid points to produce an accurate approximation. This is a consequence of the fact that for highly oscillatory functions, in order to resolve them well, it is necessary to have a fine mesh/grid  Therefore, for a highly oscillatory integrand, the computational cost associated with brute force numerical quadrature quickly becomes intractable. \n",
    "    </details>\n",
    "    </li>\n",
    "        <li> Investigate the convergence properties of the stationary phase approximation. To do so, vary the magnitude of the oscillations of the integrand (move the \"alpha\" slider in Figure 2) along with the limits of the integral (the \"integral_lims\" slider). \n",
    "    <details style=\"color: red\">\n",
    "    <summary>Solution</summary>\n",
    "       One may observe that as the amplitude of the fluctuations ($\\alpha$) is increased, the stationary phase approximation becomes an increasinly accurate estimate of the original integral. This is due to the fact that if we increase the rapidity of the fluctuations in the integrand, as we move away from the point of stationary phase, the contributions to the integral shall become increasingly negligible. Thus, an approximation of the integral by its behaviour in the vicinity of the point where its integrand's phase is stationary becomes a progressively more faithful representation of the integral as a whole. In addition, one can note that by increasing the integration limits, the value obtained through numerical quadrature and via the stationary phase approximation become better and better aligned. Note: it may be necessary to increase the number of grid points used in the numerical quadrature scheme as the integration limits are increased in order to achieve sufficient convergence of this estimate. This latter fact again illustrates the difficulty in computing such oscillatory integrals numerically; the breakdown at large integration limits is due to the increased strength of oscillations in the integrand far\n",
    "    </details>\n",
    "    </li>\n",
    "      <li> Derive the stationary phase approximations to the integrand $e^{i\\alpha( e ^{x^2}-1)}$ and use it to approximate the integral $\\int \\limits_{-L}^{L} dx e^{i\\alpha( e ^{x^2}-1)}$, as illustrated in the notebook.\n",
    "    <details style=\"color: red\">\n",
    "    <summary>Solution</summary>\n",
    "       Following the procedure outlined in the accompanying theory notebook we Taylor expand (up to second order) around the stationary point the argument of the exponential $f(x) = e^{F(x)}$, with $F(x) =  e ^{x^2}-1)$. This yields $f(x) \\approx e^{-i \\alpha x^2}$ (noting that the stationary point of the integrand is located at $x=0$). Consequently, we may then write $\\int \\limits_{-L}^{L} dx e^{i\\alpha( e ^{x^2}-1)} \\approx \\int \\limits_{-\\infty}^{\\infty} dx e^{i \\alpha x^2}$. This is just a generalized Gaussian integral that we may evaluate through the use of contour integration as $\\int \\limits_{-\\infty}^{\\infty} dx e^{i \\alpha x^2} = \\sqrt{\\frac{i\\pi}{\\alpha}} = \\sqrt{\\frac{\\pi}{\\alpha}} e ^{i \\pi/4}$.\n",
    "    </details>\n",
    "    </li>\n",
    "</ol>\n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "549bf9c5",
   "metadata": {},
   "source": [
    "# ** Difficulty of converging a highly oscillatory integral with naive numerical integration.**\n",
    "<p style=\"text-align: justify;font-size:15px;width: 90%\">\n",
    "\n",
    "Let us first look at how easy (or rather how difficult) it is to evaluate, using standard brute force numerical integration, integrals with such highly oscillatory integrands.\n",
    "    </p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4590215a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f6da1a7295ed421ba4690184eed02d37",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3bbfef5040ec46c19a42699336134592",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(interactive(children=(FloatSlider(value=2.0, continuous_update=False, description='integral_lim…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib widget\n",
    "\n",
    "from ipywidgets import interactive,Dropdown,Layout,VBox,HBox,FloatSlider\n",
    "import matplotlib.pyplot as plt1\n",
    "import numpy as np\n",
    "from IPython.display import display, clear_output\n",
    "\n",
    "x=np.linspace(-2,2,2000)\n",
    "\n",
    "def f1_x(x,alpha,a):\n",
    "    return np.exp(1.j*alpha*(x**2-a))\n",
    "\n",
    "def int_f1(x,alpha,a,Ngrid):\n",
    "\n",
    "    dx=(x[-1]-x[0])/Ngrid\n",
    "    I=0.+0.j\n",
    "    for n in range(Ngrid):\n",
    "        I += f1_x(x[0]+n*dx,alpha,a)*dx\n",
    "\n",
    "    return I\n",
    "    \n",
    "def f1_int_x_interactive(integral_lims=2.0,alpha=10.0,a=0.0, Ngrid=1000,Real_or_Imaginary_part='real'):\n",
    "    \n",
    "    \n",
    "    axs1.clear()\n",
    "    x = np.linspace(-2, 2, num=2000)\n",
    "\n",
    "    f = f1_x(x,alpha,a)\n",
    "    if(Real_or_Imaginary_part=='real'):\n",
    "        re_f=np.real(f)\n",
    "        axs1.plot(x, re_f,label=\"Re(f(x))\" )\n",
    "        axs1.fill_between(x, re_f, 0,where=re_f>0.,facecolor='yellow', alpha=0.5)\n",
    "        axs1.fill_between(x, re_f, 0,where=re_f<0.,facecolor='green', alpha=0.5)\n",
    "    else:\n",
    "        im_f = np.imag(f)\n",
    "        axs1.plot(x, np.imag(f),label=\"Im(f(x))\",color='red')\n",
    "        axs1.fill_between(x, im_f, 0,where=im_f>0.,facecolor='yellow', alpha=0.5)\n",
    "        axs1.fill_between(x, im_f, 0,where=im_f<0.,facecolor='green', alpha=0.5)\n",
    "\n",
    "    axs2.clear()\n",
    "    axs2.set_xlim(0,Ngrid)\n",
    "\n",
    "    x = np.linspace(-1.*integral_lims,integral_lims,num=Ngrid)\n",
    "    If_lst1=[]\n",
    "    grid_list=np.arange(1,Ngrid)\n",
    "    for n_points in grid_list:\n",
    "        If_lst1.append(int_f1(x, alpha,a, n_points))  \n",
    "    If_lst1=np.array(If_lst1)\n",
    "\n",
    "    if(Real_or_Imaginary_part=='real'):\n",
    "        axs2.plot(grid_list, np.real(If_lst1),label=r\"Re($\\int _{L}^{L} f(x)$)\" )\n",
    "\n",
    "    else:\n",
    "\n",
    "        axs2.plot(grid_list, np.imag(If_lst1),label=r\"Im($\\int_{L}^{L} f(x)$)\",color='red' )\n",
    "\n",
    "    axs1.set_xlabel('x')\n",
    "    axs2.set_xlabel('Ngrid')\n",
    "\n",
    "    axs1.set_title(r\"Plot of integrand $\\exp(i \\alpha (x^2-a))$.\")        \n",
    "    axs2.set_title(\"Plot of numerical integration vs number of grid points.\")\n",
    "    \n",
    "    plt1.legend()\n",
    "    #plt1.show()\n",
    "    display(plt1)\n",
    "        \n",
    "    \n",
    "fig1=plt1.figure(1)    \n",
    "fig1.set_size_inches(8,6.5)\n",
    "plt1.subplots_adjust(hspace = .5)   \n",
    "axs1=fig1.add_subplot(211)\n",
    "line1, = axs1.plot([],[])\n",
    "axs2=fig1.add_subplot(212)\n",
    "\n",
    "line2, = axs2.plot([],[])\n",
    "\n",
    "interactive_plot1 = interactive(f1_int_x_interactive,\\\n",
    "                                integral_lims=FloatSlider(value=2.0,min=1.0,max=10.0,step=0.1,continuous_update=False),\\\n",
    "                                Ngrid=1000, \\\n",
    "                                alpha=FloatSlider(value=10.0, min=1.0,max=1000.0,step=10.,continuous_update=False),\\\n",
    "                                a=FloatSlider(value=0.0, min=-1.0,max=1.0,step=.1,continuous_update=False),\\\n",
    "                                Real_or_Imaginary_part=['real', 'imaginary'])\n",
    "\n",
    "# interactive_plot\n",
    "box1 = VBox( [interactive_plot1], layout=Layout(width='800px'))\n",
    "box1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01d9a0ce",
   "metadata": {},
   "source": [
    "# ** Convergence of the stationary phase approximation to the exact result.**\n",
    "<p style=\"text-align: justify;font-size:15px;width: 90%\">\n",
    "We now look at how good a job the stationary phase approximation does of estimating the true value of the integral depending on how large $\\alpha$ (the strength of the oscillations) is. Furthermore, we see how much error is introduced by extending the integration limits to $\\pm \\infty$ from some finite interval $(x_i,x_f)$.\n",
    "    </p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "364a5366",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "417917d2a4914b5bb8b5404e23814e2f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(interactive(children=(FloatSlider(value=2.0, description='integral_lims', max=10.0, min=1.0), F…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\n",
    "\n",
    "from ipywidgets import interactive,Dropdown,Layout,VBox,HBox\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "def f_x(x,alpha):\n",
    "    return np.exp(1.j*alpha*(np.exp(1.*x**2)-1))\n",
    "    \n",
    "def int_f(x,alpha,Ngrid):\n",
    "    \n",
    "    If_lst=[]\n",
    "    dx=(x[-1]-x[0])/Ngrid\n",
    "    I=0.+0.j\n",
    "    for n in range(int(Ngrid/2)):\n",
    "        I += 2.*f_x(n*dx,alpha)*dx\n",
    "        If_lst.append(I)\n",
    "        \n",
    "    return If_lst\n",
    "    \n",
    "def f_int_x_interactive(integral_lims=2.0,alpha=10.0, Ngrid=1000,Real_or_Imaginary_part='real'):\n",
    "    \n",
    "    ax3.clear()\n",
    "    x = np.linspace(-2, 2, num=2000)\n",
    "\n",
    "    f = f_x(x,alpha)\n",
    "    f_sp = np.exp(1.j*alpha*(x**2))\n",
    "    if(Real_or_Imaginary_part=='real'):\n",
    "\n",
    "        ax3.plot(x, np.real(f),label=\"Re\" )\n",
    "        ax3.plot(x, np.real(f_sp),label=\"Re(sp approx)\" )\n",
    "\n",
    "    else:\n",
    "\n",
    "        ax3.plot(x, np.imag(f),label=\"Im\", color=\"red\")\n",
    "        ax3.plot(x, np.imag(f_sp),label=\"Im(SP approx)\", color=\"green\")\n",
    "\n",
    "    ax4.clear()\n",
    "    ax4.set_xlim(0,integral_lims)\n",
    "    x = np.linspace(0, integral_lims, num=int(Ngrid/2))    \n",
    "\n",
    "    If_lst=np.array(int_f(x, alpha, Ngrid)) \n",
    "\n",
    "    If_sp = np.sqrt(np.pi/alpha)*np.exp(1.j*np.pi/4.)*np.ones(len(x))\n",
    "    if(Real_or_Imaginary_part=='real'):\n",
    "        ax4.plot(x, np.real(If_lst),label=r\"Re($\\int _{L}^{L} f(x)$)\" )\n",
    "        ax4.plot(x, np.real(If_sp),label=\"Re(SP approx)\" )\n",
    "\n",
    "    else:\n",
    "\n",
    "        ax4.plot(x, np.imag(If_lst),label=r\"Im($\\int_{L}^{L} f(x)$)\", color=\"red\" )\n",
    "        ax4.plot(x, np.imag(If_sp),label=\"Im(SP approx)\", color=\"green\")\n",
    "\n",
    "    ax3.set_xlabel('x')\n",
    "    ax4.set_xlabel('L')\n",
    "\n",
    "    ax3.set_title(r\"Plot of integrand $\\exp(i \\alpha (e^{x^2}-1))$.\")        \n",
    "    ax4.set_title(r\"Plot of numerical estimate and stationary phase approx of $\\int_{-L}^{L}dx e^{i\\alpha( e ^{x^2}-1)}$.\")\n",
    "    ax3.legend(loc='upper right')\n",
    "    ax4.legend()\n",
    "\n",
    "    plt.show()\n",
    "          \n",
    "fig=plt.figure(2)    \n",
    "fig.set_size_inches(8,6.5)\n",
    "plt.subplots_adjust(hspace = .5)   \n",
    "ax3=fig.add_subplot(211)\n",
    "line3, = ax3.plot([],[])\n",
    "ax4=fig.add_subplot(212)\n",
    "line4, = ax4.plot([],[])\n",
    "interactive_plot = interactive(f_int_x_interactive,integral_lims=(1.0,10.0),Ngrid=1000, alpha=(0.0,100,1.), Real_or_Imaginary_part=['real', 'imaginary'])\n",
    "\n",
    "# interactive_plot\n",
    "box2 = VBox( [interactive_plot], layout=Layout(width='800px'))\n",
    "box2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73aa3481-da01-4bce-b620-3d7ad0e8ddc6",
   "metadata": {},
   "source": [
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />\n",
    "\n",
    "# Legend (How to use the interactive visualization)\n",
    "\n",
    "## Visualization\n",
    "\n",
    "The visualization consists of 4 subplots (from top to bottom):\n",
    "\n",
    "1.  A plot of the integrand $F(x)=e^{i\\alpha(x^2-a)}$, with positive and negative contributions indicated with different colours of shading.\n",
    "\n",
    "2. Plot of a numerical estimate of $\\int  \\limits_{-L}^{L}dx e^{i\\alpha(x^2-a)}$ using a very basic numerical quadrature approach wherein the integral is approximated by evaluating the sum $\\sum_n F(n dx) \\cdot dx$ where $dx=\\frac{x_f-x_i}{N_{\\text{grid}}}$ with $x_f=L$,$x_i=-L$ and $N_{\\text{grid}}$ is the number of gridpoints used in the discretization of the integration range. The value of the the integral is shown on the y-axis whilst the number of grid points used is plotted on the x-axis.\n",
    "\n",
    "3. Plot of the integrand $F(x)=e^{i\\alpha( e ^{x^2}-1)}$ and its stationary phase approximation (denoted \"sp approx\"). \n",
    "\n",
    "4. Plot of the numerical estimate of the integral $\\int \\limits_{-L}^{L}dx e^{i\\alpha( e ^{x^2}-1)}$ (control its convergence through the choice of the number of grid points) alongside the stationary phase approximation to the integral for a given value of $\\alpha$. On the y-axis is plotted the values of these estimates whilst on the x-axis the value of the integration limits is shown.\n",
    "\n",
    "## Interactivity controls\n",
    "\n",
    "* Common to all plots is the ability to choose:\n",
    "    * The limits of the integral (which are symmetric from $-L$ to $L$) via the \"integral_lims\" slider.\n",
    "    * The strength of the oscillations of the integrand, $\\alpha$, via the \"alpha\" slider.\n",
    "    * The number of grid points used in the numerical estimate of the integrals via rectangular quadrature.\n",
    "    * Observation of either the real or imaginary part of the integrands and integrals via the \"Real_or_Imaginary\" dropdown selector.\n",
    "* In the first plot (of the integrand $F(x)=e^{i\\alpha(x^2-a)}$) one may also select the value of $\\alpha$ to investigate the consequences of having more than one stationary point of the integrand.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e076feb5",
   "metadata": {},
   "source": [
    "<!-- # ** Convergence of the stationary phase approximation to the exact, numerically integrated result with increasing $\\alpha$.  **\n",
    "\n",
    "<p style=\"text-align: justify;font-size:15px;width: 90%\">\n",
    "    We now look at how the error introduced by the stationary phase approximation becomes less and less significant as we increase the value $\\alpha$ which controls the strength of the oscillations in the integrand.\n",
    "</p> -->"
   ]
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}