484 lines
14 KiB
Plaintext
484 lines
14 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 = \"Koen Vendrig, Kees van Kempen\""
|
|
]
|
|
},
|
|
{
|
|
"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": 1,
|
|
"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 numpy as np"
|
|
]
|
|
},
|
|
{
|
|
"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",
|
|
" # TODO: Docstring.\n",
|
|
" \n",
|
|
" n = len(t)\n",
|
|
" m = len(x)\n",
|
|
" \n",
|
|
" # TODO: The stepsize is defined by fixed h (k), but the input x (t)\n",
|
|
" # could allow variable step size. What should it be?\n",
|
|
" #h = x[1] - x[0]\n",
|
|
" l = x[-1]\n",
|
|
" h = l/(m - 1)\n",
|
|
" k = t[1] - t[0]\n",
|
|
" λ = a*k/h\n",
|
|
" \n",
|
|
" A = np.eye(m - 2)*2*(1 - λ**2) + ( np.eye(m - 2, m - 2, 1) + np.eye(m - 2, m - 2, -1) )*λ**2\n",
|
|
" \n",
|
|
" w = np.zeros((n, m))\n",
|
|
" \n",
|
|
" # Set initial values w[i, 0] = f[x[i]] and w[i, 1] = ...\n",
|
|
" # TODO: Fix previous comment.\n",
|
|
" w[:, 0] = f(x)\n",
|
|
" # TODO: What does the length of x[1:]?\n",
|
|
" w[:, 1] = (1 - λ**2)*f[x[1:m - 1]] + λ**2/2*f(x[2:m]) - λ**2/2*f(x[2:m]) + k*g(x[1:m - 1])\n",
|
|
" \n",
|
|
" for j in range(2, n):\n",
|
|
" w[:, j] = "
|
|
]
|
|
},
|
|
{
|
|
"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",
|
|
"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
|
|
}
|