{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Generalization: reflecting boundaries\n",
"<div id=\"wave:pde2:Neumann\"></div>\n",
"\n",
"The boundary condition $u=0$ in a wave equation reflects the wave, but\n",
"$u$ changes sign at the boundary, while the condition $u_x=0$ reflects\n",
"the wave as a mirror and preserves the sign, see a [web page](mov-wave/demo_BC_gaussian/index.html) or a\n",
"[movie file](mov-wave/demo_BC_gaussian/movie.flv) for\n",
"demonstration.\n",
"\n",
"\n",
"Our next task is to explain how to implement the boundary\n",
"condition $u_x=0$, which is\n",
"more complicated to express numerically and also to implement than\n",
"a given value of $u$.\n",
"We shall present two methods for implementing $u_x=0$\n",
"in a finite difference scheme, one based on deriving a modified\n",
"stencil at the boundary, and another one based on extending the mesh\n",
"with ghost cells and ghost points.\n",
"\n",
"\n",
"## Neumann boundary condition\n",
"<div id=\"wave:pde2:Neumann:bc\"></div>\n",
"\n",
"\n",
"When a wave hits a boundary and is to be reflected back, one applies\n",
"the condition"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde1:Neumann:0\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" \\frac{\\partial u}{\\partial n} \\equiv \\boldsymbol{n}\\cdot\\nabla u = 0\n",
"\\label{wave:pde1:Neumann:0} \\tag{1}\n",
"\\thinspace .\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The derivative $\\partial /\\partial n$ is in the\n",
"outward normal direction from a general boundary.\n",
"For a 1D domain $[0,L]$,\n",
"we have that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\left.\\frac{\\partial}{\\partial n}\\right\\vert_{x=L} =\n",
"\\left.\\frac{\\partial}{\\partial x}\\right\\vert_{x=L},\\quad\n",
"\\left.\\frac{\\partial}{\\partial n}\\right\\vert_{x=0} = -\n",
"\\left.\\frac{\\partial}{\\partial x}\\right\\vert_{x=0}\\thinspace .\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Boundary condition terminology.**\n",
"\n",
"Boundary conditions\n",
"that specify the value of $\\partial u/\\partial n$\n",
"(or shorter $u_n$) are known as\n",
"[Neumann](http://en.wikipedia.org/wiki/Neumann_boundary_condition) conditions, while [Dirichlet conditions](http://en.wikipedia.org/wiki/Dirichlet_conditions)\n",
"refer to specifications of $u$.\n",
"When the values are zero ($\\partial u/\\partial n=0$ or $u=0$) we speak\n",
"about *homogeneous* Neumann or Dirichlet conditions.\n",
"\n",
"\n",
"\n",
"## Discretization of derivatives at the boundary\n",
"<div id=\"wave:pde2:Neumann:discr\"></div>\n",
"\n",
"\n",
"How can we incorporate the condition ([1](#wave:pde1:Neumann:0))\n",
"in the finite difference scheme? Since we have used central\n",
"differences in all the other approximations to derivatives in the\n",
"scheme, it is tempting to implement ([1](#wave:pde1:Neumann:0)) at\n",
"$x=0$ and $t=t_n$ by the difference"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde1:Neumann:0:cd\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"[D_{2x} u]^n_0 = \\frac{u_{-1}^n - u_1^n}{2\\Delta x} = 0\n",
"\\thinspace .\n",
"\\label{wave:pde1:Neumann:0:cd} \\tag{2}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The problem is that $u_{-1}^n$ is not a $u$ value that is being\n",
"computed since the point is outside the mesh. However, if we combine\n",
"([2](#wave:pde1:Neumann:0:cd)) with the scheme\n",
"<!-- ([wave:pde1:step4](#wave:pde1:step4)) -->"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde1:Neumann:0:scheme\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u^{n+1}_i = -u^{n-1}_i + 2u^n_i + C^2\n",
"\\left(u^{n}_{i+1}-2u^{n}_{i} + u^{n}_{i-1}\\right),\n",
"\\label{wave:pde1:Neumann:0:scheme} \\tag{3}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"for $i=0$, we can eliminate the fictitious value $u_{-1}^n$. We see that\n",
"$u_{-1}^n=u_1^n$ from ([2](#wave:pde1:Neumann:0:cd)), which\n",
"can be used in ([3](#wave:pde1:Neumann:0:scheme)) to\n",
"arrive at a modified scheme for the boundary point $u_0^{n+1}$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto1\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u^{n+1}_i = -u^{n-1}_i + 2u^n_i + 2C^2\n",
"\\left(u^{n}_{i+1}-u^{n}_{i}\\right),\\quad i=0 \\thinspace . \\label{_auto1} \\tag{4}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Figure](#wave:pde1:fig:Neumann:stencil) visualizes this equation\n",
"for computing $u^3_0$ in terms of $u^2_0$, $u^1_0$, and\n",
"$u^2_1$.\n",
"\n",
"<!-- dom:FIGURE: [mov-wave/N_stencil_gpl/stencil_n_left.png, width=500] Modified stencil at a boundary with a Neumann condition. <div id=\"wave:pde1:fig:Neumann:stencil\"></div> -->\n",
"<!-- begin figure -->\n",
"<div id=\"wave:pde1:fig:Neumann:stencil\"></div>\n",
"\n",
"<p>Modified stencil at a boundary with a Neumann condition.</p>\n",
"<img src=\"mov-wave/N_stencil_gpl/stencil_n_left.png\" width=500>\n",
"\n",
"<!-- end figure -->\n",
"\n",
"\n",
"Similarly, ([1](#wave:pde1:Neumann:0)) applied at $x=L$\n",
"is discretized by a central difference"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde1:Neumann:0:cd2\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{u_{N_x+1}^n - u_{N_x-1}^n}{2\\Delta x} = 0\n",
"\\thinspace .\n",
"\\label{wave:pde1:Neumann:0:cd2} \\tag{5}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Combined with the scheme for $i=N_x$ we get a modified scheme for\n",
"the boundary value $u_{N_x}^{n+1}$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto2\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u^{n+1}_i = -u^{n-1}_i + 2u^n_i + 2C^2\n",
"\\left(u^{n}_{i-1}-u^{n}_{i}\\right),\\quad i=N_x \\thinspace . \\label{_auto2} \\tag{6}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The modification of the scheme at the boundary is also required for\n",
"the special formula for the first time step. How the stencil moves\n",
"through the mesh and is modified at the boundary can be illustrated by\n",
"an animation in a [web page](${doc_notes}/book/html/mov-wave/N_stencil_gpl/index.html)\n",
"or a [movie file](${docraw}/mov-wave/N_stencil_gpl/movie.ogg).\n",
"\n",
"\n",
"\n",
"\n",
"## Implementation of Neumann conditions\n",
"<div id=\"wave:pde2:Neumann:impl\"></div>\n",
"\n",
"We have seen in the preceding section\n",
"that the special formulas for the boundary points\n",
"arise from replacing $u_{i-1}^n$ by $u_{i+1}^n$ when computing\n",
"$u_i^{n+1}$ from the stencil formula for $i=0$. Similarly, we\n",
"replace $u_{i+1}^n$ by $u_{i-1}^n$ in the stencil formula\n",
"for $i=N_x$. This observation can conveniently\n",
"be used in the coding: we just work with the general stencil formula,\n",
"but write the code such that it is easy to replace `u[i-1]` by\n",
"`u[i+1]` and vice versa. This is achieved by\n",
"having the indices `i+1` and `i-1` as variables `ip1` (`i` plus 1)\n",
"and `im1` (`i` minus 1), respectively.\n",
"At the boundary we can easily define `im1=i+1` while we use\n",
"`im1=i-1` in the internal parts of the mesh. Here are the details\n",
"of the implementation (note that the updating formula for `u[i]`\n",
"is the general stencil formula):"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" i = 0\n",
" ip1 = i+1\n",
" im1 = ip1 # i-1 -> i+1\n",
" u[i] = u_n[i] + C2*(u_n[im1] - 2*u_n[i] + u_n[ip1])\n",
" \n",
" i = Nx\n",
" im1 = i-1\n",
" ip1 = im1 # i+1 -> i-1\n",
" u[i] = u_n[i] + C2*(u_n[im1] - 2*u_n[i] + u_n[ip1])\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can in fact create one loop over both the internal and boundary\n",
"points and use only one updating formula:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" for i in range(0, Nx+1):\n",
" ip1 = i+1 if i < Nx else i-1\n",
" im1 = i-1 if i > 0 else i+1\n",
" u[i] = u_n[i] + C2*(u_n[im1] - 2*u_n[i] + u_n[ip1])\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The program [`wave1D_n0.py`](${src_wave}/wave1D/wave1D_n0.py)\n",
"contains a complete implementation of the 1D wave equation with\n",
"boundary conditions $u_x = 0$ at $x=0$ and $x=L$.\n",
"\n",
"It would be nice to modify the `test_quadratic` test case from the\n",
"`wave1D_u0.py` with Dirichlet conditions, described in the section [wave:pde1:impl:vec:verify:quadratic](#wave:pde1:impl:vec:verify:quadratic). However, the Neumann\n",
"conditions require the polynomial variation in the $x$ direction to\n",
"be of third degree, which causes challenging problems when\n",
"designing a test where the numerical solution is known exactly.\n",
"[Exercise 9: Verification by a cubic polynomial in space](#wave:fd2:exer:verify:cubic) outlines ideas and code\n",
"for this purpose. The only test in `wave1D_n0.py` is to start\n",
"with a plug wave at rest and see that the initial condition is\n",
"reached again perfectly after one period of motion, but such\n",
"a test requires $C=1$ (so the numerical solution coincides with\n",
"the exact solution of the PDE, see the section [Numerical dispersion relation](wave_analysis.ipynb#wave:pde1:num:dispersion)).\n",
"\n",
"\n",
"## Index set notation\n",
"<div id=\"wave:indexset\"></div>\n",
"\n",
"\n",
"To improve our mathematical writing and our implementations,\n",
"it is wise to introduce a special notation for index sets. This means\n",
"that we write\n",
"$x_i$, followed by $i\\in\\mathcal{I}_x$, instead of $i=0,\\ldots,N_x$.\n",
"Obviously, $\\mathcal{I}_x$ must be the index set $\\mathcal{I}_x =\\{0,\\ldots,N_x\\}$, but it\n",
"is often advantageous to have a symbol for this set rather than\n",
"specifying all its elements (all the time, as we have done up to\n",
"now). This new notation saves writing and makes\n",
"specifications of algorithms and their implementation as computer code\n",
"simpler.\n",
"\n",
"The first index in the set will be denoted $\\mathcal{I}_x^0$\n",
"and the last $\\mathcal{I}_x^{-1}$. When we need to skip the first element of\n",
"the set, we use $\\mathcal{I}_x^{+}$ for the remaining subset\n",
"$\\mathcal{I}_x^{+}=\\{1,\\ldots,N_x\\}$. Similarly, if the last element is\n",
"to be dropped, we write $\\mathcal{I}_x^{-}=\\{0,\\ldots,N_x-1\\}$ for the\n",
"remaining indices.\n",
"All the\n",
"indices corresponding to inner grid points are specified by\n",
"$\\mathcal{I}_x^i=\\{1,\\ldots,N_x-1\\}$. For the time domain we find it\n",
"natural to explicitly use 0 as the first index, so we will usually\n",
"write $n=0$ and $t_0$ rather than $n=\\mathcal{I}_t^0$. We also avoid notation\n",
"like $x_{\\mathcal{I}_x^{-1}}$ and will instead use $x_i$, $i={\\mathcal{I}_x^{-1}}$.\n",
"\n",
"The Python code associated with index sets applies the following\n",
"conventions:\n",
"\n",
"\n",
"<table border=\"1\">\n",
"<thead>\n",
"<tr><th align=\"center\"> Notation </th> <th align=\"center\"> Python </th> </tr>\n",
"</thead>\n",
"<tbody>\n",
"<tr><td align=\"left\"> $\\mathcal{I}_x$ </td> <td align=\"left\"> <code>Ix</code> </td> </tr>\n",
"<tr><td align=\"left\"> $\\mathcal{I}_x^0$ </td> <td align=\"left\"> <code>Ix[0]</code> </td> </tr>\n",
"<tr><td align=\"left\"> $\\mathcal{I}_x^{-1}$ </td> <td align=\"left\"> <code>Ix[-1]</code> </td> </tr>\n",
"<tr><td align=\"left\"> $\\mathcal{I}_x^{-}$ </td> <td align=\"left\"> <code>Ix[:-1]</code> </td> </tr>\n",
"<tr><td align=\"left\"> $\\mathcal{I}_x^{+}$ </td> <td align=\"left\"> <code>Ix[1:]</code> </td> </tr>\n",
"<tr><td align=\"left\"> $\\mathcal{I}_x^i$ </td> <td align=\"left\"> <code>Ix[1:-1]</code> </td> </tr>\n",
"</tbody>\n",
"</table>\n",
"**Why index sets are useful.**\n",
"\n",
"An important feature of the index set notation is that it\n",
"keeps our formulas and code independent of how\n",
"we count mesh points. For example, the notation $i\\in\\mathcal{I}_x$ or $i=\\mathcal{I}_x^0$\n",
"remains the same whether $\\mathcal{I}_x$ is defined as above or as starting at 1,\n",
"i.e., $\\mathcal{I}_x=\\{1,\\ldots,Q\\}$. Similarly, we can in the code define\n",
"`Ix=range(Nx+1)` or `Ix=range(1,Q)`, and expressions\n",
"like `Ix[0]` and `Ix[1:-1]` remain correct. One application where\n",
"the index set notation is convenient is\n",
"conversion of code from a language where arrays has base index 0 (e.g.,\n",
"Python and C) to languages where the base index is 1 (e.g., MATLAB and\n",
"Fortran). Another important application is implementation of\n",
"Neumann conditions via ghost points (see next section).\n",
"\n",
"\n",
"\n",
"For the current problem setting in the $x,t$ plane, we work with\n",
"the index sets"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto3\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\mathcal{I}_x = \\{0,\\ldots,N_x\\},\\quad \\mathcal{I}_t = \\{0,\\ldots,N_t\\},\n",
"\\label{_auto3} \\tag{7}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"defined in Python as"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"mathcal{I}_x = range(0, Nx+1)\n",
"mathcal{I}_t = range(0, Nt+1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A finite difference scheme can with the index set notation be specified as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"u_i^{n+1} &= u^n_i - \\frac{1}{2}\n",
"C^2\\left(u^{n}_{i+1}-2u^{n}_{i} + u^{n}_{i-1}\\right),\\quad,\n",
"i\\in\\mathcal{I}_x^i,\\ n=0,\\\\ \n",
"u^{n+1}_i &= -u^{n-1}_i + 2u^n_i + C^2\n",
"\\left(u^{n}_{i+1}-2u^{n}_{i}+u^{n}_{i-1}\\right),\n",
"\\quad i\\in\\mathcal{I}_x^i,\\ n\\in\\mathcal{I}_t^i,\\\\ \n",
"u_i^{n+1} &= 0,\n",
"\\quad i=\\mathcal{I}_x^0,\\ n\\in\\mathcal{I}_t^{-},\\\\ \n",
"u_i^{n+1} &= 0,\n",
"\\quad i=\\mathcal{I}_x^{-1},\\ n\\in\\mathcal{I}_t^{-}\\thinspace .\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The corresponding implementation becomes"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Initial condition\n",
"for i in mathcal{I}_x[1:-1]:\n",
" u[i] = u_n[i] - 0.5*C2*(u_n[i-1] - 2*u_n[i] + u_n[i+1])\n",
"\n",
"# Time loop\n",
"for n in mathcal{I}_t[1:-1]:\n",
" # Compute internal points\n",
" for i in mathcal{I}_x[1:-1]:\n",
" u[i] = - u_nm1[i] + 2*u_n[i] + \\\n",
" C2*(u_n[i-1] - 2*u_n[i] + u_n[i+1])\n",
" # Compute boundary conditions\n",
" i = mathcal{I}_x[0]; u[i] = 0\n",
" i = mathcal{I}_x[-1]; u[i] = 0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Notice.**\n",
"\n",
"The program [`wave1D_dn.py`](src-wave/wave1D/python/wave1D_dn.py)\n",
"applies the index set notation and\n",
"solves the 1D wave equation $u_{tt}=c^2u_{xx}+f(x,t)$ with\n",
"quite general boundary and initial conditions:\n",
"\n",
" * $x=0$: $u=U_0(t)$ or $u_x=0$\n",
"\n",
" * $x=L$: $u=U_L(t)$ or $u_x=0$\n",
"\n",
" * $t=0$: $u=I(x)$\n",
"\n",
" * $t=0$: $u_t=V(x)$\n",
"\n",
"The program combines Dirichlet and Neumann conditions, scalar and vectorized\n",
"implementation of schemes, and the index set notation into one piece of code.\n",
"A lot of test examples are also included in the program:\n",
"\n",
" * A rectangular plug-shaped initial condition. (For $C=1$ the solution\n",
" will be a rectangle that jumps one cell per time step, making the case\n",
" well suited for verification.)\n",
"\n",
" * A Gaussian function as initial condition.\n",
"\n",
" * A triangular profile as initial condition, which resembles the\n",
" typical initial shape of a guitar string.\n",
"\n",
" * A sinusoidal variation of $u$ at $x=0$ and either $u=0$ or\n",
" $u_x=0$ at $x=L$.\n",
"\n",
" * An analytical solution $u(x,t)=\\cos(m\\pi t/L)\\sin({\\frac{1}{2}}m\\pi x/L)$, which can be used for convergence rate tests.\n",
"\n",
"\n",
"\n",
"[hpl 1: Should include some experiments here or make exercises. Qualitative\n",
"behavior of the wave equation can be exemplified.]\n",
"\n",
"## Verifying the implementation of Neumann conditions\n",
"<div id=\"wave:pde1:verify\"></div>\n",
"\n",
"\n",
"How can we test that the Neumann conditions are correctly implemented?\n",
"The `solver` function in the `wave1D_dn.py` program described in the\n",
"box above accepts Dirichlet or Neumann conditions at $x=0$ and $x=L$.\n",
"mathcal{I}_t is tempting to apply a quadratic solution as described in\n",
"the sections [wave:pde2:fd](#wave:pde2:fd) and [wave:pde1:impl:verify:quadratic](#wave:pde1:impl:verify:quadratic),\n",
"but it turns out that this solution is no longer an exact solution\n",
"of the discrete equations if a Neumann condition is implemented on\n",
"the boundary. A linear solution does not help since we only have\n",
"homogeneous Neumann conditions in `wave1D_dn.py`, and we are\n",
"consequently left with testing just a constant solution: $u=\\hbox{const}$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def test_constant():\n",
" \"\"\"\n",
" Check the scalar and vectorized versions for\n",
" a constant u(x,t). We simulate in [0, L] and apply\n",
" Neumann and Dirichlet conditions at both ends.\n",
" \"\"\"\n",
" u_const = 0.45\n",
" u_exact = lambda x, t: u_const\n",
" I = lambda x: u_exact(x, 0)\n",
" V = lambda x: 0\n",
" f = lambda x, t: 0\n",
"\n",
" def assert_no_error(u, x, t, n):\n",
" u_e = u_exact(x, t[n])\n",
" diff = np.abs(u - u_e).max()\n",
" msg = 'diff=%E, t_%d=%g' % (diff, n, t[n])\n",
" tol = 1E-13\n",
" assert diff < tol, msg\n",
"\n",
" for U_0 in (None, lambda t: u_const):\n",
" for U_L in (None, lambda t: u_const):\n",
" L = 2.5\n",
" c = 1.5\n",
" C = 0.75\n",
" Nx = 3 # Very coarse mesh for this exact test\n",
" dt = C*(L/Nx)/c\n",
" T = 18 # long time integration\n",
"\n",
" solver(I, V, f, c, U_0, U_L, L, dt, C, T,\n",
" user_action=assert_no_error,\n",
" version='scalar')\n",
" solver(I, V, f, c, U_0, U_L, L, dt, C, T,\n",
" user_action=assert_no_error,\n",
" version='vectorized')\n",
" print U_0, U_L"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The quadratic solution is very useful for testing, but it requires\n",
"Dirichlet conditions at both ends.\n",
"\n",
"Another test may utilize the fact that the approximation error vanishes\n",
"when the Courant number is unity. We can, for example, start with a\n",
"plug profile as initial condition, let this wave split into two plug waves,\n",
"one in each direction, and check that the two plug waves come back and\n",
"form the initial condition again after \"one period\" of the solution\n",
"process. Neumann conditions can be applied at both ends. A proper\n",
"test function reads"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def test_plug():\n",
" \"\"\"Check that an initial plug is correct back after one period.\"\"\"\n",
" L = 1.0\n",
" c = 0.5\n",
" dt = (L/10)/c # Nx=10\n",
" I = lambda x: 0 if abs(x-L/2.0) > 0.1 else 1\n",
"\n",
" u_s, x, t, cpu = solver(\n",
" I=I,\n",
" V=None, f=None, c=0.5, U_0=None, U_L=None, L=L,\n",
" dt=dt, C=1, T=4, user_action=None, version='scalar')\n",
" u_v, x, t, cpu = solver(\n",
" I=I,\n",
" V=None, f=None, c=0.5, U_0=None, U_L=None, L=L,\n",
" dt=dt, C=1, T=4, user_action=None, version='vectorized')\n",
" tol = 1E-13\n",
" diff = abs(u_s - u_v).max()\n",
" assert diff < tol\n",
" u_0 = np.array([I(x_) for x_ in x])\n",
" diff = np.abs(u_s - u_0).max()\n",
" assert diff < tol"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Other tests must rely on an unknown approximation error, so effectively\n",
"we are left with tests on the convergence rate.\n",
"\n",
"## Alternative implementation via ghost cells\n",
"<div id=\"wave:pde1:Neumann:ghost\"></div>\n",
"\n",
"### Idea\n",
"\n",
"Instead of modifying the scheme at the boundary, we can introduce\n",
"extra points outside the domain such that the fictitious values\n",
"$u_{-1}^n$ and $u_{N_x+1}^n$ are defined in the mesh. Adding the\n",
"intervals $[-\\Delta x,0]$ and $[L, L+\\Delta x]$, known as *ghost\n",
"cells*, to the mesh gives us all the needed mesh points, corresponding\n",
"to $i=-1,0,\\ldots,N_x,N_x+1$. The extra points with $i=-1$ and\n",
"$i=N_x+1$ are known as *ghost points*, and values at these points,\n",
"$u_{-1}^n$ and $u_{N_x+1}^n$, are called *ghost values*.\n",
"\n",
"The important idea is\n",
"to ensure that we always have"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u_{-1}^n = u_{1}^n\\hbox{ and } u_{N_x+1}^n = u_{N_x-1}^n,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"because then\n",
"the application of the standard scheme at a boundary point $i=0$ or $i=N_x$\n",
"will be correct and guarantee that the solution is compatible with the\n",
"boundary condition $u_x=0$.\n",
"\n",
"Some readers may find it strange to just extend the domain with ghost\n",
"cells as a general technique, because in some problems there is a\n",
"completely different medium with different physics and equations right\n",
"outside of a boundary. Nevertheless, one should view the ghost cell\n",
"technique as a purely mathematical technique, which is valid in the\n",
"limit $\\Delta x \\rightarrow 0$ and helps us to implement derivatives.\n",
"\n",
"\n",
"### Implementation\n",
"\n",
"The `u` array now needs extra elements corresponding to the ghost\n",
"points. Two new point values are needed:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"u = zeros(Nx+3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The arrays `u_n` and `u_nm1` must be defined accordingly.\n",
"\n",
"Unfortunately, a major indexing problem arises with ghost cells.\n",
"The reason is that Python indices *must* start\n",
"at 0 and `u[-1]` will always mean the last element in `u`.\n",
"This fact gives, apparently, a mismatch between the mathematical\n",
"indices $i=-1,0,\\ldots,N_x+1$ and the Python indices running over\n",
"`u`: `0,..,Nx+2`. One remedy is to change the mathematical indexing\n",
"of $i$ in the scheme and write"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u^{n+1}_i = \\cdots,\\quad i=1,\\ldots,N_x+1,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"instead of $i=0,\\ldots,N_x$ as we have previously used. The ghost\n",
"points now correspond to $i=0$ and $i=N_x+1$.\n",
"A better solution is to use the ideas of the section [Index set notation](#wave:indexset):\n",
"we hide the specific index value in an index set and operate with\n",
"inner and boundary points using the index set notation.\n",
"\n",
"To this end, we define `u` with proper length and `mathcal{I}_x` to be the corresponding\n",
"indices for the real physical mesh points ($1,2,\\ldots,N_x+1$):"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" u = zeros(Nx+3)\n",
" mathcal{I}_x = range(1, u.shape[0]-1)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That is, the boundary points have indices `mathcal{I}_x[0]` and `mathcal{I}_x[-1]` (as before).\n",
"We first update the solution at all physical mesh points (i.e., interior\n",
"points in the mesh):"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"for i in mathcal{I}_x:\n",
" u[i] = - u_nm1[i] + 2*u_n[i] + \\\n",
" C2*(u_n[i-1] - 2*u_n[i] + u_n[i+1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The indexing becomes a bit more complicated when we call functions like\n",
"`V(x)` and `f(x, t)`, as we must remember that the appropriate\n",
"$x$ coordinate is given as `x[i-mathcal{I}_x[0]]`:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"for i in mathcal{I}_x:\n",
" u[i] = u_n[i] + dt*V(x[i-mathcal{I}_x[0]]) + \\\n",
" 0.5*C2*(u_n[i-1] - 2*u_n[i] + u_n[i+1]) + \\\n",
" 0.5*dt2*f(x[i-mathcal{I}_x[0]], t[0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"mathcal{I}_t remains to update the solution at ghost points, i.e., `u[0]`\n",
"and `u[-1]` (or `u[Nx+2]`). For a boundary condition $u_x=0$,\n",
"the ghost value must equal the value at the associated inner mesh\n",
"point. Computer code makes this statement precise:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"i = mathcal{I}_x[0] # x=0 boundary\n",
"u[i-1] = u[i+1]\n",
"i = mathcal{I}_x[-1] # x=L boundary\n",
"u[i+1] = u[i-1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The physical solution to be plotted is now in `u[1:-1]`, or\n",
"equivalently `u[mathcal{I}_x[0]:mathcal{I}_x[-1]+1]`, so this slice is\n",
"the quantity to be returned from a solver function.\n",
"A complete implementation appears in the program\n",
"[`wave1D_n0_ghost.py`](${src_wave}/wave1D/wave1D_n0_ghost.py).\n",
"\n",
"**Warning.**\n",
"\n",
"We have to be careful with how the spatial and temporal mesh\n",
"points are stored. Say we let `x` be the physical mesh points,"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"x = linspace(0, L, Nx+1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\"Standard coding\" of the initial condition,"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"for i in mathcal{I}_x:\n",
" u_n[i] = I(x[i])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"becomes wrong, since `u_n` and `x` have different lengths and the index `i`\n",
"corresponds to two different mesh points. In fact, `x[i]` corresponds\n",
"to `u[1+i]`. A correct implementation is"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"for i in mathcal{I}_x:\n",
" u_n[i] = I(x[i-mathcal{I}_x[0]])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Similarly, a source term usually coded as `f(x[i], t[n])` is incorrect\n",
"if `x` is defined to be the physical points, so `x[i]` must be\n",
"replaced by `x[i-mathcal{I}_x[0]]`.\n",
"\n",
"An alternative remedy is to let `x` also cover the ghost points such that\n",
"`u[i]` is the value at `x[i]`.\n",
"\n",
"\n",
"\n",
"The ghost cell is only added to the boundary where we have a Neumann\n",
"condition. Suppose we have a Dirichlet condition at $x=L$ and\n",
"a homogeneous Neumann condition at $x=0$. One ghost cell $[-\\Delta x,0]$\n",
"is added to the mesh, so the index set for the physical points\n",
"becomes $\\{1,\\ldots,N_x+1\\}$. A relevant implementation\n",
"is"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"u = zeros(Nx+2)\n",
"mathcal{I}_x = range(1, u.shape[0])\n",
"...\n",
"for i in mathcal{I}_x[:-1]:\n",
" u[i] = - u_nm1[i] + 2*u_n[i] + \\\n",
" C2*(u_n[i-1] - 2*u_n[i] + u_n[i+1]) + \\\n",
" dt2*f(x[i-mathcal{I}_x[0]], t[n])\n",
"i = mathcal{I}_x[-1]\n",
"u[i] = U_0 # set Dirichlet value\n",
"i = mathcal{I}_x[0]\n",
"u[i-1] = u[i+1] # update ghost value"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The physical solution to be plotted is now in `u[1:]`\n",
"or (as always) `u[mathcal{I}_x[0]:mathcal{I}_x[-1]+1]`.\n",
"\n",
"\n",
"# Generalization: variable wave velocity\n",
"<div id=\"wave:pde2:var:c\"></div>\n",
"\n",
"\n",
"Our next generalization of the 1D wave equation ([wave:pde1](#wave:pde1)) or\n",
"([wave:pde2](#wave:pde2)) is to allow for a variable wave velocity $c$:\n",
"$c=c(x)$, usually motivated by wave motion in a domain composed of\n",
"different physical media. When the media differ in physical properties\n",
"like density or porosity, the wave velocity $c$ is affected and\n",
"will depend on the position in space.\n",
"[Figure](#wave:pde1:fig:pulse1:two:media) shows a wave\n",
"propagating in one medium $[0, 0.7]\\cup [0.9,1]$ with wave\n",
"velocity $c_1$ (left) before it enters a second medium $(0.7,0.9)$\n",
"with wave velocity $c_2$ (right). When the wave meets the boundary\n",
"where $c$ jumps from $c_1$ to $c_2$, a part of the wave is reflected back\n",
"into the first medium (the *reflected* wave), while one part is\n",
"transmitted through the second medium (the *transmitted* wave).\n",
"\n",
"\n",
"<!-- dom:FIGURE: [fig-wave/pulse1_in_two_media.png, width=800] Left: wave entering another medium; right: transmitted and reflected wave. <div id=\"wave:pde1:fig:pulse1:two:media\"></div> -->\n",
"<!-- begin figure -->\n",
"<div id=\"wave:pde1:fig:pulse1:two:media\"></div>\n",
"\n",
"<p>Left: wave entering another medium; right: transmitted and reflected wave.</p>\n",
"<img src=\"fig-wave/pulse1_in_two_media.png\" width=800>\n",
"\n",
"<!-- end figure -->\n",
"\n",
"\n",
"\n",
"## The model PDE with a variable coefficient\n",
"\n",
"Instead of working with the squared quantity $c^2(x)$, we\n",
"shall for notational convenience introduce $q(x) = c^2(x)$.\n",
"A 1D wave equation with variable wave velocity often takes the form"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:var:c:pde\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial^2 u}{\\partial t^2} =\n",
"\\frac{\\partial}{\\partial x}\\left( q(x)\n",
"\\frac{\\partial u}{\\partial x}\\right) + f(x,t)\n",
"\\label{wave:pde2:var:c:pde} \\tag{8}\n",
"\\thinspace .\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the most frequent form of a wave\n",
"equation with variable wave velocity,\n",
"but other forms also appear, see the section [wave:app:string](#wave:app:string)\n",
"and equation ([wave:app:string:model2](#wave:app:string:model2)).\n",
"\n",
"As usual, we sample ([8](#wave:pde2:var:c:pde)) at a mesh point,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{\\partial^2 }{\\partial t^2} u(x_i,t_n) =\n",
"\\frac{\\partial}{\\partial x}\\left( q(x_i)\n",
"\\frac{\\partial}{\\partial x} u(x_i,t_n)\\right) + f(x_i,t_n),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where the only new term to discretize is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{\\partial}{\\partial x}\\left( q(x_i)\n",
"\\frac{\\partial}{\\partial x} u(x_i,t_n)\\right) = \\left[\n",
"\\frac{\\partial}{\\partial x}\\left( q(x)\n",
"\\frac{\\partial u}{\\partial x}\\right)\\right]^n_i\n",
"\\thinspace .\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Discretizing the variable coefficient\n",
"<div id=\"wave:pde2:var:c:ideas\"></div>\n",
"\n",
"The principal idea is to first discretize the outer derivative.\n",
"Define"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\phi = q(x)\n",
"\\frac{\\partial u}{\\partial x},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and use a centered derivative around $x=x_i$ for the derivative of $\\phi$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\left[\\frac{\\partial\\phi}{\\partial x}\\right]^n_i\n",
"\\approx \\frac{\\phi_{i+\\frac{1}{2}} - \\phi_{i-\\frac{1}{2}}}{\\Delta x}\n",
"= [D_x\\phi]^n_i\n",
"\\thinspace .\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then discretize"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\phi_{i+\\frac{1}{2}} = q_{i+\\frac{1}{2}}\n",
"\\left[\\frac{\\partial u}{\\partial x}\\right]^n_{i+\\frac{1}{2}}\n",
"\\approx q_{i+\\frac{1}{2}} \\frac{u^n_{i+1} - u^n_{i}}{\\Delta x}\n",
"= [q D_x u]_{i+\\frac{1}{2}}^n\n",
"\\thinspace .\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Similarly,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\phi_{i-\\frac{1}{2}} = q_{i-\\frac{1}{2}}\n",
"\\left[\\frac{\\partial u}{\\partial x}\\right]^n_{i-\\frac{1}{2}}\n",
"\\approx q_{i-\\frac{1}{2}} \\frac{u^n_{i} - u^n_{i-1}}{\\Delta x}\n",
"= [q D_x u]_{i-\\frac{1}{2}}^n\n",
"\\thinspace .\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These intermediate results are now combined to"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:var:c:formula\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\left[\n",
"\\frac{\\partial}{\\partial x}\\left( q(x)\n",
"\\frac{\\partial u}{\\partial x}\\right)\\right]^n_i\n",
"\\approx \\frac{1}{\\Delta x^2}\n",
"\\left( q_{i+\\frac{1}{2}} \\left({u^n_{i+1} - u^n_{i}}\\right)\n",
"- q_{i-\\frac{1}{2}} \\left({u^n_{i} - u^n_{i-1}}\\right)\\right)\n",
"\\label{wave:pde2:var:c:formula} \\tag{9}\n",
"\\thinspace .\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With operator notation we can write the discretization as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:var:c:formula:op\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\left[\n",
"\\frac{\\partial}{\\partial x}\\left( q(x)\n",
"\\frac{\\partial u}{\\partial x}\\right)\\right]^n_i\n",
"\\approx [D_x (\\overline{q}^{x} D_x u)]^n_i\n",
"\\label{wave:pde2:var:c:formula:op} \\tag{10}\n",
"\\thinspace .\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Do not use the chain rule on the spatial derivative term!**\n",
"\n",
"Many are tempted to use the chain rule on the\n",
"term $\\frac{\\partial}{\\partial x}\\left( q(x)\n",
"\\frac{\\partial u}{\\partial x}\\right)$, but this is not a good idea\n",
"when discretizing such a term.\n",
"\n",
"The term with a variable coefficient expresses the net flux\n",
"$qu_x$ into a small volume (i.e., interval in 1D):"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{\\partial}{\\partial x}\\left( q(x)\n",
"\\frac{\\partial u}{\\partial x}\\right) \\approx\n",
"\\frac{1}{\\Delta x}(q(x+\\Delta x)u_x(x+\\Delta x) - q(x)u_x(x))\\thinspace .\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our discretization reflects this\n",
"principle directly: $qu_x$ at the right end of the cell minus $qu_x$\n",
"at the left end, because this follows from the formula\n",
"([9](#wave:pde2:var:c:formula)) or $[D_x(q D_x u)]^n_i$.\n",
"\n",
"When using the chain rule, we get two\n",
"terms $qu_{xx} + q_xu_x$. The typical discretization is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:var:c:chainrule_scheme\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"[D_x q D_x u + D_{2x}q D_{2x} u]_i^n,\n",
"\\label{wave:pde2:var:c:chainrule_scheme} \\tag{11}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Writing this out shows that it is different from\n",
"$[D_x(q D_x u)]^n_i$ and lacks the physical interpretation of\n",
"net flux into a cell. With a smooth and slowly varying $q(x)$ the\n",
"differences between the two discretizations are not substantial.\n",
"However, when $q$ exhibits (potentially large) jumps,\n",
"$[D_x(q D_x u)]^n_i$ with harmonic averaging of $q$ yields\n",
"a better solution than arithmetic averaging or\n",
"([11](#wave:pde2:var:c:chainrule_scheme)).\n",
"In the literature, the discretization $[D_x(q D_x u)]^n_i$ totally\n",
"dominates and very few mention the alternative in\n",
"([11](#wave:pde2:var:c:chainrule_scheme)).\n",
"\n",
"\n",
"\n",
"<!-- Needs some better explanation here - maybe the exact solution of a -->\n",
"<!-- poisson type problem (piecewise linear solution) failes if we use -->\n",
"<!-- the chain rule? Wesserling has an example, but it is tedious to -->\n",
"<!-- work out. -->\n",
"\n",
"\n",
"## Computing the coefficient between mesh points\n",
"<div id=\"wave:pde2:var:c:means\"></div>\n",
"\n",
"\n",
"If $q$ is a known function of $x$, we can easily evaluate\n",
"$q_{i+\\frac{1}{2}}$ simply as $q(x_{i+\\frac{1}{2}})$ with $x_{i+\\frac{1}{2}} = x_i +\n",
"\\frac{1}{2}\\Delta x$. However, in many cases $c$, and hence $q$, is only\n",
"known as a discrete function, often at the mesh points $x_i$.\n",
"Evaluating $q$ between two mesh points $x_i$ and $x_{i+1}$ must then\n",
"be done by *interpolation* techniques, of which three are of\n",
"particular interest in this context:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:var:c:mean:arithmetic\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"q_{i+\\frac{1}{2}} \\approx\n",
"\\frac{1}{2}\\left( q_{i} + q_{i+1}\\right) =\n",
"[\\overline{q}^{x}]_i\n",
"\\quad \\hbox{(arithmetic mean)}\n",
"\\label{wave:pde2:var:c:mean:arithmetic} \\tag{12}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:var:c:mean:harmonic\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"q_{i+\\frac{1}{2}} \\approx\n",
"2\\left( \\frac{1}{q_{i}} + \\frac{1}{q_{i+1}}\\right)^{-1}\n",
"\\quad \\hbox{(harmonic mean)}\n",
"\\label{wave:pde2:var:c:mean:harmonic} \\tag{13}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:var:c:mean:geometric\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"q_{i+\\frac{1}{2}} \\approx\n",
"\\left(q_{i}q_{i+1}\\right)^{1/2}\n",
"\\quad \\hbox{(geometric mean)}\n",
"\\label{wave:pde2:var:c:mean:geometric} \\tag{14}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The arithmetic mean in ([12](#wave:pde2:var:c:mean:arithmetic)) is by\n",
"far the most commonly used averaging technique and is well suited\n",
"for smooth $q(x)$ functions.\n",
"The harmonic mean is often preferred when $q(x)$ exhibits large\n",
"jumps (which is typical for geological media).\n",
"The geometric mean is less used, but popular in\n",
"discretizations to linearize quadratic\n",
"% if BOOK == \"book\":\n",
"nonlinearities (see the section [vib:ode2:fdm:fquad](#vib:ode2:fdm:fquad) for an example).\n",
"% else:\n",
"nonlinearities.\n",
"% endif\n",
"\n",
"With the operator notation from ([12](#wave:pde2:var:c:mean:arithmetic))\n",
"we can specify the discretization of the complete variable-coefficient\n",
"wave equation in a compact way:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:var:c:scheme:op\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\lbrack D_tD_t u = D_x\\overline{q}^{x}D_x u + f\\rbrack^{n}_i\n",
"\\thinspace .\n",
"\\label{wave:pde2:var:c:scheme:op} \\tag{15}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Strictly speaking, $\\lbrack D_x\\overline{q}^{x}D_x u\\rbrack^{n}_i\n",
"= \\lbrack D_x (\\overline{q}^{x}D_x u)\\rbrack^{n}_i$.\n",
"\n",
"From the compact difference notation we immediately see what kind of differences that\n",
"each term is approximated with. The notation $\\overline{q}^{x}$\n",
"also specifies that the variable coefficient is approximated by\n",
"an arithmetic mean, the definition being\n",
"$[\\overline{q}^{x}]_{i+\\frac{1}{2}}=(q_i+q_{i+1})/2$.\n",
"\n",
"Before implementing, it remains to solve\n",
"([15](#wave:pde2:var:c:scheme:op)) with respect to $u_i^{n+1}$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u^{n+1}_i = - u_i^{n-1} + 2u_i^n + \\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\quad \\left(\\frac{\\Delta t}{\\Delta x}\\right)^2 \\left(\n",
"\\frac{1}{2}(q_{i} + q_{i+1})(u_{i+1}^n - u_{i}^n) -\n",
"\\frac{1}{2}(q_{i} + q_{i-1})(u_{i}^n - u_{i-1}^n)\\right)\n",
"+ \\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:var:c:scheme:impl\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
" \\quad \\Delta t^2 f^n_i\n",
"\\thinspace .\n",
"\\label{wave:pde2:var:c:scheme:impl} \\tag{16}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## How a variable coefficient affects the stability\n",
"<div id=\"wave:pde2:var:c:stability\"></div>\n",
"\n",
"\n",
"The stability criterion derived later (the section [wave:pde1:stability](#wave:pde1:stability))\n",
"reads $\\Delta t\\leq \\Delta x/c$. If $c=c(x)$, the criterion will depend\n",
"on the spatial location. We must therefore choose a $\\Delta t$ that\n",
"is small enough such that no mesh cell has $\\Delta t > \\Delta x/c(x)$.\n",
"That is, we must use the largest $c$ value in the criterion:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto4\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\Delta t \\leq \\beta \\frac{\\Delta x}{\\max_{x\\in [0,L]}c(x)}\n",
"\\thinspace .\n",
"\\label{_auto4} \\tag{17}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The parameter $\\beta$ is included as a safety factor: in some problems with a\n",
"significantly varying $c$ it turns out that one must choose $\\beta <1$ to\n",
"have stable solutions ($\\beta =0.9$ may act as an all-round value).\n",
"\n",
"A different strategy to handle the stability criterion with variable\n",
"wave velocity is to use a spatially varying $\\Delta t$. While the idea\n",
"is mathematically attractive at first sight, the implementation\n",
"quickly becomes very complicated, so we stick to a constant $\\Delta t$\n",
"and a worst case value of $c(x)$ (with a safety factor $\\beta$).\n",
"\n",
"## Neumann condition and a variable coefficient\n",
"<div id=\"wave:pde2:var:c:Neumann\"></div>\n",
"\n",
"Consider a Neumann condition $\\partial u/\\partial x=0$ at $x=L=N_x\\Delta x$,\n",
"discretized as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"[D_{2x} u]^n_i =\n",
"\\frac{u_{i+1}^{n} - u_{i-1}^n}{2\\Delta x} = 0\\quad\\Rightarrow\\quad\n",
"u_{i+1}^n = u_{i-1}^n,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"for $i=N_x$. Using the scheme ([16](#wave:pde2:var:c:scheme:impl))\n",
"at the end point $i=N_x$ with $u_{i+1}^n=u_{i-1}^n$ results in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u^{n+1}_i = - u_i^{n-1} + 2u_i^n + \\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto5\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"\\quad \\left(\\frac{\\Delta t}{\\Delta x}\\right)^2 \\left(\n",
"q_{i+\\frac{1}{2}}(u_{i-1}^n - u_{i}^n) -\n",
"q_{i-\\frac{1}{2}}(u_{i}^n - u_{i-1}^n)\\right)\n",
"+ \\Delta t^2 f^n_i\n",
"\\label{_auto5} \\tag{18}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:var:c:scheme:impl:Neumann0\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"= - u_i^{n-1} + 2u_i^n + \\left(\\frac{\\Delta t}{\\Delta x}\\right)^2\n",
"(q_{i+\\frac{1}{2}} + q_{i-\\frac{1}{2}})(u_{i-1}^n - u_{i}^n) +\n",
"\\Delta t^2 f^n_i\n",
"\\label{wave:pde2:var:c:scheme:impl:Neumann0} \\tag{19}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:var:c:scheme:impl:Neumann\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"\\approx - u_i^{n-1} + 2u_i^n + \\left(\\frac{\\Delta t}{\\Delta x}\\right)^2\n",
"2q_{i}(u_{i-1}^n - u_{i}^n) + \\Delta t^2 f^n_i\n",
"\\thinspace .\n",
"\\label{wave:pde2:var:c:scheme:impl:Neumann} \\tag{20}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we used the approximation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"q_{i+\\frac{1}{2}} + q_{i-\\frac{1}{2}} =\n",
"q_i + \\left(\\frac{dq}{dx}\\right)_i \\Delta x\n",
"+ \\left(\\frac{d^2q}{dx^2}\\right)_i \\Delta x^2 + \\cdots\n",
"+\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\quad q_i - \\left(\\frac{dq}{dx}\\right)_i \\Delta x\n",
"+ \\left(\\frac{d^2q}{dx^2}\\right)_i \\Delta x^2 + \\cdots\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"= 2q_i + 2\\left(\\frac{d^2q}{dx^2}\\right)_i \\Delta x^2 + {\\cal O}(\\Delta x^4)\n",
"\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto6\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"\\approx 2q_i\n",
"\\thinspace .\n",
"\\label{_auto6} \\tag{21}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An alternative derivation may apply the arithmetic mean of\n",
"$q_{n-\\frac{1}{2}}$ and $q_{n+\\frac{1}{2}}$ in\n",
"([19](#wave:pde2:var:c:scheme:impl:Neumann0)), leading to the term"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(q_i + \\frac{1}{2}(q_{i+1}+q_{i-1}))(u_{i-1}^n-u_i^n)\\thinspace .\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since $\\frac{1}{2}(q_{i+1}+q_{i-1}) = q_i + {\\cal O}(\\Delta x^2)$,\n",
"we can approximate with $2q_i(u_{i-1}^n-u_i^n)$ for $i=N_x$ and\n",
"get the same term as we did above.\n",
"\n",
"A common technique when implementing $\\partial u/\\partial x=0$\n",
"boundary conditions, is to assume $dq/dx=0$ as well. This implies\n",
"$q_{i+1}=q_{i-1}$ and $q_{i+1/2}=q_{i-1/2}$ for $i=N_x$.\n",
"The implications for the scheme are"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u^{n+1}_i = - u_i^{n-1} + 2u_i^n + \\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\quad \\left(\\frac{\\Delta t}{\\Delta x}\\right)^2 \\left(\n",
"q_{i+\\frac{1}{2}}(u_{i-1}^n - u_{i}^n) -\n",
"q_{i-\\frac{1}{2}}(u_{i}^n - u_{i-1}^n)\\right)\n",
"+ \\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto7\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
" \\quad \\Delta t^2 f^n_i\n",
"\\label{_auto7} \\tag{22}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:var:c:scheme:impl:Neumann2\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"= - u_i^{n-1} + 2u_i^n + \\left(\\frac{\\Delta t}{\\Delta x}\\right)^2\n",
"2q_{i-\\frac{1}{2}}(u_{i-1}^n - u_{i}^n) +\n",
"\\Delta t^2 f^n_i\n",
"\\thinspace .\n",
"\\label{wave:pde2:var:c:scheme:impl:Neumann2} \\tag{23}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Implementation of variable coefficients\n",
"<div id=\"wave:pde2:var:c:impl\"></div>\n",
"\n",
"The implementation of the scheme with a variable wave velocity $q(x)=c^2(x)$\n",
"may assume that $q$ is available as an array `q[i]` at\n",
"the spatial mesh points. The following loop is a straightforward\n",
"implementation of the scheme ([16](#wave:pde2:var:c:scheme:impl)):"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"for i in range(1, Nx):\n",
" u[i] = - u_nm1[i] + 2*u_n[i] + \\\n",
" C2*(0.5*(q[i] + q[i+1])*(u_n[i+1] - u_n[i]) - \\\n",
" 0.5*(q[i] + q[i-1])*(u_n[i] - u_n[i-1])) + \\\n",
" dt2*f(x[i], t[n])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The coefficient `C2` is now defined as `(dt/dx)**2`, i.e., *not* as the\n",
"squared Courant number, since the wave velocity is variable and appears\n",
"inside the parenthesis.\n",
"\n",
"With Neumann conditions $u_x=0$ at the\n",
"boundary, we need to combine this scheme with the discrete\n",
"version of the boundary condition, as shown in the section [Neumann condition and a variable coefficient](#wave:pde2:var:c:Neumann).\n",
"Nevertheless, it would be convenient to reuse the formula for the\n",
"interior points and just modify the indices `ip1=i+1` and `im1=i-1`\n",
"as we did in the section [Implementation of Neumann conditions](#wave:pde2:Neumann:impl). Assuming\n",
"$dq/dx=0$ at the boundaries, we can implement the scheme at\n",
"the boundary with the following code."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"i = 0\n",
"ip1 = i+1\n",
"im1 = ip1\n",
"u[i] = - u_nm1[i] + 2*u_n[i] + \\\n",
" C2*(0.5*(q[i] + q[ip1])*(u_n[ip1] - u_n[i]) - \\\n",
" 0.5*(q[i] + q[im1])*(u_n[i] - u_n[im1])) + \\\n",
" dt2*f(x[i], t[n])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With ghost cells we can just reuse the formula for the interior\n",
"points also at the boundary, provided that the ghost values of both\n",
"$u$ and $q$ are correctly updated to ensure $u_x=0$ and $q_x=0$.\n",
"\n",
"A vectorized version of the scheme with a variable coefficient\n",
"at internal mesh points becomes"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"u[1:-1] = - u_nm1[1:-1] + 2*u_n[1:-1] + \\\n",
" C2*(0.5*(q[1:-1] + q[2:])*(u_n[2:] - u_n[1:-1]) -\n",
" 0.5*(q[1:-1] + q[:-2])*(u_n[1:-1] - u_n[:-2])) + \\\n",
" dt2*f(x[1:-1], t[n])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## A more general PDE model with variable coefficients\n",
"\n",
"\n",
"Sometimes a wave PDE has a variable coefficient in front of\n",
"the time-derivative term:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:var:c:pde2\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\varrho(x)\\frac{\\partial^2 u}{\\partial t^2} =\n",
"\\frac{\\partial}{\\partial x}\\left( q(x)\n",
"\\frac{\\partial u}{\\partial x}\\right) + f(x,t)\n",
"\\label{wave:pde2:var:c:pde2} \\tag{24}\n",
"\\thinspace .\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One example appears when modeling elastic waves in a rod\n",
"with varying density, cf. ([wave:app:string](#wave:app:string)) with $\\varrho (x)$.\n",
"\n",
"A natural scheme for ([24](#wave:pde2:var:c:pde2)) is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto8\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"[\\varrho D_tD_t u = D_x\\overline{q}^xD_x u + f]^n_i\n",
"\\thinspace .\n",
"\\label{_auto8} \\tag{25}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We realize that the $\\varrho$ coefficient poses no particular\n",
"difficulty, since $\\varrho$ enters the formula just as a simple factor\n",
"in front of a derivative. There is hence no need for any averaging\n",
"of $\\varrho$. Often, $\\varrho$ will be moved to the right-hand side,\n",
"also without any difficulty:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto9\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"[D_tD_t u = \\varrho^{-1}D_x\\overline{q}^xD_x u + f]^n_i\n",
"\\thinspace .\n",
"\\label{_auto9} \\tag{26}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generalization: damping\n",
"\n",
"\n",
"Waves die out by two mechanisms. In 2D and 3D the energy of the wave\n",
"spreads out in space, and energy conservation then requires\n",
"the amplitude to decrease. This effect is not present in 1D.\n",
"Damping is another cause of amplitude reduction. For example,\n",
"the vibrations of a string die out because of damping due to\n",
"air resistance and non-elastic effects in the string.\n",
"\n",
"The simplest way of including damping is to add a first-order derivative\n",
"to the equation (in the same way as friction forces enter a vibrating\n",
"mechanical system):"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde3\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial^2 u}{\\partial t^2} + b\\frac{\\partial u}{\\partial t} =\n",
"c^2\\frac{\\partial^2 u}{\\partial x^2}\n",
" + f(x,t),\n",
"\\label{wave:pde3} \\tag{27}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $b \\geq 0$ is a prescribed damping coefficient.\n",
"\n",
"A typical discretization of ([27](#wave:pde3)) in terms of centered\n",
"differences reads"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde3:fd\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"[D_tD_t u + bD_{2t}u = c^2D_xD_x u + f]^n_i\n",
"\\thinspace .\n",
"\\label{wave:pde3:fd} \\tag{28}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Writing out the equation and solving for the unknown $u^{n+1}_i$\n",
"gives the scheme"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde3:fd2\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u^{n+1}_i = (1 + {\\frac{1}{2}}b\\Delta t)^{-1}(({\\frac{1}{2}}b\\Delta t -1)\n",
"u^{n-1}_i + 2u^n_i + C^2\n",
"\\left(u^{n}_{i+1}-2u^{n}_{i} + u^{n}_{i-1}\\right) + \\Delta t^2 f^n_i),\n",
"\\label{wave:pde3:fd2} \\tag{29}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"for $i\\in\\mathcal{I}_x^i$ and $n\\geq 1$.\n",
"New equations must be derived for $u^1_i$, and for boundary points in case\n",
"of Neumann conditions.\n",
"\n",
"The damping is very small in many wave phenomena and thus only evident\n",
"for very long time simulations. This makes the standard wave equation\n",
"without damping relevant for a lot of applications.\n",
"\n",
"\n",
"# Building a general 1D wave equation solver\n",
"<div id=\"wave:pde2:software\"></div>\n",
"\n",
"\n",
"The program [`wave1D_dn_vc.py`](${src_wave}/wave1D/wave1D_dn_vc.py)\n",
"is a fairly general code for 1D wave propagation problems that\n",
"targets the following initial-boundary value problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:software:ueq\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u_{tt} = (c^2(x)u_x)_x + f(x,t),\\quad x\\in (0,L),\\ t\\in (0,T]\n",
"\\label{wave:pde2:software:ueq} \\tag{30}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto10\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u(x,0) = I(x),\\quad x\\in [0,L]\n",
"\\label{_auto10} \\tag{31}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto11\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u_t(x,0) = V(t),\\quad x\\in [0,L]\n",
"\\label{_auto11} \\tag{32}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto12\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u(0,t) = U_0(t)\\hbox{ or } u_x(0,t)=0,\\quad t\\in (0,T]\n",
"\\label{_auto12} \\tag{33}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:pde2:software:bcL\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u(L,t) = U_L(t)\\hbox{ or } u_x(L,t)=0,\\quad t\\in (0,T]\n",
"\\label{wave:pde2:software:bcL} \\tag{34}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The only new feature here is the time-dependent Dirichlet conditions, but\n",
"they are trivial to implement:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"i = mathcal{I}_x[0] # x=0\n",
"u[i] = U_0(t[n+1])\n",
"\n",
"i = mathcal{I}_x[-1] # x=L\n",
"u[i] = U_L(t[n+1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `solver` function is a natural extension of the simplest\n",
"`solver` function in the initial `wave1D_u0.py` program,\n",
"extended with Neumann boundary conditions ($u_x=0$),\n",
"time-varying Dirichlet conditions, as well as\n",
"a variable wave velocity. The different code segments needed\n",
"to make these extensions have been shown and commented upon in the\n",
"preceding text. We refer to the `solver` function in the\n",
"`wave1D_dn_vc.py` file for all the details. Note in that\n",
" `solver` function, however, that the technique of \"hashing\" is\n",
"used to check whether a certain simulation has been run before, or not.\n",
"% if BOOK == 'book':\n",
"This technique is further explained in the section [softeng2:wave1D:filestorage:hash](#softeng2:wave1D:filestorage:hash).\n",
"% endif\n",
"\n",
"The vectorization is only applied inside the time loop, not for the\n",
"initial condition or the first time steps, since this initial work\n",
"is negligible for long time simulations in 1D problems.\n",
"\n",
"The following sections explain various more advanced programming\n",
"techniques applied in the general 1D wave equation solver.\n",
"\n",
"## User action function as a class\n",
"\n",
"A useful feature in the `wave1D_dn_vc.py` program is the specification\n",
"of the `user_action` function as a class. This part of the program may\n",
"need some motivation and explanation. Although the `plot_u_st`\n",
"function (and the `PlotMatplotlib` class) in the `wave1D_u0.viz`\n",
"function remembers the local variables in the `viz` function, it is a\n",
"cleaner solution to store the needed variables together with the\n",
"function, which is exactly what a class offers.\n",
"\n",
"### The code\n",
"\n",
"A class for flexible plotting, cleaning up files, making movie\n",
"files, like the function `wave1D_u0.viz` did, can be coded as follows:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"class PlotAndStoreSolution:\n",
" \"\"\"\n",
" Class for the user_action function in solver.\n",
" Visualizes the solution only.\n",
" \"\"\"\n",
" def __init__(\n",
" self,\n",
" casename='tmp', # Prefix in filenames\n",
" umin=-1, umax=1, # Fixed range of y axis\n",
" pause_between_frames=None, # Movie speed\n",
" backend='matplotlib', # or 'gnuplot' or None\n",
" screen_movie=True, # Show movie on screen?\n",
" title='', # Extra message in title\n",
" skip_frame=1, # Skip every skip_frame frame\n",
" filename=None): # Name of file with solutions\n",
" self.casename = casename\n",
" self.yaxis = [umin, umax]\n",
" self.pause = pause_between_frames\n",
" self.backend = backend\n",
" if backend is None:\n",
" # Use native matplotlib\n",
" import matplotlib.pyplot as plt\n",
" elif backend in ('matplotlib', 'gnuplot'):\n",
" module = 'scitools.easyviz.' + backend + '_'\n",
" exec('import %s as plt' % module)\n",
" self.plt = plt\n",
" self.screen_movie = screen_movie\n",
" self.title = title\n",
" self.skip_frame = skip_frame\n",
" self.filename = filename\n",
" if filename is not None:\n",
" # Store time points when u is written to file\n",
" self.t = []\n",
" filenames = glob.glob('.' + self.filename + '*.dat.npz')\n",
" for filename in filenames:\n",
" os.remove(filename)\n",
"\n",
" # Clean up old movie frames\n",
" for filename in glob.glob('frame_*.png'):\n",
" os.remove(filename)\n",
"\n",
" def __call__(self, u, x, t, n):\n",
" \"\"\"\n",
" Callback function user_action, call by solver:\n",
" Store solution, plot on screen and save to file.\n",
" \"\"\"\n",
" # Save solution u to a file using numpy.savez\n",
" if self.filename is not None:\n",
" name = 'u%04d' % n # array name\n",
" kwargs = {name: u}\n",
" fname = '.' + self.filename + '_' + name + '.dat'\n",
" np.savez(fname, **kwargs)\n",
" self.t.append(t[n]) # store corresponding time value\n",
" if n == 0: # save x once\n",
" np.savez('.' + self.filename + '_x.dat', x=x)\n",
"\n",
" # Animate\n",
" if n % self.skip_frame != 0:\n",
" return\n",
" title = 't=%.3f' % t[n]\n",
" if self.title:\n",
" title = self.title + ' ' + title\n",
" if self.backend is None:\n",
" # native matplotlib animation\n",
" if n == 0:\n",
" self.plt.ion()\n",
" self.lines = self.plt.plot(x, u, 'r-')\n",
" self.plt.axis([x[0], x[-1],\n",
" self.yaxis[0], self.yaxis[1]])\n",
" self.plt.xlabel('x')\n",
" self.plt.ylabel('u')\n",
" self.plt.title(title)\n",
" self.plt.legend(['t=%.3f' % t[n]])\n",
" else:\n",
" # Update new solution\n",
" self.lines[0].set_ydata(u)\n",
" self.plt.legend(['t=%.3f' % t[n]])\n",
" self.plt.draw()\n",
" else:\n",
" # scitools.easyviz animation\n",
" self.plt.plot(x, u, 'r-',\n",
" xlabel='x', ylabel='u',\n",
" axis=[x[0], x[-1],\n",
" self.yaxis[0], self.yaxis[1]],\n",
" title=title,\n",
" show=self.screen_movie)\n",
" # pause\n",
" if t[n] == 0:\n",
" time.sleep(2) # let initial condition stay 2 s\n",
" else:\n",
" if self.pause is None:\n",
" pause = 0.2 if u.size < 100 else 0\n",
" time.sleep(pause)\n",
"\n",
" self.plt.savefig('frame_%04d.png' % (n))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Dissection\n",
"\n",
"Understanding this class requires quite some familiarity with Python\n",
"in general and class programming in particular.\n",
"The class supports plotting with Matplotlib (`backend=None`) or\n",
"SciTools (`backend=matplotlib` or `backend=gnuplot`) for maximum\n",
"flexibility.\n",
"\n",
"<!-- Since all the plot frames are to be collected in a separate subdirectory, -->\n",
"<!-- we demand a (logical) \"casename\" from the user that is used as -->\n",
"<!-- subdirectory name in the `make_movie_file` method. The statements -->\n",
"<!-- in this method perform actions normally done in the operating -->\n",
"<!-- system, but the Python interface via `shutil.rmtree`, `os.mkdir`, -->\n",
"<!-- `os.chdir`, etc., works on all platforms where Python works. -->\n",
"\n",
"The constructor shows how we can flexibly import the plotting engine\n",
"as (typically) `scitools.easyviz.gnuplot_` or\n",
"`scitools.easyviz.matplotlib_` (note the trailing underscore - it is required).\n",
"With the `screen_movie` parameter\n",
"we can suppress displaying each movie frame on the screen.\n",
"Alternatively, for slow movies associated with\n",
"fine meshes, one can set\n",
"`skip_frame=10`, causing every 10 frames to be shown.\n",
"\n",
"The `__call__` method makes `PlotAndStoreSolution` instances behave like\n",
"functions, so we can just pass an instance, say `p`, as the\n",
"`user_action` argument in the `solver` function, and any call to\n",
"`user_action` will be a call to `p.__call__`. The `__call__`\n",
"method plots the solution on the screen,\n",
"saves the plot to file, and stores the solution in a file for\n",
"later retrieval.\n",
"\n",
"More details on storing the solution in files appear in\n",
"in\n",
"the document\n",
"[Scientific software engineering; wave equation case](http://tinyurl.com/k3sdbuv/pub/softeng2)\n",
"[[Langtangen_deqbook_softeng2]](#Langtangen_deqbook_softeng2).\n",
"\n",
"## Pulse propagation in two media\n",
"\n",
"\n",
"The function `pulse` in `wave1D_dn_vc.py` demonstrates wave motion in\n",
"heterogeneous media where $c$ varies. One can specify an interval\n",
"where the wave velocity is decreased by a factor `slowness_factor`\n",
"(or increased by making this factor less than one).\n",
"[Figure](#wave:pde1:fig:pulse1:two:media) shows a typical simulation\n",
"scenario.\n",
"\n",
"Four types of initial conditions are available:\n",
"\n",
"1. a rectangular pulse (`plug`),\n",
"\n",
"2. a Gaussian function (`gaussian`),\n",
"\n",
"3. a \"cosine hat\" consisting of one period of the cosine function\n",
" (`cosinehat`),\n",
"\n",
"4. frac{1}{2} a period of a \"cosine hat\" (`frac{1}{2}-cosinehat`)\n",
"\n",
"These peak-shaped initial conditions can be placed in the middle\n",
"(`loc='center'`) or at the left end (`loc='left'`) of the domain.\n",
"With the pulse in the middle, it splits in two parts, each with frac{1}{2}\n",
"the initial amplitude, traveling in opposite directions. With the\n",
"pulse at the left end, centered at $x=0$, and using the symmetry\n",
"condition $\\partial u/\\partial x=0$, only a right-going pulse is\n",
"generated. There is also a left-going pulse, but it travels from $x=0$\n",
"in negative $x$ direction and is not visible in the domain $[0,L]$.\n",
"\n",
"The `pulse` function is a flexible tool for playing around with\n",
"various wave shapes and jumps in the wave velocity (i.e.,\n",
"discontinuous media). The code is shown to demonstrate how easy it is\n",
"to reach this flexibility with the building blocks we have already\n",
"developed:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"def pulse(\n",
" C=1, # Maximum Courant number\n",
" Nx=200, # spatial resolution\n",
" animate=True,\n",
" version='vectorized',\n",
" T=2, # end time\n",
" loc='left', # location of initial condition\n",
" pulse_thinspace .='gaussian', # pulse/init.cond. type\n",
" slowness_factor=2, # inverse of wave vel. in right medium\n",
" medium=[0.7, 0.9], # interval for right medium\n",
" skip_frame=1, # skip frames in animations\n",
" sigma=0.05 # width measure of the pulse\n",
" ):\n",
" \"\"\"\n",
" Various peaked-shaped initial conditions on [0,1].\n",
" Wave velocity is decreased by the slowness_factor inside\n",
" medium. The loc parameter can be 'center' or 'left',\n",
" depending on where the initial pulse is to be located.\n",
" The sigma parameter governs the width of the pulse.\n",
" \"\"\"\n",
" # Use scaled parameters: L=1 for domain length, c_0=1\n",
" # for wave velocity outside the domain.\n",
" L = 1.0\n",
" c_0 = 1.0\n",
" if loc == 'center':\n",
" xc = L/2\n",
" elif loc == 'left':\n",
" xc = 0\n",
"\n",
" if pulse_thinspace . in ('gaussian','Gaussian'):\n",
" def I(x):\n",
" return np.exp(-0.5*((x-xc)/sigma)**2)\n",
" elif pulse_thinspace . == 'plug':\n",
" def I(x):\n",
" return 0 if abs(x-xc) > sigma else 1\n",
" elif pulse_thinspace . == 'cosinehat':\n",
" def I(x):\n",
" # One period of a cosine\n",
" w = 2\n",
" a = w*sigma\n",
" return 0.5*(1 + np.cos(np.pi*(x-xc)/a)) \\\n",
" if xc - a <= x <= xc + a else 0\n",
"\n",
" elif pulse_thinspace . == 'frac{1}{2}-cosinehat':\n",
" def I(x):\n",
" # Half a period of a cosine\n",
" w = 4\n",
" a = w*sigma\n",
" return np.cos(np.pi*(x-xc)/a) \\\n",
" if xc - 0.5*a <= x <= xc + 0.5*a else 0\n",
" else:\n",
" raise ValueError('Wrong pulse_thinspace .=\"%s\"' % pulse_thinspace .)\n",
"\n",
" def c(x):\n",
" return c_0/slowness_factor \\\n",
" if medium[0] <= x <= medium[1] else c_0\n",
"\n",
" umin=-0.5; umax=1.5*I(xc)\n",
" casename = '%s_Nx%s_sf%s' % \\\n",
" (pulse_thinspace ., Nx, slowness_factor)\n",
" action = PlotMediumAndSolution(\n",
" medium, casename=casename, umin=umin, umax=umax,\n",
" skip_frame=skip_frame, screen_movie=animate,\n",
" backend=None, filename='tmpdata')\n",
"\n",
" # Choose the stability limit with given Nx, worst case c\n",
" # (lower C will then use this dt, but smaller Nx)\n",
" dt = (L/Nx)/c_0\n",
" cpu, hashed_input = solver(\n",
" I=I, V=None, f=None, c=c,\n",
" U_0=None, U_L=None,\n",
" L=L, dt=dt, C=C, T=T,\n",
" user_action=action,\n",
" version=version,\n",
" stability_safety_factor=1)\n",
"\n",
" if cpu > 0: # did we generate new data?\n",
" action.close_file(hashed_input)\n",
" action.make_movie_file()\n",
" print 'cpu (-1 means no new data generated):', cpu\n",
"\n",
"def convergence_rates(\n",
" u_exact,\n",
" I, V, f, c, U_0, U_L, L,\n",
" dt0, num_meshes,\n",
" C, T, version='scalar',\n",
" stability_safety_factor=1.0):\n",
" \"\"\"\n",
" Half the time step and estimate convergence rates for\n",
" for num_meshes simulations.\n",
" \"\"\"\n",
" class ComputeError:\n",
" def __init__(self, norm_type):\n",
" self.error = 0\n",
"\n",
" def __call__(self, u, x, t, n):\n",
" \"\"\"Store norm of the error in self.E.\"\"\"\n",
" error = np.abs(u - u_exact(x, t[n])).max()\n",
" self.error = max(self.error, error)\n",
"\n",
" E = []\n",
" h = [] # dt, solver adjusts dx such that C=dt*c/dx\n",
" dt = dt0\n",
" for i in range(num_meshes):\n",
" error_calculator = ComputeError('Linf')\n",
" solver(I, V, f, c, U_0, U_L, L, dt, C, T,\n",
" user_action=error_calculator,\n",
" version='scalar',\n",
" stability_safety_factor=1.0)\n",
" E.append(error_calculator.error)\n",
" h.append(dt)\n",
" dt /= 2 # halve the time step for next simulation\n",
" print 'E:', E\n",
" print 'h:', h\n",
" r = [np.log(E[i]/E[i-1])/np.log(h[i]/h[i-1])\n",
" for i in range(1,num_meshes)]\n",
" return r\n",
"\n",
"def test_convrate_sincos():\n",
" n = m = 2\n",
" L = 1.0\n",
" u_exact = lambda x, t: np.cos(m*np.pi/L*t)*np.sin(m*np.pi/L*x)\n",
"\n",
" r = convergence_rates(\n",
" u_exact=u_exact,\n",
" I=lambda x: u_exact(x, 0),\n",
" V=lambda x: 0,\n",
" f=0,\n",
" c=1,\n",
" U_0=0,\n",
" U_L=0,\n",
" L=L,\n",
" dt0=0.1,\n",
" num_meshes=6,\n",
" C=0.9,\n",
" T=1,\n",
" version='scalar',\n",
" stability_safety_factor=1.0)\n",
" print 'rates sin(x)*cos(t) solution:', \\\n",
" [round(r_,2) for r_ in r]\n",
" assert abs(r[-1] - 2) < 0.002"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `PlotMediumAndSolution` class used here is a subclass of\n",
"`PlotAndStoreSolution` where the medium with reduced $c$ value,\n",
"as specified by the `medium` interval,\n",
"is visualized in the plots.\n",
"\n",
"**Comment on the choices of discretization parameters.**\n",
"\n",
"The argument $N_x$ in the `pulse` function does not correspond to\n",
"the actual spatial resolution of $C<1$, since the `solver`\n",
"function takes a fixed $\\Delta t$ and $C$, and adjusts $\\Delta x$\n",
"accordingly. As seen in the `pulse` function,\n",
"the specified $\\Delta t$ is chosen according to the\n",
"limit $C=1$, so if $C<1$, $\\Delta t$ remains the same, but the\n",
"`solver` function operates with a larger $\\Delta x$ and smaller\n",
"$N_x$ than was specified in the call to `pulse`. The practical reason\n",
"is that we always want to keep $\\Delta t$ fixed such that\n",
"plot frames and movies are synchronized in time regardless of the\n",
"value of $C$ (i.e., $\\Delta x$ is varied when the\n",
"Courant number varies).\n",
"\n",
"\n",
"\n",
"The reader is encouraged to play around with the `pulse` function:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To easily kill the graphics by Ctrl-C and restart a new simulation it might be\n",
"easier to run the above two statements from the command line\n",
"with"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" Terminal> python -c 'import wave1D_dn_vc as w; w.pulse(...)'\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Exercises\n",
"\n",
"\n",
"\n",
"<!-- --- begin exercise --- -->\n",
"\n",
"## Exercise 1: Find the analytical solution to a damped wave equation\n",
"<div id=\"wave:exer:standingwave:damped:uex\"></div>\n",
"\n",
"Consider the wave equation with damping ([27](#wave:pde3)).\n",
"The goal is to find an exact solution to a wave problem with damping and zero source term.\n",
"A starting point is the standing wave solution from\n",
"[wave:exer:standingwave](#wave:exer:standingwave). mathcal{I}_t becomes necessary to\n",
"include a damping term $e^{-\\beta t}$ and also have both a sine and cosine\n",
"component in time:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\uex(x,t) = e^{-\\beta t}\n",
"\\sin kx \\left( A\\cos\\omega t\n",
"+ B\\sin\\omega t\\right)\n",
"\\thinspace .\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Find $k$ from the boundary conditions\n",
"$u(0,t)=u(L,t)=0$. Then use the PDE to find constraints on\n",
"$\\beta$, $\\omega$, $A$, and $B$.\n",
"Set up a complete initial-boundary value problem\n",
"and its solution.\n",
"\n",
"\n",
"<!-- --- begin solution of exercise --- -->\n",
"**Solution.**\n",
"Mathematical model:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{\\partial^2 u}{\\partial t^2} + b\\frac{\\partial u}{\\partial t} =\n",
"c^2\\frac{\\partial^2 u}{\\partial x^2},\n",
"\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$b \\geq 0$ is a prescribed damping coefficient.\n",
"\n",
"Ansatz:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(x,t) = e^{-\\beta t}\n",
"\\sin kx \\left( A\\cos\\omega t\n",
"+ B\\sin\\omega t\\right)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Boundary condition: $u=0$ for $x=0,L$. Fulfilled for $x=0$. Requirement\n",
"at $x=L$ gives"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"kL = m\\pi,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"for an arbitrary integer $m$. Hence, $k=m\\pi/L$.\n",
"\n",
"Inserting the ansatz in the PDE and dividing by $e^{-\\beta t}$ results in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"(\\beta^2 sin kx -\\omega^2 sin kx - b\\beta sin kx) (A\\cos\\omega t + B\\sin\\omega t) &+ \\nonumber \\\\ \n",
"(b\\omega sin kx - 2\\beta\\omega sin kx) (-A\\sin\\omega t + B\\cos\\omega t) &= -(A\\cos\\omega t + B\\sin\\omega t)k^2c^2 \\nonumber\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This gives us two requirements:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\beta^2 - \\omega^2 + b\\beta + k^2c^2 = 0\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"-2\\beta\\omega + b\\omega = 0\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since $b$, $c$ and $k$ are to be given in advance, we may solve these two equations to get"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\beta &= \\frac{b}{2} \\nonumber \\\\ \n",
"\\omega &= \\sqrt{c^2k^2 - \\frac{b^2}{4}} \\nonumber\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From the initial condition on the derivative, i.e. $\\frac{\\partial u_e}{\\partial t} = 0$, we find that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"B\\omega = \\beta A\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Inserting the expression for $\\omega$, we find that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"B = \\frac{b}{2\\sqrt{c^2k^2 - \\frac{b^2}{4}}} A\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"for $A$ prescribed.\n",
"\n",
"Using $t = 0$ in the expression for $u_e$ gives us the initial condition as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"I(x) = A sin kx\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Summarizing, the PDE problem can then be states as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{\\partial^2 u}{\\partial t^2} + b\\frac{\\partial u}{\\partial t} =\n",
"c^2 \\frac{\\partial^2 u}{\\partial x^2}, \\quad x\\in (0,L),\\ t\\in (0,T]\n",
"\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(x,0) = I(x), \\quad x\\in [0,L]\n",
"\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{\\partial}{\\partial t}u(x,0) = 0, \\quad x\\in [0,L]\n",
"\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(0,t) = 0, \\quad t\\in (0,T]\n",
"\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(L,t) = 0, \\quad t\\in (0,T]\n",
"\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where constants $c$, $A$, $b$ and $k$, as well as $I(x)$, are prescribed.\n",
"\n",
"The solution to the problem is then given as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\uex(x,t) = e^{-\\beta t}\n",
"\\sin kx \\left( A\\cos\\omega t\n",
"+ B\\sin\\omega t\\right)\n",
"\\thinspace .\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"with $k=m\\pi/L$ for arbitrary integer $m$, $\\beta = \\frac{b}{2}$,\n",
"$\\omega = \\sqrt{c^2k^2 - \\frac{b^2}{4}}$, $B = \\frac{b}{2\\sqrt{c^2k^2 - \\frac{b^2}{4}}} A$\n",
"and $I(x) = A sin kx$.\n",
"\n",
"<!-- --- end solution of exercise --- -->\n",
"Filename: `damped_waves`.\n",
"\n",
"<!-- --- end exercise --- -->\n",
"\n",
"\n",
"\n",
"\n",
"<!-- --- begin exercise --- -->\n",
"\n",
"## Problem 2: Explore symmetry boundary conditions\n",
"<div id=\"wave:exer:symmetry:bc\"></div>\n",
"\n",
"Consider the simple \"plug\" wave where $\\Omega = [-L,L]$ and"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"I(x) = \\left\\lbrace\\begin{array}{ll}\n",
"1, & x\\in [-\\delta, \\delta],\\\\ \n",
"0, & \\hbox{otherwise}\n",
"\\end{array}\\right.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"for some number $0 < \\delta < L$. The other initial condition is\n",
"$u_t(x,0)=0$ and there is no source term $f$.\n",
"The boundary conditions can be set to $u=0$.\n",
"The solution to this problem is symmetric around $x=0$.\n",
"This means that we can simulate the wave process in only frac{1}{2}\n",
"of the domain $[0,L]$.\n",
"\n",
"\n",
"**a)**\n",
"Argue why the symmetry boundary condition\n",
"is $u_x=0$ at $x=0$.\n",
"\n",
"<!-- --- begin hint in exercise --- -->\n",
"\n",
"**Hint.**\n",
"Symmetry of a function about $x=x_0$ means that\n",
"$f(x_0+h) = f(x_0-h)$.\n",
"\n",
"<!-- --- end hint in exercise --- -->\n",
"\n",
"\n",
"<!-- --- begin solution of exercise --- -->\n",
"**Solution.**\n",
"A symmetric $u$ around $x=0$ means that $u(-x,t)=u(x,t)$.\n",
"Let $x_0=0$ and $x=x_0+h$. Then we can use a *centered* finite difference\n",
"definition of the derivative:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{\\partial}{\\partial x}u(x_0,t) =\n",
"\\lim_{h\\rightarrow 0}\\frac{u(x_0+h,t)- u(x_0-h)}{2h} =\n",
"\\lim_{h\\rightarrow 0}\\frac{u(h,t)- u(-h,t)}{2h} = 0,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"since $u(h,t)=u(-h,t)$ for any $h$. Symmetry around a point $x=x_0$\n",
"therefore always implies $u_x(x_0,t)=0$.\n",
"\n",
"<!-- --- end solution of exercise --- -->\n",
"\n",
"**b)**\n",
"Perform simulations of the complete wave problem\n",
"on $[-L,L]$. Thereafter, utilize the\n",
"symmetry of the solution and run a simulation\n",
"in frac{1}{2} of the domain $[0,L]$, using a boundary condition\n",
"at $x=0$. Compare plots from the two solutions and\n",
"confirm that they are the same.\n",
"\n",
"\n",
"<!-- --- begin solution of exercise --- -->\n",
"**Solution.**\n",
"We can utilize the `wave1D_dn.py` code which allows Dirichlet and\n",
"Neumann conditions. The `solver` and `viz` functions must take $x_0$\n",
"and $x_L$ as parameters instead of just $L$ such that we can solve the\n",
"wave equation in $[x_0, x_L]$. The we can call up `solver` for the two\n",
"problems on $[-L,L]$ and $[0,L]$ with boundary conditions\n",
"$u(-L,t)=u(L,t)=0$ and $u_x(0,t)=u(L,t)=0$, respectively.\n",
"\n",
"The original `wave1D_dn.py` code makes a movie by playing all the\n",
"`.png` files in a browser. mathcal{I}_t can then be wise to let the `viz`\n",
"function create a movie directory and place all the frames and HTML\n",
"player file in that directory. Alternatively, one can just make\n",
"some ordinary movie file (Ogg, WebM, MP4, Flash) with `ffmpeg` or\n",
"`ffmpeg` and give it a name. mathcal{I}_t is a point that the name is\n",
"transferred to `viz` so it is easy to call `viz` twice and get two\n",
"separate movie files or movie directories.\n",
"\n",
"The plots produced by the code (below) shows that the solutions indeed\n",
"are the same.\n",
"\n",
"<!-- --- end solution of exercise --- -->\n",
"\n",
"**c)**\n",
"Prove the symmetry property of the solution\n",
"by setting up the complete initial-boundary value problem\n",
"and showing that if $u(x,t)$ is a solution, then also $u(-x,t)$\n",
"is a solution.\n",
"\n",
"\n",
"<!-- --- begin solution of exercise --- -->\n",
"**Solution.**\n",
"The plan in this proof is to introduce $v(x,t)=u(-x,t)$\n",
"and show that $v$ fulfills the same\n",
"initial-boundary value problem as $u$. If the problem has a unique\n",
"solution, then $v=u$. Or, in other words, the solution is\n",
"symmetric: $u(-x,t)=u(x,t)$.\n",
"\n",
"We can work with a general initial-boundary value problem on the form"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto13\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u_tt(x,t) = c^2u_{xx}(x,t) + f(x,t)\n",
"\\label{_auto13} \\tag{35}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto14\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u(x,0) = I(x)\n",
"\\label{_auto14} \\tag{36}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto15\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u_t(x,0) = V(x)\n",
"\\label{_auto15} \\tag{37}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto16\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u(-L,0) = 0\n",
"\\label{_auto16} \\tag{38}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto17\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u(L,0) = 0\n",
"\\label{_auto17} \\tag{39}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Introduce a new coordinate $\\bar x = -x$. We have that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{\\partial^2 u}{\\partial x^2} = \\frac{\\partial}{\\partial x}\n",
"\\left(\n",
"\\frac{\\partial u}{\\partial\\bar x}\n",
"\\frac{\\partial\\bar x}{\\partial x}\n",
"\\right)\n",
"= \\frac{\\partial}{\\partial x}\n",
"\\left(\n",
"\\frac{\\partial u}{\\partial\\bar x} (-1)\\right)\n",
"= (-1)^2 \\frac{\\partial^2 u}{\\partial \\bar x^2}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The derivatives in time are unchanged.\n",
"\n",
"Substituting $x$ by $-\\bar x$ leads to"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto18\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u_{tt}(-\\bar x,t) = c^2u_{\\bar x\\bar x}(-\\bar x,t) + f(-\\bar x,t)\n",
"\\label{_auto18} \\tag{40}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto19\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u(-\\bar x,0) = I(-\\bar x)\n",
"\\label{_auto19} \\tag{41}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto20\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u_t(-\\bar x,0) = V(-\\bar x)\n",
"\\label{_auto20} \\tag{42}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto21\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u(L,0) = 0\n",
"\\label{_auto21} \\tag{43}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto22\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u(-L,0) = 0\n",
"\\label{_auto22} \\tag{44}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, dropping the bars and introducing $v(x,t)=u(-x,t)$, we find that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto23\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"v_{tt}(x,t) = c^2v_{xx}(x,t) + f(-x,t)\n",
"\\label{_auto23} \\tag{45}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto24\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"v(x,0) = I(-x)\n",
"\\label{_auto24} \\tag{46}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto25\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"v_t(x ,0) = V(-x)\n",
"\\label{_auto25} \\tag{47}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto26\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"v(-L,0) = 0\n",
"\\label{_auto26} \\tag{48}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto27\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"v(L,0) = 0\n",
"\\label{_auto27} \\tag{49}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Provided that $I$, $f$, and $V$ are all symmetric* around $x=0$\n",
"such that $I(x)=I(-x)$, $V(x)=V(-x)$, and $f(x,t)=f(-x,t)$, we\n",
"can express the initial-boundary value problem as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto28\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"v_{tt}(x,t) = c^2v_{xx}(x,t) + f(x,t)\n",
"\\label{_auto28} \\tag{50}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto29\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"v(x,0) = I(x)\n",
"\\label{_auto29} \\tag{51}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto30\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"v_t(x ,0) = V(x)\n",
"\\label{_auto30} \\tag{52}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto31\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"v(-L,0) = 0\n",
"\\label{_auto31} \\tag{53}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto32\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"v(L,0) = 0\n",
"\\label{_auto32} \\tag{54}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the same problem as the one that $u$ fulfills. If the solution\n",
"is unique, which can be proven, then $v=u$, and $u(-x,t)=u(x,t)$.\n",
"\n",
"To summarize, the necessary conditions for symmetry are that\n",
"\n",
" * all involved functions $I$, $V$, and $f$ must be symmetric, and\n",
"\n",
" * the boundary conditions are symmetric in the sense that they\n",
" can be flipped (the condition at $x=-L$ can be applied\n",
" at $x=L$ and vice versa).\n",
"\n",
"<!-- --- end solution of exercise --- -->\n",
"\n",
"**d)**\n",
"If the code works correctly, the solution $u(x,t) = x(L-x)(1+\\frac{t}{2})$\n",
"should be reproduced exactly. Write a test function `test_quadratic` that\n",
"checks whether this is the case. Simulate for $x$ in $[0, \\frac{L}{2}]$ with\n",
"a symmetry condition at the end $x = \\frac{L}{2}$.\n",
"\n",
"\n",
"<!-- --- begin solution of exercise --- -->\n",
"**Solution.**\n",
"Running the code below, shows that the test case indeed is reproduced exactly."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"#!/usr/bin/env python\n",
"from scitools.std import *\n",
"\n",
"# Add an x0 coordinate for solving the wave equation on [x0, xL]\n",
"\n",
"def solver(I, V, f, c, U_0, U_L, x0, xL, Nx, C, T,\n",
" user_action=None, version='scalar'):\n",
" \"\"\"\n",
" Solve u_tt=c^2*u_xx + f on (0,L)x(0,T].\n",
" u(0,t)=U_0(t) or du/dn=0 (U_0=None), u(L,t)=U_L(t) or du/dn=0 (u_L=None).\n",
" \"\"\"\n",
" x = linspace(x0, xL, Nx+1) # Mesh points in space\n",
" dx = x[1] - x[0]\n",
" dt = C*dx/c\n",
" Nt = int(round(T/dt))\n",
" t = linspace(0, Nt*dt, Nt+1) # Mesh points in time\n",
" C2 = C**2; dt2 = dt*dt # Help variables in the scheme\n",
"\n",
" # Wrap user-given f, V, U_0, U_L\n",
" if f is None or f == 0:\n",
" f = (lambda x, t: 0) if version == 'scalar' else \\\n",
" lambda x, t: zeros(x.shape)\n",
" if V is None or V == 0:\n",
" V = (lambda x: 0) if version == 'scalar' else \\\n",
" lambda x: zeros(x.shape)\n",
" if U_0 is not None:\n",
" if isinstance(U_0, (float,int)) and U_0 == 0:\n",
" U_0 = lambda t: 0\n",
" if U_L is not None:\n",
" if isinstance(U_L, (float,int)) and U_L == 0:\n",
" U_L = lambda t: 0\n",
"\n",
" u = zeros(Nx+1) # Solution array at new time level\n",
" u_1 = zeros(Nx+1) # Solution at 1 time level back\n",
" u_2 = zeros(Nx+1) # Solution at 2 time levels back\n",
"\n",
" mathcal{I}_x = range(0, Nx+1)\n",
" mathcal{I}_t = range(0, Nt+1)\n",
"\n",
" import time; t0 = time.clock() # CPU time measurement\n",
"\n",
" # Load initial condition into u_1\n",
" for i in mathcal{I}_x:\n",
" u_1[i] = I(x[i])\n",
"\n",
" if user_action is not None:\n",
" user_action(u_1, x, t, 0)\n",
"\n",
" # Special formula for the first step\n",
" for i in mathcal{I}_x[1:-1]:\n",
" u[i] = u_1[i] + dt*V(x[i]) + \\\n",
" 0.5*C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1]) + \\\n",
" 0.5*dt2*f(x[i], t[0])\n",
"\n",
" i = mathcal{I}_x[0]\n",
" if U_0 is None:\n",
" # Set boundary values du/dn = 0\n",
" # x=0: i-1 -> i+1 since u[i-1]=u[i+1]\n",
" # x=L: i+1 -> i-1 since u[i+1]=u[i-1])\n",
" ip1 = i+1\n",
" im1 = ip1 # i-1 -> i+1\n",
" u[i] = u_1[i] + dt*V(x[i]) + \\\n",
" 0.5*C2*(u_1[im1] - 2*u_1[i] + u_1[ip1]) + \\\n",
" 0.5*dt2*f(x[i], t[0])\n",
" else:\n",
" u[0] = U_0(dt)\n",
"\n",
" i = mathcal{I}_x[-1]\n",
" if U_L is None:\n",
" im1 = i-1\n",
" ip1 = im1 # i+1 -> i-1\n",
" u[i] = u_1[i] + dt*V(x[i]) + \\\n",
" 0.5*C2*(u_1[im1] - 2*u_1[i] + u_1[ip1]) + \\\n",
" 0.5*dt2*f(x[i], t[0])\n",
" else:\n",
" u[i] = U_L(dt)\n",
"\n",
" if user_action is not None:\n",
" user_action(u, x, t, 1)\n",
"\n",
" # Update data structures for next step\n",
" u_2[:], u_1[:] = u_1, u\n",
"\n",
" for n in mathcal{I}_t[1:-1]:\n",
" # Update all inner points\n",
" if version == 'scalar':\n",
" for i in mathcal{I}_x[1:-1]:\n",
" u[i] = - u_2[i] + 2*u_1[i] + \\\n",
" C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1]) + \\\n",
" dt2*f(x[i], t[n])\n",
"\n",
" elif version == 'vectorized':\n",
" u[1:-1] = - u_2[1:-1] + 2*u_1[1:-1] + \\\n",
" C2*(u_1[0:-2] - 2*u_1[1:-1] + u_1[2:]) + \\\n",
" dt2*f(x[1:-1], t[n])\n",
" else:\n",
" raise ValueError('version=%s' % version)\n",
"\n",
" # Insert boundary conditions\n",
" i = mathcal{I}_x[0]\n",
" if U_0 is None:\n",
" # Set boundary values\n",
" # x=0: i-1 -> i+1 since u[i-1]=u[i+1] when du/dn=0\n",
" # x=L: i+1 -> i-1 since u[i+1]=u[i-1] when du/dn=0\n",
" ip1 = i+1\n",
" im1 = ip1\n",
" u[i] = - u_2[i] + 2*u_1[i] + \\\n",
" C2*(u_1[im1] - 2*u_1[i] + u_1[ip1]) + \\\n",
" dt2*f(x[i], t[n])\n",
" else:\n",
" u[0] = U_0(t[n+1])\n",
"\n",
" i = mathcal{I}_x[-1]\n",
" if U_L is None:\n",
" im1 = i-1\n",
" ip1 = im1\n",
" u[i] = - u_2[i] + 2*u_1[i] + \\\n",
" C2*(u_1[im1] - 2*u_1[i] + u_1[ip1]) + \\\n",
" dt2*f(x[i], t[n])\n",
" else:\n",
" u[i] = U_L(t[n+1])\n",
"\n",
" if user_action is not None:\n",
" if user_action(u, x, t, n+1):\n",
" break\n",
"\n",
" # Update data structures for next step\n",
" u_2[:], u_1[:] = u_1, u\n",
"\n",
" cpu_time = t0 - time.clock()\n",
" return u, x, t, cpu_time\n",
"\n",
"\n",
"def viz(I, V, f, c, U_0, U_L, x0, xL, Nx, C, T, umin, umax,\n",
" version='scalar', animate=True,\n",
" movie_dir='tmp'):\n",
" \"\"\"Run solver and visualize u at each time level.\"\"\"\n",
" import scitools.std as plt, time, glob, os\n",
"\n",
" def plot_u(u, x, t, n):\n",
" \"\"\"user_action function for solver.\"\"\"\n",
" plt.plot(x, u, 'r-',\n",
" xlabel='x', ylabel='u',\n",
" axis=[x0, xL, umin, umax],\n",
" title='t=%f' % t[n])\n",
" # Let the initial condition stay on the screen for 2\n",
" # seconds, else insert a pause of 0.2 s between each plot\n",
" time.sleep(2) if t[n] == 0 else time.sleep(0.2)\n",
" plt.savefig('frame_%04d.png' % n) # for movie making\n",
"\n",
" # Clean up old movie frames\n",
" for filename in glob.glob('frame_*.png'):\n",
" os.remove(filename)\n",
"\n",
" user_action = plot_u if animate else None\n",
" u, x, t, cpu = solver(I, V, f, c, U_0, U_L, L, Nx, C, T,\n",
" user_action, version)\n",
" if animate:\n",
" # Make a directory with the frames\n",
" if os.path.isdir(movie_dir):\n",
" shutil.rmtree(movie_dir)\n",
" os.mkdir(movie_dir)\n",
" os.chdir(movie_dir)\n",
" # Move all frame_*.png files to this subdirectory\n",
" for filename in glob.glob(os.path.join(os.pardir, 'frame_*.png')):\n",
" os.renamve(os.path.join(os.pardir, filename), filename)\n",
" plt.movie('frame_*.png', encoder='html', fps=4,\n",
" output_file='movie.html')\n",
" # Invoke movie.html in a browser to steer the movie\n",
"\n",
" return cpu\n",
"\n",
"import nose.tools as nt\n",
"\n",
"def test_quadratic():\n",
" \"\"\"\n",
" Check the scalar and vectorized versions work for\n",
" a quadratic u(x,t)=x(L-x)(1+t/2) that is exactly reproduced.\n",
" We simulate in [0, L/2] and apply a symmetry condition\n",
" at the end x=L/2.\n",
" \"\"\"\n",
" exact_solution = lambda x, t: x*(L-x)*(1+0.5*t)\n",
" I = lambda x: exact_solution(x, 0)\n",
" V = lambda x: 0.5*exact_solution(x, 0)\n",
" f = lambda x, t: 2*(1+0.5*t)*c**2\n",
" U_0 = lambda t: exact_solution(0, t)\n",
" U_L = None\n",
" L = 2.5\n",
" c = 1.5\n",
" Nx = 3 # very coarse mesh\n",
" C = 1\n",
" T = 18 # long time integration\n",
"\n",
" def assert_no_error(u, x, t, n):\n",
" u_e = exact_solution(x, t[n])\n",
" diff = abs(u - u_e).max()\n",
" nt.assert_almost_equal(diff, 0, places=13)\n",
"\n",
" solver(I, V, f, c, U_0, U_L, 0, L/2, Nx, C, T,\n",
" user_action=assert_no_error, version='scalar')\n",
" solver(I, V, f, c, U_0, U_L, 0, L/2, Nx, C, T,\n",
" user_action=assert_no_error, version='vectorized')\n",
"\n",
"\n",
"def plug(C=1, Nx=50, animate=True, version='scalar', T=2):\n",
" \"\"\"Plug profile as initial condition.\"\"\"\n",
" L = 1.\n",
" c = 1\n",
" delta = 0.1\n",
"\n",
" def I(x):\n",
" if abs(x) > delta:\n",
" return 0\n",
" else:\n",
" return 1\n",
"\n",
" # Solution on [-L,L]\n",
" cpu = viz(I=I, V=0, f=0, c, U_0=0, U_L=0, -L, L, 2*Nx, C, T,\n",
" umin=-1.1, umax=1.1, version=version, animate=animate,\n",
" movie_dir='full')\n",
"\n",
" # Solution on [0,L]\n",
" cpu = viz(I=I, V=0, f=0, c, U_0=None, U_L=0, 0, L, Nx, C, T,\n",
" umin=-1.1, umax=1.1, version=version, animate=animate,\n",
" movie_dir='frac{1}{2}')\n",
"\n",
"\n",
"if __name__ == '__main__':\n",
" plug()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- --- end solution of exercise --- -->\n",
"\n",
"\n",
"Filename: `wave1D_symmetric`.\n",
"\n",
"<!-- --- end exercise --- -->\n",
"\n",
"\n",
"\n",
"\n",
"<!-- --- begin exercise --- -->\n",
"\n",
"## Exercise 3: Send pulse waves through a layered medium\n",
"<div id=\"wave:app:exer:pulse1D\"></div>\n",
"\n",
"Use the `pulse` function in `wave1D_dn_vc.py` to investigate\n",
"sending a pulse, located with its peak at $x=0$, through two\n",
"media with different wave velocities. The (scaled) velocity in\n",
"the left medium is 1 while it is $\\frac{1}{s_f}$ in the right medium.\n",
"Report what happens with a Gaussian pulse, a \"cosine hat\" pulse,\n",
"frac{1}{2} a \"cosine hat\" pulse, and a plug pulse for resolutions\n",
"$N_x=40,80,160$, and $s_f=2,4$. Simulate until $T=2$.\n",
"\n",
"\n",
"<!-- --- begin solution of exercise --- -->\n",
"**Solution.**\n",
"In all cases, the change in velocity causes some of the wave to\n",
"be reflected back (while the rest is let through). When the waves\n",
"go from higher to lower velocity, the amplitude builds, and vice versa."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"import wave1D_dn_vc as wave\n",
"import os, sys, shutil, glob\n",
"\n",
"for pulse_thinspace . in 'gaussian', 'cosinehat', 'frac{1}{2}-cosinehat', 'plug':\n",
" for Nx in 40, 80, 160:\n",
" for sf in 2, 4:\n",
" if sf == 1 and Nx > 40:\n",
" continue # homogeneous medium with C=1: Nx=40 enough\n",
" print 'wave1D.pulse:', pulse_thinspace ., Nx, sf\n",
"\n",
" wave.pulse(C=1, Nx=Nx, animate=False, # just hardcopies\n",
" version='vectorized',\n",
" T=2, loc='left', pulse_thinspace .=pulse_thinspace .,\n",
" slowness_factor=sf, medium=[0.7, 0.9],\n",
" skip_frame = 1,\n",
" sigma=0.05)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- --- end solution of exercise --- -->\n",
"Filename: `pulse1D`.\n",
"\n",
"<!-- --- end exercise --- -->\n",
"\n",
"\n",
"\n",
"\n",
"<!-- --- begin exercise --- -->\n",
"\n",
"## Exercise 4: Explain why numerical noise occurs\n",
"<div id=\"wave:app:exer:pulse1D:analysis\"></div>\n",
"\n",
"The experiments performed in [Exercise 3: Send pulse waves through a layered medium](#wave:app:exer:pulse1D) shows\n",
"considerable numerical noise in the form of non-physical waves,\n",
"especially for $s_f=4$ and the plug pulse or the frac{1}{2} a \"cosinehat\"\n",
"pulse. The noise is much less visible for a Gaussian pulse. Run the\n",
"case with the plug and frac{1}{2} a \"cosinehat\" pulse for $s_f=1$, $C=0.9,\n",
"0.25$, and $N_x=40,80,160$. Use the numerical dispersion relation to\n",
"explain the observations.\n",
"Filename: `pulse1D_analysis`.\n",
"\n",
"<!-- --- end exercise --- -->\n",
"\n",
"\n",
"\n",
"\n",
"<!-- --- begin exercise --- -->\n",
"\n",
"## Exercise 5: Investigate harmonic averaging in a 1D model\n",
"<div id=\"wave:app:exer:pulse1D:harmonic\"></div>\n",
"\n",
"Harmonic means are often used if the wave velocity is non-smooth or\n",
"discontinuous. Will harmonic averaging of the wave velocity give less\n",
"numerical noise for the case $s_f=4$ in [Exercise 3: Send pulse waves through a layered medium](#wave:app:exer:pulse1D)?\n",
"Filename: `pulse1D_harmonic`.\n",
"\n",
"<!-- --- end exercise --- -->\n",
"\n",
"\n",
"\n",
"\n",
"<!-- --- begin exercise --- -->\n",
"\n",
"## Problem 6: Implement open boundary conditions\n",
"<div id=\"wave:app:exer:radiationBC\"></div>\n",
"\n",
"<!-- Solution file is actually periodic.py from Exer [Exercise 7: Implement periodic boundary conditions](#wave:exer:periodic), -->\n",
"<!-- just remove the periodc stuff ;-) -->\n",
"\n",
"\n",
"To enable a wave to leave the computational domain and travel\n",
"undisturbed through\n",
"the boundary $x=L$, one can in a one-dimensional problem impose the\n",
"following condition, called a *radiation condition* or\n",
"*open boundary condition*:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:app:exer:radiationBC:eq\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial u}{\\partial t} + c\\frac{\\partial u}{\\partial x} = 0\\thinspace .\n",
"\\label{wave:app:exer:radiationBC:eq} \\tag{55}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The parameter $c$ is the wave velocity.\n",
"\n",
"Show that ([55](#wave:app:exer:radiationBC:eq)) accepts\n",
"a solution $u = g_R(x-ct)$ (right-going wave),\n",
"but not $u = g_L(x+ct)$ (left-going wave). This means\n",
"that ([55](#wave:app:exer:radiationBC:eq)) will allow any\n",
"right-going wave $g_R(x-ct)$ to pass through the boundary undisturbed.\n",
"\n",
"A corresponding open boundary condition for a left-going wave\n",
"through $x=0$ is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:app:exer:radiationBC:eqL\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial u}{\\partial t} - c\\frac{\\partial u}{\\partial x} = 0\\thinspace .\n",
"\\label{wave:app:exer:radiationBC:eqL} \\tag{56}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**a)**\n",
"A natural idea for discretizing\n",
"the condition ([55](#wave:app:exer:radiationBC:eq))\n",
"at the spatial end point $i=N_x$ is to apply\n",
"centered differences in time and space:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:app:exer:radiationBC:eq:op\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"[D_{2t}u + cD_{2x}u =0]^n_{i},\\quad i=N_x\\thinspace .\n",
"\\label{wave:app:exer:radiationBC:eq:op} \\tag{57}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Eliminate the fictitious value $u_{N_x+1}^n$ by using\n",
"the discrete equation at the same point.\n",
"\n",
"The equation for the first step, $u_i^1$, is in principle also affected,\n",
"but we can then use the condition $u_{N_x}=0$ since the wave\n",
"has not yet reached the right boundary.\n",
"\n",
"**b)**\n",
"A much more convenient implementation of the open boundary condition\n",
"at $x=L$ can be based on an explicit discretization"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:app:exer:radiationBC:eq:op:1storder\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"[D^+_tu + cD_x^- u = 0]_i^n,\\quad i=N_x\\thinspace .\n",
"\\label{wave:app:exer:radiationBC:eq:op:1storder} \\tag{58}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From this equation, one can solve for $u^{n+1}_{N_x}$ and apply the\n",
"formula as a Dirichlet condition at the boundary point.\n",
"However, the finite difference approximations involved are of\n",
"first order.\n",
"\n",
"Implement this scheme for a wave equation\n",
"$u_{tt}=c^2u_{xx}$ in a domain $[0,L]$,\n",
"where you have $u_x=0$ at $x=0$, the condition ([55](#wave:app:exer:radiationBC:eq))\n",
"at $x=L$, and an initial disturbance in the middle\n",
"of the domain, e.g., a plug profile like"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(x,0) = \\left\\lbrace\\begin{array}{ll} 1,& L/2-\\ell \\leq x \\leq L/2+\\ell,\\\\ \n",
"0,\\hbox{otherwise}\\end{array}\\right.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Observe that the initial wave is split in two, the left-going wave\n",
"is reflected at $x=0$, and both waves travel out of $x=L$,\n",
"leaving the solution as $u=0$ in $[0,L]$. Use a unit Courant number\n",
"such that the numerical solution is exact.\n",
"Make a movie to illustrate what happens.\n",
"\n",
"Because this simplified\n",
"implementation of the open boundary condition works, there is no\n",
"need to pursue the more complicated discretization in a).\n",
"\n",
"<!-- --- begin hint in exercise --- -->\n",
"\n",
"**Hint.**\n",
"Modify the solver function in\n",
"[`wave1D_dn.py`](${src_wave}/wave1D/wave1D_dn.py).\n",
"\n",
"<!-- --- end hint in exercise --- -->\n",
"\n",
"**c)**\n",
"Add the possibility to have either $u_x=0$ or an open boundary\n",
"condition at the left boundary. The latter condition is discretized\n",
"as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:app:exer:radiationBC:eq:op:1storder2\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"[D^+_tu - cD_x^+ u = 0]_i^n,\\quad i=0,\n",
"\\label{wave:app:exer:radiationBC:eq:op:1storder2} \\tag{59}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"leading to an explicit update of the boundary value $u^{n+1}_0$.\n",
"\n",
"The implementation can be tested with a Gaussian function as initial condition:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"g(x;m,s) = \\frac{1}{\\sqrt{2\\pi}s}e^{-\\frac{(x-m)^2}{2s^2}}\\thinspace .\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Run two tests:\n",
"\n",
"1. Disturbance in the middle of the domain, $I(x)=g(x;L/2,s)$, and\n",
" open boundary condition at the left end.\n",
"\n",
"2. Disturbance at the left end, $I(x)=g(x;0,s)$, and $u_x=0$\n",
" as symmetry boundary condition at this end.\n",
"\n",
"Make test functions for both cases, testing that the solution is zero\n",
"after the waves have left the domain.\n",
"\n",
"**d)**\n",
"In 2D and 3D it is difficult to compute the correct wave velocity\n",
"normal to the boundary, which is needed in generalizations of\n",
"the open boundary conditions in higher dimensions. Test the effect\n",
"of having a slightly wrong wave velocity in\n",
"([58](#wave:app:exer:radiationBC:eq:op:1storder)).\n",
"Make movies to illustrate what happens.\n",
"\n",
"\n",
"\n",
"Filename: `wave1D_open_BC`.\n",
"\n",
"<!-- Closing remarks for this Problem -->\n",
"\n",
"### Remarks\n",
"\n",
"The condition ([55](#wave:app:exer:radiationBC:eq))\n",
"works perfectly in 1D when $c$ is known. In 2D and 3D, however, the\n",
"condition reads $u_t + c_x u_x + c_y u_y=0$, where $c_x$ and\n",
"$c_y$ are the wave speeds in the $x$ and $y$ directions. Estimating\n",
"these components (i.e., the direction of the wave) is often\n",
"challenging. Other methods are normally used in 2D and 3D to\n",
"let waves move out of a computational domain.\n",
"\n",
"\n",
"<!-- --- end exercise --- -->\n",
"\n",
"\n",
"\n",
"\n",
"<!-- --- begin exercise --- -->\n",
"\n",
"## Exercise 7: Implement periodic boundary conditions\n",
"<div id=\"wave:exer:periodic\"></div>\n",
"\n",
"\n",
"mathcal{I}_t is frequently of interest to follow wave motion over large\n",
"distances and long times. A straightforward approach is to\n",
"work with a very large domain, but that might lead to a lot of\n",
"computations in areas of the domain where the waves cannot\n",
"be noticed. A more efficient approach is to let a right-going\n",
"wave out of the domain and at the same time let it enter\n",
"the domain on the left. This is called a *periodic boundary\n",
"condition*.\n",
"\n",
"The boundary condition at the right end $x=L$ is an open boundary\n",
"condition (see [Problem 6: Implement open boundary conditions](#wave:app:exer:radiationBC)) to let a\n",
"right-going wave out of the domain. At the left end, $x=0$, we apply,\n",
"in the beginning of the simulation, either a symmetry boundary\n",
"condition (see [Problem 2: Explore symmetry boundary conditions](#wave:exer:symmetry:bc)) $u_x=0$, or an\n",
"open boundary condition.\n",
"\n",
"This initial wave will split in two and either be reflected or\n",
"transported out of the domain at $x=0$. The purpose of the exercise is\n",
"to follow the right-going wave. We can do that with a *periodic\n",
"boundary condition*. This means that when the right-going wave hits\n",
"the boundary $x=L$, the open boundary condition lets the wave out of\n",
"the domain, but at the same time we use a boundary condition on the\n",
"left end $x=0$ that feeds the outgoing wave into the domain\n",
"again. This periodic condition is simply $u(0)=u(L)$. The switch from\n",
"$u_x=0$ or an open boundary condition at the left end to a periodic\n",
"condition can happen when $u(L,t)>\\epsilon$, where $\\epsilon =10^{-4}$\n",
"might be an appropriate value for determining when the right-going\n",
"wave hits the boundary $x=L$.\n",
"\n",
"The open boundary conditions can conveniently be discretized as\n",
"explained in [Problem 6: Implement open boundary conditions](#wave:app:exer:radiationBC). Implement the\n",
"described type of boundary conditions and test them on two different\n",
"initial shapes: a plug $u(x,0)=1$ for $x\\leq 0.1$, $u(x,0)=0$ for\n",
"$x>0.1$, and a Gaussian function in the middle of the domain:\n",
"$u(x,0)=\\exp{(-\\frac{1}{2}(x-0.5)^2/0.05)}$. The domain is the unit\n",
"interval $[0,1]$. Run these two shapes for Courant numbers 1 and\n",
"0.5. Assume constant wave velocity. Make movies of the four cases.\n",
"Reason why the solutions are correct.\n",
"Filename: `periodic`.\n",
"\n",
"<!-- --- end exercise --- -->\n",
"\n",
"\n",
"\n",
"\n",
"<!-- --- begin exercise --- -->\n",
"\n",
"## Exercise 8: Compare discretizations of a Neumann condition\n",
"\n",
"We have a 1D wave equation with variable wave velocity:\n",
"$u_{tt}=(qu_x)_x$.\n",
"A Neumann condition $u_x$ at $x=0, L$ can be\n",
"discretized as shown in ([20](#wave:pde2:var:c:scheme:impl:Neumann))\n",
"and ([23](#wave:pde2:var:c:scheme:impl:Neumann2)).\n",
"\n",
"The aim of this exercise is to examine the rate of the numerical\n",
"error when using different ways of discretizing the Neumann condition.\n",
"\n",
"\n",
"**a)**\n",
"As a test problem, $q=1+(x-L/2)^4$ can be used, with $f(x,t)$\n",
"adapted such that the solution has a simple form, say\n",
"$u(x,t)=\\cos (\\pi x/L)\\cos (\\omega t)$ for, e.g., $\\omega = 1$.\n",
"Perform numerical experiments and find the convergence rate of the\n",
"error using the approximation\n",
"([20](#wave:pde2:var:c:scheme:impl:Neumann)).\n",
"\n",
"**b)**\n",
"Switch to $q(x)=1+\\cos(\\pi x/L)$, which is symmetric at $x=0,L$,\n",
"and check the convergence rate\n",
"of the scheme\n",
"([23](#wave:pde2:var:c:scheme:impl:Neumann2)). Now,\n",
"$q_{i-1/2}$ is a 2nd-order approximation to $q_i$,\n",
"$q_{i-1/2}=q_i + 0.25q_i''\\Delta x^2 + \\cdots$, because $q_i'=0$\n",
"for $i=N_x$ (a similar argument can be applied to the case $i=0$).\n",
"\n",
"**c)**\n",
"A third discretization can be based on a simple and convenient,\n",
"but less accurate, one-sided difference:\n",
"$u_{i}-u_{i-1}=0$ at $i=N_x$ and $u_{i+1}-u_i=0$ at $i=0$.\n",
"Derive the resulting scheme in detail and implement it.\n",
"Run experiments with $q$ from a) or b) to establish the rate of convergence\n",
"of the scheme.\n",
"\n",
"**d)**\n",
"A fourth technique is to view the scheme as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"[D_tD_tu]^n_i = \\frac{1}{\\Delta x}\\left(\n",
"[qD_xu]_{i+\\frac{1}{2}}^n - [qD_xu]_{i-\\frac{1}{2}}^n\\right)\n",
"+ [f]_i^n,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and place the boundary at $x_{i+\\frac{1}{2}}$, $i=N_x$, instead of\n",
"exactly at the physical boundary. With this idea of approximating (moving) the\n",
"boundary,\n",
"we can just set $[qD_xu]_{i+\\frac{1}{2}}^n=0$.\n",
"Derive the complete scheme\n",
"using this technique. The implementation of the boundary condition at\n",
"$L-\\Delta x/2$ is $\\Oof{\\Delta x^2}$ accurate, but the interesting question\n",
"is what impact the movement of the boundary has on the convergence\n",
"rate. Compute the errors as usual over the entire mesh and use $q$ from\n",
"a) or b).\n",
"\n",
"\n",
"Filename: `Neumann_discr`.\n",
"\n",
"<!-- --- end exercise --- -->\n",
"\n",
"\n",
"\n",
"\n",
"<!-- --- begin exercise --- -->\n",
"\n",
"## Exercise 9: Verification by a cubic polynomial in space\n",
"<div id=\"wave:fd2:exer:verify:cubic\"></div>\n",
"\n",
"The purpose of this exercise is to verify the implementation of the\n",
"`solver` function in the program [`wave1D_n0.py`](#src_wave/wave1D/wave1D_n0.py) by using an exact numerical solution\n",
"for the wave equation $u_{tt}=c^2u_{xx} + f$ with Neumann boundary\n",
"conditions $u_x(0,t)=u_x(L,t)=0$.\n",
"\n",
"A similar verification is used in the file [`wave1D_u0.py`](#src_wave}/wave1D/wave1D_u0.py), which solves the same PDE, but with\n",
"Dirichlet boundary conditions $u(0,t)=u(L,t)=0$. The idea of the\n",
"verification test in function `test_quadratic` in `wave1D_u0.py` is to\n",
"produce a solution that is a lower-order polynomial such that both the\n",
"PDE problem, the boundary conditions, and all the discrete equations\n",
"are exactly fulfilled. Then the `solver` function should reproduce\n",
"this exact solution to machine precision. More precisely, we seek\n",
"$u=X(x)T(t)$, with $T(t)$ as a linear function and $X(x)$ as a\n",
"parabola that fulfills the boundary conditions. Inserting this $u$ in\n",
"the PDE determines $f$. mathcal{I}_t turns out that $u$ also fulfills the\n",
"discrete equations, because the truncation error of the discretized\n",
"PDE has derivatives in $x$ and $t$ of order four and higher. These\n",
"derivatives all vanish for a quadratic $X(x)$ and linear $T(t)$.\n",
"\n",
"mathcal{I}_t would be attractive to use a similar approach in the case of\n",
"Neumann conditions. We set $u=X(x)T(t)$ and seek lower-order\n",
"polynomials $X$ and $T$.\n",
"To force $u_x$ to vanish at the boundary, we let $X_x$ be\n",
"a parabola. Then $X$ is a cubic polynomial. The fourth-order\n",
"derivative of a cubic polynomial vanishes, so $u=X(x)T(t)$\n",
"will fulfill the discretized PDE also in this case, if $f$\n",
"is adjusted such that $u$ fulfills the PDE.\n",
"\n",
"However, the discrete boundary condition is not exactly fulfilled\n",
"by this choice of $u$. The reason is that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"wave:fd2:exer:verify:cubic:D2x\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"[D_{2x}u]^n_i = u_{x}(x_i,t_n) + \\frac{1}{6}u_{xxx}(x_i,t_n)\\Delta x^2\n",
"+ \\Oof{\\Delta x^4}\\thinspace .\n",
"\\label{wave:fd2:exer:verify:cubic:D2x} \\tag{60}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"At the two boundary points, we must demand that\n",
"the derivative $X_x(x)=0$ such that $u_x=0$.\n",
"However, $u_{xxx}$ is a constant and not zero\n",
"when $X(x)$ is a cubic polynomial.\n",
"Therefore, our $u=X(x)T(t)$ fulfills"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"[D_{2x}u]^n_i = \\frac{1}{6}u_{xxx}(x_i,t_n)\\Delta x^2,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and not"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"[D_{2x}u]^n_i =0, \\quad i=0,N_x,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"as it should. (Note that all the higher-order terms $\\Oof{\\Delta x^4}$\n",
"also have higher-order derivatives that vanish for a cubic polynomial.)\n",
"So to summarize, the fundamental problem is that $u$ as a product of\n",
"a cubic polynomial and a linear or quadratic polynomial in time\n",
"is not an exact solution of the discrete boundary conditions.\n",
"\n",
"To make progress,\n",
"we assume that $u=X(x)T(t)$, where $T$ for simplicity is taken as a\n",
"prescribed linear function $1+\\frac{1}{2}t$, and $X(x)$ is taken\n",
"as an *unknown* cubic polynomial $\\sum_{j=0}^3 a_jx^j$.\n",
"There are two different ways of determining the coefficients\n",
"$a_0,\\ldots,a_3$ such that both the discretized PDE and the\n",
"discretized boundary conditions are fulfilled, under the\n",
"constraint that we can specify a function $f(x,t)$ for the PDE to feed\n",
"to the `solver` function in `wave1D_n0.py`. Both approaches\n",
"are explained in the subexercises.\n",
"\n",
"\n",
"<!-- {wave:fd2:exer:verify:cubic:D2x} -->\n",
"\n",
"\n",
"**a)**\n",
"One can insert $u$ in the discretized PDE and find the corresponding $f$.\n",
"Then one can insert $u$ in the discretized boundary conditions.\n",
"This yields two equations for the four coefficients $a_0,\\ldots,a_3$.\n",
"To find the coefficients, one can set $a_0=0$ and $a_1=1$ for\n",
"simplicity and then determine $a_2$ and $a_3$. This approach will make\n",
"$a_2$ and $a_3$ depend on $\\Delta x$ and $f$ will depend on both\n",
"$\\Delta x$ and $\\Delta t$.\n",
"\n",
"Use `sympy` to perform analytical computations.\n",
"A starting point is to define $u$ as follows:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"def test_cubic1():\n",
" import sympy as sm\n",
" x, t, c, L, dx, dt = sm.symbols('x t c L dx dt')\n",
" i, n = sm.symbols('i n', integer=True)\n",
"\n",
" # Assume discrete solution is a polynomial of degree 3 in x\n",
" T = lambda t: 1 + sm.Rational(1,2)*t # Temporal term\n",
" a = sm.symbols('a_0 a_1 a_2 a_3')\n",
" X = lambda x: sum(a[q]*x**q for q in range(4)) # Spatial term\n",
" u = lambda x, t: X(x)*T(t)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The symbolic expression for $u$ is reached by calling `u(x,t)`\n",
"with `x` and `t` as `sympy` symbols.\n",
"\n",
"Define `DxDx(u, i, n)`, `DtDt(u, i, n)`, and `D2x(u, i, n)`\n",
"as Python functions for returning the difference\n",
"approximations $[D_xD_x u]^n_i$, $[D_tD_t u]^n_i$, and\n",
"$[D_{2x}u]^n_i$. The next step is to set up the residuals\n",
"for the equations $[D_{2x}u]^n_0=0$ and $[D_{2x}u]^n_{N_x}=0$,\n",
"where $N_x=L/\\Delta x$. Call the residuals `R_0` and `R_L`.\n",
"Substitute $a_0$ and $a_1$ by 0 and 1, respectively, in\n",
"`R_0`, `R_L`, and `a`:"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"R_0 = R_0.subs(a[0], 0).subs(a[1], 1)\n",
"R_L = R_L.subs(a[0], 0).subs(a[1], 1)\n",
"a = list(a) # enable in-place assignment\n",
"a[0:2] = 0, 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Determining $a_2$ and $a_3$ from the discretized boundary conditions\n",
"is then about solving two equations with respect to $a_2$ and $a_3$,\n",
"i.e., `a[2:]`:"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"s = sm.solve([R_0, R_L], a[2:])\n",
"# s is dictionary with the unknowns a[2] and a[3] as keys\n",
"a[2:] = s[a[2]], s[a[3]]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, `a` contains computed values and `u` will automatically use\n",
"these new values since `X` accesses `a`.\n",
"\n",
"Compute the source term $f$ from the discretized PDE:\n",
"$f^n_i = [D_tD_t u - c^2D_xD_x u]^n_i$. Turn $u$, the time\n",
"derivative $u_t$ (needed for the initial condition $V(x)$),\n",
"and $f$ into Python functions. Set numerical values for\n",
"$L$, $N_x$, $C$, and $c$. Prescribe the time interval as\n",
"$\\Delta t = CL/(N_xc)$, which imply $\\Delta x = c\\Delta t/C = L/N_x$.\n",
"Define new functions `I(x)`, `V(x)`, and `f(x,t)` as wrappers of the ones\n",
"made above, where fixed values of $L$, $c$, $\\Delta x$, and $\\Delta t$\n",
"are inserted, such that `I`, `V`, and `f` can be passed on to the\n",
"`solver` function. Finally, call `solver` with a `user_action`\n",
"function that compares the numerical solution to this exact\n",
"solution $u$ of the discrete PDE problem.\n",
"\n",
"<!-- --- begin hint in exercise --- -->\n",
"\n",
"**Hint.**\n",
"To turn a `sympy` expression `e`, depending on a series of\n",
"symbols, say `x`, `t`, `dx`, `dt`, `L`, and `c`, into a plain\n",
"Python function `e_exact(x,t,L,dx,dt,c)`, one can write"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"e_exact = sm.lambdify([x,t,L,dx,dt,c], e, 'numpy')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `'numpy'` argument is a good habit as the `e_exact` function\n",
"will then work with array arguments if it contains mathematical\n",
"functions (but here we only do plain arithmetics, which automatically\n",
"work with arrays).\n",
"\n",
"<!-- --- end hint in exercise --- -->\n",
"\n",
"**b)**\n",
"An alternative way of determining $a_0,\\ldots,a_3$ is to reason as\n",
"follows. We first construct $X(x)$ such that the boundary conditions\n",
"are fulfilled: $X=x(L-x)$. However, to compensate for the fact\n",
"that this choice of $X$ does not fulfill the discrete boundary\n",
"condition, we seek $u$ such that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u_x = \\frac{\\partial}{\\partial x}x(L-x)T(t) - \\frac{1}{6}u_{xxx}\\Delta x^2,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"since this $u$ will fit the discrete boundary condition.\n",
"Assuming $u=T(t)\\sum_{j=0}^3a_jx^j$, we can use the above equation to\n",
"determine the coefficients $a_1,a_2,a_3$. A value, e.g., 1 can be used for\n",
"$a_0$. The following `sympy` code computes this $u$:"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"def test_cubic2():\n",
" import sympy as sm\n",
" x, t, c, L, dx = sm.symbols('x t c L dx')\n",
" T = lambda t: 1 + sm.Rational(1,2)*t # Temporal term\n",
" # Set u as a 3rd-degree polynomial in space\n",
" X = lambda x: sum(a[i]*x**i for i in range(4))\n",
" a = sm.symbols('a_0 a_1 a_2 a_3')\n",
" u = lambda x, t: X(x)*T(t)\n",
" # Force discrete boundary condition to be zero by adding\n",
" # a correction term the analytical suggestion x*(L-x)*T\n",
" # u_x = x*(L-x)*T(t) - 1/6*u_xxx*dx**2\n",
" R = sm.diff(u(x,t), x) - (\n",
" x*(L-x) - sm.Rational(1,6)*sm.diff(u(x,t), x, x, x)*dx**2)\n",
" # R is a polynomial: force all coefficients to vanish.\n",
" # Turn R to Poly to extract coefficients:\n",
" R = sm.poly(R, x)\n",
" coeff = R.all_coeffs()\n",
" s = sm.solve(coeff, a[1:]) # a[0] is not present in R\n",
" # s is dictionary with a[i] as keys\n",
" # Fix a[0] as 1\n",
" s[a[0]] = 1\n",
" X = lambda x: sm.simplify(sum(s[a[i]]*x**i for i in range(4)))\n",
" u = lambda x, t: X(x)*T(t)\n",
" print 'u:', u(x,t)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The next step is to find the source term `f_e` by inserting `u_e`\n",
"in the PDE. Thereafter, turn `u`, `f`, and the time derivative of `u`\n",
"into plain Python functions as in a), and then wrap these functions\n",
"in new functions `I`, `V`, and `f`, with the right signature as\n",
"required by the `solver` function. Set parameters as in a) and\n",
"check that the solution is exact to machine precision at each\n",
"time level using an appropriate `user_action` function.\n",
"\n",
"Filename: `wave1D_n0_test_cubic`.\n",
"\n",
"<!-- --- end exercise --- -->"
]
}
],
"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.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}