\n",
+ " \n",
+ " \n",
+ " \n",
+ "
\n",
+ "Modified stencil at a boundary with a Neumann condition.
\n", + "![](\"mov-wave/N_stencil_gpl/stencil_n_left.png\")
\n",
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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). Using $i = 0$ as an example, the steps above can be implemented in Devito by"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from devito import *\n",
+ "\n",
+ "# du/dt(x,0) = V(x)\n",
+ "# We take the original PDE and apply central finite differences\n",
+ "# But if we let time = 0, we will require a value with time = -1\n",
+ "# We can also apply central finite differences to du/dt(x,0) = V(x),\n",
+ "# and rewrite the term with time = -1 in terms of real terms, and do substitution to get the equation below\n",
+ "\n",
+ "# Solving for u(t+dt, x) in the actual PDE\n",
+ "eq_u = Eq(u.dt2, c**2*u.dx2 + F)\n",
+ "stencil_u = solve(eq_u.evaluate, u.forward)\n",
+ "update_u = Eq(u.forward, stencil_u)\n",
+ "\n",
+ "# Applying central differences to du/dt(x,0) = V(x) and solving for u(-dt, x) in terms of points with physical meaning\n",
+ "Neumann_eq = Eq(u.dtc, V)\n",
+ "u_negative_time = solve(Neumann_eq.evaluate, u.backward)._subs(u.time_dim, 0)\n",
+ "\n",
+ "# Get the stencil at t=0 and replace u(-dt, x) with {u_negative_time} that we computed above\n",
+ "stencil_at_t0 = Eq(u.forward._subs(u.time_dim, 0), stencil_u._subs(u.backward.evaluate, u_negative_time.evaluate))\n",
+ "stencil_at_t0 = solve(stencil_at_t0.evaluate, u.forward._subs(u.time_dim, 0))._subs(u.time_dim, 0)\n",
+ "lower_Neumann = Eq(u.forward._subs(u.time_dim, 0), stencil_at_t0)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We would replace `._subs(u.time_dim, 0)` with `._sub(u.time_dim, Nx)` for the $i = N_x$ case"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Implementation of Neumann conditions\n",
+ "\n",
+ "\n",
+ "To deal with potential out of bounds problem, 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."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Dealing with boundary points is made very easy by Devito as Devito automatically adds ghost points (called halo) to defined functions. The size of the `halo` is equal to its `space_order` (refer to `examples/compiler/01_data_region` example in the main Devito Github repo). The code below implements said function, an additional test case can be found in [`wave1D_n0.py`](${src_wave}/wave1D/wave1D_n0.py) along with the implementation."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def solver(i, v, f, c, L, dt, C, T, user_action=None):\n",
+ " \"\"\"\n",
+ " Solve u_tt=c^2*u_xx + f on (0,L)x(0,T].\n",
+ "\n",
+ " Boundary Conditions:\n",
+ " du/dn = 0 at both U(0,t) and U(L,t)\n",
+ " u(x,0) = I(x)\n",
+ " du/dt(x,0) = V(x)\n",
+ " \"\"\"\n",
+ "\n",
+ " Nt = int(round(T/dt))\n",
+ " dx = dt*c/float(C)\n",
+ " Nx = int(round(L/dx))\n",
+ " \n",
+ " # Wrap user-given f, V\n",
+ " if f is None or f == 0:\n",
+ " f = (lambda x, t: 0)\n",
+ " if v is None or v == 0:\n",
+ " v = (lambda x: 0)\n",
+ " \n",
+ " # A padding, called halo in this case, of size {space_order} is automatically added when the grid is initialized\n",
+ " grid = Grid(shape=(Nx+1,), dtype=np.float64)\n",
+ "\n",
+ " # We need to have u(t,x) instead of u(x,t) so that the produced code loops through time {t} on the outside\n",
+ " # This is necessary because our boundary conditions are on the space-axis\n",
+ " u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)\n",
+ "\n",
+ " # Devito-ize external functions\n",
+ " t = grid.stepping_dim\n",
+ " x = grid.dimensions[0]\n",
+ " F = TimeFunction(name='F', dimensions=(t, x), shape=(Nt+1, Nx+1), space_order=2)\n",
+ " for time in range(0, Nt+1):\n",
+ " for space in range(0, Nx+1):\n",
+ " # f is given with arguments (x,t) not (t,x)\n",
+ " F.data[time][space] = f(space * dx, time * dt)\n",
+ "\n",
+ " I = Function(name='I', dimensions=(x,), shape=(Nx+1,), space_order=2)\n",
+ " for space in range(0, Nx+1):\n",
+ " I.data[space] = i(space*dx)\n",
+ "\n",
+ " V = Function(name='V', dimensions=(x,), shape=(Nx+1,), space_order=2)\n",
+ " for space in range(0, Nx+1):\n",
+ " V.data[space] = v(space*dx)\n",
+ "\n",
+ " # du/dt(x,0) = V(x)\n",
+ " # We take the original PDE and apply central finite differences\n",
+ " # Let time = 0, we will require a value with time = -1\n",
+ " # But we can also apply central finite differences to du/dt(x,0) = V(x),\n",
+ " # We write the term with time = -1 in terms of real terms, and do substitution to get the equation below\n",
+ "\n",
+ " # Solving for u(t+dt, x) in the actual PDE\n",
+ " eq_u = Eq(u.dt2, c**2*u.dx2 + F)\n",
+ " stencil_u = solve(eq_u.evaluate, u.forward)\n",
+ " update_u = Eq(u.forward, stencil_u)\n",
+ "\n",
+ " # Boundary conditions\n",
+ " boundary_eqns = []\n",
+ " bottom_Dirichlet = Eq(u[0, x], I)\n",
+ " boundary_eqns.append(bottom_Dirichlet)\n",
+ " \n",
+ " # Applying central differences to du/dt(x,0) = V(x) and solving for u(t-dt, x)\n",
+ " Neumann_eq = Eq(u.dtc, V)\n",
+ " u_negative_time = solve(Neumann_eq.evaluate, u.backward)._subs(u.time_dim, 0)\n",
+ "\n",
+ " # Get the stencil at t=0 and replace u(t-dt, x) with {u_negative_time} that we computed above\n",
+ " stencil_at_t0 = Eq(u.forward._subs(u.time_dim, 0), stencil_u._subs(u.backward.evaluate, u_negative_time.evaluate))\n",
+ " stencil_at_t0 = solve(stencil_at_t0.evaluate, u.forward._subs(u.time_dim, 0))._subs(u.time_dim, 0)\n",
+ " lower_Neumann = Eq(u.forward._subs(u.time_dim, 0), stencil_at_t0)\n",
+ " boundary_eqns.append(lower_Neumann)\n",
+ "\n",
+ " left_Neumann = Eq(u[t, -1], u[t,1])\n",
+ " boundary_eqns.append(left_Neumann)\n",
+ " right_Neumann = Eq(u[t, Nx+1], u[t, Nx-1])\n",
+ " boundary_eqns.append(right_Neumann)\n",
+ "\n",
+ " op = Operator([bottom_Dirichlet, lower_Neumann, left_Neumann, right_Neumann, update_u])\n",
+ "\n",
+ " import time; t0 = time.time()\n",
+ " op.apply(dt=dt, t_m=1, t_M=Nt, x_m=0, x_M=Nx)\n",
+ " cpu_time = time.time() - t0\n",
+ "\n",
+ " # Nt%3 as u.data is a 3 by Nx array. The 3 rows are used to store values at T=t, T=t-1, and T=t-2\n",
+ " # The exact array you return might have be checked case by case\n",
+ " return u.data[Nt%3], np.linspace(0, L, Nx+1), np.linspace(0, Nt*dt, Nt+1), cpu_time\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The program `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))."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "**Notice.**\n",
+ "\n",
+ "The program [`wave1D_dn.py`](src-wave/wave1D/python/wave1D_dn.py)\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 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",
+ "## Verifying the implementation of Neumann conditions\n",
+ "\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",
+ "It 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": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "Left: wave entering another medium; right: transmitted and reflected wave.
\n", + "![](\"fig-wave/pulse1_in_two_media.png\")
\n",
+ "\n",
+ "\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": [
+ "\n",
+ "\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",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Computing the coefficient between mesh points\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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",
+ "\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": [
+ "\n",
+ "\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",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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",
+ "\n",
+ "\n",
+ "To implement this using Devito we can just take the code from the previous ghost points example and replace the equation accordingly. If necessary we can replace values of `q` between mesh points with approximations that we derived in the section right above.\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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",
+ "\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": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "!WHAT DO I DO WITH THIS SECTION?"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\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",
+ "\n",
+ "\n",
+ "## Exercise 1: Find the analytical solution to a damped wave equation\n",
+ "\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",
+ "\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",
+ "\n",
+ "Filename: `damped_waves`.\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Problem 2: Explore symmetry boundary conditions\n",
+ "\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",
+ "\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",
+ "\n",
+ "\n",
+ "\n",
+ "\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",
+ "\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",
+ "\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",
+ "\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",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation} \n",
+ "u(-L,0) = 0\n",
+ "\\label{_auto16} \\tag{38}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation} \n",
+ "u(L,0) = 0\n",
+ "\\label{_auto21} \\tag{43}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation} \n",
+ "v(-L,0) = 0\n",
+ "\\label{_auto26} \\tag{48}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation} \n",
+ "v(-L,0) = 0\n",
+ "\\label{_auto31} \\tag{53}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\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",
+ "\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",
+ "\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": [
+ "\n",
+ "\n",
+ "\n",
+ "Filename: `wave1D_symmetric`.\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Exercise 3: Send pulse waves through a layered medium\n",
+ "\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, half 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",
+ "\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": 9,
+ "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": [
+ "\n",
+ "Filename: `pulse1D`.\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Exercise 4: Explain why numerical noise occurs\n",
+ "\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 half a \"cosinehat\"\n",
+ "pulse. The noise is much less visible for a Gaussian pulse. Run the\n",
+ "case with the plug and half 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",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Exercise 5: Investigate harmonic averaging in a 1D model\n",
+ "\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",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Problem 6: Implement open boundary conditions\n",
+ "\n",
+ "\n",
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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": [
+ "\n",
+ "\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",
+ "\n",
+ "\n",
+ "**Hint.**\n",
+ "Modify the solver function in\n",
+ "[`wave1D_dn.py`](${src_wave}/wave1D/wave1D_dn.py).\n",
+ "\n",
+ "\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": [
+ "\n",
+ "\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",
+ "\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",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Exercise 7: Implement periodic boundary conditions\n",
+ "\n",
+ "\n",
+ "\n",
+ "It 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",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\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",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Exercise 9: Verification by a cubic polynomial in space\n",
+ "\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": [
+ "\n",
+ "\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",
+ "+ \\mathcal{O}(\\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 $\\mathcal{O}(\\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",
+ "\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",
+ "\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",
+ "\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",
+ ""
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.6"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/fdm-jupyter-book/notebooks/02_wave/wave1D_prog.ipynb b/fdm-jupyter-book/notebooks/02_wave/wave1D_prog.ipynb
index 07b646ec..deb0d231 100644
--- a/fdm-jupyter-book/notebooks/02_wave/wave1D_prog.ipynb
+++ b/fdm-jupyter-book/notebooks/02_wave/wave1D_prog.ipynb
@@ -1366,7 +1366,7 @@
" u_1 = np.zeros(Nx+1) # Solution at 1 time level back\n",
" u_2 = np.zeros(Nx+1) # Solution at 2 time levels back\n",
"\n",
- " import time; t0 = time.clock() # for measuring CPU time\n",
+ " import time; t0 = time.perf_counter() # for measuring CPU time\n",
"\n",
" # Load initial condition into u_1\n",
" for i in range(0,Nx+1):\n",
@@ -1405,7 +1405,7 @@
" # Switch variables before next step\n",
" u_2[:] = u_1; u_1[:] = u\n",
"\n",
- " cpu_time = t0 - time.clock()\n",
+ " cpu_time = t0 - time.perf_counter()\n",
" return u, x, t, cpu_time\n",
"\n",
"def test_quadratic():\n",
@@ -1642,7 +1642,7 @@
" u_1 = np.zeros(Nx+1) # Solution at 1 time level back\n",
" u_2 = np.zeros(Nx+1) # Solution at 2 time levels back\n",
"\n",
- " import time; t0 = time.clock() # for measuring CPU time\n",
+ " import time; t0 = time.perf_counter() # for measuring CPU time\n",
"\n",
" # Load initial condition into u_1\n",
" for i in range(0,Nx+1):\n",
@@ -1681,7 +1681,7 @@
" # Switch variables before next step\n",
" u_2[:] = u_1; u_1[:] = u\n",
"\n",
- " cpu_time = t0 - time.clock()\n",
+ " cpu_time = t0 - time.perf_counter()\n",
" return u, x, t, cpu_time\n",
"\n",
"def test_quadratic():\n",
@@ -1928,7 +1928,7 @@
" u_1 = np.zeros(Nx+1) # Solution at 1 time level back\n",
" u_2 = np.zeros(Nx+1) # Solution at 2 time levels back\n",
"\n",
- " import time; t0 = time.clock() # for measuring CPU time\n",
+ " import time; t0 = time.perf_counter() # for measuring CPU time\n",
"\n",
" # Load initial condition into u_1\n",
" for i in range(0,Nx+1):\n",
@@ -1967,7 +1967,7 @@
" # Switch variables before next step\n",
" u_2[:] = u_1; u_1[:] = u\n",
"\n",
- " cpu_time = time.clock() - t0\n",
+ " cpu_time = time.perf_counter() - t0\n",
" return u, x, t, cpu_time\n",
"\n",
"def viz(\n",
@@ -2572,7 +2572,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Python 3",
+ "display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
@@ -2586,7 +2586,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.8.2"
+ "version": "3.10.6"
}
},
"nbformat": 4,
diff --git a/fdm-jupyter-book/notebooks/02_wave/wave2D_fd.ipynb b/fdm-jupyter-book/notebooks/02_wave/wave2D_fd.ipynb
new file mode 100644
index 00000000..6762cc45
--- /dev/null
+++ b/fdm-jupyter-book/notebooks/02_wave/wave2D_fd.ipynb
@@ -0,0 +1,587 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Finite difference methods for 2D and 3D wave equations\n",
+ "\n",
+ "\n",
+ "A natural next step is to consider extensions of the methods for\n",
+ "various\n",
+ "variants of the one-dimensional wave equation to two-dimensional (2D) and\n",
+ "three-dimensional (3D) versions of the wave equation.\n",
+ "\n",
+ "## Multi-dimensional wave equations\n",
+ "\n",
+ "\n",
+ "The general wave equation in $d$ space dimensions, with constant\n",
+ "wave velocity $c$,\n",
+ "can be written in the compact form"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\frac{\\partial^2 u}{\\partial t^2} = c^2\\nabla^2 u\\hbox{ for } \\textbf{x}\\in\\Omega\\subset\\mathbb{R}^d,\\ t\\in (0,T] ,\n",
+ "\\label{wave:2D3D:model1} \\tag{1}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\nabla^2 u = \\frac{\\partial^2 u}{\\partial x^2} +\n",
+ "\\frac{\\partial^2 u}{\\partial y^2} ,\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "in a 2D problem ($d=2$) and"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\nabla^2 u = \\frac{\\partial^2 u}{\\partial x^2} +\n",
+ "\\frac{\\partial^2 u}{\\partial y^2} + \\frac{\\partial^2 u}{\\partial z^2},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "in three space dimensions ($d=3$).\n",
+ "\n",
+ "Many applications involve variable coefficients, and the general\n",
+ "wave equation in $d$ dimensions is in this case written as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\varrho\\frac{\\partial^2 u}{\\partial t^2} = \\nabla\\cdot (q\\nabla u) + f\\hbox{ for } \\textbf{x} \\in\\Omega\\subset\\mathbb{R}^d,\\ t\\in (0,T],\n",
+ "\\label{wave:2D3D:model2} \\tag{2}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "which in, e.g., 2D becomes"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\varrho(x,y)\n",
+ "\\frac{\\partial^2 u}{\\partial t^2} =\n",
+ "\\frac{\\partial}{\\partial x}\\left( q(x,y)\n",
+ "\\frac{\\partial u}{\\partial x}\\right)\n",
+ "+\n",
+ "\\frac{\\partial}{\\partial y}\\left( q(x,y)\n",
+ "\\frac{\\partial u}{\\partial y}\\right)\n",
+ "+ f(x,y,t)\n",
+ "\\thinspace .\n",
+ "\\label{_auto1} \\tag{3}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "To save some writing and space we may use the index notation, where\n",
+ "subscript $t$, $x$, or $y$ means differentiation with respect\n",
+ "to that coordinate. For example,"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\begin{align*}\n",
+ "\\frac{\\partial^2 u}{\\partial t^2} &= u_{tt},\\\\\n",
+ "\\frac{\\partial}{\\partial y}\\left( q(x,y)\n",
+ "\\frac{\\partial u}{\\partial y}\\right) &= (q u_y)_y\n",
+ "\\thinspace .\n",
+ "\\end{align*}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "These comments extend straightforwardly to 3D, which means that\n",
+ "the 3D versions of the\n",
+ "two wave PDEs, with and without variable coefficients,\n",
+ "can be stated as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "u_{tt} = c^2(u_{xx} + u_{yy} + u_{zz}) + f,\n",
+ "\\label{wave:2D3D:model1:v2} \\tag{4}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation} \n",
+ "\\varrho u_{tt} = (q u_x)_x + (q u_y)_y + (q u_z)_z + f\\thinspace .\n",
+ "\\label{wave:2D3D:model2:v2} \\tag{5}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "At *each point* of the boundary $\\partial\\Omega$ (of $\\Omega$) we need\n",
+ "*one* boundary condition involving the unknown $u$.\n",
+ "The boundary conditions are of three principal types:\n",
+ "\n",
+ "1. $u$ is prescribed ($u=0$ or a known time variation of $u$ at\n",
+ " the boundary points, e.g.,\n",
+ " modeling an incoming wave),\n",
+ "\n",
+ "2. $\\partial u/\\partial n = \\boldsymbol{n}\\cdot\\nabla u$ is prescribed\n",
+ " (zero for reflecting boundaries),\n",
+ "\n",
+ "3. an open boundary condition (also called radiation condition)\n",
+ " is specified to let waves travel undisturbed out of the domain,\n",
+ " see [wave:app:exer:radiationBC](#wave:app:exer:radiationBC) for details.\n",
+ "\n",
+ "All the listed wave equations with *second-order* derivatives in\n",
+ "time need *two* initial conditions:\n",
+ "\n",
+ "1. $u=I$,\n",
+ "\n",
+ "2. $u_t = V$.\n",
+ "\n",
+ "## Mesh\n",
+ "\n",
+ "\n",
+ "We introduce a mesh in time and in space. The mesh in time consists\n",
+ "of time points"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "t_0=0 < t_1 < \\cdots < t_{N_t},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "normally, for wave equation problems, with a constant spacing $\\Delta\n",
+ "t= t_{n+1}-t_{n}$, $n\\in\\setl{\\mathcal{I}_t}$.\n",
+ "\n",
+ "Finite difference methods are easy to implement on simple rectangle-\n",
+ "or box-shaped *spatial domains*. More complicated shapes of the\n",
+ "spatial domain require substantially more advanced techniques and\n",
+ "implementational efforts (and a finite element method is usually a more\n",
+ "convenient approach). On a rectangle- or box-shaped domain, mesh\n",
+ "points are introduced separately in the various space directions:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\begin{align*}\n",
+ "&x_0 < x_1 < \\cdots < x_{N_x} \\hbox{ in the }x \\hbox{ direction},\\\\\n",
+ "&y_0 < y_1 < \\cdots < y_{N_y} \\hbox{ in the }y \\hbox{ direction},\\\\\n",
+ "&z_0 < z_1 < \\cdots < z_{N_z} \\hbox{ in the }z \\hbox{ direction}\\thinspace .\n",
+ "\\end{align*}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can write a general mesh point as $(x_i,y_j,z_k,t_n)$, with\n",
+ "$i\\in\\mathcal{I}_x$, $j\\in\\Iy$, $k\\in\\Iz$, and $n\\in\\mathcal{I}_t$.\n",
+ "\n",
+ "It is a very common choice to use constant mesh spacings:\n",
+ "$\\Delta x = x_{i+1}-x_{i}$, $i\\in\\setl{\\mathcal{I}_x}$,\n",
+ "$\\Delta y = y_{j+1}-y_{j}$, $j\\in\\setl{\\Iy}$, and\n",
+ "$\\Delta z = z_{k+1}-z_{k}$, $k\\in\\setl{\\Iz}$.\n",
+ "With equal mesh spacings one often introduces\n",
+ "$h = \\Delta x = \\Delta y =\\Delta z$.\n",
+ "\n",
+ "The unknown $u$ at mesh point $(x_i,y_j,z_k,t_n)$ is denoted by\n",
+ "$u^{n}_{i,j,k}$. In 2D problems we just skip the $z$ coordinate\n",
+ "(by assuming no variation in that direction: $\\partial/\\partial z=0$)\n",
+ "and write $u^n_{i,j}$.\n",
+ "\n",
+ "\n",
+ "## Discretization\n",
+ "\n",
+ "\n",
+ "Two- and three-dimensional wave equations are easily discretized by\n",
+ "assembling building blocks for discretization of\n",
+ "1D wave equations, because the multi-dimensional versions just contain\n",
+ "terms of the same type as those in 1D.\n",
+ "\n",
+ "### Discretizing the PDEs\n",
+ "\n",
+ "Equation ([4](#wave:2D3D:model1:v2)) can be discretized as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "[D_tD_t u = c^2(D_xD_x u + D_yD_yu + D_zD_z u) + f]^n_{i,j,k}\n",
+ "\\thinspace .\n",
+ "\\label{_auto2} \\tag{6}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "A 2D version might be instructive to write out in detail:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "[D_tD_t u = c^2(D_xD_x u + D_yD_yu) + f]^n_{i,j},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "which becomes"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\frac{u^{n+1}_{i,j} - 2u^{n}_{i,j} + u^{n-1}_{i,j}}{\\Delta t^2}\n",
+ "= c^2\n",
+ "\\frac{u^{n}_{i+1,j} - 2u^{n}_{i,j} + u^{n}_{i-1,j}}{\\Delta x^2}\n",
+ "+ c^2\n",
+ "\\frac{u^{n}_{i,j+1} - 2u^{n}_{i,j} + u^{n}_{i,j-1}}{\\Delta y^2}\n",
+ "+ f^n_{i,j},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Assuming, as usual, that all values at time levels $n$ and $n-1$\n",
+ "are known, we can solve for the only unknown $u^{n+1}_{i,j}$. The\n",
+ "result can be compactly written as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "u^{n+1}_{i,j} = 2u^n_{i,j} + u^{n-1}_{i,j} + c^2\\Delta t^2[D_xD_x u + D_yD_y u]^n_{i,j}\\thinspace .\n",
+ "\\label{wave:2D3D:models:unp1} \\tag{7}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "As in the 1D case, we need to develop a special formula for $u^1_{i,j}$\n",
+ "where we combine the general scheme for $u^{n+1}_{i,j}$, when $n=0$,\n",
+ "with the discretization of the initial condition:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "[D_{2t}u = V]^0_{i,j}\\quad\\Rightarrow\\quad u^{-1}_{i,j} = u^1_{i,j} - 2\\Delta t V_{i,j}\n",
+ "\\thinspace .\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The result becomes, in compact form,"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "u^{1}_{i,j} = u^0_{i,j} -2\\Delta V_{i,j} + {\\frac{1}{2}}\n",
+ "c^2\\Delta t^2[D_xD_x u + D_yD_y u]^0_{i,j}\\thinspace .\n",
+ "\\label{wave:2D3D:models:u1} \\tag{8}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The PDE ([5](#wave:2D3D:model2:v2))\n",
+ "with variable coefficients is discretized term by term using\n",
+ "the corresponding elements from the 1D case:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "[\\varrho D_tD_t u = (D_x\\overline{q}^x D_x u +\n",
+ "D_y\\overline{q}^y D_yu + D_z\\overline{q}^z D_z u) + f]^n_{i,j,k}\n",
+ "\\thinspace .\n",
+ "\\label{_auto3} \\tag{9}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "When written out and solved for the unknown $u^{n+1}_{i,j,k}$, one gets the\n",
+ "scheme"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\begin{align*}\n",
+ "u^{n+1}_{i,j,k} &= - u^{n-1}_{i,j,k} + 2u^{n}_{i,j,k} + \\\\\n",
+ "&\\quad \\frac{1}{\\varrho_{i,j,k}}\\frac{1}{\\Delta x^2} ( \\frac{1}{2}(q_{i,j,k} + q_{i+1,j,k})(u^{n}_{i+1,j,k} - u^{n}_{i,j,k}) - \\\\\n",
+ "&\\qquad\\qquad \\frac{1}{2}(q_{i-1,j,k} + q_{i,j,k})(u^{n}_{i,j,k} - u^{n}_{i-1,j,k})) + \\\\\n",
+ "&\\quad \\frac{1}{\\varrho_{i,j,k}}\\frac{1}{\\Delta y^2} ( \\frac{1}{2}(q_{i,j,k} + q_{i,j+1,k})(u^{n}_{i,j+1,k} - u^{n}_{i,j,k}) - \\\\\n",
+ "&\\qquad\\qquad\\frac{1}{2}(q_{i,j-1,k} + q_{i,j,k})(u^{n}_{i,j,k} - u^{n}_{i,j-1,k})) + \\\\\n",
+ "&\\quad \\frac{1}{\\varrho_{i,j,k}}\\frac{1}{\\Delta z^2} ( \\frac{1}{2}(q_{i,j,k} + q_{i,j,k+1})(u^{n}_{i,j,k+1} - u^{n}_{i,j,k}) -\\\\\n",
+ "&\\qquad\\qquad \\frac{1}{2}(q_{i,j,k-1} + q_{i,j,k})(u^{n}_{i,j,k} - u^{n}_{i,j,k-1})) + \\\\\n",
+ "&\\quad \\Delta t^2 f^n_{i,j,k}\n",
+ "\\thinspace .\n",
+ "\\end{align*}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Also here we need to develop a special formula for $u^1_{i,j,k}$\n",
+ "by combining the scheme for $n=0$ with the discrete initial condition,\n",
+ "which is just a matter of inserting\n",
+ "$u^{-1}_{i,j,k}=u^1_{i,j,k} - 2\\Delta tV_{i,j,k}$ in the scheme\n",
+ "and solving for $u^1_{i,j,k}$.\n",
+ "\n",
+ "### Handling boundary conditions where $u$ is known\n",
+ "\n",
+ "The schemes listed above are valid for the internal points in the mesh.\n",
+ "After updating these, we need to visit all the mesh points at the\n",
+ "boundaries and set the prescribed $u$ value.\n",
+ "\n",
+ "### Discretizing the Neumann condition\n",
+ "\n",
+ "The condition $\\partial u/\\partial n = 0$ was implemented in 1D by\n",
+ "discretizing it with a $D_{2x}u$ centered difference, followed by\n",
+ "eliminating the fictitious $u$ point outside the mesh by using the\n",
+ "general scheme at the boundary point. Alternatively, one can introduce\n",
+ "ghost cells and update a ghost value for use in the Neumann\n",
+ "condition. Exactly the same ideas are reused in multiple dimensions.\n",
+ "\n",
+ "Consider the condition $\\partial u/\\partial n = 0$\n",
+ "at a boundary $y=0$ of a rectangular domain $[0, L_x]\\times [0,L_y]$ in 2D.\n",
+ "The normal direction is then in $-y$ direction,\n",
+ "so"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\frac{\\partial u}{\\partial n} = -\\frac{\\partial u}{\\partial y},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "and we set"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "[-D_{2y} u = 0]^n_{i,0}\\quad\\Rightarrow\\quad \\frac{u^n_{i,1}-u^n_{i,-1}}{2\\Delta y} = 0\n",
+ "\\thinspace .\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "From this it follows that $u^n_{i,-1}=u^n_{i,1}$.\n",
+ "The discretized PDE at the boundary point $(i,0)$ reads"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\frac{u^{n+1}_{i,0} - 2u^{n}_{i,0} + u^{n-1}_{i,0}}{\\Delta t^2}\n",
+ "= c^2\n",
+ "\\frac{u^{n}_{i+1,0} - 2u^{n}_{i,0} + u^{n}_{i-1,0}}{\\Delta x^2}\n",
+ "+ c^2\n",
+ "\\frac{u^{n}_{i,1} - 2u^{n}_{i,0} + u^{n}_{i,-1}}{\\Delta y^2}\n",
+ "+ f^n_{i,j},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can then just insert $u^n_{i,1}$ for $u^n_{i,-1}$ in this equation\n",
+ "and solve for the boundary value $u^{n+1}_{i,0}$, just as was done in 1D.\n",
+ "\n",
+ "From these calculations, we see a pattern:\n",
+ "the general scheme applies at the boundary $j=0$ too if we just\n",
+ "replace $j-1$ by $j+1$. Such a pattern is particularly useful for\n",
+ "implementations. The details follow from the explained 1D case\n",
+ "in the section [wave:pde2:Neumann:impl](#wave:pde2:Neumann:impl).\n",
+ "\n",
+ "The alternative approach to eliminating fictitious values outside the\n",
+ "mesh is to have $u^n_{i,-1}$ available as a ghost value. The mesh is\n",
+ "extended with one extra line (2D) or plane (3D) of ghost cells at a\n",
+ "Neumann boundary. In the present example it means that we need a line\n",
+ "with ghost cells below the $y$ axis. The ghost values must be updated\n",
+ "according to $u^{n+1}_{i,-1}=u^{n+1}_{i,1}$."
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.6"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/fdm-jupyter-book/notebooks/02_wave/wave_analysis.ipynb b/fdm-jupyter-book/notebooks/02_wave/wave_analysis.ipynb
new file mode 100644
index 00000000..50e943b3
--- /dev/null
+++ b/fdm-jupyter-book/notebooks/02_wave/wave_analysis.ipynb
@@ -0,0 +1,1567 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "# Analysis of the difference equations\n",
+ "\n",
+ "\n",
+ "## Properties of the solution of the wave equation\n",
+ "\n",
+ "\n",
+ "The wave equation"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\frac{\\partial^2 u}{\\partial t^2} =\n",
+ "c^2 \\frac{\\partial^2 u}{\\partial x^2}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "has solutions of the form"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "u(x,t) = g_R(x-ct) + g_L(x+ct),\n",
+ "\\label{wave:pde1:gensol} \\tag{1}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "for any functions $g_R$ and $g_L$ sufficiently smooth to be differentiated\n",
+ "twice. The result follows from inserting ([1](#wave:pde1:gensol))\n",
+ "in the wave equation. A function of the form $g_R(x-ct)$ represents a\n",
+ "signal\n",
+ "moving to the right in time with constant velocity $c$.\n",
+ "This feature can be explained as follows.\n",
+ "At time $t=0$ the signal looks like $g_R(x)$. Introducing a\n",
+ "moving horizontal coordinate $\\xi = x-ct$, we see the function\n",
+ "$g_R(\\xi)$ is \"at rest\"\n",
+ "in the $\\xi$ coordinate system, and the shape is always\n",
+ "the same. Say the $g_R(\\xi)$ function has a peak at $\\xi=0$. This peak\n",
+ "is located at $x=ct$, which means that it moves with the velocity\n",
+ "$dx/dt=c$ in the $x$ coordinate system. Similarly, $g_L(x+ct)$\n",
+ "is a function, initially with shape $g_L(x)$, that moves in the negative\n",
+ "$x$ direction with constant velocity $c$ (introduce $\\xi=x+ct$,\n",
+ "look at the point $\\xi=0$, $x=-ct$, which has velocity $dx/dt=-c$).\n",
+ "\n",
+ "With the particular initial conditions"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "u(x,0)=I(x),\\quad \\frac{\\partial}{\\partial t}u(x,0) =0,\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "we get, with $u$ as in ([1](#wave:pde1:gensol)),"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "g_R(x) + g_L(x) = I(x),\\quad -cg_R'(x) + cg_L'(x) = 0\\thinspace .\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The former suggests $g_R=g_L$, and the former then leads to\n",
+ "$g_R=g_L=I/2$. Consequently,"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "u(x,t) = \\frac{1}{2} I(x-ct) + \\frac{1}{2} I(x+ct) \\thinspace .\n",
+ "\\label{wave:pde1:gensol2} \\tag{2}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The interpretation of ([2](#wave:pde1:gensol2)) is that\n",
+ "the initial shape of $u$ is split into two parts, each with the same\n",
+ "shape as $I$ but frac{1}{2}\n",
+ "of the initial amplitude. One part is traveling to the left and the\n",
+ "other one to the right.\n",
+ "\n",
+ "\n",
+ "The solution has two important physical features: constant amplitude\n",
+ "of the left and right wave, and constant velocity of these two waves.\n",
+ "mathcal{I}_t turns out that the numerical solution will also preserve the\n",
+ "constant amplitude, but the velocity depends on the mesh parameters\n",
+ "$\\Delta t$ and $\\Delta x$.\n",
+ "\n",
+ "The solution ([2](#wave:pde1:gensol2)) will be influenced by\n",
+ "boundary conditions when the parts\n",
+ "$\\frac{1}{2} I(x-ct)$ and $\\frac{1}{2} I(x+ct)$ hit the boundaries and get, e.g.,\n",
+ "reflected back into the domain. However, when $I(x)$ is nonzero\n",
+ "only in a small part in the middle\n",
+ "of the spatial domain $[0,L]$, which means that the\n",
+ "boundaries are placed far away from the initial disturbance of $u$,\n",
+ "the solution ([2](#wave:pde1:gensol2)) is very clearly observed\n",
+ "in a simulation.\n",
+ "\n",
+ "\n",
+ "\n",
+ "A useful representation of solutions of wave equations is a linear\n",
+ "combination of sine and/or cosine waves. Such a sum of waves is a\n",
+ "solution if the governing PDE is linear and each sine or cosine\n",
+ "wave fulfills the\n",
+ "equation. To ease analytical calculations by hand we shall work with\n",
+ "complex exponential functions instead of real-valued sine or cosine\n",
+ "functions. The real part of complex expressions will typically be\n",
+ "taken as the physical relevant quantity (whenever a physical relevant\n",
+ "quantity is strictly needed).\n",
+ "The idea now is to build $I(x)$ of complex wave components\n",
+ "$e^{ikx}$:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation} I(x) \\approx \\sum_{k\\in K} b_k e^{ikx} \\thinspace .\n",
+ "\\label{wave:Fourier:I} \\tag{3}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Here, $k$ is the frequency of a component,\n",
+ "$K$ is some set of all the discrete\n",
+ "$k$ values needed to approximate $I(x)$ well,\n",
+ "and $b_k$ are\n",
+ "constants that must be determined. We will very seldom\n",
+ "need to compute the $b_k$ coefficients: most of the insight\n",
+ "we look for, and the understanding of the numerical methods we want to\n",
+ "establish, come from\n",
+ "investigating how the PDE and the scheme treat a single\n",
+ "component $e^{ikx}$ wave.\n",
+ "\n",
+ "Letting the number of $k$ values in $K$ tend to infinity, makes the sum\n",
+ "([3](#wave:Fourier:I)) converge to $I(x)$. This sum is known as a\n",
+ "*Fourier series* representation of $I(x)$. Looking at\n",
+ "([2](#wave:pde1:gensol2)), we see that the solution $u(x,t)$, when\n",
+ "$I(x)$ is represented as in ([3](#wave:Fourier:I)), is also built of\n",
+ "basic complex exponential wave components of the form $e^{ik(x\\pm\n",
+ "ct)}$ according to"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "u(x,t) = \\frac{1}{2} \\sum_{k\\in K} b_k e^{ik(x - ct)}\n",
+ "+ \\frac{1}{2} \\sum_{k\\in K} b_k e^{ik(x + ct)} \\thinspace .\n",
+ "\\label{wave:Fourier:u1} \\tag{4}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "mathcal{I}_t is common to introduce the frequency in time $\\omega = kc$ and\n",
+ "assume that $u(x,t)$ is a sum of basic wave components\n",
+ "written as $e^{ikx -\\omega t}$.\n",
+ "(Observe that inserting such a wave component in the governing PDE reveals that\n",
+ "$\\omega^2 = k^2c^2$, or $\\omega =\\pm kc$, reflecting the\n",
+ "two solutions: one ($+kc$) traveling to the right and the other ($-kc$)\n",
+ "traveling to the left.)\n",
+ "\n",
+ "## More precise definition of Fourier representations\n",
+ "\n",
+ "\n",
+ "The above introduction to function representation by sine and cosine\n",
+ "waves was quick and intuitive, but will suffice as background\n",
+ "knowledge for the following material of single wave component\n",
+ "analysis.\n",
+ "However, to understand\n",
+ "all details of how different wave components sum up to the analytical\n",
+ "and numerical solutions, a more precise mathematical treatment is helpful\n",
+ "and therefore summarized below.\n",
+ "\n",
+ "mathcal{I}_t is well known that periodic functions can be represented by\n",
+ "Fourier series. A generalization of the Fourier series idea to\n",
+ "non-periodic functions defined on the real line is the *Fourier transform*:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "I(x) = \\int_{-\\infty}^\\infty A(k)e^{ikx}dk,\n",
+ "\\label{wave:pde1:Fourier:I} \\tag{5} \n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation} \n",
+ "A(k) = \\int_{-\\infty}^\\infty I(x)e^{-ikx}dx\\thinspace .\n",
+ "\\label{wave:pde1:Fourier:A} \\tag{6}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The function $A(k)$ reflects the weight of each wave component $e^{ikx}$\n",
+ "in an infinite sum of such wave components. That is, $A(k)$\n",
+ "reflects the frequency content in the function $I(x)$. Fourier transforms\n",
+ "are particularly fundamental for analyzing and understanding time-varying\n",
+ "signals.\n",
+ "\n",
+ "The solution of the linear 1D wave PDE can be expressed as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "u(x,t) = \\int_{-\\infty}^\\infty A(k)e^{i(kx-\\omega(k)t)}dx\\thinspace .\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In a finite difference method, we represent $u$ by a mesh function\n",
+ "$u^n_q$, where $n$ counts temporal mesh points and $q$ counts\n",
+ "the spatial ones (the usual counter for spatial points, $i$, is\n",
+ "here already used as imaginary unit). Similarly, $I(x)$ is approximated by\n",
+ "the mesh function $I_q$, $q=0,\\ldots,N_x$.\n",
+ "On a mesh, it does not make sense to work with wave\n",
+ "components $e^{ikx}$ for very large $k$, because the shortest possible\n",
+ "sine or cosine wave that can be represented uniquely\n",
+ "on a mesh with spacing $\\Delta x$\n",
+ "is the wave with wavelength $2\\Delta x$. This wave has its peaks\n",
+ "and throughs at every two mesh points. That is, the wave \"jumps up and down\"\n",
+ "between the mesh points.\n",
+ "\n",
+ "The corresponding $k$ value for the shortest possible wave in the mesh is\n",
+ "$k=2\\pi /(2\\Delta x) = \\pi/\\Delta x$. This maximum frequency is\n",
+ "known as the *Nyquist frequency*.\n",
+ "Within the range of\n",
+ "relevant frequencies $(0,\\pi/\\Delta x]$ one defines\n",
+ "the [discrete Fourier transform](http://en.wikipedia.org/wiki/Discrete_Fourier_transform), using $N_x+1$ discrete frequencies:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "I_q = \\frac{1}{N_x+1}\\sum_{k=0}^{N_x} A_k e^{i2\\pi k q/(N_x+1)},\\quad\n",
+ "q=0,\\ldots,N_x,\n",
+ "\\label{_auto1} \\tag{7}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation} \n",
+ "A_k = \\sum_{q=0}^{N_x} I_q e^{-i2\\pi k q/(N_x+1)},\n",
+ "\\quad k=0,\\ldots,N_x\\thinspace .\n",
+ "\\label{_auto2} \\tag{8}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The $A_k$ values represent the discrete Fourier transform of the $I_q$ values,\n",
+ "which themselves are the inverse discrete Fourier transform of the $A_k$\n",
+ "values.\n",
+ "\n",
+ "The discrete Fourier transform is efficiently computed by the\n",
+ "*Fast Fourier transform* algorithm. For a real function $I(x)$,\n",
+ "the relevant Python code for computing and plotting\n",
+ "the discrete Fourier transform appears in the example below."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%matplotlib inline\n",
+ "\n",
+ "import numpy as np\n",
+ "from numpy import sin, pi\n",
+ "\n",
+ "def I(x):\n",
+ " return sin(2*pi*x) + 0.5*sin(4*pi*x) + 0.1*sin(6*pi*x)\n",
+ "\n",
+ "# Mesh\n",
+ "L = 10; Nx = 100\n",
+ "x = np.linspace(0, L, Nx+1)\n",
+ "dx = L/float(Nx)\n",
+ "\n",
+ "# Discrete Fourier transform\n",
+ "A = np.fft.rfft(I(x))\n",
+ "A_amplitude = np.abs(A)\n",
+ "\n",
+ "# Compute the corresponding frequencies\n",
+ "freqs = np.linspace(0, pi/dx, A_amplitude.size)\n",
+ "\n",
+ "import matplotlib.pyplot as plt\n",
+ "plt.plot(freqs, A_amplitude)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Stability\n",
+ "\n",
+ "\n",
+ "\n",
+ "The scheme"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "[D_tD_t u = c^2 D_xD_x u]^n_q\n",
+ "\\label{wave:pde1:analysis:scheme} \\tag{9}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "for the wave equation $u_{tt} = c^2u_{xx}$ allows basic wave components"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "u^n_q=e^{i(kx_q - \\tilde\\omega t_n)}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "as solution, but it turns out that\n",
+ "the frequency in time, $\\tilde\\omega$, is not equal to\n",
+ "the exact frequency $\\omega = kc$. The goal now is to\n",
+ "find exactly what $\\tilde \\omega$ is. We ask two key\n",
+ "questions:\n",
+ "\n",
+ " * How accurate is $\\tilde\\omega$\n",
+ " compared to $\\omega$?\n",
+ "\n",
+ " * Does the amplitude of such a wave component\n",
+ " preserve its (unit) amplitude, as it should,\n",
+ " or does it get amplified or damped in time (because of a complex $\\tilde\\omega$)?\n",
+ "\n",
+ "The following analysis will answer these questions. We shall\n",
+ "continue using $q$ as an identifier for a certain mesh point in\n",
+ "the $x$ direction.\n",
+ "\n",
+ "\n",
+ "### Preliminary results\n",
+ "\n",
+ "A key result needed in the investigations is the finite difference\n",
+ "approximation of a second-order derivative acting on a complex\n",
+ "wave component:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "[D_tD_t e^{i\\omega t}]^n = -\\frac{4}{\\Delta t^2}\\sin^2\\left(\n",
+ "\\frac{\\omega\\Delta t}{2}\\right)e^{i\\omega n\\Delta t}\n",
+ "\\thinspace .\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "By just changing symbols ($\\omega\\rightarrow k$,\n",
+ "$t\\rightarrow x$, $n\\rightarrow q$) it follows that"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "[D_xD_x e^{ikx}]_q = -\\frac{4}{\\Delta x^2}\\sin^2\\left(\n",
+ "\\frac{k\\Delta x}{2}\\right)e^{ikq\\Delta x} \\thinspace .\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Numerical wave propagation\n",
+ "\n",
+ "Inserting a basic wave component $u^n_q=e^{i(kx_q-\\tilde\\omega t_n)}$ in\n",
+ "([9](#wave:pde1:analysis:scheme)) results in the need to\n",
+ "evaluate two expressions:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\lbrack D_tD_t e^{ikx}e^{-i\\tilde\\omega t}\\rbrack^n_q = \\lbrack D_tD_t e^{-i\\tilde\\omega t}\\rbrack^ne^{ikq\\Delta x}\\nonumber\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation} = -\\frac{4}{\\Delta t^2}\\sin^2\\left(\n",
+ "\\frac{\\tilde\\omega\\Delta t}{2}\\right)e^{-i\\tilde\\omega n\\Delta t}e^{ikq\\Delta x}\n",
+ "\\label{_auto3} \\tag{10}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\lbrack D_xD_x e^{ikx}e^{-i\\tilde\\omega t}\\rbrack^n_q = \\lbrack D_xD_x e^{ikx}\\rbrack_q e^{-i\\tilde\\omega n\\Delta t}\\nonumber\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation} = -\\frac{4}{\\Delta x^2}\\sin^2\\left(\n",
+ "\\frac{k\\Delta x}{2}\\right)e^{ikq\\Delta x}e^{-i\\tilde\\omega n\\Delta t} \\thinspace . \\label{_auto4} \\tag{11}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Then the complete scheme,"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\lbrack D_tD_t e^{ikx}e^{-i\\tilde\\omega t} = c^2D_xD_x e^{ikx}e^{-i\\tilde\\omega t}\\rbrack^n_q\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "leads to the following equation for the unknown numerical\n",
+ "frequency $\\tilde\\omega$\n",
+ "(after dividing by $-e^{ikx}e^{-i\\tilde\\omega t}$):"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\frac{4}{\\Delta t^2}\\sin^2\\left(\\frac{\\tilde\\omega\\Delta t}{2}\\right)\n",
+ "= c^2 \\frac{4}{\\Delta x^2}\\sin^2\\left(\\frac{k\\Delta x}{2}\\right),\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "or"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\sin^2\\left(\\frac{\\tilde\\omega\\Delta t}{2}\\right)\n",
+ "= C^2\\sin^2\\left(\\frac{k\\Delta x}{2}\\right),\n",
+ "\\label{wave:pde1:analysis:sineq1} \\tag{12}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "C = \\frac{c\\Delta t}{\\Delta x}\n",
+ "\\label{_auto5} \\tag{13}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "is the Courant number.\n",
+ "Taking the square root of ([12](#wave:pde1:analysis:sineq1)) yields"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\sin\\left(\\frac{\\tilde\\omega\\Delta t}{2}\\right)\n",
+ "= C\\sin\\left(\\frac{k\\Delta x}{2}\\right),\n",
+ "\\label{wave:pde1:analysis:sineq2} \\tag{14}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Since the exact $\\omega$ is real it is reasonable to look for a real\n",
+ "solution $\\tilde\\omega$ of ([14](#wave:pde1:analysis:sineq2)).\n",
+ "The right-hand side of\n",
+ "([14](#wave:pde1:analysis:sineq2)) must then be in $[-1,1]$ because\n",
+ "the sine function on the left-hand side has values in $[-1,1]$\n",
+ "for real $\\tilde\\omega$. The magnitude of the sine function on\n",
+ "the right-hand side attains the value 1 when"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\frac{k\\Delta x}{2} = \\frac{\\pi}{2} + m\\pi,\\quad m \\in\\mathbb{Z}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "With $m=0$ we have $k\\Delta x = \\pi$, which means that\n",
+ "the wavelength $\\lambda = 2\\pi/k$ becomes $2\\Delta x$. This is\n",
+ "the absolutely shortest wavelength that can be represented on the mesh:\n",
+ "the wave jumps up and down between each mesh point. Larger values of $|m|$\n",
+ "are irrelevant since these correspond to $k$ values whose\n",
+ "waves are too short to be represented\n",
+ "on a mesh with spacing $\\Delta x$.\n",
+ "For the shortest possible wave in the mesh, $\\sin\\left(k\\Delta x/2\\right)=1$,\n",
+ "and we must require"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "C\\leq 1 \\thinspace .\n",
+ "\\label{wave:pde1:stability:crit} \\tag{15}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Consider a right-hand side in ([14](#wave:pde1:analysis:sineq2)) of\n",
+ "magnitude larger\n",
+ "than unity. The solution $\\tilde\\omega$ of ([14](#wave:pde1:analysis:sineq2))\n",
+ "must then be a complex number\n",
+ "$\\tilde\\omega = \\tilde\\omega_r + i\\tilde\\omega_i$ because\n",
+ "the sine function is larger than unity for a complex argument.\n",
+ "One can show that for any $\\omega_i$ there will also be a\n",
+ "corresponding solution with $-\\omega_i$. The component with $\\omega_i>0$\n",
+ "gives an amplification factor $e^{\\omega_it}$ that grows exponentially\n",
+ "in time. We cannot allow this and must therefore require $C\\leq 1$\n",
+ "as a *stability criterion*.\n",
+ "\n",
+ "**Remark on the stability requirement.**\n",
+ "\n",
+ "For smoother wave components with longer wave lengths per length $\\Delta x$,\n",
+ "([15](#wave:pde1:stability:crit)) can in theory be relaxed. However,\n",
+ "small round-off errors are always present in a numerical solution and these\n",
+ "vary arbitrarily from mesh point to mesh point and can be viewed as\n",
+ "unavoidable noise with wavelength $2\\Delta x$. As explained, $C>1$\n",
+ "will for this very small noise lead to exponential growth of\n",
+ "the shortest possible wave component in the mesh. This noise will\n",
+ "therefore grow with time and destroy the whole solution.\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Numerical dispersion relation\n",
+ "\n",
+ "\n",
+ "\n",
+ "Equation ([14](#wave:pde1:analysis:sineq2)) can be solved with respect\n",
+ "to $\\tilde\\omega$:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\tilde\\omega = \\frac{2}{\\Delta t}\n",
+ "\\sin^{-1}\\left( C\\sin\\left(\\frac{k\\Delta x}{2}\\right)\\right) \\thinspace .\n",
+ "\\label{wave:pde1:disprel} \\tag{16}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The relation between the numerical frequency $\\tilde\\omega$ and\n",
+ "the other parameters $k$, $c$, $\\Delta x$, and $\\Delta t$ is called\n",
+ "a *numerical dispersion relation*. Correspondingly,\n",
+ "$\\omega =kc$ is the *analytical dispersion relation*.\n",
+ "In general, dispersion refers to the phenomenon where the wave\n",
+ "velocity depends on the spatial frequency ($k$, or the\n",
+ "wave length $\\lambda = 2\\pi/k$) of the wave.\n",
+ "Since the wave velocity is $\\omega/k =c$, we realize that the\n",
+ "analytical dispersion relation reflects the fact that there is no\n",
+ "dispersion. However, in a numerical scheme we have dispersive waves\n",
+ "where the wave velocity depends on $k$.\n",
+ "\n",
+ "The special case $C=1$ deserves attention since then the right-hand side\n",
+ "of ([16](#wave:pde1:disprel)) reduces to"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\frac{2}{\\Delta t}\\frac{k\\Delta x}{2} = \\frac{1}{\\Delta t}\n",
+ "\\frac{\\omega\\Delta x}{c} = \\frac{\\omega}{C} = \\omega \\thinspace .\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "That is, $\\tilde\\omega = \\omega$ and the numerical solution is exact\n",
+ "at all mesh points regardless of $\\Delta x$ and $\\Delta t$!\n",
+ "This implies that the numerical solution method is also an analytical\n",
+ "solution method, at least for computing $u$ at discrete points (the\n",
+ "numerical method says nothing about the\n",
+ "variation of $u$ *between* the mesh points, and employing the\n",
+ "common linear interpolation for extending the discrete solution\n",
+ "gives a curve that in general deviates from the exact one).\n",
+ "\n",
+ "For a closer examination of the error in the numerical dispersion\n",
+ "relation when $C<1$, we can study\n",
+ "$\\tilde\\omega -\\omega$, $\\tilde\\omega/\\omega$, or the similar\n",
+ "error measures in wave velocity: $\\tilde c - c$ and $\\tilde c/c$,\n",
+ "where $c=\\omega /k$ and $\\tilde c = \\tilde\\omega /k$.\n",
+ "mathcal{I}_t appears that the most convenient expression to work with is $\\tilde c/c$,\n",
+ "since it can be written as a function of just two parameters:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\frac{\\tilde c}{c} = \\frac{1}{Cp}{\\sin}^{-1}\\left(C\\sin p\\right),\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "with $p=k\\Delta x/2$ as a non-dimensional measure of the spatial frequency.\n",
+ "In essence, $p$ tells how many spatial mesh points we have per\n",
+ "wave length in space for the wave component with frequency $k$ (recall\n",
+ "that the wave\n",
+ "length is $2\\pi/k$). That is, $p$ reflects how well the\n",
+ "spatial variation of the wave component\n",
+ "is resolved in the mesh. Wave components with wave length\n",
+ "less than $2\\Delta x$ ($2\\pi/k < 2\\Delta x$) are not visible in the mesh,\n",
+ "so it does not make sense to have $p>\\pi/2$.\n",
+ "\n",
+ "We may introduce the function $r(C, p)=\\tilde c/c$ for further investigation\n",
+ "of numerical errors in the wave velocity:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "r(C, p) = \\frac{1}{Cp}{\\sin}^{-1}\\left(C\\sin p\\right), \\quad C\\in (0,1],\\ p\\in (0,\\pi/2] \\thinspace .\n",
+ "\\label{wave:pde1:disprel2} \\tag{17}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This function is very well suited for plotting since it combines several\n",
+ "parameters in the problem into a dependence on two dimensionless\n",
+ "numbers, $C$ and $p$.\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "The fractional error in the wave velocity for different Courant numbers.
\n", + "![](\"fig-wave/disprel.png\")
\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "Defining"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ " def r(C, p):\n",
+ " return 2/(C*p)*asin(C*sin(p))\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "we can plot $r(C,p)$ as a function of $p$ for various values of\n",
+ "$C$, see [Figure](#wave:pde1:fig:disprel). Note that the shortest\n",
+ "waves have the most erroneous velocity, and that short waves move\n",
+ "more slowly than they should.\n",
+ "\n",
+ "\n",
+ "We can also easily make a Taylor series expansion in the\n",
+ "discretization parameter $p$:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import sympy as sym\n",
+ "C, p = sym.symbols('C p')\n",
+ "# Compute the 7 first terms around p=0 with no O() term\n",
+ "rs = r(C, p).series(p, 0, 7).removeO()\n",
+ "rs"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Pick out the leading order term, but drop the constant 1\n",
+ "rs_error_leading_order = (rs - 1).extract_leading_order(p)\n",
+ "rs_error_leading_order"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Turn the series expansion into a Python function\n",
+ "rs_pyfunc = lambdify([C, p], rs, modules='numpy')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Check: rs_pyfunc is exact (=1) for C=1\n",
+ "rs_pyfunc(1, 0.1)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Note that without the `.removeO()` call the series gets an `O(x**7)` term\n",
+ "that makes it impossible to convert the series to a Python function\n",
+ "(for, e.g., plotting).\n",
+ "\n",
+ "From the `rs_error_leading_order` expression above, we see that the leading\n",
+ "order term in the error of this series expansion is"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\frac{1}{6}\\left(\\frac{k\\Delta x}{2}\\right)^2(C^2-1)\n",
+ "= \\frac{k^2}{24}\\left( c^2\\Delta t^2 - \\Delta x^2\\right),\n",
+ "\\label{_auto6} \\tag{18}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "pointing to an error $\\Oof{\\Delta t^2, \\Delta x^2}$, which is\n",
+ "compatible with the errors in the difference approximations ($D_tD_tu$\n",
+ "and $D_xD_xu$).\n",
+ "\n",
+ "We can do more with a series expansion, e.g., factor it to see how\n",
+ "the factor $C-1$ plays a significant role.\n",
+ "To this end, we make a list of the terms, factor each term,\n",
+ "and then sum the terms:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "rs = r(C, p).series(p, 0, 4).removeO().as_ordered_terms()\n",
+ "rs"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "rs = [factor(t) for t in rs]\n",
+ "rs"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "rs = sum(rs) # Python's sum function sums the list\n",
+ "rs"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We see from the last expression\n",
+ "that $C=1$ makes all the terms in `rs` vanish.\n",
+ "Since we already know that the numerical solution is exact for $C=1$, the\n",
+ "remaining terms in the Taylor series expansion\n",
+ "will also contain factors of $C-1$ and cancel for $C=1$.\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Extending the analysis to 2D and 3D\n",
+ "\n",
+ "\n",
+ "The typical analytical solution of a 2D wave equation"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "u_{tt} = c^2(u_{xx} + u_{yy}),\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "is a wave traveling in the direction of $\\kk = k_x\\ii + k_y\\jj$, where\n",
+ "$\\ii$ and $\\jj$ are unit vectors in the $x$ and $y$ directions, respectively\n",
+ "($\\ii$ should not be confused with $i=\\sqrt{-1}$ here).\n",
+ "Such a wave can be expressed by"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "u(x,y,t) = g(k_xx + k_yy - kct)\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "for some twice differentiable function $g$, or with $\\omega =kc$, $k=|\\kk|$:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "u(x,y,t) = g(k_xx + k_yy - \\omega t)\\thinspace .\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can, in particular, build a solution by adding complex Fourier components\n",
+ "of the form"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "e^{(i(k_xx + k_yy - \\omega t))}\n",
+ "\\thinspace .\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "A discrete 2D wave equation can be written as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\lbrack D_tD_t u = c^2(D_xD_x u + D_yD_y u)\\rbrack^n_{q,r}\n",
+ "\\thinspace .\n",
+ "\\label{wave:pde1:analysis:scheme2D} \\tag{19}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This equation admits a Fourier component"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "u^n_{q,r} = e^{\\left( i(k_x q\\Delta x + k_y r\\Delta y -\n",
+ "\\tilde\\omega n\\Delta t)\\right)},\n",
+ "\\label{wave:pde1:analysis:numsol2D} \\tag{20}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "as solution. Letting the operators $D_tD_t$, $D_xD_x$, and $D_yD_y$\n",
+ "act on $u^n_{q,r}$ from ([20](#wave:pde1:analysis:numsol2D)) transforms\n",
+ "([19](#wave:pde1:analysis:scheme2D)) to"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\frac{4}{\\Delta t^2}\\sin^2\\left(\\frac{\\tilde\\omega\\Delta t}{2}\\right)\n",
+ "= c^2 \\frac{4}{\\Delta x^2}\\sin^2\\left(\\frac{k_x\\Delta x}{2}\\right)\n",
+ "+ c^2 \\frac{4}{\\Delta y^2}\\sin^2\\left(\\frac{k_y\\Delta y}{2}\\right) \\thinspace . \\label{_auto7} \\tag{21}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "or"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\sin^2\\left(\\frac{\\tilde\\omega\\Delta t}{2}\\right)\n",
+ "= C_x^2\\sin^2 p_x\n",
+ "+ C_y^2\\sin^2 p_y, \\label{_auto8} \\tag{22}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where we have eliminated the factor 4 and introduced the symbols"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "C_x = \\frac{c\\Delta t}{\\Delta x},\\quad\n",
+ "C_y = \\frac{c\\Delta t}{\\Delta y}, \\quad\n",
+ "p_x = \\frac{k_x\\Delta x}{2},\\quad\n",
+ "p_y = \\frac{k_y\\Delta y}{2}\\thinspace .\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "For a real-valued $\\tilde\\omega$ the right-hand side\n",
+ "must be less than or equal to unity in absolute value, requiring in general\n",
+ "that"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "C_x^2 + C_y^2 \\leq 1 \\thinspace .\n",
+ "\\label{wave:pde1:analysis:2DstabC} \\tag{23}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This gives the stability criterion, more commonly expressed directly\n",
+ "in an inequality for the time step:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\Delta t \\leq \\frac{1}{c} \\left( \\frac{1}{\\Delta x^2} +\n",
+ "\\frac{1}{\\Delta y^2}\\right)^{-\\frac{1}{2}i}\n",
+ "\\label{wave:pde1:analysis:2Dstab} \\tag{24}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "A similar, straightforward analysis for the 3D case leads to"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\Delta t \\leq \\frac{1}{c}\\left( \\frac{1}{\\Delta x^2} +\n",
+ "\\frac{1}{\\Delta y^2} + \\frac{1}{\\Delta z^2}\\right)^{-\\frac{1}{2}i}\n",
+ "\\label{_auto9} \\tag{25}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In the case of a variable coefficient $c^2=c^2(\\xpoint)$, we must use\n",
+ "the worst-case value"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\bar c = \\sqrt{\\max_{\\xpoint\\in\\Omega} c^2(\\xpoint)}\n",
+ "\\label{_auto10} \\tag{26}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "in the stability criteria. Often, especially in the variable wave\n",
+ "velocity case, it is wise to introduce a safety factor $\\beta\\in (0,1]$ too:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\Delta t \\leq \\beta \\frac{1}{\\bar c}\n",
+ "\\left( \\frac{1}{\\Delta x^2} +\n",
+ "\\frac{1}{\\Delta y^2} + \\frac{1}{\\Delta z^2}\\right)^{-\\frac{1}{2}i}\n",
+ "\\label{_auto11} \\tag{27}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The exact numerical dispersion relations in 2D and 3D becomes, for constant $c$,"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation}\n",
+ "\\tilde\\omega = \\frac{2}{\\Delta t}\\sin^{-1}\\left(\n",
+ "\\left( C_x^2\\sin^2 p_x + C_y^2\\sin^2 p_y\\right)^\\frac{1}{2}\\right),\n",
+ "\\label{_auto12} \\tag{28}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation} \n",
+ "\\tilde\\omega = \\frac{2}{\\Delta t}\\sin^{-1}\\left(\n",
+ "\\left( C_x^2\\sin^2 p_x + C_y^2\\sin^2 p_y + C_z^2\\sin^2 p_z\\right)^\\frac{1}{2}\\right)\\thinspace .\n",
+ "\\label{_auto13} \\tag{29}\n",
+ "\\end{equation}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can visualize the numerical dispersion error in 2D much like we did\n",
+ "in 1D. To this end, we need to reduce the number of parameters in\n",
+ "$\\tilde\\omega$. The direction of the wave is parameterized by the\n",
+ "polar angle $\\theta$, which means that"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "k_x = k\\sin\\theta,\\quad k_y=k\\cos\\theta\\thinspace .\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "A simplification is to set $\\Delta x=\\Delta y=h$.\n",
+ "Then $C_x=C_y=c\\Delta t/h$, which we call $C$. Also,"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "p_x=\\frac{1}{2} kh\\cos\\theta,\\quad p_y=\\frac{1}{2} kh\\sin\\theta\\thinspace .\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The numerical frequency $\\tilde\\omega$\n",
+ "is now a function of three parameters:\n",
+ "\n",
+ " * $C$, reflecting the number of cells a wave is displaced during a time step,\n",
+ "\n",
+ " * $p=\\frac{1}{2} kh$, reflecting the number of cells per wave length in space,\n",
+ "\n",
+ " * $\\theta$, expressing the direction of the wave.\n",
+ "\n",
+ "We want to visualize the error in the numerical frequency. To avoid having\n",
+ "$\\Delta t$ as a free parameter in $\\tilde\\omega$, we work with\n",
+ "$\\tilde c/c = \\tilde\\omega/(kc)$. The coefficient in front of the\n",
+ "$\\sin^{-1}$ factor is then"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\frac{2}{kc\\Delta t} = \\frac{2}{2kc\\Delta t h/h} =\n",
+ "\\frac{1}{Ckh} = \\frac{2}{Cp},\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "and"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\frac{\\tilde c}{c} = \\frac{2}{Cp}\n",
+ "\\sin^{-1}\\left(C\\left(\\sin^2 (p\\cos\\theta)\n",
+ "+ \\sin^2(p\\sin\\theta) \\right)^\\frac{1}{2}\\right)\\thinspace .\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We want to visualize this quantity as a function of\n",
+ "$p$ and $\\theta$ for some values of $C\\leq 1$. mathcal{I}_t is\n",
+ "instructive\n",
+ "to make color contour plots of $1-\\tilde c/c$ in\n",
+ "*polar coordinates* with $\\theta$ as the angular coordinate and\n",
+ "$p$ as the radial coordinate.\n",
+ "\n",
+ "The stability criterion ([23](#wave:pde1:analysis:2DstabC))\n",
+ "becomes $C\\leq C_{\\max} = 1/\\sqrt{2}$ in the present 2D case with the\n",
+ "$C$ defined above. Let us plot $1-\\tilde c/c$ in polar coordinates\n",
+ "for $C_{\\max}, 0.9C_{\\max}, 0.5C_{\\max}, 0.2C_{\\max}$.\n",
+ "The program below does the somewhat tricky\n",
+ "work in Matplotlib, and the result appears\n",
+ "in [Figure](#wave:pde1:fig:disprel2D). From the figure we clearly\n",
+ "see that the maximum $C$ value gives the best results, and that\n",
+ "waves whose propagation direction makes an angle of 45 degrees with\n",
+ "an axis are the most accurate."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def dispersion_relation_2D(p, theta, C):\n",
+ " arg = C*sqrt(sin(p*cos(theta))**2 +\n",
+ " sin(p*sin(theta))**2)\n",
+ " c_frac = 2./(C*p)*arcsin(arg)\n",
+ "\n",
+ " return c_frac\n",
+ "\n",
+ "import numpy as np\n",
+ "from numpy import \\\n",
+ " cos, sin, arcsin, sqrt, pi # for nicer math formulas\n",
+ "\n",
+ "r = p = np.linspace(0.001, pi/2, 101)\n",
+ "theta = np.linspace(0, 2*pi, 51)\n",
+ "r, theta = np.meshgrid(r, theta)\n",
+ "\n",
+ "# Make 2x2 filled contour plots for 4 values of C\n",
+ "import matplotlib.pyplot as plt\n",
+ "C_max = 1/sqrt(2)\n",
+ "C = [[C_max, 0.9*C_max], [0.5*C_max, 0.2*C_max]]\n",
+ "fix, axes = plt.subplots(2, 2, subplot_kw=dict(polar=True))\n",
+ "for row in range(2):\n",
+ " for column in range(2):\n",
+ " error = 1 - dispersion_relation_2D(\n",
+ " p, theta, C[row][column])\n",
+ " print error.min(), error.max()\n",
+ " # use vmin=error.min(), vmax=error.max()\n",
+ " cax = axes[row][column].contourf(\n",
+ " theta, r, error, 50, vmin=-1, vmax=-0.28)\n",
+ " axes[row][column].set_xticks([])\n",
+ " axes[row][column].set_yticks([])\n",
+ "\n",
+ "# Add colorbar to the last plot\n",
+ "cbar = plt.colorbar(cax)\n",
+ "cbar.ax.set_ylabel('error in wave velocity')\n",
+ "plt.savefig('disprel2D.png'); plt.savefig('disprel2D.pdf')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "Error in numerical dispersion in 2D.
\n", + "![](\"fig-wave/disprel2D.png\")
\n",
+ "\n",
+ ""
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.6"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}