Files
cds-numerical-methods/Week 6/10 Hyperbolic PDEs.ipynb
2022-03-11 09:12:38 +01:00

464 lines
13 KiB
Plaintext

{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "4ec40081b048ce2f34f3f4fedbb0be10",
"grade": false,
"grade_id": "cell-98f724ece1aacb67",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"# CDS: Numerical Methods Assignments\n",
"\n",
"- See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.\n",
"\n",
"- Solutions must be submitted via the Jupyter Hub.\n",
"\n",
"- Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\".\n",
"\n",
"## Submission\n",
"\n",
"1. Name all team members in the the cell below\n",
"2. make sure everything runs as expected\n",
"3. **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart)\n",
"4. **run all cells** (in the menubar, select Cell$\\rightarrow$Run All)\n",
"5. Check all outputs (Out[\\*]) for errors and **resolve them if necessary**\n",
"6. submit your solutions **in time (before the deadline)**"
]
},
{
"cell_type": "raw",
"metadata": {},
"source": [
"team_members = \"\""
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "62fbcf7aaa9e165abd27319c5bd35555",
"grade": false,
"grade_id": "cell-e9dc85dd9ad77a5e",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"## Dynamic Behavior from Hyperbolic PDEs\n",
"\n",
"In the following you will derive and implement finite-difference methods to study generalized hyperbolic partial differential equations."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "910692c2f4b93f251a3a128caf2b0c37",
"grade": true,
"grade_id": "cell-acb87410eaef9fe1",
"locked": false,
"points": 0,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# Import packages here ...\n",
"\n",
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "7238a56701aaf868838ffdb706f5eb8b",
"grade": false,
"grade_id": "cell-a27c54a734bf2232",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"### Task 1: finite-difference hyperbolic PDE solver\n",
"\n",
"Our aim is to implement a Python function to find the solution of the following PDE:\n",
"\n",
"\\begin{align*}\n",
" \\frac{\\partial^2}{\\partial t^2} u(x,t)\n",
" - \\alpha^2 \\frac{\\partial^2}{\\partial x^2} u(x,t) = 0,\n",
" \\qquad\n",
" 0 \\leq x \\leq l,\n",
" \\qquad\n",
" 0 \\leq t\n",
"\\end{align*}\n",
"\n",
"with the boundary conditions\n",
"\n",
"\\begin{align*}\n",
" u(0,t) = u(l,t) &= 0, \n",
" &&\\text{for } t > 0 \\\\\n",
" u(x,0) &= f(x) \\\\\n",
" \\frac{\\partial}{\\partial t} u(x,0) &= g(x)\n",
" &&\\text{for } 0 \\leq x \\leq l.\n",
"\\end{align*}\n",
" \n",
"By approximating the partial derivatives with finite differences we can recast the problem into the following form\n",
"\t\t\n",
"\\begin{align*}\n",
" \\vec{w}_{j+1} = \\mathbf{A} \\vec{w}_{j}\n",
" - \\vec{w}_{j-1}.\n",
"\\end{align*}\n",
"\n",
"Here $\\vec{w}_j$ is a vector of length $m$ in the discretized spatial coordinate $x_i$ at time step $t_j$. The spatial coordinates are defined as $x_i = i h$, where $i=0,1,\\dots,m-1$ and $h = l/(m-1)$. The time steps are defined as $t_j = j k$, where $j = 0, 1, \\dots, n-1$.\n",
"\n",
"The tri-diagonal matrix $\\mathbf{A}$ has size $(m-2)\\times(m-2)$ and is defined by\n",
"\n",
"\\begin{align*}\n",
" \\mathbf{A} =\n",
" \\left( \\begin{array}{cccc}\n",
" 2(1-\\lambda^2) & \\lambda^2 & & 0\\\\\n",
" \\lambda^2 & \\ddots & \\ddots & \\\\\n",
" & \\ddots & \\ddots & \\lambda^2 \\\\\n",
" 0 & & \\lambda^2 & 2(1-\\lambda^2) \n",
" \\end{array} \\right),\n",
"\\end{align*}\n",
" \n",
"where $\\lambda = \\alpha k / h$. This $(m-2)\\times(m-2)$ structure accounts for the first set of boundary conditions. Note that the product $\\mathbf{A} \\vec{w}_{j}$ is thus only performed over the $m-2$ subset, i.e. $i=1,2,\\dots,m-2$. The other boundary conditions are accounted for by initializing the first two time steps with\n",
"\n",
"\\begin{align*}\n",
" w_{i,j=0} &= f(x_i) \\\\\n",
" w_{i,j=1} &= (1-\\lambda^2) f(x_i)\n",
" + \\frac{\\lambda^2}{2} f(x_{i+1})\n",
" + \\frac{\\lambda^2}{2} f(x_{i-1})\n",
" + kg(x_i).\n",
"\\end{align*}\n",
" \n",
"Implement a Python function of the form $\\text{pdeHyperbolic(a, x, t, f, g)}$, where $\\text{a}$ represents the PDE parameter $\\alpha$, $\\text{x}$ and $\\text{t}$ are the discretized spatial and time grids, and $\\text{f}$ and $\\text{g}$ are the functions defining the boundary conditions. This function should return a two-dimensional array $\\text{w[:,:]}$, which stores the spatial vector $\\vec{w}_j$ at each time coordinate $t_j$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "95d463713a216d57cbf814a0f40a482d",
"grade": true,
"grade_id": "cell-1f9655333f2dc3cb",
"locked": false,
"points": 3,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def pdeHyperbolic(a, x, t, f, g):\n",
" # YOUR CODE HERE\n",
" raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "d5ee65aba16eef98296a6a66af3c0888",
"grade": false,
"grade_id": "cell-047f2dfab5489a88",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"### Task 2 \n",
"\n",
"Use your implementation to solve the following problems. Compare in the first problem your numerical solution to the analytic one and use it to debug your code.\n",
"\n",
"#### Problem 1:\n",
"\\begin{align*}\n",
" \\frac{\\partial^2}{\\partial t^2} u(x,t)\n",
" - \\frac{\\partial^2}{\\partial x^2} u(x,t) &= 0,\n",
" &&\\text{for }0 \\leq x \\leq 1 \\text{ and } 0 \\leq t \\leq 1\\\\\n",
" u(0,t) = u(l,t) &= 0, \n",
" &&\\text{for } t > 0 \\\\\n",
" u(x,0) &= \\operatorname{sin}\\left(2 \\pi x \\right), \\\\\n",
" \\frac{\\partial}{\\partial t} u(x,0) &= 2 \\pi \\operatorname{sin}\\left(2 \\pi x \\right)\n",
" &&\\text{for } 0 \\leq x \\leq 1.\n",
"\\end{align*}\n",
"\n",
"Compare your numerical solution to the analytic one and use it to debug your code. The corresponding analytic solution is\n",
"\n",
"\\begin{equation}\n",
" u(x,t) = \\operatorname{sin}(2 \\pi x) (\\operatorname{cos}(2 \\pi t) + \\operatorname{sin}(2 \\pi t)).\n",
"\\end{equation}\n",
"\n",
"Then write a unit test for your function based on this problem."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "17a4067a37c5edfaffd6ebffa47942d6",
"grade": true,
"grade_id": "cell-8b8cb282dfe95a03",
"locked": false,
"points": 0,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# Do your own testing here ...\n",
"\n",
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "fbb85c41f7de70517abcdc482a875fd5",
"grade": true,
"grade_id": "cell-e83b24067743e433",
"locked": false,
"points": 3,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def test_pdeHyperbolic():\n",
" # YOUR CODE HERE\n",
" raise NotImplementedError()\n",
" \n",
"test_pdeHyperbolic()"
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "42880f6025ac92dc08f1c44a5bbf4312",
"grade": false,
"grade_id": "cell-99399435c164c96d",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"#### Problem 2:\n",
"\n",
"\\begin{align*}\n",
" \\frac{\\partial^2}{\\partial t^2} u(x,t)\n",
" - \\frac{\\partial^2}{\\partial x^2} u(x,t) &= 0,\n",
" &&\\text{for }0 \\leq x \\leq 1 \\text{ and } 0 \\leq t \\leq 2\\\\\n",
" u(0,t) = u(l,t) &= 0, \n",
" &&\\text{for } t > 0 \\\\\n",
" u(x,0) &= \\left\\{ \\begin{array}{cc} +1 & \\text{for } x<0.5 \\\\ -1 & \\text{for } x \\geq 0.5 \\end{array} \\right., \\\\\n",
" \\frac{\\partial}{\\partial t} u(x,0) &= 0\n",
" &&\\text{for } 0 \\leq x \\leq 1.\n",
"\\end{align*}\n",
"\n",
"Use $m=200$ and $n=400$ to discretize the spatial and time grids, respectively."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "1b6e620f32a96a23c736b827e282bc6b",
"grade": true,
"grade_id": "cell-ea865d2efbc574d6",
"locked": false,
"points": 2,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "3bb27614fe36d18302f4ace6442c67d9",
"grade": false,
"grade_id": "cell-a783908a0db32a24",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"### Task 3\n",
"\n",
"Animate your solutions! To this end you can use the following code:\n",
"\n",
"```python\n",
"\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",
"# 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",
"\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( ... ) # update plot\n",
"\n",
"# Show the animation\n",
"frames = np.arange(1, np.size(t)) # t is the time grid here\n",
"myAnimation = animation.FuncAnimation(fig, animate, frames, interval = 20)\n",
"myAnimation\n",
"\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "5fa181b7e8df93bd0dd35613bc553701",
"grade": true,
"grade_id": "cell-0bd05d93759bd652",
"locked": false,
"points": 3,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# Animate problem 1 here ...\n",
"\n",
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "d4f93f8615bd3576e0248a27d5c1d664",
"grade": true,
"grade_id": "cell-43b3222743c428f5",
"locked": false,
"points": 0,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# Animate problem 2 here ...\n",
"\n",
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
}
],
"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.10"
}
},
"nbformat": 4,
"nbformat_minor": 4
}