diff --git a/Week 6/10 Hyperbolic PDEs.ipynb b/Week 6/10 Hyperbolic PDEs.ipynb index ea72bb0..7559a98 100644 --- a/Week 6/10 Hyperbolic PDEs.ipynb +++ b/Week 6/10 Hyperbolic PDEs.ipynb @@ -318,36 +318,36 @@ ], "source": [ "# There is an 𝑙 in the boundary conditions we assume should be a 1.\n", - "a = 1\n", - "l = 1\n", - "x = np.linspace(0, l, 1200)\n", - "t = np.linspace(0, 1, 1200)\n", - "f = lambda x: np.sin(2*np.pi*x)\n", - "g = lambda x: 2*np.pi*np.sin(2*np.pi*x)\n", + "a1 = 1\n", + "l1 = 1\n", + "x1 = np.linspace(0, l1, 1200)\n", + "t1 = np.linspace(0, 1, 1200)\n", + "f1 = lambda x: np.sin(2*np.pi*x)\n", + "g1 = lambda x: 2*np.pi*np.sin(2*np.pi*x)\n", "\n", - "w = pdeHyperbolic(a, x, t, f, g)\n", - "u = lambda x, t: np.sin(2*np.pi*x)*(np.cos(2*np.pi*t) + np.sin(2*np.pi*t))\n", + "w1 = pdeHyperbolic(a1, x1, t1, f1, g1)\n", + "u1 = lambda x, t: np.sin(2*np.pi*x)*(np.cos(2*np.pi*t) + np.sin(2*np.pi*t))\n", "\n", "# Plot the whole.\n", "fig = plt.figure(figsize=(24, 8))\n", "\n", - "grid = np.meshgrid(x, t)\n", + "grid1 = np.meshgrid(x1, t1)\n", "\n", "# Plot w\n", "ax0 = fig.add_subplot(1, 3, 1, projection='3d')\n", - "ax0.plot_surface(*grid, w, cmap=\"turbo\")\n", + "ax0.plot_surface(*grid1, w1, cmap=\"turbo\")\n", "ax0.set_title('numerically found solution $w[t_j, x_i]$')\n", "ax0.set_zlabel('$w[t_j, x_i]$')\n", "\n", "# Plot the difference u - w\n", "ax1 = fig.add_subplot(1, 3, 2, projection='3d')\n", - "ax1.plot_surface(*grid, u(*grid) - w, cmap=\"turbo\")\n", + "ax1.plot_surface(*grid1, u1(*grid1) - w1, cmap=\"turbo\")\n", "ax1.set_title('difference $u(x_i,t_j) - w[t_j, x_i]$')\n", "ax1.set_zlabel('$u(x_i,t_j) - w[t_j, x_i]$')\n", "\n", "# Plot u\n", "ax2 = fig.add_subplot(1, 3, 3, projection='3d')\n", - "ax2.plot_surface(*grid, u(*grid), cmap=\"turbo\")\n", + "ax2.plot_surface(*grid1, u1(*grid1), cmap=\"turbo\")\n", "ax2.set_title('exact solution $u(x_i,t_j)$')\n", "ax2.set_zlabel('$u(x_i, t_j)$')\n", "\n", @@ -387,6 +387,7 @@ "outputs": [], "source": [ "def test_pdeHyperbolic():\n", + " # TODO: Docstring\n", " a = 1\n", " l = 1\n", " x = np.linspace(0, l, 1200)\n", @@ -470,25 +471,25 @@ ], "source": [ "# There is an 𝑙 in the boundary conditions we assume should be a 1.\n", - "a = 1\n", - "l = 1\n", - "m = 200\n", - "n = 400\n", - "x = np.linspace(0, l, m)\n", - "t = np.linspace(0, 2, n)\n", - "f = lambda x: -1 + 2.*(x < .5)\n", - "g = lambda x: 0\n", + "a2 = 1\n", + "l2 = 1\n", + "m2 = 200\n", + "n2 = 400\n", + "x2 = np.linspace(0, l2, m2)\n", + "t2 = np.linspace(0, 2, n2)\n", + "f2 = lambda x: -1 + 2.*(x < .5)\n", + "g2 = lambda x: 0\n", "\n", - "w = pdeHyperbolic(a, x, t, f, g)\n", + "w2 = pdeHyperbolic(a2, x2, t2, f2, g2)\n", "\n", "# Plot the whole.\n", "fig = plt.figure(figsize=(16, 16))\n", "\n", - "grid = np.meshgrid(x, t)\n", + "grid2 = np.meshgrid(x2, t2)\n", "\n", "# Plot w\n", "ax = fig.add_subplot(1, 1, 1, projection='3d')\n", - "ax.plot_surface(*grid, w, cmap=\"turbo\")\n", + "ax.plot_surface(*grid2, w2, cmap=\"turbo\")\n", "ax.set_title('numerically found solution $w[t_j, x_i]$')\n", "ax.set_zlabel('$w[t_j, x_i]$')\n", "ax.set_xlabel('$x_i$')\n", @@ -553,7 +554,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": { "deletable": false, "nbgrader": { @@ -568,7 +569,235985 @@ "task": false } }, - "outputs": [], + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
\n", + " \n", + "
\n", + " \n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
\n", + "
\n", + "
\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "# use matplotlib's animation package\n", "import matplotlib.pylab as plt\n", @@ -580,8 +236559,10 @@ "# create a figure for the animation\n", "fig = plt.figure()\n", "plt.grid(True)\n", - "plt.xlim( ... ) # fix x limits\n", - "plt.ylim( ... ) # fix y limits\n", + "plt.xlim(0, 1) # fix x limits\n", + "plt.ylim(-1.5, 1.5) # fix w limits\n", + "plt.xlabel('$x_i$')\n", + "plt.ylabel('$w[t_j, x_i]$')\n", "\n", "# Create an empty plot object and prevent its showing (we will fill it each frame)\n", "myPlot, = plt.plot([0], [0])\n", @@ -589,17 +236570,17 @@ "\n", "# This function is called each frame to generate the animation (f is the frame number)\n", "def animate(f): \n", - " myPlot.set_data( ... ) # update plot\n", + " myPlot.set_data(x1, w1[f]) # update plot\n", "\n", "# Show the animation\n", - "frames = np.arange(1, np.size(t)) # t is the time grid here\n", + "frames = np.arange(1, np.size(t1)) # t is the time grid here\n", "myAnimation = animation.FuncAnimation(fig, animate, frames, interval = 20)\n", "myAnimation" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": { "deletable": false, "nbgrader": { @@ -614,12 +236595,96259 @@ "task": false } }, - "outputs": [], + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
\n", + " \n", + "
\n", + " \n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
\n", + "
\n", + "
\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "# Animate problem 2 here ...\n", + "# use matplotlib's animation package\n", + "import matplotlib.pylab as plt\n", + "import matplotlib\n", + "import matplotlib.animation as animation\n", + "# set the animation style to \"jshtml\" (for the use in Jupyter)\n", + "matplotlib.rcParams['animation.html'] = 'jshtml'\n", "\n", - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "# create a figure for the animation\n", + "fig = plt.figure()\n", + "plt.grid(True)\n", + "plt.xlim(0, 1) # fix x limits\n", + "plt.ylim(-1.5, 1.5) # fix w limits\n", + "plt.xlabel('$x_i$')\n", + "plt.ylabel('$w[t_j, x_i]$')\n", + "\n", + "# Create an empty plot object and prevent its showing (we will fill it each frame)\n", + "myPlot, = plt.plot([0], [0])\n", + "plt.close()\n", + "\n", + "# This function is called each frame to generate the animation (f is the frame number)\n", + "def animate(f): \n", + " myPlot.set_data(x2, w2[f]) # update plot\n", + "\n", + "# Show the animation\n", + "frames = np.arange(1, np.size(t2)) # t is the time grid here\n", + "myAnimation = animation.FuncAnimation(fig, animate, frames, interval = 20)\n", + "myAnimation" ] } ],