Explicit Stencil Computation Schemes Generated by Poisson's Formula for the 2D Wave Equation

TL;DR

This paper introduces a novel explicit stencil computation method for the 2D wave equation using Poisson's formula, polynomial interpolation, and exact integration, resulting in schemes with improved accuracy and stability.

Contribution

It presents a new approach to generate explicit stencil schemes for the 2D wave equation based on Poisson's formula, with demonstrated stability and accuracy improvements.

Findings

Constructed three explicit schemes with 5, 9, and 13 points.

Schemes show larger stability regions than conventional methods.

Simulation results confirm higher accuracy of the new schemes.

Abstract

A new approach to building explicit time-marching stencil computation schemes for the transient 2D acoustic wave equation is implemented. It is based on using Poisson's formula and its three time level modification combined with polynomial stencil interpolation of the solution at each time-step and exact integration. The time-stepping algorithm consists of two explicit stencil computation procedures: a first time-step procedure incorporating the initial conditions and a two-step scheme for the second and next time-steps. Three particular explicit stencil schemes (with five, nine and 13 space points) are constructed using this approach. Their stability regions are presented. Accuracy advantages of the new schemes in comparison with conventional finite-difference schemes are demonstrated by simulation using an exact benchmark solution.

Equations

$$\frac{{\partial }^{2}u}{\partial t}-c^2\left(\frac{{\partial }^{2}u}{\partial x_1^2}+\frac{{\partial }^{2}u}{\partial x_2^2}\right)=0,u=u(\mathbf{x},t),\mathbf{x}=(x_1,x_2)\text{mbox}in{\mathbb{R}}^2,$$

$$u{|}_{t=0}=u_0(\mathbf{x}),\frac{\partial u}{\partial t}{|}_{t=0}=v_0(\mathbf{x}).$$

$$u(\mathbf{x},t)=\frac{1}{2\pi c{t}^2}{\int }_{|\mathbf{y}-\mathbf{x}{|}^2<t^2}\frac{t{u}_0(\mathbf{y})+t\nabla u_0(\mathbf{y})\cdot (\mathbf{y}-\mathbf{x})+t^2{v}_0(\mathbf{y})}{\sqrt{c^2{t}^2-|\mathbf{y}-\mathbf{x}{|}^2}}d{y}_{1}d{y}_2$$

$$u(\mathbf{x},t)=\frac{1}{2\pi }{\int }_{|\mathbf{z}|<1}\frac{u_0(\mathbf{x}+ct\mathbf{z})+ct\nabla u_0(\mathbf{x}+ct\mathbf{z})\cdot \mathbf{z}}{\sqrt{1-|\mathbf{z}{|}^2}}d{z}_{1}d{z}_2+\frac{t}{2\pi }{\int }_{|\mathbf{z}|<1}\frac{v_0(\mathbf{x}+ct\mathbf{z})}{\sqrt{1-|\mathbf{z}{|}^2}}d{z}_{1}d{z}_2$$

$$\begin{gathered} u(\mathbf{x},-t)=\frac{1}{2\pi }{\int }_{|\mathbf{z}|<1}\frac{u_0(\mathbf{x}+ct\mathbf{z})+ct\nabla u_0(\mathbf{x}+ct\mathbf{z})\cdot \mathbf{z}}{\sqrt{1-|\mathbf{z}{|}^2}}d{z}_{1}d{z}_2 \\ -\frac{t}{2\pi }{\int }_{|\mathbf{z}|<1}\frac{v_0(\mathbf{x}+ct\mathbf{z})}{\sqrt{1-|\mathbf{z}{|}^2}}d{z}_{1}d{z}_2. \end{gathered}$$

$$u(\mathbf{x},t_{\ast }+\tau )+u(\mathbf{x},t_{\ast }-\tau )=\frac{1}{\pi }{\int }_{|\mathbf{z}|<1}\frac{u(\mathbf{x}+ct\mathbf{z},t_{\ast })+ct\nabla u(\mathbf{x}+ct\mathbf{z},t_{\ast })\cdot \mathbf{z}}{\sqrt{1-|\mathbf{z}{|}^2}}d{z}_{1}d{z}_2$$

$$A(\mathbf{x},\tau )f(\cdot )=\frac{1}{2\pi }{\int }_{|\mathbf{z}|<1}\frac{f(\mathbf{x}+c\tau \mathbf{z})+c\tau \nabla f(\mathbf{x}+c\tau \mathbf{z})\cdot \mathbf{z}}{\sqrt{1-|\mathbf{z}{|}^2}}d{z}_{1}d{z}_2,$$

$$B(\mathbf{x},\tau )f(\cdot )=\frac{\tau }{2\pi }{\int }_{|\mathbf{z}|<1}\frac{f(\mathbf{x}+c\tau \mathbf{z})}{\sqrt{1-|\mathbf{z}{|}^2}}d{z}_{1}d{z}_2$$

$$u_1(\mathbf{x})=A(\mathbf{x},\tau )u_0(\cdot )+B(\mathbf{x},\tau )v_0(\cdot ),\mathbf{x}\in {\mathbb{R}}^2$$

$$u_{k+1}(\mathbf{x})=2A(\mathbf{x},\tau )u_k(\cdot )-u_{k-1}(\mathbf{x}),\mathbf{x}\in {\mathbb{R}}^2,k=1,2,...$$

$${\mathbf{x}}_{i,j}+{\mathbf{x}}_{q_1,q_2}={\mathbf{x}}_{i+q_1,j+q_2},(q_1,q_2)\in Q_m.$$

$${\mathcal{M}}_m=\{x_1^{\alpha (q_1)}{x}_2^{\alpha (q_2)},(q_1,q_2)\in Q_m\}$$

$$\alpha (q)=\{\begin{array}{ll}2|q|-1& \text{if }q<0\\ 2q& \text{if }q\ge 0.\end{array}$$

$$q(\alpha )=(-1{)}^{\alpha }\left[\frac{\alpha +1}{2}\right]$$

$$\{{\mu }_1(\mathbf{x}),\dots ,{\mu }_m(\mathbf{x})\}.$$

$$\{{\mathbf{x}}_{(1)},...,{\mathbf{x}}_{(m)}\}.$$

$$D={\left[{\mu }_s({\mathbf{x}}_{(r)})\right]}_{m\times m}$$

$$C={\left[c_{sr}\right]}_{m\times m}=D^{-1}$$

$$L_s(\mathbf{x})=\sum _{r=1}^m{c}_{sr}{\mu }_r(\mathbf{x}),s=1,...,m.$$

$$s=\gamma (q_1,q_2)=g(\alpha (q_1),\alpha (q_2)).$$

$${\varphi }_{q_1,q_2}(\mathbf{x})=L_{\gamma (q_1,q_2)}(\mathbf{x}),$$

$$\tilde{f}(\mathbf{x})=\sum _{(q_1,q_2)\in Q_m}{f}_{i+q_1,j+q_2}{\varphi }_{q_1{q}_2}(\mathbf{x})$$

$$s=g({\alpha }_1,{\alpha }_2)=\frac{({\alpha }_1+{\alpha }_2)({\alpha }_1+{\alpha }_2+1)}{2}+\{\begin{array}{ll}{\alpha }_1-{\alpha }_2& \text{if }{\alpha }_2<{\alpha }_1,\\ {\alpha }_2-{\alpha }_1+1& \text{if }{\alpha }_2\ge {\alpha }_1.\end{array}$$

$$\mu (\mathbf{y})={\left(\frac{y_1}{h}\right)}^{{\alpha }_1}{\left(\frac{y_2}{h}\right)}^{{\alpha }_2},{\alpha }_1,{\alpha }_2=0,1,\dots .$$

$$\lambda =\frac{c\tau }{h}$$

$$A({\bf x}_{ij},\tau)\mu(\cdot)=0,\,B({\bf x}_{ij},\tau)\mu(\cdot)=0\,\mbox{ if $\alpha_{1}$ or $\alpha_{2}$ are non-negative odd integers},$$

$$A({\mathbf{x}}_{ij},\tau )\mu (\cdot )=\frac{({\alpha }_1-1)!!({\alpha }_2-1)!!}{({\alpha }_1+{\alpha }_2-1)!!}{\lambda }^{{\alpha }_1+{\alpha }_2},$$

$$B({\mathbf{x}}_{ij},\tau )\mu (\cdot )=\frac{\tau }{{\alpha }_1+{\alpha }_2+1}A({\mathbf{x}}_{ij},\tau )\mu (\cdot )$$

$$\displaystyle\mbox{if $\alpha_{1}$ and $\alpha_{2}$ are both non-negative even integers},$$

$$u_{ij}^0=u({\mathbf{x}}_{ij},0),v_{ij}^0=v({\mathbf{x}}_{ij},0),$$

$$u_{ij}^k=u({\mathbf{x}}_{ij},k\tau ),k=1,2,\dots .$$

$${\mathcal{M}}_{\le 6}^{\ast }=\{1,x_1,x_2,x_1{x}_2,x_1^2,x_2^2\}$$

Full text

1 Introduction

Stencil computations are widely implemented in many numerical algorithms that involve structured grids. In acoustic field simulation based on two- or three-dimensional transient wave equations, the finite difference time domain method is a standard approach leading to stencil operations [1]. In an explicit time-stepping finite difference scheme, a solution value at each point in a time-space grid is calculated using a linear combination of values at its spatial neighbors from previous time steps. Among such schemes the main attention in the literature has been given to two-step schemes (which operate over three time levels $t_{k+1}=(k+1)\tau$, $t_{k}=k\tau$ and $t_{k-1}=(k-1)\tau$ where $\tau$ is a fixed time increment). They have been intensively studied and reviewed in many articles and books (see, e.g., 2-7).

Explicit two-step numerical schemes for the scalar wave equation can also be devised based on spherical means representations as it is done in [8] where an integral evolution formula with three time levels is derived. For the 2D case, the evolution formula has been implemented in [9] and [10] using piecewise polynomial interpolation in 2D mesh cells and numerical integration.

The two-step schemes pose some challenges when imposing the initial conditions. To calculate the value of a sought solution $u$ at the first time level $t_{1}$, one needs values of $u$ from time levels $t_{0}$ and $t_{-1}$. The initial condition for $u$ provides the required values at $t_{0}$. However, $u({\bf x},t_{-1})$ must be inferred from the other initial condition. So, there is a need for transforming the two-step scheme at the first time step to a one-step form which has been discussed, e.g., in [6] and [7]. A conventional approach for the 2D case (described, e.g., in [7]) uses the central difference in time for approximating the initial condition for $v=\partial u/\partial t$. For any point ${\bf x}$ in $R^{2}$, the value of $u({\bf x},t_{-1})$ is inferred as $u({\bf x},t_{1})-2\tau v({\bf x},0)$. Even though this approach is attractive due to its simplicity, it is worthwhile to consider the numerical schemes for both the first time step and the next steps derived using a unified approach.

The present work focuses on building explicit time-stepping stencil computation schemes for the transient 2D acoustic (scalar) wave equation using spherical means formulas including Poisson’s formula [11] and a similar integral formula involving three time levels [8] that are described in section 2. A general form of the implemented time-stepping algorithm is presented in section 3 where two different integral expressions are given for the first time-step and for the next steps. Sections 4 and 5 consider some known results regarding polynomial interpolation on stencils and exact integration needed for deriving explicit stencil computation schemes from these integral expressions. Particular numerical schemes for five-, nine- and 13-point stencils are obtained in section 6. For each scheme, two separate expressions are derived: 1) a one-step expression for the first time-step through the stencil interpolated initial conditions; 2) a two-step expression for the second and next time-steps. All the obtained stencil expressions for the first time-step have not been previously presented in the literature. A derived two-step expression for the nine-point stencil is also new. In contrast, the obtained two-step expressions for five- and 13-point stencils coincide with the corresponding finite-difference stencil expressions. It is shown by simulation that the derived numerical schemes can significantly improve accuracy of stencil computations in comparison with conventional approaches.

2 Representation Formulas

Consider a Cauchy problem for the transient 2D scalar wave equation

[TABLE]

Its solution is given by the representation formula which is often named Poisson’s formula for the 2D wave equation [11]:

[TABLE]

where $\nabla$ denotes the gradient operator and ${\bf y}=(y_{1},y_{2})\in\mathbb{R}^{2}$ is a variable of integration.

Rewriting the previous integral on the unit disk, we get

[TABLE]

where ${\bf z}=(z_{1},z_{2})\in\mathbb{R}^{2}$ is a new variable of integration.

Formulas (3) and (4) are also valid for negative $t$ which can be proved based on the time reversal property of the wave equation. By substituting $-t$ for $t$ in formula (4) and changing ${\bf z}$ to $-{\bf z}$ inside its integrals, we obtain the following analog of that formula for negative time:

[TABLE]

The only difference between the right-hand sides of (4) and (5) is the opposite signs of the second term. So, one can eliminate this term by summing (4) and (5). Shifting the initial moment in the resulting formula from $t=0$ to $t=t_{*}$, we obtain the following expression involving three time points with a time increment $\tau$:

[TABLE]

where the right-hand side does not include the time derivative $v$. The same representation formula (in a different form) has been derived in [8] without using Poisson’s formula.

Both formulas (4) and (6) will be used below to build a time-marching stencil computation algorithm for the wave equation.

3 An Integral Time-Stepping Algorithm

Consider a uniform time grid $\{t_{0}=0,t_{1}=\tau,\ldots,t_{k}=k\tau,\ldots\}$ where $\tau$ is a fixed time-step. Denote by $u_{k}({\bf x})$ the restriction of $u({\bf x},t)$ to a moment $t=t_{k}$. Next, denote by $A({\bf x},\tau)$ and $B({\bf x},\tau)$ the following integral operators acting on continuous functions defined in $\mathbb{R}^{2}$:

[TABLE]

[TABLE]

The time-stepping algorithm proposed here consists of two procedures based on the representation formulas (4) and (6):

1. The procedure for the first time-step which, according to (4), calculates $u_{1}({\bf x})$ as

[TABLE]

2. The procedure for the second and next time-steps which, using (6) for $t_{*}=t_{k}$, calculates $u_{k+1}({\bf x})$ as

[TABLE]

While formula (10) involves three time levels (two time-steps), formula (9) for the first time-step incorporates only two time levels (one time-step) without using any finite difference approximation of the time derivative.

4 Using Polynomial Interpolation on Stencils

Consider a two-dimensional uniform Cartesian grid $\{(i_{1}h,i_{2}h)\}$ where $i_{1}$ and $i_{2}$ are integers, and $h$ is the grid spacing in both directions $x_{1}$ and $x_{2}$. Suppose that the evaluation point ${\bf x}$ in formulas (9) and (10) is a grid point ${\bf x}_{ij}=(ih,jh)$. Our intention is to choose a stencil in the Cartesian grid and reduce the integral formulas (9) and (10) to linear combinations of the stencil node values of $u_{k}({\bf x})\,(k=0,1,2,\ldots)$ and $v_{0}({\bf x})$. Such a reduction will be done by using polynomial interpolation.

Assume that a particular stencil with $m$ nodes is chosen for polynomial interpolation. The corresponding index set $\{(q_{1},q_{2})\}$ is denoted by $Q_{m}$. The stencil index components $q_{1}$ and $q_{2}$ are numbered relative to the referencing point located at the evaluation point ${\bf x}_{ij}$. So, polynomial interpolation in a neighborhood of the evaluation point ${\bf x}_{ij}$ will be carried out using interpolation points

[TABLE]

Following [12] we associate with $Q_{m}$ a set of $m$ distinct bivariate monomials

[TABLE]

where

[TABLE]

The above function has a unique inverse function

[TABLE]

where $[\cdot]$ is the whole part function. So, for each index value $(q_{1},q_{2}))\in Q_{m}$ there exists a unique monomial from $\mathcal{M}_{m}$ and vice versa.

Consider a polynomial space $\mathcal{P}_{m}$ spanned by $\mathcal{M}_{m}$. We will use only those stencils for which the Lagrange interpolation problem is unisolvent in $\mathcal{P}_{m}$ (see [13]). In this case, there exists a Lagrange basis for $\mathcal{P}_{m}$ that can be built as described below.

Suppose that there is an ordering imposed on the monomials in $\mathcal{M}_{m}$

[TABLE]

The corresponding stencil nodes are numbered accordingly using (12)-(14):

[TABLE]

Thus, one can compute the following matrix [13]:

[TABLE]

If this matrix is non-singular, i.e., $\det(D)\neq 0$, which means that the Lagrange basis exists, then the inverse matrix

[TABLE]

can be calculated. Its components are instrumental in expressing the Lagrange basis functions through the chosen set of monomials:

[TABLE]

Once the Lagrange basis is obtained, we again need two indexes to denote the Lagrange basis functions in accordance with the two index notation for grid points. The sequence (15) defines a relationship $s=g(\alpha_{1},\alpha_{2})$ between the ordinal number $s$ and monomial exponents $(\alpha_{1},\alpha_{2})$. Then, the resulting relationship between an index pair $(q_{1},q_{2})$ and the corresponding ordinal number $s$ is given by the following expression:

[TABLE]

where $\alpha(\cdot)$ is specified in (13). Therefore, by introducing a new (two index) notation for the Lagrange basis functions

[TABLE]

we get the interpolation formula for a continuous function $f({\bf x})$ in the form

[TABLE]

Even though different sets of monomials and their sequences can be employed for building Lagrange bases, we will use a particular method of monomial ordering that is described below.

Denote by $\mathcal{M}$ the set of all monomials $x_{1}^{\alpha_{1}}x_{2}^{\alpha_{2}}$ where $\alpha_{1}$ and $\alpha_{2}$ are natural numbers including [math]. For each monomial $x_{1}^{\alpha_{1}}x_{2}^{\alpha_{2}}$, the corresponding ordinal number $s$ will be assigned using the following function:

[TABLE]

It is easy to prove that the function (23) provides a one-to-one correspondence between $\mathcal{M}$ and the set of all positive natural numbers with the usual ordering. The order induced by (23) uses the total degree as the first sorting parameter (similarly to the more common graded lexicographic order) while the difference between individual degrees with the same total degree is used as the next sorting variable.

The set $\mathcal{M}$ endowed with the order induced by (23) will be denoted by $\mathcal{M}^{*}$. The initial segment of $\mathcal{M}^{*}$ with $m$ members is denoted below by $\mathcal{M}_{\leq m}^{*}$. It will be shown in section 6 that such ordered monomial sets play a useful role in building particular numerical schemes.

5 Calculating the Integrals

Let ${\bf x}={\bf x}_{ij}$ in the integral operators (7) and (8). Without loss of generality, assume that the origin of the 2D coordinate system is located at ${\bf x}_{ij}$ which can be achieved by a parallel translation of coordinates. Let $f({\bf y})$ appearing in (7) and (8) be a monomial in local scaled variables $y_{1}/h$ and $y_{2}/h$:

[TABLE]

Then the integrals in (7) and (8) can be exactly calculated and expressed through the Courant number

[TABLE]

Indeed, using the table of integrals of Gradshteyn & Ryzhik [14], expressions (7) and (8) are reduced to the following exact values:

[TABLE]

[TABLE]

where $(\cdot)!!$ is the double factorial. It is assumed that $(-1)!!=1,0!!=1$.

The above formulas allow one to exactly calculate integrals (7) and (8) when $f(\cdot)$ is a polynomial from the Lagrange basis (see section 4).

6 Particular Explicit Two-Step Schemes

Now we can start building some numerical schemes by transforming the procedures (9) and (10) into algebraic expressions. All the functions $u_{0}({\bf x}),v_{0}({\bf x})$ and $u_{k}({\bf x}),k=1,2,...$ included in these procedures will be interpolated in a stencil’s center neighborhood using the same stencil nodes. The following standard notations for grid values of the solution and initial conditions will be used in the computation schemes:

[TABLE]

[TABLE]

6.1 The five-point stencil

Consider building the Lagrange basis for the space of complete second degree polynomials. In this case, the monomial basis sequence ordered according to (23) is as follows:

[TABLE]

The corresponding interpolation stencil node sequence will be written according to (14) as

[TABLE]

The matrix $D$ for this stencil is non-singular with $\det(D)=4h^{8}$. As a result, we get the Lagrange basis as

[TABLE]

One can see that the Lagrange basis (31) includes monomials in the scaled variables $x_{1}/h$ and $x_{2}/h$ with coefficients independent of $h$.

Using formulas from section 5 and notation (21) we get

[TABLE]

[TABLE]

So, all coefficients for node $(-h,-h)$ disappear and the corresponding computational scheme contains only 5 spatial points which is shown in Figure 1.

6.1.1 A new first time-step expression

As a result, the proposed five-point numerical scheme is as follows:

1. for the first time-step

[TABLE]

2. for the second and next time-steps

[TABLE]

A conventional approach for the five-point stencil (see, e.g., [7]) uses the same procedure (35) for the second and next time-steps. However, for the first time-step, this approach uses another formula (rather than (34)) based on the central difference in time for approximating the initial condition for $v$ combined with (35) for $k=0$:

[TABLE]

Comparing (34) and (36) one can see that the first three terms of the right-hand parts of these equations coincide, but the forth term present in the right-hand side of (34) is absent in (36). Thus, the difference between these two first time-step expressions depends on properties of $v_{0}({\bf x})$. If, for example, $v_{0}({\bf x})$ is a linear function of spatial coordinate, then there is no difference between (34) and (36). However, for more general cases encountered in practice, the difference may exist and can be influential from the accuracy point of view which is shown in subsection 6.1.2.

The stability condition for numerical schemes utilizing (35) is well known from von Neumann stability analysis [6]:

[TABLE]

6.1.2 Numerical simulation

To make a numerical comparison of both approaches, we use the initial and boundary conditions corresponding to a standing wave exact solution:

[TABLE]

This solution creates the following pair of initial conditions for the numerical simulation:

[TABLE]

We consider the unit square $\Omega=[0,1]^{2}$ as a space region for the numerical solution and apply boundary condition $u=0$ on $\partial\Omega$ generated by (38). In addition, we assume that $c=1$. Let $n_{t}$ be the number of time-steps and $n\geq 2$ be the space discretization number related to $h$ by equation $nh=1$. The above numerical schemes have been employed over a spatial index set $\{(i,j),0<i,j<n\}$ using grid boundary conditions

[TABLE]

The accuracy of both numerical schemes have been estimated using the relative $L^{2}$ error defined as

[TABLE]

The relative $L^{2}$ errors for the proposed new five-point scheme (34)$\&$(35) (based on Poisson’s formula) and the conventional one (36)$\&$(35) are denoted by $E_{P5}(n,n_{t})$ and $E_{C5}(n,n_{t})$, respectively. Calculated values of $E_{P5}(n,n_{t})$, $E_{C5}(n,n_{t})$ and their ratio for different combinations of $n$, $n_{t}$ and $\lambda=0.707$ are presented in table 1 below. For every value of $n$ used, three values of $n_{t}$ are considered: $n_{t}=1$ (to check errors after the first step), $n_{t}=n$ and $n_{t}=2n$ (to demonstrate the error accumulation process).

Table 1 demonstrated a much higher accuracy of the new scheme in comparison with the conventional one (with the error ratio $E_{C5}/E_{P5}$ exceeding $10^{3}$ for more dense grids). Observed increases of the relative $L^{2}$ errors from $n_{t}=n$ to $n_{t}=2n$ are negligible (less than $10^{-6}$).

6.2 A nine-point square stencil

Consider the initial segment $\mathcal{M}_{\leq 11}^{*}$ of $\mathcal{M}^{*}$ with the last member equal to $x_{1}^{2}x_{2}^{2}$:

[TABLE]

This set includes all the monomials of total degree $\leq 3$ and one monomial of total degree $=4$, and is the minimal initial segment of $\mathcal{M}^{*}$ that includes the bivariate tensor-product of the monomial bases for the second degree polynomials in each spatial coordinate. The corresponding index set $Q_{11}$ is shown in Figure 2 where both solid dots and empty circles denote interpolation points.

The corresponding Lagrange basis can be easily calculated using matrix (17) and presented similar to (31). However, it is preferable to avoid presenting long expressions that include 11 monomials and display only terms that will be used later. Denote by $\phi^{e}_{q_{1},q_{2}}$ a part of $\phi_{q_{1},q_{2}}$ that includes all the monomial terms with even exponents in both coordinates. Then we get

[TABLE]

6.2.1 A new nine-point time-stepping scheme

It follows from the results (26) and (27) of section 5 that functions $\phi^{e}_{q_{1},q_{2}}$ rather than the complete Lagrange basis will be used in building numerical schemes. Therefore, two points (-2,0) and (0.-2) (presented by empty circles in Figure 2) will disappear in the corresponding numerical scheme. The remaining nodes (solid circles in Figure 2) create the nine-point square-shaped computational stencil. Using (26)-(27) we obtain:

[TABLE]

[TABLE]

Now we can use (9) and (10) to build a new nine-point numerical scheme similar to (34)-(35). However, to avoid writing long expressions, some additional notations will be needed:

[TABLE]

Thus, the following time-stepping numerical scheme is derived:

1. for the first time-step

[TABLE]

2. for the second and next time-steps

[TABLE]

Using the von Neumann stability analysis method (see, e.g., [6]), we obtain the stability condition for the above numerical scheme as

[TABLE]

To our knowledge this numerical scheme has not been presented previously in the literature. A conventional explicit nine-point square-shaped scheme (dubbed the isotropic scheme (see, e.g., [15]) has a different form obtained using a nine-point finite difference approximation of the Laplace operator in the two-dimensional space [16]:

[TABLE]

Comparing the coefficients in (47) and (49) it is easy to see that the term $\delta^{k}_{ij}(1,1)$ is less influential in the new nine-point scheme than in the conventional nine-point scheme (with a coefficient ratio equal to $\lambda^{2}$). The stability condition for the conventional scheme is less restrictive than for the new one: $\lambda\leq\lambda_{\max}=\sqrt{3}/{2}\approx 0.866$. On the other hand, the new scheme has some accuracy advantages in comparison with the conventional nine-point approach for $\lambda\leq 0.796$ which is shown in the next subsection.

6.2.2 Numerical simulation

We consider simulation results for the proposed nine-point scheme (46)-(47) using initial and boundary conditions generated by the exact solutions of the previous subsection. Let $c=1$. We assume that the unit square $\Omega=[0,1]^{2}$ is used as a space grid region, the time grid interval is $[0,\lambda]$ with $n_{t}=n$. A comparison is made with simulated results on the same grids for the conventional scheme (49). Even though no expression for the first time-step corresponding to the conventional nine-point scheme is presented in [15], the usual approach based on the central difference for approximating the initial condition for $v$ combined with (49) provides the corresponding expression

[TABLE]

which is used for numerical simulation.

The simulation results are presented in the next table for standing wave simulations. We have used two values of $\lambda$ for the simulation: $\lambda=0.707$, as in subsection 6.1, and $\lambda=0.796$, according to the stability condition (48). Relative $L^{2}$ errors for the schemes (46)$\&$(47) and (50)$\&$(49) are denoted by $E_{P9}$ and $E_{C9}$, respectively. According to the simulated data, the new scheme (46)$\&$(47) appears to be more accurate than the other scheme for both $\lambda=0.707$ and $\lambda=0.796$.

The obtained results show a loss in accuracy of the both nine-point schemes for the simulated problem in comparison with the new five-point scheme results presented in Table 1.

6.3 A 13-point stencil

The next numerical stencil is based on complete bivariate interpolation polynomials of the fourth degree. Consider the initial segment $\mathcal{M}_{\leq 15}^{*}$ of $\mathcal{M}^{*}$:

[TABLE]

This set includes all the monomials of total degree $\leq 4$. The corresponding index set $Q_{15}$ is presented in Figure 3 where both thirteen solid dots and two empty circles denote interpolation points.

After calculating the corresponding Lagrange basis (see section 4) and using (26)-(27), one can determine that all the coefficients in expressions (9)-(10) related to the two empty circles will disappear. As a result, we get a numerical scheme that involves the 13-point stencil (the solid dots in Figure 3). The notations (45) will be used in presenting the scheme to make the expressions more compact:

1. for the first time-step

[TABLE]

2. for the second and next time-steps

[TABLE]

Formula (53) for the second and next time-steps completely coincides with that obtained previously [4, 5] using the finite-difference method.

However, expression (52) for the first time-step has not been presented in the literature before. An advantage of using this expression rather than the conventional one (based on the central difference for approximating the initial condition for $v$) is shown in Table 3 where simulation results for the standing wave (38) solution are presented. The relative $L^{2}$ errors for the conventional and new (Poisson’s) approaches are denoted by $E_{C13}$ and $E_{P13}$, respectively. For simplicity sake, the periodic boundary conditions have been incorporated in the simulation taking into account that the solution (38) is periodic in both spatial directions. Since the maximal Courant number needed for stability of this scheme is $1/\sqrt{2}$, a value of $\lambda=0.707$ was used.

A large advantage in accuracy for the new scheme demonstrated in Table 3 can be attributed to a higher accuracy of the new first time-step expression (52) in comparison with the conventional one. On the other hand, it is worth to notice that the relative $L^{2}$ error $E_{C13}$ (corresponding to the conventional first time-step approach) in Table 3 has almost the same values as the error $E_{C5}$ in Table 1 despite using a higher degree interpolation stencil in Table 3. That is, an error introduced at the first time-step probably suppresses advantages of using a higher degree interpolation at later time-steps.

7 Summary

A new method is implemented to build explicit time-stepping stencil computation schemes for the transient 2D acoustic wave equation. It is based on using Poisson’s formula and a similar three time level expression combined with polynomial stencil interpolation of the solution at each time-step and exact integration. As a result, for each chosen 2D stencil and a set of monomials, a unified time-marching scheme is created that includes two explicit computation procedures: for the first time-step and for the next steps.

Particular explicit stencil computation schemes (with five, nine and 13 space points) are derived. All of the obtained first time-step computation expressions are different from those used in conventional finite-difference methods. The obtained two-step stencil expressions for the five- and 13-point stencils (where the complete interpolation polynomials have been used) coincide with the corresponding finite difference schemes. The obtained two-step expression for the nine-point stencil is new. Its stability region is determined by the von Neumann analysis.

Simulation comparison results are presented for a benchmark problem with an exact solution. It is demonstrated by simulation that the proposed stencil computation approach maintains an accuracy advantage in comparison with conventional finite difference schemes which is mostly attributed to the new first time-step computation expressions.