{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "A very wide range of physical processes lead to wave motion, where\n", "signals are propagated through a medium in space and time, normally\n", "with little or no permanent movement of the medium itself.\n", "The shape of the signals may undergo changes as they travel through\n", "matter, but usually not so much that the signals cannot be recognized\n", "at some later point in space and time.\n", "Many types of wave motion can be described by the equation\n", "$u_{tt}=\\nabla\\cdot (c^2\\nabla u) + f$, which we will solve\n", "in the forthcoming text by finite difference methods.\n", "\n", "# Simulation of waves on a string\n", "
\n", "\n", "We begin our study of wave equations by simulating one-dimensional\n", "waves on a string, say on a guitar or violin.\n", "Let the string in the undeformed state\n", "coincide with the interval\n", "$[0,L]$ on the $x$ axis, and let $u(x,t)$ be the displacement at\n", "time $t$ in the $y$ direction of a point initially at $x$.\n", "The displacement function $u$ is governed by the mathematical model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial^2 u}{\\partial t^2} =\n", "c^2 \\frac{\\partial^2 u}{\\partial x^2}, \\quad x\\in (0,L),\\ t\\in (0,T]\n", "\\label{wave:pde1} \\tag{1}\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{wave:pde1:ic:u} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\frac{\\partial}{\\partial t}u(x,0) = 0, \\quad x\\in [0,L]\n", "\\label{wave:pde1:ic:ut} \\tag{3}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(0,t) = 0, \\quad t\\in (0,T]\n", "\\label{wave:pde1:bc:0} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(L,t) = 0, \\quad t\\in (0,T]\n", "\\label{wave:pde1:bc:L} \\tag{5}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The constant $c$ and the function $I(x)$ must be prescribed.\n", "\n", "Equation ([1](#wave:pde1)) is known as the one-dimensional\n", "*wave equation*. Since this PDE contains a second-order derivative\n", "in time, we need *two initial conditions*. The condition\n", "([2](#wave:pde1:ic:u)) specifies\n", "the initial shape of the string, $I(x)$, and\n", "([3](#wave:pde1:ic:ut)) expresses that the initial velocity of the\n", "string is zero. In addition, PDEs need *boundary conditions*, given here as\n", "([4](#wave:pde1:bc:0)) and ([5](#wave:pde1:bc:L)). These two\n", "conditions specify that\n", "the string is fixed at the ends, i.e., that the displacement $u$ is zero.\n", "\n", "The solution $u(x,t)$ varies in space and time and describes waves that\n", "move with velocity $c$ to the left and right.\n", "\n", "\n", "Sometimes we will use a more compact notation for the partial derivatives\n", "to save space:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u_t = \\frac{\\partial u}{\\partial t}, \\quad\n", "u_{tt} = \\frac{\\partial^2 u}{\\partial t^2},\n", "\\label{_auto1} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and similar expressions\n", "for derivatives with respect to other variables. Then the\n", "wave equation can be written compactly as $u_{tt} = c^2u_{xx}$.\n", "\n", "\n", "\n", "The PDE problem ([1](#wave:pde1))-([5](#wave:pde1:bc:L)) will now be\n", "discretized in space and time by a finite difference method.\n", "\n", "\n", "## Discretizing the domain\n", "
\n", "\n", "\n", "The temporal domain $[0,T]$ is represented by a finite number of mesh points" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "0 = t_0 < t_1 < t_2 < \\cdots < t_{N_t-1} < t_{N_t} = T \\label{_auto2} \\tag{7}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly, the spatial domain $[0,L]$ is replaced by a set of mesh points" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "0 = x_0 < x_1 < x_2 < \\cdots < x_{N_x-1} < x_{N_x} = L \\label{_auto3} \\tag{8}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One may view the mesh as two-dimensional in the $x,t$ plane, consisting\n", "of points $(x_i, t_n)$, with $i=0,\\ldots,N_x$ and $n=0,\\ldots,N_t$.\n", "\n", "\n", "\n", "### Uniform meshes\n", "\n", "For uniformly distributed mesh points we can introduce the constant\n", "mesh spacings $\\Delta t$ and $\\Delta x$. We have that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "x_i = i\\Delta x,\\ i=0,\\ldots,N_x,\\quad\n", "t_n = n\\Delta t,\\ n=0,\\ldots,N_t\n", "\\label{_auto4} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We also have that $\\Delta x = x_i-x_{i-1}$, $i=1,\\ldots,N_x$, and\n", "$\\Delta t = t_n - t_{n-1}$, $n=1,\\ldots,N_t$. [Figure](#wave:pde1:fig:mesh)\n", "displays a mesh in the $x,t$ plane with $N_t=5$, $N_x=5$, and constant\n", "mesh spacings.\n", "\n", "## The discrete solution\n", "
\n", "\n", "\n", "The solution $u(x,t)$ is sought at the mesh points. We introduce\n", "the mesh function $u_i^n$, which approximates the exact\n", "solution at the\n", "mesh point $(x_i,t_n)$ for $i=0,\\ldots,N_x$ and $n=0,\\ldots,N_t$.\n", "Using the finite difference method, we shall\n", "develop algebraic equations for computing the mesh function.\n", "\n", "## Fulfilling the equation at the mesh points\n", "
\n", "\n", "\n", "In the finite difference method, we relax\n", "the condition that ([1](#wave:pde1)) holds at all points in\n", "the space-time domain $(0,L)\\times (0,T]$ to the requirement that the PDE is\n", "fulfilled at the *interior* mesh points only:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial^2}{\\partial t^2} u(x_i, t_n) =\n", "c^2\\frac{\\partial^2}{\\partial x^2} u(x_i, t_n),\n", "\\label{wave:pde1:step2} \\tag{10}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "for $i=1,\\ldots,N_x-1$ and $n=1,\\ldots,N_t-1$. For $n=0$ we have\n", "the initial conditions $u=I(x)$ and $u_t=0$,\n", "and at the boundaries $i=0,N_x$ we\n", "have the boundary condition $u=0$.\n", "\n", "## Replacing derivatives by finite differences\n", "
\n", "\n", "The second-order derivatives can be replaced by central\n", "differences. The most widely used difference approximation of\n", "the second-order derivative is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial^2}{\\partial t^2}u(x_i,t_n)\\approx\n", "\\frac{u_i^{n+1} - 2u_i^n + u^{n-1}_i}{\\Delta t^2}\\\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "mathcal{I}_t is convenient to introduce the finite difference operator notation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "[D_tD_t u]^n_i = \\frac{u_i^{n+1} - 2u_i^n + u^{n-1}_i}{\\Delta t^2}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A similar approximation of the second-order derivative in the $x$\n", "direction reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial^2}{\\partial x^2}u(x_i,t_n)\\approx\n", "\\frac{u_{i+1}^{n} - 2u_i^n + u^{n}_{i-1}}{\\Delta x^2} = [D_xD_x u]^n_i\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Algebraic version of the PDE\n", "\n", "We can now replace the derivatives in ([10](#wave:pde1:step2))\n", "and get" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{u_i^{n+1} - 2u_i^n + u^{n-1}_i}{\\Delta t^2} =\n", "c^2\\frac{u_{i+1}^{n} - 2u_i^n + u^{n}_{i-1}}{\\Delta x^2},\n", "\\label{wave:pde1:step3b} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "or written more compactly using the operator notation:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "[D_tD_t u = c^2 D_xD_x]^{n}_i\n", "\\label{wave:pde1:step3a} \\tag{12}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interpretation of the equation as a stencil\n", "\n", "A characteristic feature of ([11](#wave:pde1:step3b)) is that it\n", "involves $u$ values from neighboring points only: $u_i^{n+1}$,\n", "$u^n_{i\\pm 1}$, $u^n_i$, and $u^{n-1}_i$. The circles in [Figure](#wave:pde1:fig:mesh) illustrate such neighboring mesh points that\n", "contribute to an algebraic equation. In this particular case, we have\n", "sampled the PDE at the point $(2,2)$ and constructed\n", "([11](#wave:pde1:step3b)), which then involves a coupling of $u_1^2$,\n", "$u_2^3$, $u_2^2$, $u_2^1$, and $u_3^2$. The term *stencil* is often\n", "used about the algebraic equation at a mesh point, and the geometry of\n", "a typical stencil is illustrated in [Figure](#wave:pde1:fig:mesh). One also often refers to the algebraic\n", "equations as *discrete equations*, *(finite) difference equations* or\n", "a *finite difference scheme*.\n", "\n", "\n", "\n", "
\n", "\n", "

Mesh in space and time. The circles show points connected in a finite difference equation.

\n", "\n", "\n", "\n", "\n", "\n", "\n", "### Algebraic version of the initial conditions\n", "\n", "We also need to replace the derivative in the initial condition\n", "([3](#wave:pde1:ic:ut)) by a finite difference approximation.\n", "A centered difference of the type" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial}{\\partial t} u(x_i,t_0)\\approx\n", "\\frac{u^1_i - u^{-1}_i}{2\\Delta t} = [D_{2t} u]^0_i,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "seems appropriate. Writing out this equation and ordering the terms give" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u^{-1}_i=u^{1}_i,\\quad i=0,\\ldots,N_x\n", "\\label{wave:pde1:step3c} \\tag{13}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The other initial condition can be computed by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u_i^0 = I(x_i),\\quad i=0,\\ldots,N_x\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Formulating a recursive algorithm\n", "
\n", "\n", "\n", "We assume that $u^n_i$ and $u^{n-1}_i$ are available for\n", "$i=0,\\ldots,N_x$. The only unknown quantity in\n", "([11](#wave:pde1:step3b)) is therefore $u^{n+1}_i$, which we now can\n", "solve for:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u^{n+1}_i = -u^{n-1}_i + 2u^n_i + C^2\n", "\\left(u^{n}_{i+1}-2u^{n}_{i} + u^{n}_{i-1}\\right)\n", "\\label{wave:pde1:step4} \\tag{14}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have here introduced the parameter" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "C = c\\frac{\\Delta t}{\\Delta x},\n", "\\label{_auto5} \\tag{15}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "known as the *Courant number*.\n", "\n", "\n", "**$C$ is the key parameter in the discrete wave equation.**\n", "\n", "\n", "We see that the discrete version of the PDE features only one\n", "parameter, $C$, which is therefore the key parameter, together with\n", "$N_x$, that governs the quality of the numerical solution (see the section [Analysis of the difference equations](wave_analysis.ipynb) for details). Both the primary physical\n", "parameter $c$ and the numerical parameters $\\Delta x$ and $\\Delta t$\n", "are lumped together in $C$. Note that $C$ is a dimensionless\n", "parameter.\n", "\n", "\n", "\n", "Given that $u^{n-1}_i$ and $u^n_i$ are known for $i=0,\\ldots,N_x$,\n", "we find new values at the next time level by applying the formula\n", "([14](#wave:pde1:step4)) for $i=1,\\ldots,N_x-1$. [Figure](#wave:pde1:fig:mesh) illustrates the points that are used to\n", "compute $u^3_2$. For the boundary points, $i=0$ and $i=N_x$, we apply\n", "the boundary conditions $u_i^{n+1}=0$.\n", "\n", "\n", "Even though sound reasoning leads up to\n", "([14](#wave:pde1:step4)), there is still a minor challenge with it that needs\n", "to be resolved. Think of the very first computational step to be made.\n", "The scheme ([14](#wave:pde1:step4)) is supposed to start at $n=1$, which means\n", "that we compute $u^2$ from $u^1$ and $u^0$. Unfortunately, we do not know the\n", "value of $u^1$, so how to proceed? A standard procedure in such cases is to\n", "apply ([14](#wave:pde1:step4)) also for $n=0$. This immediately seems strange,\n", "since it involves $u^{-1}_i$, which is an undefined\n", "quantity outside the time mesh (and the time domain). However, we can\n", "use the initial condition ([13](#wave:pde1:step3c)) in combination with\n", "([14](#wave:pde1:step4)) when $n=0$ to eliminate $u^{-1}_i$ and\n", "arrive at a special formula for $u_i^1$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u_i^1 = u^0_i - \\frac{1}{2}\n", "C^2\\left(u^{0}_{i+1}-2u^{0}_{i} + u^{0}_{i-1}\\right)\n", "\\label{wave:pde1:step4:1} \\tag{16}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Figure](#wave:pde1:fig:stencil:u1) illustrates how ([16](#wave:pde1:step4:1))\n", "connects four instead of five points: $u^1_2$, $u_1^0$, $u_2^0$, and $u_3^0$.\n", "\n", "\n", "\n", "
\n", "\n", "

Modified stencil for the first time step.

\n", "\n", "\n", "\n", "\n", "\n", "We can now summarize the computational algorithm:\n", "\n", "1. Compute $u^0_i=I(x_i)$ for $i=0,\\ldots,N_x$\n", "\n", "2. Compute $u^1_i$ by ([16](#wave:pde1:step4:1)) for $i=1,2,\\ldots,N_x-1$ and set $u_i^1=0$\n", " for the boundary points given by $i=0$ and $i=N_x$,\n", "\n", "3. For each time level $n=1,2,\\ldots,N_t-1$\n", "\n", " a. apply ([14](#wave:pde1:step4)) to find $u^{n+1}_i$ for $i=1,\\ldots,N_x-1$\n", "\n", " b. set $u^{n+1}_i=0$ for the boundary points having $i=0$, $i=N_x$.\n", "\n", "\n", "The algorithm essentially consists of moving a finite difference\n", "stencil through all the mesh points, which can be seen as an animation\n", "in a [web page](mov-wave/D_stencil_gpl/index.html)\n", "or a [movie file](mov-wave/D_stencil_gpl/movie.ogg)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Sketch of an implementation\n", "
\n", "\n", "We start by defining some constants that will be used throughout our Devito code." ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "# Given mesh points as arrays x and t (x[i], t[n]),\n", "# constant c and function I for initial condition\n", "x = np.linspace(0, 2, 101)\n", "t = np.linspace(0, 2, 101)\n", "c = 1\n", "I = lambda x: np.sin(x)\n", "\n", "dx = x[1] - x[0]\n", "dt = t[1] - t[0]\n", "C = c*dt/dx # Courant number\n", "Nx = len(x)-1\n", "Nt = len(t)-1\n", "C2 = C**2 # Help variable in the scheme\n", "L = 2." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we define our 1D computational grid and create a function `u` as a symbolic `devito.TimeFunction`. We need to specify the `space_order` as 2 since our wave equation involves second-order derivatives with respect to $x$. Similarly, we specify the `time_order` as 2, as our equation involves second-order derivatives with respect to $t$. Setting these parameters allows us to use `u.dx2` and `u.dt2`." ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [], "source": [ "from devito import Grid, TimeFunction\n", "\n", "# Initialise `u` for space and time order 2, using initialisation function I\n", "grid = Grid(shape=(Nx+1), extent=(L))\n", "u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)\n", "u.data[:,:] = I(x[:])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have initialised `u`, we can solve our wave equation for the unknown quantity $u^{n+1}_i$ using forward and backward differences in space and time." ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "LHS: u(t + dt, x)\n", "RHS: 1.0*dt**2*(-2.0*u(t, x)/h_x**2 + u(t, x - h_x)/h_x**2 + u(t, x + h_x)/h_x**2 + 2.0*u(t, x)/dt**2 - 1.0*u(t - dt, x)/dt**2)\n" ] } ], "source": [ "from devito import Constant, Eq, solve\n", "\n", "# Set up wave equation and solve for forward stencil point in time\n", "pde = (1/c**2)*u.dt2-u.dx2\n", "stencil = Eq(u.forward, solve(pde, u.forward))\n", "\n", "print(\"LHS: %s\" % stencil.lhs)\n", "print(\"RHS: %s\" % stencil.rhs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Great! From these print statements, we can see that Devito has taken the wave equation in ([1](#wave:pde1)) and solved it for $u^{n+1}_i$, giving us equation ([14](#wave:pde1:step4)). Note that `dx` is denoted as `h_x`, while `u(t, x)`, `u(t, x - h_x)` and `u(t, x + h_x)` denote the equivalent of $u^{n}_{i}$, $u^{n}_{i-1}$ and $u^{n}_{i+1}$ respectively. \n", "\n", "We also need to create a separate stencil for the first timestep, where we substitute $u^{1}_i$ for $u^{-1}_i$, as given in ([13](#wave:pde1:step3c))." ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "stencil_init = stencil.subs(u.backward, u.forward)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can create expressions for our boundary conditions and build the operator. The results are plotted below." ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Data type float64 of runtime value `dt` does not match the Constant data type \n", "Operator `Kernel` run in 0.01 s\n", "Data type float64 of runtime value `dt` does not match the Constant data type \n", "Operator `Kernel` run in 0.01 s\n" ] }, { "data": { "text/plain": [ "PerformanceSummary([(PerfKey(name='section0', rank=None),\n", " PerfEntry(time=5.899999999999998e-05, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#NBVAL_IGNORE_OUTPUT\n", "from devito import Operator\n", "\n", "t_s = grid.stepping_dim \n", "\n", "# Boundary conditions\n", "bc = [Eq(u[t_s+1, 0], 0)]\n", "bc += [Eq(u[t_s+1, Nx], 0)]\n", "\n", "# Defining one Operator for initial timestep and one for the rest\n", "op_init = Operator([stencil_init]+bc)\n", "op = Operator([stencil]+bc)\n", " \n", "op_init.apply(time_M=1, dt=dt)\n", "op.apply(time_m=1,time_M=Nt, dt=dt)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can plot our results using `matplotlib`:" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU9b3/8ddnZrKQEMKSsO+7IIgQENS6YuuOrV1EFG1RlFavrb3eX3ttba+9t612u7fVVnGpG2u1VaooVsW6IJCw7xL2AELYwhKyTb6/P2aC0wgSNGfOZOb9fDzycObMmczbk2E+c87nnO/XnHOIiEjqCvgdQERE/KVCICKS4lQIRERSnAqBiEiKUyEQEUlxIb8DnKq8vDzXvXt3v2OIiDQpixYt2uOcyz/eY02uEHTv3p2ioiK/Y4iINClmtuVEj+nQkIhIilMhEBFJcSoEIiIpToVARCTFqRCIiKQ4zwqBmT1pZrvNbOUJHjcz+72ZFZvZcjMb6lUWERE5MS/3CJ4CLv2Uxy8D+kR/JgJ/8jCLiIicgGfXETjn3jGz7p+yyhjgGRcZB3u+mbU0sw7OuZ1eZRKRxFFVU8uB8ir2lVex70gVhypqOFJZw+HKGiqra6kK11JVU0vsUPnBQID0UIC0oJGVHiI7I0jzjBAts9JomZVO66x0cpulEQiYj/9nTY+fF5R1ArbF3C+JLvtEITCziUT2GujatWtcwonI51NVU8vWfUfYUHqEbfvKIz/7j/JRWQW7Dlaw90hVg36PRT/TGzp1SlrQaJuTSdsWGXRs2YzOrZrRpVUWPfKy6ZmfTfsWmZipUMRqElcWO+cmA5MBCgoKNJOOSIL5qKyCFdvLWLPzIGt2HmTdR4fYsq+ccO3H/1xzMkJ0aZ1Fh9xMhnRtSbucTNo0T6d1djots9JokZlG84wQ2RkhMtMi3/zTg4F/+dAO1zqqaiJ7CuXVNRypDHO4soayo9XsP1LF3iNV7Dlcya6DkWKzesdB/rFqF1Xh2mO/Izs9SO+2zenfvgX9O+QwoEMLBnbKpXlGk/g49ISf/+fbgS4x9ztHl4lIAqusCbOipIzCzftZvHU/y7YdYPehSiDy7b1b6yz6tc/h8kEd6NU2m555zenWJovcZmmf+5t4MGA0Sw/SLD1ILmkNek5trWPXoQo27YnsnWwsPcz6XYd5Y80uZhRtO5a7d35zzujSkuHdWzGsW2t65WenzJ6Dn4VgFnCHmU0HzgLK1B8QSTzV4VqWbTvAvA17mbdhD4u3HqCqJvINu2deNuf0zmNw51wGd86lf/sWZCfYN+tAwOiQ24wOuc04u1feseXOOUoPVbJyRxnLSyI/b67ZxfOLSgBok53OyJ5tGNWrDef0zqN7m6ykLQye/cXMbBpwAZBnZiXATyBSwp1zjwCzgcuBYqAc+KZXWUTk1Ow+WMHcdbt5e10p763fw6HKGsxgQIcWjB/ZjeE9WlPQrRVtmmf4HfUzMzPatsjkohaZXNS/HRApDhtKj7Boyz4WbNrHvOK9vLIi8v20a+ssLuiXzwX98jm7Vx6ZaUE/4zcqa2qT1xcUFDiNPirS+DaWHua1VR/x+qpdLN12AID2LTK5oF8+5/XNZ1TPNrTKTvc5ZXw559i05wjvF+/h7XWlzNuwl6PVYZqlBTmvbx5fHNCe0QPakdusYYep/GRmi5xzBcd9TIVAJHVt2XuEvy/bwSsrPmLNzoMADO6cyyWntWP0gHb0b5+TtIdDPovKmjDzN+7jjdW7eGPNLnaWVZAWNM7tnccVgzvypYHtyMlMzKKgQiAix+w/UsXLy3fwtyXbWbw18s1/WLdWXDGoA5cNak+H3GY+J2wanHMsKylj9oqdvLJ8J9sPHCUjFGD0gHZ8eUgnzu+XT1owcUbxUSEQSXHhWsd7xXuYWbTt2OmU/drl8OWhnbjqjI50aqkP/8/DOcfirQd4ael2Xl6+k31HqsjPyeArQzvxtWFd6N22ud8RVQhEUlXpoUpmFm1j6oKtbD9wlFZZaVxzZie+OqwzAzvm+h0vKVWHa5m7djd/WVTCW2t3E651nNWjNeNGduNLA9uREfKnyaxCIJJiFm/dz1Pvb+bVlTupDjvO7tWG68/qyiUD/PsgSkW7D1XwwqLtTFu4la37yslrns7YEV25YWQ32rXIjGsWFQKRFFATrmX2yo944r1NLNt2gJyMEF8r6MK4kV3ple//oYlUVlvreLd4D89+sJk31+4maMYVgztw6xd6cnqn+OyZqRCIJLHyqhpmFm7j8fc2UbL/KD3ysrn57O5cO6xzSg+bkKg27znCMx9sYWbRNg5X1nB2rzbcdn4vzu+b7+nrqhCIJKGyo9U8M28zT76/if3l1RR0a8Vt5/fi4v5tNfpmE3CwopppC7by5Pub2HWwkhcmjWJYt9aevd6nFQJ9XRBpYg6UV/HYuxt5Zt4WDlXWcHH/tnz7wl6efohI42uRmcZt5/firJ5tuObh9yk7Wu1bFhUCkSai7Gg1T7y3iSff28ThyhouH9Se71zYW2f/NHGh6N5bTdi/ozMqBCIJ7mhVmD/P28Qjb2/gYEUNl53enrtG96F/+xZ+R5NGEApGCkHskN1xz+DbK4vIp6oJ1zK9cBu/f3M9uw9VcnH/ttz9xb7aA0gydXsE1SoEIlLHOceba3bzi1fXsKH0CAXdWvHwuKEM764eQDIKBSLDUIRra0+ypocZfHtlEfmE1TsO8rOXV/PBxr30zMtm8o3DuGRAOw38lsSC6hGICMDew5X85h8fMn3hVnKbpXH/mIGMHdE1oQYtE2/U9QhqdGhIJDXVhGuZsmArv359HUerwtx8dg/uurgPuVmJOZSxNL5jewQqBCKpZ9GWffz4xVWs3nmQc3vn8dOrB9C7bY7fsSTO0up6BGH1CERSxoHyKh54bS3TFm6jfYtMHr5+KJcPaq8+QIoK6tCQSOpwzjFr2Q5+9vJq9pdXc+sXevDd0X0TbrJ3ia+QDg2JpIaS/eXc+7eV/PPDUs7o0pKnv3W6rgcQIPb0URUCkaRUW+t4dv4WHnhtLQA/uWoA40d1P9YgFNEQEyJJbMveI9zz/HIWbtrHeX3z+Z9rTqdL6yy/Y0mCCQQMM6jRBWUiyaNuL+CXr64lFDAe/Opgvjass5rBckKhgKlHIJIsdhw4yj3PL+P94r2c3zefX147iA65mhhePl0oEFCPQKSpc87x4tLt3PfSKsK1jl98ZRDXDe+ivQBpkFDAqNZ1BCJNV9nRau792wpeXr6Tgm6t+M3Xz6Bbm2y/Y0kTEgya9ghEmqoFG/dy98xl7DpYwT1f6sft5/fSGUFyykKBgHoEIk1NTbiW37+5nofmFtO1dRbPTzqbIV1a+h1LmqhQwAjr9FGRpmPHgaPcNX0JhZv3c+3Qztw/ZqCuDpbPJRgwqnX6qEjT8OaaXdw9cxk14Vr+9xtDuObMTn5HkiSQph6BSOKrDtfy6znrePSdjQzs2IKHrh9Kjzw1hKVxBH2+jsDTWS/M7FIzW2dmxWb2g+M83tXM5prZEjNbbmaXe5lH5LPYWXaUsZPn8+g7G7lhZFdemHS2ioA0qlAgQE0ynj5qZkHgYeASoAQoNLNZzrnVMav9CJjpnPuTmQ0AZgPdvcokcqrmFe/hzmlLqKgO8/uxZ3L1GR39jiRJKBhI3kNDI4Bi59xGADObDowBYguBA1pEb+cCOzzMI9JgtbWOR97ZwK/nrKNnfnMeuWEYvds29zuWJKm0YPIOMdEJ2BZzvwQ4q946PwVeN7M7gWxg9PF+kZlNBCYCdO3atdGDisQ6XFnD92cuZc6qXVw5uAMPXDtYZwWJp/zeI/B7ZuyxwFPOuc7A5cCzZvaJTM65yc65AudcQX5+ftxDSurYUHqYax5+nzfW7OZHV5zGH8aeqSIgngsFAkk7xMR2oEvM/c7RZbEmAJcCOOc+MLNMIA/Y7WEukeN6a+0u7pq2lLRQgGcnjODsXnl+R5IUEQr6O9aQl3sEhUAfM+thZunAdcCseutsBS4GMLPTgEyg1MNMIp/gnOOPbxcz4ekiurbJYtYd56gISFz5ffqoZ3sEzrkaM7sDmAMEgSedc6vM7H6gyDk3C/g+8JiZfY9I4/hm55x/W0NSztGqMP/vheXMWraDq87oyIPXDqZZetDvWJJiQgFL3hnKnHOziZwSGrvsvpjbq4FzvMwgciK7DlZw6zNFrNhexn9c2o9J5/fSsNHii1BQg86JxN2KkjJueaaQQxU1TL6xgEsGtPM7kqSwUMAIa6whkfh5beVOvjtjKW2yM3hh0tmc1qHFyZ8k4qFgMh8aEkkkzjkmv7ORX7y6liFdWvLY+ALyczL8jiWiOYtF4qE6XMt9L61i2sKtXDG4A7/52hlkpqkpLIkhFNScxSKeOlRRzbenLObd9Xv4zoW9+P4l/QhoFjFJIJE9AvUIRDyx62AFN/+5kA93HeKBawfxjeEaokQSj3oEIh5Z99EhvvnnhZQdrebJm4dzfl8NTyKJKU2nj4o0vgUb93LLM0U0Swsy47ZRnN4p1+9IIifk96BzKgSSdF5dsZO7ZiylS6tmPP2tEXRuleV3JJFPFQr4O9aQCoEklWc/2Mx9s1YxtGsrHh9fQKvsdL8jiZxUSHMWi3x+zjl+98Z6fv/mekaf1o6Hrj9Tp4dKkxEMRHoEzjlfhjlRIZAmL1zruO+llUxZsJWvF3Tm518eRCjo91QbIg0Xip7OHK51hIIqBCKnpLImzN0zlvHKip1MuqAX//Glfho4TpqcYLQQ1NQ6Qj7syKoQSJNVXlXDbc8u4t31e7j38tO49byefkcS+UzSgh/vEfhBhUCapLLyar71dCFLtu7nwWsH8/XhXU7+JJEEFQxEDmX6dS2BCoE0OXsOV3LjEwsp3n2Ih68fymWDOvgdSeRzqesR1Ph0CqkKgTQpO8uOMu7xBew4cJQnbhrOebpaWJJASIeGRBpmy94jjHt8AWXl1Tw74SyGd2/tdySRRhGKaRb78vq+vKrIKSrefZhxj8+nsqaWqbeOZFBnDRkhyeNYj8CngedUCCThrf3oIDc8vgAwZkwcRb/2OX5HEmlUH+8R+NMj0FU3ktBWbi/jusnzCQaMGbeNVBGQpKQegcgJLNt2gBufWEBOZhpTbz2Lbm2y/Y4k4gn1CESOY/HW/dz0xEJaZqcx7daRGkFUkpp6BCL1LNqyj5ueLKRN83Sm3TqSji2b+R1JxFN1h4b86hGoEEhCWbRlH+OfWEjbFplMu3Uk7XMz/Y4k4rnYQed8eX1fXlXkOOqKQLsWmUybOJJ2LVQEJDXUDTpX7dOhIZ01JAlBRUBSWVp02HS/9ghUCMR3i7fu56YnCyOHg1QEJAUFdR2BpLJl2w5w0xMLyYs2hlUEJBV9POic9ggkxazcXsaNTyygZXYaU9UYlhQW9Pk6Ak8LgZldambrzKzYzH5wgnW+bmarzWyVmU31Mo8kjtU7DnJD9GIxnSIqqc7vHoFnZw2ZWRB4GLgEKAEKzWyWc251zDp9gB8C5zjn9ptZW6/ySOJYv+sQNz6xgGZpQV0sJkJy9whGAMXOuY3OuSpgOjCm3jq3Ag875/YDOOd2e5hHEsDmPZGhpAMBY8otZ9G1jYqASDL3CDoB22Lul0SXxeoL9DWz981svplderxfZGYTzazIzIpKS0s9iiteK9lfzvWPzaem1jHllrPomd/c70giCSGU4qePhoA+wAXAWOAxM2tZfyXn3GTnXIFzriA/XzNSNUW7D1Yw7vEFHK6s4dkJI+jbTqOIitTxe9A5LwvBdiB2RvHO0WWxSoBZzrlq59wm4EMihUGSyL4jVYx7fAF7DlXy1LdGMLCjJpURiZXMPYJCoI+Z9TCzdOA6YFa9dV4ksjeAmeUROVS00cNMEmcHK6oZ/+QCtu4r5/GbhjO0ayu/I4kknDSfRx/1rBA452qAO4A5wBpgpnNulZndb2ZXR1ebA+w1s9XAXOAe59xerzJJfB2tCjPhqULWfXSIR24cxqhebfyOJJKQgsk8MY1zbjYwu96y+2JuO+Du6I8kkcqaMLc9t4hFW/bzh7FDubCfzgwWOZG6HkG1hqGWZBGudXxvxlLe+bCUB68dzBWDO/gdSSSh1fUIwsl2aEhSk3OO//zrCmav+IgfXXEaXx/e5eRPEklxyXzWkKSgX766lhlF27jzot7c8oWefscRaRLMjGDAUvY6Akkif3p7A4++s5EbR3bj7kv6+h1HpEkJBsy3HoEKgTSKGYVbeeC1tVx9Rkf+6+qBmJnfkUSalLSAqUcgTdecVR/xw7+u4Py++fz6a2cQCKgIiJyqYMDUI5Cmaf7Gvdw5bQlndGnJn24YSnpIbymRzyIUDCTllcWS5FbtKOPWp4vo1jqLP988nKx0nY0s8lmF1CyWpmbbvnJu/nMhOZkhnpkwgpZZ6X5HEmnSQgFLviEmJHntPVzJ+CcXUlVTyzMTRtAhV7OLiXxewaD2CKSJOFJZw7eeKmRn2VGevLmA3m01nLRIYwgFAlSrEEiiqw7X8u0pi1mxvYyHxg5lWLfWfkcSSRqRHoHGGpIE5pzjh39dwT8/LOWXXxnE6AHt/I4kklSCPvYIGlQIzOy+4y13zt3fuHEkUf3m9Q95flEJ3x3dh+tGdPU7jkjSCQX9u46goXsER2JuZwJXEpljQFLAc/O38NDcYsaO6MJdF2sCOREvhAKBxC4EzrnfxN43s18TmVRGktwbq3dx30sruah/W3425nQNHSHiET97BJ+1WZxFZA5iSWJLtx3gjmmLOb1TLg9dfyahoM4tEPFKMGBUJ3iPYAVQlzAI5APqDySxLXuPMOGpQtrmZPLETbpqWMRroaBRUZ3YZw1dGXO7BtgVnZNYktC+I1Xc/OdCws7x1DeHk5+T4XckkaQX6RGE/XnthqzknNvidRBJDBXVYW59pojtB44y7daz6Jnf3O9IIimhKfYIJAnV1jq+P3MZi7fu53+/MUQXjInEkZ/XEagQyDEPvLaWV1bs5D8vO43LB2nCeZF4Sgv6d/qoCoEAkWsFHn1nI+NHdeOWL/TwO45IytGcxeKruet2H7tW4L4rB+haAREfhAJGdVg9AvHB6h0HuWPKYvq3b8EfxupaARG/hDQMtfhh18EKJjxdSE5mGk/ePJzsDF0rIOKXYKIPMSHJp7yqhglPF1J2tJq/3D6K9rmZfkcSSWmRGcp0aEjiJFzruGv6UlbvOMhD15/JwI65fkcSSXnBgH+jj6oQpKBfvrqGf6zexX1XDuCi/ppXQCQRpKlHIPEydcFWHnt3EzeN6sbN5+g0UZFE4WePQIUghby3fg8/fmklF/TL58dXDvA7jojESNoegZldambrzKzYzH7wKetda2bOzAq8zJPKincfZtKURfTOb67TREUSUCho1LrIUC/x5tmngZkFgYeBy4ABwFgz+8TXUDPLAe4CFniVJdXtP1LFhKcLyQgFeOLmAnIy0/yOJCL1hAKRCznDLokKATACKHbObXTOVQHTgTHHWe9nwANAhYdZUlZVTS23PbeInWUVTB5fQOdWWX5HEpHjCAYiH8d+DDznZSHoBGyLuV8SXXaMmQ0FujjnXvm0X2RmE82syMyKSktLGz9pknLO8aMXV7Bw0z5+9dXBDO3ayu9IInICacHIHkGND0NR+3ag2MwCwG+B759sXefcZOdcgXOuID8/3/twSeKxdzcys6iEOy/qzZghnU7+BBHxTbDu0FAy9QiA7UCXmPudo8vq5ACnA2+b2WZgJDBLDePG8eaaXfzi1bVcPqg93xvd1+84InISdT0CP+Yt9rIQFAJ9zKyHmaUD1wGz6h50zpU55/Kcc92dc92B+cDVzrkiDzOlhHUfHeLfpi1hYMcW/OZrQwgENJqoSKKr6xEk1R5BdE7jO4A5wBpgpnNulZndb2ZXe/W6qW7v4UomPF1IdkaIx8YX0Cw96HckEWmAkI89Ak8HnXPOzQZm11t23wnWvcDLLKmgqqaWSc8tpvRQJTNuG0WH3GZ+RxKRBgr52CPQ6KNJwjnHfS+tZOHmffzfdUMY0qWl35FE5BQEk7RHIHH01LzNTC/cxh0X6gwhkaYoLZiEPQKJn3fXl/Kzl1fzxQHtuPsSnSEk0hTV7RGk1HUE0jg27TnCd6Yspm+7HH73DZ0hJNJU1fUIku3KYvHYwYpqbnm6kFAwwGPjCzTVpEgTVjcQpB9DUasQNFHhWsdd05awZW85fxw3lC6tNYaQSFOms4bklD04Zy1z15Xy39eczsiebfyOIyKfk3oEckpeWrqdR/+5kXFndeWGkd38jiMijUA9AmmwFSVl/MfzyxnRozU/uWqg33FEpJGEdPqoNETpoUomPltEXvMM/jhuKOkh/flEksWxPQL1COREqmpq+faURewvr+KFSWeT1zzD70gi0oiO9Qh8mLdYhaCJ+OnfV1G4eT9/GHsmAzvm+h1HRBrZxxPT6NCQHMdz87cwdcFWJl3Qi6vO6Oh3HBHxQFIOQy2NY+Gmffx01iou6JfPv3+xn99xRMQjH09Mo9NHJcaOA0f59pRFdGmdxf9dd+axY4giknz8nKpSPYIEVVEd5vbnFlFRXcv0icPIbZbmdyQR8VDIxx6BCkECcs5x799WsrykjMfGF9C7bY7fkUTEYyH1CCTWU/M288LiEr47ug+XDGjndxwRiYOgegRS54MNe/nvV9ZwyYB2/NtFffyOIyJxUnf6qPYIUtz2A0e5Y+piurfJ4rdfP0NzC4ikkKCPVxarECSIiuowtz+7iKqaWiaPLyAnU81hkVRS1yPwY9A5NYsTQF1zeMX2Mh4fX0Cv/OZ+RxKROAsGDDMIaxjq1PTMB1uONYdHqzkskrJCAdOhoVS0cNM+fvbyakaf1lbNYZEUF1QhSD0flVUcu3L4t5p4XiTlhQIBTUyTSiprwkyasojyqjCP3jiMFmoOi6S8UNB86RGoWeyT//r7apZsPcAfxw2lbztdOSwi6hGklBmFW5m6YCu3n9+Lywd18DuOiCSIYMB0aCgVLN12gB+/uIpze+dxz5c0rLSIfCwUCGiPINntOVzJpOcWkZ+Twe/HalhpEflXfvUIPC0EZnapma0zs2Iz+8FxHr/bzFab2XIze9PMunmZx0814VrunLqEfUeqePTGYbTOTvc7kogkmGDAqE6mPQIzCwIPA5cBA4CxZjag3mpLgALn3GDgeeBBr/L47Vdz1vHBxr38z5cHcXonzTksIp+UFggQTrIewQig2Dm30TlXBUwHxsSu4Jyb65wrj96dD3T2MI9vXlm+k0ff2ciNI7vx1WFJ+b8oIo0gGS8o6wRsi7lfEl12IhOAV4/3gJlNNLMiMysqLS1txIjeW7/rEPc8v4yhXVvy4yvr7xCJiHwsKXsEDWVmNwAFwK+O97hzbrJzrsA5V5Cfnx/fcJ/DoYpqbnt2EVnpQf44bhjpoYTY3CKSoPzaI/DygrLtQJeY+52jy/6FmY0G7gXOd85Vepgnrpxz/PtflrFlXzlTbjmL9rmZfkcSkQSXloRDTBQCfcysh5mlA9cBs2JXMLMzgUeBq51zuz3MEnd/+ucG5qzaxQ8v68/Inm38jiMiTUAwYMk1Q5lzrga4A5gDrAFmOudWmdn9ZnZ1dLVfAc2Bv5jZUjObdYJf16S8t34Pv56zjisGd2DCuT38jiMiTUQoaFQn21hDzrnZwOx6y+6LuT3ay9f3w/YDR/m36Uvold+cB68djJkuGhORhgkl2x5BKqqsCfPt5yLTTT5y4zCyMzSmn4g0XNCnHoE+qRrR/X9fzbKSMh65YZimmxSRUxYZfTRFTx9NBn8p2saUBVu57fyeXHp6e7/jiEgTFAom3wVlKWPVjjJ+9OJKRvVswz1f1IiiIvLZqEfQRJWVV3P7c4tolZXOH64/k1BQm1REPhv1CJqg2lrH3TOX8lFZBdMnjiKveYbfkUSkCVOPoAl6eG4xb67dzY+vHMCwbq38jiMiTVxkrCEdGmoy3vmwlN++8SHXDOnIjSOTdhoFEYkjzVnchJTsL+eu6Uvo2zaHn39lkC4aE5FG4VePQIXgFFXWhPnOlMXUhB1/umEoWelqs4hI40gL+tMj0KfYKYq9aKynLhoTkUaUdIPOJaMXFpXoojER8UwoYFSHHc7FtxioEDTQmp0HuffFFZzVo7UuGhMRTwQDkY/keO8UqBA0wMGKaiY9t4gWmWm6aExEPBMKRk48iXefQD2Ck3DO8e8zl1Gy/yjTJo6kbY5mGhMRb4QCkUIQ7z6BvtqexOR3NvL66l384LL+DO/e2u84IpLEgtFCUB3nU0hVCD7FBxv28sBra7l8UHvNNCYinkuLHnbWHkGC2HWwgjunLaF7XjYPaKYxEYmDuj0C9QgSQHW4ljumLuZIZQ1Tbz2LnMw0vyOJSAqo6xHE++piFYLjePC1tRRu3s//XTeEvu1y/I4jIikipENDieHVFTt57N1NjB/VjTFDOvkdR0RSyLE9AhUC/2wsPcw9zy9nSJeW3HvFaX7HEZEUEzx2+mh8ewQqBFHlVTVMem4xaUHj4XFDyQgF/Y4kIikm5NPpo+oRELlo7Ed/W8mHuw/x9DdH0KllM78jiUgKUo/AR1MXbuWvS7bz3Yv7cl7ffL/jiEiKUo/AJ8tLDvBfs1ZzQb987ryot99xRCSFHbuOIKweQdzsP1LFpOcWk5+Twe++PoRAQBeNiYh/Ph50Tj2CuKitdXxv5lJKD1Xyl9tH0So73e9IIpLiQgH1COLqobnFvL2ulB9fNYAzurT0O46ISMygczo05Ll315fyuzc+5JohHbnhrK5+xxERASJzFoP2CDy348BR7pq+lD5tm/PzrwzSYHIikjDqDg29uvIj9h2pitvreloIzOxSM1tnZsVm9oPjPJ5hZjOijy8ws+5e5qmqqeU7UxdTWR3mTzcMIys9ZVskIpKAerXN5stnduKFxSV84YG3+NWctRwo974geFYIzCwIPAxcBgwAxprZgHqrTQD2O+d6A78DHvAqD8DPZ69hydYDPPjVM+iV39zLlxIROWUZoSC/+8YQXv/ueVzYvy1/fHsD5z4wl9+8vo6y8mrPXtfLPYIRQLFzbqNzru5ScAYAAAehSURBVAqYDoypt84Y4Ono7eeBi82jYzV/X7aDp+Zt5lvn9OCKwR28eAkRkUbRp10OD10/lNfuOo/z+ubxh7eKOfeBt5i1bIcnr+dlIegEbIu5XxJddtx1nHM1QBnQpv4vMrOJZlZkZkWlpaWfKUzr7HQuGdCOH17e/zM9X0Qk3vq1z+GP44bx6l1f4OzebejWOsuT12kSB8mdc5OByQAFBQWfqZ1+Tu88zumd16i5RETi4bQOLXj0xgLPfr+XewTbgS4x9ztHlx13HTMLAbnAXg8ziYhIPV4WgkKgj5n1MLN04DpgVr11ZgE3RW9/FXjLORffE2hFRFKcZ4eGnHM1ZnYHMAcIAk8651aZ2f1AkXNuFvAE8KyZFQP7iBQLERGJI097BM652cDsesvui7ldAXzNywwiIvLpUu7KYhER+VcqBCIiKU6FQEQkxakQiIikOGtqZ2uaWSmw5TM+PQ/Y04hxGotynRrlOnWJmk25Ts3nydXNOXfcSdmbXCH4PMysyDnn3eV5n5FynRrlOnWJmk25To1XuXRoSEQkxakQiIikuFQrBJP9DnACynVqlOvUJWo25To1nuRKqR6BiIh8UqrtEYiISD0qBCIiKS5pCoGZXWpm68ys2Mx+cJzHM8xsRvTxBWbWPeaxH0aXrzOzL8U5191mttrMlpvZm2bWLeaxsJktjf7UH8Lb61w3m1lpzOvfEvPYTWa2PvpzU/3nepzrdzGZPjSzAzGPebm9njSz3Wa28gSPm5n9Ppp7uZkNjXnMk+3VgEzjollWmNk8Mzsj5rHN0eVLzayosTKdQrYLzKws5u91X8xjn/oe8DjXPTGZVkbfU62jj3myzcysi5nNjX4OrDKzu46zjrfvL+dck/8hMsz1BqAnkA4sAwbUW+fbwCPR29cBM6K3B0TXzwB6RH9PMI65LgSyorcn1eWK3j/s4/a6GXjoOM9tDWyM/rdV9HareOWqt/6dRIY393R7RX/3ecBQYOUJHr8ceBUwYCSwIA7b62SZzq57LeCyukzR+5uBPB+31wXAy5/3PdDYueqtexWROVI83WZAB2Bo9HYO8OFx/j16+v5Klj2CEUCxc26jc64KmA6MqbfOGODp6O3ngYvNzKLLpzvnKp1zm4Di6O+LSy7n3FznXHn07nwiM7l5rSHb60S+BPzDObfPObcf+AdwqU+5xgLTGum1P5Vz7h0ic2acyBjgGRcxH2hpZh3wcHudLJNzbl70NSF+76261z7Z9jqRz/PebOxccXl/Oed2OucWR28fAtbwyfndPX1/JUsh6ARsi7lfwic35LF1nHM1QBnQpoHP9TJXrAlEqn6dTDMrMrP5ZnZNI2U6lVzXRndDnzezumlHE2J7RQ+h9QDeilns1fZqiBNl93J7nYr67y0HvG5mi8xsog95AEaZ2TIze9XMBkaXJcT2MrMsIh+oL8Qs9nybWeSQ9ZnAgnoPefr+ahKT16cCM7sBKADOj1nczTm33cx6Am+Z2Qrn3IY4Rfo7MM05V2lmtxHZm7ooTq/dENcBzzvnwjHL/NxeCcvMLiRSCM6NWXxudFu1Bf5hZmuj35bjZTGRv9dhM7sceBHoE8fXP5mrgPedc7F7D55uMzNrTqTwfNc5d7Cxfm9DJMsewXagS8z9ztFlx13HzEJALrC3gc/1MhdmNhq4F7jaOVdZt9w5tz36343A20S+KcQll3Nub0yWx4FhDX2ul7liXEe93XYPt1dDnCi7l9vrpMxsMJG/3xjn3N665THbajfwNxrvcGiDOOcOOucOR2/PBtLMLA+ft1eMT3t/Nfo2M7M0IkVginPur8dZxdv3V2M3Pvz4IbJns5HIoYK6BtPAeut8h39tFs+M3h7IvzaLN9J4zeKG5DqTSHOsT73lrYCM6O08YD2N1DRrYK4OMbe/DMx3HzenNkXztYrebh2vXNH1+hNp3Fk8tlfMa3TnxM3PK/jXZt5Cr7dXAzJ1JdLzOrve8mwgJ+b2PODSxtxWDcjWvu7vR+QDdWt02zXoPeBVrujjuUT6CNnx2GbR/+9ngP/9lHU8fX816h/ezx8iXfUPiXyo3htddj+Rb9kAmcBfov8wFgI9Y557b/R564DL4pzrDWAXsDT6Myu6/GxgRfQfwgpgQpxz/QJYFX39uUD/mOd+K7odi4FvxjNX9P5PgV/We57X22sasBOoJnIcdgJwO3B79HEDHo7mXgEUeL29GpDpcWB/zHurKLq8Z3Q7LYv+je9tzG3VwGx3xLy/5hNTrI73HohXrug6NxM5gST2eZ5tMyKH7BywPOZvdXk8318aYkJEJMUlS49AREQ+IxUCEZEUp0IgIpLiVAhERFKcCoGISIpTIRARSXEqBCIiKU6FQORzMrPh0cH5Ms0sOzqm/Ol+5xJpKF1QJtIIzOy/iVy93gwocc79wudIIg2mQiDSCMwsHSgEKogMlxA+yVNEEoYODYk0jjZAcyIzTGX6nEXklGiPQKQRROdInk5k1MwOzrk7fI4k0mCamEbkczKz8UC1c26qmQWBeWZ2kXPurZM9VyQRaI9ARCTFqUcgIpLiVAhERFKcCoGISIpTIRARSXEqBCIiKU6FQEQkxakQiIikuP8PMe82xrtZbxcAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "plt.plot(x, u.data[-1])\n", "plt.xlabel('x')\n", "plt.ylabel('u')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Verification\n", "\n", "\n", "Before implementing the algorithm, it is convenient to add a source\n", "term to the PDE ([1](#wave:pde1)), since that gives us more freedom in\n", "finding test problems for verification. Physically, a source term acts\n", "as a generator for waves in the interior of the domain.\n", "\n", "## A slightly generalized model problem\n", "
\n", "\n", "We now address the following extended initial-boundary value problem\n", "for one-dimensional wave phenomena:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u_{tt} = c^2 u_{xx} + f(x,t), \\quad x\\in (0,L),\\ t\\in (0,T]\n", "\\label{wave:pde2} \\tag{17}\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{wave:pde2:ic:u} \\tag{18}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u_t(x,0) = V(x), \\quad x\\in [0,L]\n", "\\label{wave:pde2:ic:ut} \\tag{19}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(0,t) = 0, \\quad t>0\n", "\\label{wave:pde2:bc:0} \\tag{20}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(L,t) = 0, \\quad t>0\n", "\\label{wave:pde2:bc:L} \\tag{21}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sampling the PDE at $(x_i,t_n)$ and using the same finite difference\n", "approximations as above, yields" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "[D_tD_t u = c^2 D_xD_x u + f]^{n}_i\n", "\\label{wave:pde2:fdop} \\tag{22}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Writing this out and solving for the unknown $u^{n+1}_i$ results in" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u^{n+1}_i = -u^{n-1}_i + 2u^n_i + C^2\n", "(u^{n}_{i+1}-2u^{n}_{i} + u^{n}_{i-1}) + \\Delta t^2 f^n_i\n", "\\label{wave:pde2:step3b} \\tag{23}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The equation for the first time step must be rederived. The discretization\n", "of the initial condition $u_t = V(x)$ at $t=0$\n", "becomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "[D_{2t}u = V]^0_i\\quad\\Rightarrow\\quad u^{-1}_i = u^{1}_i - 2\\Delta t V_i,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which, when inserted in ([23](#wave:pde2:step3b)) for $n=0$, gives\n", "the special formula" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u^{1}_i = u^0_i + \\Delta t V_i + {\\frac{1}{2}}\n", "C^2\n", "\\left(u^{0}_{i+1}-2u^{0}_{i} + u^{0}_{i-1}\\right) + \\frac{1}{2}\\Delta t^2 f^0_i\n", "\\label{wave:pde2:step3c} \\tag{24}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using an analytical solution of physical significance\n", "
\n", "\n", "\n", "Many wave problems feature sinusoidal oscillations in time\n", "and space. For example, the original PDE problem\n", "([1](#wave:pde1))-([5](#wave:pde1:bc:L)) allows an exact solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u_e(x,t) = A\\sin\\left(\\frac{\\pi}{L}x\\right)\n", "\\cos\\left(\\frac{\\pi}{L}ct\\right)\n", "\\label{wave:pde2:test:ue} \\tag{25}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This $u_e$ fulfills the PDE with $f=0$, boundary conditions\n", "$u_e(0,t)=u_e(L,t)=0$, as well as initial\n", "conditions $I(x)=A\\sin\\left(\\frac{\\pi}{L}x\\right)$ and $V=0$.\n", "\n", "**How to use exact solutions for verification.**\n", "\n", "It is common to use such exact solutions of physical interest\n", "to verify implementations. However, the numerical\n", "solution $u^n_i$ will only be an approximation to $u_e(x_i,t_n)$.\n", "We have no knowledge of the precise size of the error in\n", "this approximation, and therefore we can never know if discrepancies\n", "between $u^n_i$ and $u_e(x_i,t_n)$ are caused\n", "by mathematical approximations or programming errors.\n", "In particular, if plots of the computed solution $u^n_i$ and\n", "the exact one ([25](#wave:pde2:test:ue)) look similar, many\n", "are tempted to claim that the implementation works. However,\n", "even if color plots look nice and the accuracy is \"deemed good\",\n", "there can still be serious programming errors present!\n", "\n", "The only way to use exact physical solutions like\n", "([25](#wave:pde2:test:ue)) for serious and thorough verification is to\n", "run a series of simulations on finer and finer meshes, measure the\n", "integrated error in each mesh, and from this information estimate the\n", "empirical convergence rate of the method.\n", "\n", "\n", "\n", "An introduction to the computing of convergence rates is given in Section 3.1.6\n", " in [[Langtangen_decay]](#Langtangen_decay).\n", "There is also a detailed example on computing convergence rates in\n", "the [verification section](../01_vib/vib_undamped.ipynb#vib:ode1:verify) of the Vibration ODEs chapter.\n", "\n", "\n", "In the present problem, one expects the method to have a convergence rate\n", "of 2 (see the section [Analysis of the difference equations](wave_analysis.ipynb)), so if the computed rates\n", "are close to 2 on a sufficiently fine mesh, we have good evidence that\n", "the implementation is free of programming mistakes.\n", "\n", "## Manufactured solution and estimation of convergence rates\n", "
\n", "\n", "### Specifying the solution and computing corresponding data\n", "\n", "One problem with the exact solution ([25](#wave:pde2:test:ue)) is\n", "that it requires a simplification (${V}=0, f=0$) of the implemented problem\n", "([17](#wave:pde2))-([21](#wave:pde2:bc:L)). An advantage of using\n", "a *manufactured solution* is that we can test all terms in the\n", "PDE problem. The idea of this approach is to set up some chosen\n", "solution and fit the source term, boundary conditions, and initial\n", "conditions to be compatible with the chosen solution.\n", "Given that our boundary conditions in the implementation are\n", "$u(0,t)=u(L,t)=0$, we must choose a solution that fulfills these\n", "conditions. One example is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u_e(x,t) = x(L-x)\\sin t\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Inserted in the PDE $u_{tt}=c^2u_{xx}+f$ we get" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "-x(L-x)\\sin t = -c^2 2\\sin t + f\\quad\\Rightarrow f = (2c^2 - x(L-x))\\sin t\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The initial conditions become" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "u(x,0) =& I(x) = 0,\\\\ \n", "u_t(x,0) &= V(x) = x(L-x)\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Defining a single discretization parameter\n", "\n", "To verify the code, we compute the convergence rates in a series of\n", "simulations, letting each simulation use a finer mesh than the\n", "previous one. Such empirical estimation of convergence rates relies on\n", "an assumption that some measure $E$ of the numerical error is related\n", "to the discretization parameters through" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E = C_t\\Delta t^r + C_x\\Delta x^p,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $C_t$, $C_x$, $r$, and $p$ are constants. The constants $r$ and\n", "$p$ are known as the *convergence rates* in time and space,\n", "respectively. From the accuracy in the finite difference\n", "approximations, we expect $r=p=2$, since the error terms are of order\n", "$\\Delta t^2$ and $\\Delta x^2$. This is confirmed by truncation error\n", "analysis and other types of analysis.\n", "\n", "By using an exact solution of the PDE problem, we will next compute\n", "the error measure $E$ on a sequence of refined meshes and see if\n", "the rates $r=p=2$ are obtained. We will not be concerned with estimating\n", "the constants $C_t$ and $C_x$, simply because we are not interested in\n", "their values.\n", "\n", "mathcal{I}_t is advantageous to introduce a single discretization parameter\n", "$h=\\Delta t=\\hat c \\Delta x$ for some constant $\\hat c$. Since\n", "$\\Delta t$ and $\\Delta x$ are related through the Courant number,\n", "$\\Delta t = C\\Delta x/c$, we set $h=\\Delta t$, and then $\\Delta x =\n", "hc/C$. Now the expression for the error measure is greatly\n", "simplified:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E = C_t\\Delta t^r + C_x\\Delta x^r =\n", "C_t h^r + C_x\\left(\\frac{c}{C}\\right)^r h^r\n", "= Dh^r,\\quad D = C_t+C_x\\left(\\frac{c}{C}\\right)^r \n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Computing errors\n", "\n", "We choose an initial discretization parameter $h_0$ and run\n", "experiments with decreasing $h$: $h_i=2^{-i}h_0$, $i=1,2,\\ldots,m$.\n", "Halving $h$ in each experiment is not necessary, but it is a common\n", "choice. For each experiment we must record $E$ and $h$. Standard\n", "choices of error measure are the $\\ell^2$ and $\\ell^\\infty$ norms of the\n", "error mesh function $e^n_i$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "E = ||e^n_i||_{\\ell^2} = \\left( \\Delta t\\Delta x\n", "\\sum_{n=0}^{N_t}\\sum_{i=0}^{N_x}\n", "(e^n_i)^2\\right)^{\\frac{1}{2}},\\quad e^n_i = u_e(x_i,t_n)-u^n_i,\n", "\\label{wave:pde2:fd:MMS:E:l2} \\tag{26}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "E = ||e^n_i||_{\\ell^\\infty} = \\max_{i,n} |e^n_i|\n", "\\label{wave:pde2:fd:MMS:E:linf} \\tag{27}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In Python, one can compute $\\sum_{i}(e^{n}_i)^2$ at each time step\n", "and accumulate the value in some sum variable, say `e2_sum`. At the\n", "final time step one can do `sqrt(dt*dx*e2_sum)`. For the\n", "$\\ell^\\infty$ norm one must compare the maximum error at a time level\n", "(`e.max()`) with the global maximum over the time domain: `e_max = max(e_max, e.max())`.\n", "\n", "An alternative error measure is to use a spatial norm at one time step\n", "only, e.g., the end time $T$ ($n=N_t$):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "E = ||e^n_i||_{\\ell^2} = \\left( \\Delta x\\sum_{i=0}^{N_x}\n", "(e^n_i)^2\\right)^{\\frac{1}{2}},\\quad e^n_i = u_e(x_i,t_n)-u^n_i,\n", "\\label{_auto6} \\tag{28}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "E = ||e^n_i||_{\\ell^\\infty} = \\max_{0\\leq i\\leq N_x} |e^{n}_i|\n", "\\label{_auto7} \\tag{29}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The important point is that the error measure ($E$) for the simulation is represented by a single number.\n", "\n", "### Computing rates\n", "\n", "Let $E_i$ be the error measure in experiment (mesh) number $i$ \n", "(not to be confused with the spatial index $i$) and\n", "let $h_i$ be the corresponding discretization parameter ($h$).\n", "With the error model $E_i = Dh_i^r$, we can\n", "estimate $r$ by comparing two consecutive\n", "experiments:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "E_{i+1}& =D h_{i+1}^{r},\\\\ \n", "E_{i}& =D h_{i}^{r}\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dividing the two equations eliminates the (uninteresting) constant $D$.\n", "Thereafter, solving for $r$ yields" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "r = \\frac{\\ln E_{i+1}/E_{i}}{\\ln h_{i+1}/h_{i}}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since $r$ depends on $i$, i.e., which simulations we compare,\n", "we add an index to $r$: $r_i$, where $i=0,\\ldots,m-2$, if we\n", "have $m$ experiments: $(h_0,E_0),\\ldots,(h_{m-1}, E_{m-1})$.\n", "\n", "In our present discretization of the wave equation we expect $r=2$, and\n", "hence the $r_i$ values should converge to 2 as $i$ increases.\n", "\n", "## Constructing an exact solution of the discrete equations\n", "
\n", "\n", "With a manufactured or known analytical solution, as outlined above,\n", "we can estimate convergence rates and see if they have the correct\n", "asymptotic behavior. Experience shows that this is a quite good\n", "verification technique in that many common bugs will destroy the\n", "convergence rates. A significantly better test though,\n", "would be to check that the\n", "numerical solution is exactly what it should be. This will in general\n", "require exact knowledge of the numerical error, which we do not normally have\n", "(although we in the section [Analysis of the difference equations](wave_analysis.ipynb) establish such knowledge\n", "in simple cases).\n", "However, it is possible to look for solutions where we can show that\n", "the numerical error vanishes, i.e., the solution of the original continuous\n", "PDE problem is\n", "also a solution of the discrete equations. This property often arises\n", "if the exact solution of the PDE\n", "is a lower-order polynomial. (Truncation error\n", "analysis leads to error measures that involve derivatives of the\n", "exact solution. In the present problem, the truncation error involves\n", "4th-order derivatives of $u$ in space and time. Choosing $u$\n", "as a polynomial of degree three or less\n", "will therefore lead to vanishing error.)\n", "\n", "We shall now illustrate the construction of an exact solution to both the\n", "PDE itself and the discrete equations.\n", "Our chosen manufactured solution is quadratic in space\n", "and linear in time. More specifically, we set" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u_e (x,t) = x(L-x)(1+{\\frac{1}{2}}t),\n", "\\label{wave:pde2:fd:verify:quadratic:uex} \\tag{30}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which by insertion in the PDE leads to $f(x,t)=2(1+t)c^2$. This $u_e$\n", "fulfills the boundary conditions $u=0$ and demands $I(x)=x(L-x)$\n", "and $V(x)={\\frac{1}{2}}x(L-x)$.\n", "\n", "To realize that the chosen $u_e$ is also an exact\n", "solution of the discrete equations,\n", "we first remind ourselves that $t_n=n\\Delta t$ so that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\lbrack D_tD_t t^2\\rbrack^n = \\frac{t_{n+1}^2 - 2t_n^2 + t_{n-1}^2}{\\Delta t^2}\n", "= (n+1)^2 -2n^2 + (n-1)^2 = 2,\n", "\\label{_auto8} \\tag{31}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\lbrack D_tD_t t\\rbrack^n = \\frac{t_{n+1} - 2t_n + t_{n-1}}{\\Delta t^2}\n", "= \\frac{((n+1) -2n + (n-1))\\Delta t}{\\Delta t^2} = 0\n", "\\label{_auto9} \\tag{32}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hence," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "[D_tD_t u_e]^n_i = x_i(L-x_i)[D_tD_t (1+{\\frac{1}{2}}t)]^n =\n", "x_i(L-x_i){\\frac{1}{2}}[D_tD_t t]^n = 0\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly, we get that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\lbrack D_xD_x u_e\\rbrack^n_i &=\n", "(1+{\\frac{1}{2}}t_n)\\lbrack D_xD_x (xL-x^2)\\rbrack_i\\\\ \n", "& =\n", "(1+{\\frac{1}{2}}t_n)\\lbrack LD_xD_x x - D_xD_x x^2\\rbrack_i \\\\ \n", "&= -2(1+{\\frac{1}{2}}t_n)\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, $f^n_i = 2(1+{\\frac{1}{2}}t_n)c^2$, which results in" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "[D_tD_t u_e - c^2D_xD_xu_e - f]^n_i = 0 +\n", "c^2 2(1 + {\\frac{1}{2}}t_{n}) +\n", "2(1+{\\frac{1}{2}}t_n)c^2 = 0\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Moreover, $u_e(x_i,0)=I(x_i)$,\n", "$\\partial u_e/\\partial t = V(x_i)$ at $t=0$, and\n", "$u_e(x_0,t)=u_e(x_{N_x},0)=0$. Also the modified scheme for the\n", "first time step is fulfilled by $u_e(x_i,t_n)$.\n", "\n", "Therefore, the exact solution $u_e(x,t)=x(L-x)(1+t/2)$ of the PDE\n", "problem is also an exact solution of the discrete problem. This means\n", "that we know beforehand what numbers the numerical algorithm should\n", "produce. We can use this fact to check that the computed $u^n_i$\n", "values from an implementation equals $u_e(x_i,t_n)$, within machine\n", "precision. This result is valid *regardless of the mesh spacings*\n", "$\\Delta x$ and $\\Delta t$! Nevertheless, there might be stability\n", "restrictions on $\\Delta x$ and $\\Delta t$, so the test can only be run\n", "for a mesh that is compatible with the stability criterion (which in\n", "the present case is $C\\leq 1$, to be derived later).\n", "\n", "**Notice.**\n", "\n", "A product of quadratic or linear expressions in the various\n", "independent variables, as shown above, will often fulfill both the\n", "PDE problem and the discrete equations, and can therefore be very useful\n", "solutions for verifying implementations.\n", "\n", "However, for 1D wave\n", "equations of the type $u_{tt}=c^2u_{xx}$ we shall see that there is always\n", "another much more powerful way of generating exact\n", "solutions (which consists in just setting $C=1$ (!), as shown in\n", "the section [Analysis of the difference equations](wave_analysis.ipynb))." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.3" } }, "nbformat": 4, "nbformat_minor": 4 }