.. !split Wave equations (2) =========================== .. _trunc:wave:1D: Linear wave equation in 1D -------------------------- The standard, linear wave equation in 1D for a function :math:`u(x,t)` reads .. _Eq:trunc:wave:pde1D: .. math:: \tag{721} \frac{\partial^2 u}{\partial t^2} = c^2\frac{\partial^2 u}{\partial x^2} + f(x,t),\quad x\in (0, L),\ t\in (0,T], where :math:`c` is the constant wave velocity of the physical medium in :math:`[0,L]`. The equation can also be more compactly written as .. _Eq:trunc:wave:pde1D:v2: .. math:: \tag{722} u_{tt} = c^2u_{xx} + f,\quad x\in (0, L),\ t\in (0,T], Centered, second-order finite differences are a natural choice for discretizing the derivatives, leading to .. _Eq:trunc:wave:pde1D:fd: .. math:: \tag{723} [D_t D_t u = c^2 D_xD_x u + f]^n_i {\thinspace .} Inserting the exact solution :math:`{u_{\small\mbox{e}}}(x,t)` in :ref:`(723) ` makes this function fulfill the equation if we add the term :math:`R`: .. _Eq:trunc:wave:pde1D:fd:R: .. math:: \tag{724} [D_t D_t {u_{\small\mbox{e}}} = c^2 D_xD_x {u_{\small\mbox{e}}} + f + R]^n_i Our purpose is to calculate the truncation error :math:`R`. From :ref:`(670) `-:ref:`(671) ` we have that .. math:: [D_t D_t{u_{\small\mbox{e}}}]_i^n = {u_{\small\mbox{e}, tt}}(x_i,t_n) + \frac{1}{12}{u_{\small\mbox{e}, tttt}}(x_i,t_n)\Delta t^2 + {\mathcal{O}(\Delta t^4)}, when we use a notation taking into account that :math:`{u_{\small\mbox{e}}}` is a function of two variables and that derivatives must be partial derivatives. The notation :math:`{u_{\small\mbox{e}, tt}}` means :math:`\partial^2{u_{\small\mbox{e}}} /\partial t^2`. The same formula may also be applied to the :math:`x`-derivative term: .. math:: [D_xD_x{u_{\small\mbox{e}}}]_i^n = {u_{\small\mbox{e}, xx}}(x_i,t_n) + \frac{1}{12}{u_{\small\mbox{e}, xxxx}}(x_i,t_n)\Delta x^2 + {\mathcal{O}(\Delta x^4)}, Equation :ref:`(724) ` now becomes .. math:: \begin{align*} {u_{\small\mbox{e}, tt}} + \frac{1}{12}{u_{\small\mbox{e}, tttt}}(x_i,t_n)\Delta t^2 &= c^2{u_{\small\mbox{e}, xx}} + c^2\frac{1}{12}{u_{\small\mbox{e}, xxxx}}(x_i,t_n)\Delta x^2 + f(x_i,t_n) + \\ & \quad {\mathcal{O}(\Delta t^4,\Delta x^4)} + R^n_i {\thinspace .} \end{align*} Because :math:`{u_{\small\mbox{e}}}` fulfills the partial differential equation (PDE) :ref:`(722) `, the first, third, and fifth term cancel out, and we are left with .. _Eq:trunc:wave:1D:R: .. math:: \tag{725} R^n_i = \frac{1}{12}{u_{\small\mbox{e}, tttt}}(x_i,t_n)\Delta t^2 - c^2\frac{1}{12}{u_{\small\mbox{e}, xxxx}}(x_i,t_n)\Delta x^2 + {\mathcal{O}(\Delta t^4,\Delta x^4)}, showing that the scheme :ref:`(723) ` is of second order in the time and space mesh spacing. .. _trunc:wave:1D:corr: Finding correction terms ------------------------ Can we add correction terms to the PDE and increase the order of :math:`R^n_i` in :ref:`(725) `? The starting point is .. _Eq:trunc:wave:pde1D:fd:R2: .. math:: \tag{726} [D_t D_t {u_{\small\mbox{e}}} = c^2 D_xD_x {u_{\small\mbox{e}}} + f + C + R]^n_i From the previous analysis we simply get :ref:`(725) ` again, but now with :math:`C`: .. _Eq:trunc:wave:1D:R:C: .. math:: \tag{727} R^n_i + C_i^n = \frac{1}{12}{u_{\small\mbox{e}, tttt}}(x_i,t_n)\Delta t^2 - c^2\frac{1}{12}{u_{\small\mbox{e}, xxxx}}(x_i,t_n)\Delta x^2 + {\mathcal{O}(\Delta t^4,\Delta x^4)}{\thinspace .} The idea is to let :math:`C_i^n` cancel the :math:`\Delta t^2` and :math:`\Delta x^2` terms to make :math:`R^n_i = {\mathcal{O}(\Delta t^4,\Delta x^4)}`: .. math:: C_i^n = \frac{1}{12}{u_{\small\mbox{e}, tttt}}(x_i,t_n)\Delta t^2 - c^2\frac{1}{12}{u_{\small\mbox{e}, xxxx}}(x_i,t_n)\Delta x^2{\thinspace .} Essentially, it means that we add a new term .. math:: C = \frac{1}{12}\left( u_{tttt}\Delta t^2 - c^2u_{xxxx}\Delta x^2\right), to the right-hand side of the PDE. We must either discretize these 4th-order derivatives directly or rewrite them in terms of lower-order derivatives with the aid of the PDE. The latter approach is more feasible. From the PDE we have the operator equality .. math:: \frac{\partial^2}{\partial t^2} = c^2\frac{\partial^2}{\partial x^2}, so .. math:: u_{tttt} = c^2u_{xxtt},\quad u_{xxxx} = c^{-2}u_{ttxx}{\thinspace .} Assuming :math:`u` is smooth enough, so that :math:`u_{xxtt}=u_{ttxx}`, these relations lead to .. math:: C = \frac{1}{12}((c^2\Delta t^2 - \Delta x^2)u_{xx})_{tt}{\thinspace .} A natural discretization is .. math:: C^n_i = \frac{1}{12}((c^2\Delta t^2 - \Delta x^2) [D_xD_xD_tD_t u]^n_i{\thinspace .} Writing out :math:`[D_xD_xD_tD_t u]^n_i` as :math:`[D_xD_x (D_tD_t u)]^n_i` gives .. math:: \begin{align*} \frac{1}{\Delta t^2}\biggl( &\frac{u^{n+1}_{i+1} - 2u^{n}_{i+1} + u^{n-1}_{i+1}}{\Delta x^2} -2\\ &\frac{u^{n+1}_{i} - 2u^{n}_{i} + u^{n-1}_{i}}{\Delta x^2} + &\frac{u^{n+1}_{i-1} - 2u^{n}_{i-1} + u^{n-1}_{i-1}}{\Delta x^2} \biggr) \end{align*} Now the unknown values :math:`u^{n+1}_{i+1}`, :math:`u^{n+1}_{i}`, and :math:`u^{n+1}_{i-1}` are *coupled*, and we must solve a tridiagonal system to find them. This is in principle straightforward, but it results in an implicit finite difference schemes, while we had a convenient explicit scheme without the correction terms. .. _trunc:wave:1D:varcoeff: Extension to variable coefficients (2) ----------------------------------------------- Now we address the variable coefficient version of the linear 1D wave equation, .. math:: \frac{\partial^2 u}{\partial t^2} = \frac{\partial}{\partial x} \left( \lambda(x)\frac{\partial u}{\partial x}\right), or written more compactly as .. _Eq:trunc:wave:1D:varcoeff:pde: .. math:: \tag{728} u_{tt} = (\lambda u_x)_x{\thinspace .} The discrete counterpart to this equation, using arithmetic mean for :math:`\lambda` and centered differences, reads .. _Eq:trunc:wave:1D:varcoeff:fd: .. math:: \tag{729} [D_t D_t u = D_x \overline{\lambda}^{x}D_x u]^n_i{\thinspace .} The truncation error is the residual :math:`R` in the equation .. _Eq:trunc:wave:1D:varcoef:fd:R: .. math:: \tag{730} [D_t D_t {u_{\small\mbox{e}}} = D_x \overline{\lambda}^{x}D_x {u_{\small\mbox{e}}} + R]^n_i{\thinspace .} The difficulty with :ref:`(730) ` is how to compute the truncation error of the term :math:`[D_x \overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_i`. We start by writing out the outer operator: .. _Eq:trunc:wave:1D:varcoeff:outer: .. math:: \tag{731} [D_x \overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_i = \frac{1}{\Delta x}\left( [\overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_{i+\frac{1}{2}} - [\overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_{i-\frac{1}{2}} \right). With the aid of :ref:`(658) `-:ref:`(659) ` and :ref:`(674) `-:ref:`(675) ` we have .. math:: \begin{align*} \lbrack D_x {u_{\small\mbox{e}}} \rbrack^n_{i+\frac{1}{2}} & = {u_{\small\mbox{e}, x}}(x_{i+\frac{1}{2}},t_n) + \frac{1}{24}{u_{\small\mbox{e}, xxx}}(x_{i+\frac{1}{2}},t_n)\Delta x^2 + {\mathcal{O}(\Delta x^4)},\\ \lbrack\overline{\lambda}^{x}\rbrack_{i+\frac{1}{2}} &= \lambda(x_{i+\frac{1}{2}}) + \frac{1}{8}\lambda''(x_{i+\frac{1}{2}})\Delta x^2 + {\mathcal{O}(\Delta x^4)},\\ [\overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_{i+\frac{1}{2}} &= (\lambda(x_{i+\frac{1}{2}}) + \frac{1}{8}\lambda''(x_{i+\frac{1}{2}})\Delta x^2 + {\mathcal{O}(\Delta x^4)})\times\\ &\quad ({u_{\small\mbox{e}, x}}(x_{i+\frac{1}{2}},t_n) + \frac{1}{24}{u_{\small\mbox{e}, xxx}}(x_{i+\frac{1}{2}},t_n)\Delta x^2 + {\mathcal{O}(\Delta x^4)})\\ &= \lambda(x_{i+\frac{1}{2}}){u_{\small\mbox{e}, x}}(x_{i+\frac{1}{2}},t_n) + \lambda(x_{i+\frac{1}{2}}) \frac{1}{24}{u_{\small\mbox{e}, xxx}}(x_{i+\frac{1}{2}},t_n)\Delta x^2 + \\ &\quad {u_{\small\mbox{e}, x}}(x_{i+\frac{1}{2}}) \frac{1}{8}\lambda''(x_{i+\frac{1}{2}})\Delta x^2 +{\mathcal{O}(\Delta x^4)}\\ &= [\lambda {u_{\small\mbox{e}, x}}]^n_{i+\frac{1}{2}} + G^n_{i+\frac{1}{2}}\Delta x^2 +{\mathcal{O}(\Delta x^4)}, \end{align*} where we have introduced the short form .. math:: G^n_{i+\frac{1}{2}} = (\frac{1}{24}{u_{\small\mbox{e}, xxx}}(x_{i+\frac{1}{2}},t_n)\lambda((x_{i+\frac{1}{2}}) + {u_{\small\mbox{e}, x}}(x_{i+\frac{1}{2}},t_n) \frac{1}{8}\lambda''(x_{i+\frac{1}{2}}))\Delta x^2{\thinspace .} Similarly, we find that .. math:: \lbrack\overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}\rbrack^n_{i-\frac{1}{2}} = [\lambda {u_{\small\mbox{e}, x}}]^n_{i-\frac{1}{2}} + G^n_{i-\frac{1}{2}}\Delta x^2 +{\mathcal{O}(\Delta x^4)}{\thinspace .} Inserting these expressions in the outer operator :ref:`(731) ` results in .. math:: \begin{align*} \lbrack D_x \overline{\lambda}^{x}D_x {u_{\small\mbox{e}}} \rbrack^n_i &= \frac{1}{\Delta x}( [\overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_{i+\frac{1}{2}} - [\overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_{i-\frac{1}{2}} )\\ &= \frac{1}{\Delta x}( [\lambda {u_{\small\mbox{e}, x}}]^n_{i+\frac{1}{2}} + G^n_{i+\frac{1}{2}}\Delta x^2 - [\lambda {u_{\small\mbox{e}, x}}]^n_{i-\frac{1}{2}} - G^n_{i-\frac{1}{2}}\Delta x^2 + {\mathcal{O}(\Delta x^4)} )\\ &= [D_x \lambda {u_{\small\mbox{e}, x}}]^n_i + [D_x G]^n_i\Delta x^2 + {\mathcal{O}(\Delta x^4)}{\thinspace .} \end{align*} The reason for :math:`{\mathcal{O}(\Delta x^4)}` in the remainder is that there are coefficients in front of this term, say :math:`H\Delta x^4`, and the subtraction and division by :math:`\Delta x` results in :math:`[D_x H]^n_i\Delta x^4`. We can now use :ref:`(658) `-:ref:`(659) ` to express the :math:`D_x` operator in :math:`[D_x \lambda {u_{\small\mbox{e}, x}}]^n_i` as a derivative and a truncation error: .. math:: [D_x \lambda {u_{\small\mbox{e}, x}}]^n_i = \frac{\partial}{\partial x}\lambda(x_i){u_{\small\mbox{e}, x}}(x_i,t_n) + \frac{1}{24}(\lambda{u_{\small\mbox{e}, x}})_{xxx}(x_i,t_n)\Delta x^2 + {\mathcal{O}(\Delta x^4)}{\thinspace .} Expressions like :math:`[D_x G]^n_i\Delta x^2` can be treated in an identical way, .. math:: [D_x G]^n_i\Delta x^2 = G_x(x_i,t_n)\Delta x^2 + \frac{1}{24}G_{xxx}(x_i,t_n)\Delta x^4 + {\mathcal{O}(\Delta x^4)}{\thinspace .} There will be a number of terms with the :math:`\Delta x^2` factor. We lump these now into :math:`{\mathcal{O}(\Delta x^2)}`. The result of the truncation error analysis of the spatial derivative is therefore summarized as .. math:: [D_x \overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_i = \frac{\partial}{\partial x} \lambda(x_i){u_{\small\mbox{e}, x}}(x_i,t_n) + {\mathcal{O}(\Delta x^2)}{\thinspace .} After having treated the :math:`[D_tD_t{u_{\small\mbox{e}}}]^n_i` term as well, we achieve .. math:: R^n_i = {\mathcal{O}(\Delta x^2)} + \frac{1}{12}{u_{\small\mbox{e}, tttt}}(x_i,t_n)\Delta t^2 {\thinspace .} The main conclusion is that the scheme is of second-order in time and space also in this variable coefficient case. The key ingredients for second order are the centered differences and the arithmetic mean for :math:`\lambda`: all those building blocks feature second-order accuracy. 1D wave equation on a staggered mesh ------------------------------------ .. _trunc:wave:2D: Linear wave equation in 2D/3D ----------------------------- The two-dimensional extension of :ref:`(721) ` takes the form .. _Eq:trunc:wave:pde2D: .. math:: \tag{732} \frac{\partial^2 u}{\partial t^2} = c^2\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) + f(x,y,t),\quad (x,y)\in (0, L)\times (0,H),\ t\in (0,T], where now :math:`c(x,y)` is the constant wave velocity of the physical medium :math:`[0,L]\times [0,H]`. In the compact notation, the PDE :ref:`(732) ` can be written .. _Eq:trunc:wave:pde2D:v2: .. math:: \tag{733} u_{tt} = c^2(u_{xx} + u_{yy}) + f(x,y,t),\quad (x,y)\in (0, L)\times (0,H), \ t\in (0,T], in 2D, while the 3D version reads .. _Eq:trunc:wave:pde3D:v2: .. math:: \tag{734} u_{tt} = c^2(u_{xx} + u_{yy} + u_{zz}) + f(x,y,z,t), for :math:`(x,y,z)\in (0, L)\times (0,H)\times (0,B)` and :math:`t\in (0,T]`. Approximating the second-order derivatives by the standard formulas :ref:`(670) `-:ref:`(671) ` yields the scheme .. _Eq:_auto267: .. math:: \tag{735} [D_t D_t u = c^2(D_xD_x u + D_yD_y u) + f]^n_{i,j,k} {\thinspace .} The truncation error is found from .. _Eq:_auto268: .. math:: \tag{736} [D_t D_t {u_{\small\mbox{e}}} = c^2(D_xD_x {u_{\small\mbox{e}}} + D_yD_y {u_{\small\mbox{e}}}) + f + R]^n {\thinspace .} The calculations from the 1D case can be repeated to the terms in the :math:`y` and :math:`z` directions. Collecting terms that fulfill the PDE, we end up with .. _Eq:_auto269: .. math:: \tag{737} R^n_{i,j,k} = [\frac{1}{12}{u_{\small\mbox{e}, tttt}}\Delta t^2 - c^2\frac{1}{12}\left( {u_{\small\mbox{e}, xxxx}}\Delta x^2 + {u_{\small\mbox{e}, yyyy}}\Delta x^2 + {u_{\small\mbox{e}, zzzz}}\Delta z^2\right)]^n_{i,j,k} + .. math:: \quad {\mathcal{O}(\Delta t^4,\Delta x^4,\Delta y^4,\Delta z^4)}\nonumber {\thinspace .}