improvements to notebook 3

PSPython_03-FFT-aliasing-filtering.ipynb

@@ -119,7 +119,7 @@
     "\n",
     "A numerical grid has a limited resolution.  If you try to represent a rapidly-oscillating function with relatively few grid points, you will observe an effect known as **aliasing**.  This naturally limits the range of frequencies that can be modelled on a given grid.  It can also lead to instabilities in pseudospectral simulations, when generation of high frequencies leads to buildup of lower-frequency energy due to aliasing.\n",
     "\n",
-    "The code below plots a sine wave of a given wavenumber $p$, along with its representation on a grid with $m$ points.  Try changing $p$ and notice how for $  |p| \\ge m$ the function looks like a lower-frequency mode."
+    "The code below plots a sine wave of a given wavenumber $\\xi$, along with its representation on a grid with $m$ points.  Try changing $\\xi$ and notice how for $  |\\xi| \\ge m$ the function looks like a lower-frequency mode."
    ]
   },
   {
@@ -334,10 +334,11 @@
     "    line2.set_data(np.sort(xi),power_spectrum[np.argsort(xi)])\n",
     "    axes.set_title('t= %.2e' % tt[i])\n",
     "    axes.set_xlim((-np.pi,np.pi))\n",
-    "    axes.set_ylim((-100,3000))\n",
+    "    axes.set_ylim((-200,3000))\n",
     "    \n",
     "anim = matplotlib.animation.FuncAnimation(fig, plot_frame,\n",
-    "                                   frames=len(frames), interval=100)\n",
+    "                                          frames=len(frames), interval=100,\n",
+    "                                          repeat=False)\n",
     "\n",
     "HTML(anim.to_jshtml())"
    ]
@@ -346,7 +347,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In the output, we're plotting the solution (top plot) and its power spectrum ($|\\hat{u}|^2$) (bottom plot).  There are a lot of interesting things to say about the solution, but for now let's focus on the Fourier transform.  Notice how the wavenumbers present in the solution remain in the lower half of those representable on the grid (this region is delimited by the dashed red lines).  Because of this, no aliasing occurs.\n",
+    "In the output, we're plotting the solution (top plot) and its power spectrum ($|\\hat{u}|^2$) (bottom plot).  There are a lot of interesting things to say about the solution, but for now let's focus on the Fourier transform.  Notice how the amount of energy in the high wavenumbers present (those outside the dashed red lines) remains relatively small.  Because of this, no aliasing occurs.\n",
     "\n",
     "Now change the code above to use only $m=128$ grid points.  What happens?"
    ]
@@ -359,7 +360,7 @@
     "\n",
     "Here we will give a somewhat simplified explanation of the blow-up just observed.  First, this blowup has nothing to do with the absoute stability condition -- when we change $m$, the time step is automatically changed in a way that will ensure absolute stability.  If you're not convinced, try taking the time step even smaller; you will still observe the blowup.\n",
     "\n",
-    "By taking $m=128$, we cut by half the wavenumbers that can be represented on the grid.  As you can see from the plots, this means that some of the wavenumbers present in the initial data are in the upper half of the representable range (i.e., outside the dashed red lines).  That means that the highest wavenumbers generated by the quadratic term will be aliased -- and they will be aliased back into that upper-half range.  This leads to a gradual accumulation of high-wavenumber modes, easily visible in both plots.  Eventually the high-wavenumber modes dominate the numerical solution and lead to blowup.\n",
+    "By taking $m=128$, we cut by half the wavenumbers that can be represented on the grid.  As you can see from the plots, this means that after the solution begins to steepen, a significant amount of the energy present in the solution is in the upper half of the representable range of wavenumbers (i.e., outside the dashed red lines).  That means that the highest wavenumbers generated by the quadratic term will be aliased -- and they will be aliased back into that upper-half range.  This leads to a gradual and incorrect accumulation of high-wavenumber modes, easily visible in both plots.  Eventually the high-wavenumber modes dominate the numerical solution and lead to blowup.\n",
     "\n",
     "For a detailed discussion of aliasing instabilities, see Chapter 11 of John Boyd's \"Chebyshev and Fourier Spectral Methods\"."
    ]
@@ -391,7 +392,7 @@
     "def rhs(u, xi, filtr, equation='KdV'):\n",
     "    uhat = np.fft.fft(u)\n",
     "    if equation == 'Burgers': \n",
-    "        return -u*np.real(np.fft.ifft(1j*xi*uhat)) \\\n",
+    "        return -u*np.real(np.fft.ifft(1j*xi*uhat*filtr)) \\\n",
     "                + np.real(np.fft.ifft(-xi**2*uhat))\n",
     "    elif equation == 'KdV':\n",
     "        return -u*np.real(np.fft.ifft(1j*xi*uhat*filtr)) \\\n",
@@ -417,11 +418,11 @@
     "A = 25; B = 16;\n",
     "u = 3*A**2/np.cosh(0.5*(A*(x+2.)))**2 + 3*B**2/np.cosh(0.5*(B*(x+1)))**2\n",
     "#u = 1500*np.exp(-10*(x+2)**2)\n",
-    "tmax = 0.006\n",
+    "tmax = 0.012\n",
     "\n",
     "uhat2 = np.abs(np.fft.fft(u))\n",
     "\n",
-    "num_plots = 50\n",
+    "num_plots = 200\n",
     "nplt = np.floor((tmax/num_plots)/dt)\n",
     "nmax = int(round(tmax/dt))\n",
     "\n",
@@ -473,11 +474,34 @@
    "source": [
     "Notice how the solution remains stable, but small wiggles appear throughout the domain.  These are a hint that something is not sufficiently resolved."
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Exercises\n",
+    "The examples we have looked at behave in a way that agrees with the ideas that have been presented here.  However, the stability of nonlinear pseudospectral PDE discretizations is complicated, and it's easy to find examples that aren't fully explained by the discussion above.  You can explore this by experiments like the following.\n",
+    "\n",
+    "1.  Try decreasing $m$ even further in the example above.  You may find that the solution blows up for some values.  Can you prevent this by adjusting $\\Delta t$?  Can you prevent it by filtering more wavenumbers?\n",
+    "\n",
+    "2.  Try solving the inviscid Burgers' equation with a 2/3 filter.  For this problem, there is no value of $m$ large enough to resolve the gradients that appear, since the gradient blows up in finite time.  After the blowup time, traditionally one employs the notion of \"vanishing viscosity weak solutions\" in which a discontinuous solution is acceptable if it is given by the limit of solutions of\n",
+    "$$\n",
+    "u_t + uu_x = \\epsilon u_{xx}\n",
+    "$$\n",
+    "as $\\epsilon \\to 0$.  By using filtering, can you obtain such solutions?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },