Complete centered finite difference method for Helmholtz equation

Abstract

A new approach in the finite difference framework is developed, which consists of three steps: choosing the dimension of the local approximation subspace, constructing a vector basis for this subspace, and determining the coefficients of the linear combination. New schemes were developed to form the basis of the local approximation subspace, which were derived by approximating only the $k^{2}u$ term of the Helmholtz equation. The construction of a basis of the local approximation subspace allows the new approach to be able to represent any finite difference scheme that belongs to this subspace. The new method is both consistent and capable of minimizing the dispersion relation for all stencils in all dimensions. In the one-dimensional case and 3-point stencil, pollution error is eliminated. In the two-dimensional (2D) case and 5-point stencil, the Complete Centered Finite Difference Method presents a dispersion relation equivalent to Galerkin/Least-Squares Finite Element Method. In the 2D case and 9-point stencil, two versions were developed using two different bases for the local approximation space. Both versions are equivalent and exhibit a dispersion relation similar to Quasi Stabilized Finite Element Method. Additionally, the dispersion analysis revealed a connection between the coefficients of the linear system and the stencil symmetry.

Key words
dispersion analysis; finite difference method; Helmholtz equation; local truncation error; stabilization

INTRODUCTION

The Helmholtz equation has great applicability in physics and engineering because it describes the time-harmonic acoustic, elastic and electromagnetic waves. Stable and accurate numerical solutions using the finite difference method (FDM) and finite element method (FEM) have been a major challenge to date as reported in a vast literature (Harari & Hughes 1991, Ihlenburg & Babuška 1995a, Singer & Turkel 1998, Alvarez et al. 2006). It is well known that, when the wavenumber of the exact solution $k$ increases, the numerical solutions lose stability and accuracy, even if the ‘rule of thumb’ (Harari & Hughes 1991, Ihlenburg & Babuška 1995a) is satisfied. This heuristic rule prescribes a relationship between $k$ and the mesh parameter $h$ necessary for the numerical solution to adequately approximate a wavelength $kh=\frac{2\pi }{\#p}$, where a reasonable intuitive choice for the number of nodes per wavelength is $7\le \#p\le 10$. However, this rule does not guarantee stability and accuracy for medium and high $k$ values (Ihlenburg & Babuška 1995b, Singer & Turkel 1998, Ihlenburg 1998, Alvarez et al. 2006). In these cases the wavenumber of the numerical solution $k^h$ can be very different from $k$. This misbehavior, known as pollution of the finite element solution or pollution effect (Ihlenburg & Babuška 1995a, Babuška et al. 1995, Babuška & Sauter 1997), is also present in numerical solutions using FDM.

In the 1990s, FEM emerged that eliminated the pollution effect in the one-dimensional (1D) case with pioneering work conducted by Professor Thomas J.R. Hughes (Harari & Hughes 1991, Harari & Hughes 1992). The Galerkin Least-Squares Finite Element Method (GLS) introduces a free parameter $\tau$ that can be chosen so that $k^h-k=0$ (dispersion analysis) in the 1D case (Harari & Hughes 1992). However, when applied to the two-dimensional (2D) case, the GLS method presents a pollution effect similar to the classical Galerkin FEM (Thompson & Pinsky 1995). Currently, the idea is accepted that in 2D and three-dimensional (3D) cases, it is impossible to completely eliminate this pollution effect (Babuška & Sauter 1997). Therefore, the new challenge was to develop methods that minimize this pollution effect with the work led by Professor Ivo M. Babuška (Babuška et al. 1995). The Quasi Stabilized Finite Element Method (QSFEM), as presented in (Babuška et al. 1995), introduces modifications at the algebraic level, and no variational formulation is associated with this method. For that reason, there was no FEM derived from a variational formulation with a dispersion relation equivalent to QSFEM. Later, methods with a dispersion relationship equivalent to QSFEM will emerge (Alvarez et al. 2006, Loula et al. 2007, Rochinha et al. 2007, Carmo et al. 2008).

Within the framework of the FDM, some schemes that improve the quality of the numerical solution have been developed. For example, in Singer & Turkel 1998, fourth-order schemes are developed based on generalizations of Padé’s approximations. In Sutmann 2007, new compact finite difference schemes of sixth order are developed for the 3D case. In Nabavi et al. 2007, a new 9-point sixth-order accurate compact method is developed for the 1D and 2D cases. In Fernandes & Loula 2010, a quasi-optimal finite difference method on non-uniform and unstructured meshes is developed. In Wu (2017) a dispersion minimizing compact finite difference scheme of fourth order is developed for the 2D case. In Wu & Xu (2018), an optimal compact finite difference scheme of sixth order is developed for the 2D case. However, all of these schemes can be seen as isolated or incomplete attempts to improve the numerical solution. This is because none of these schemes consider a completely local approximation subspace.

In the present work, the local approximation subspace is considered in its entirety, since a basis of vector is constructed for this subspace. Thus, it is possible to develop a new formalism within the finite difference framework, which always generates the numerical solution with the dispersion relation with the lowest possible error. The new approach, called the Complete Centered Finite Difference Method, is consistent and presents the highest possible stability (Lax & Richtmyer 1956). That is, for the Helmholtz equation, the new method is convergent and exhibits the lowest possible pollution effect. In this first article, our analysis is performed for the Helmholtz scalar equation in an interior domain with Dirichlet boundary conditions, considering only the 1D and 2D cases on uniform meshes. This choice was made with the aim of presenting the new approach in a simple way. In future work, non-uniform meshes and the 3D case will be considered, in addition to Neumann boundary conditions.

The work is organized as follows. The Helmholtz equation and the classical centered scheme for the 1D and 2D cases are presented first. In next section, the new method for the 1D and 2D cases using different uniform stencils is developed, along with new centered schemes necessary to form the basis of the local approximation subspace. Some of the new schemes were developed by simply making new approximations of the term $k^{2}u$ in the Helmholtz equation. Later, the dispersion analysis of all presented schemes is carried out, along with the selection of free parameters for the new method. In next section, the local truncation error of all analyzed schemes is presented, demonstrating that all schemes are consistent. That is, as the mesh is refined, the numerical solution tends toward the exact solution. Later section presents some numerical experiments, which confirm the good performance of the new method and the theoretical analysis carried out. In the final section, some concluding remarks and suggestions for future work are presented.

Helmholtz equation and centered finite difference method

Consider the Helmholtz scalar equation in an interior domain and with Dirichlet boundary conditions.

$$\Delta u+k^{2}u=f\mathrm{\text{in}}\Omega ,$$(1)

$$u=g\mathrm{\text{on}}\partial \Omega ,$$(2)

where $k$ is the wavenumber, $\Omega$ is the interior of the bounded domain with a Lipschitz boundary $\partial \Omega$, $f$ is a source term, and $u$ is a scalar field describing the time-harmonic acoustic waves.

In uniform meshes or grids and 3-point stencils, for the 1D case, the well-known second order centered finite difference scheme approximates the second derivative by

$$\frac{d^{2}u}{d{x}^2}\approx \frac{U_{i-1}-2U_i+U_{i+1}}{h^2},$$(3)

where $U_{i-1}=U(x_i-h)$, $U_i=U\left(x_i\right)$, $U_{i+1}=U(x_i+h)$ and $h=\Delta x$ (LeVeque 2007). The second term of the Helmholtz equation is approximated as $k^{2}u\left(x\right)\approx k^2{U}_i$. This results in the equation for the node $i$ of the stencil

$$\frac{U_{i-1}-2U_i+U_{i+1}}{h^2}+k^2{U}_i=f_i\mathrm{\text{or}}\frac{1}{h^2}(A{U}_{i-1}+B{U}_i+A{U}_{i+1})=f_i,$$(4)

where, $A=1$ and $B=-2+{\left(kh\right)}^2$. If the generating vector of the stencil equation is defined as $G=\frac{1}{h^2}[A,B,A]$ and the unknown vector as $U=[U_{i-1},U_i,U_{i+1}]$, then equation (4) for every interior point $i$ can be seen as the inner product of these two vectors $<G,U>=f_i$. In the particular case of the homogeneous Helmholtz equation ($f=0$), these vectors are orthogonal $<G,U>=0$. This mathematical formalism for representing a finite difference scheme will facilitate the development of the new method in the next section. It should be noted that the dimension of the vectors $G$ and $U$ coincides with the number of points $L$ of the stencil considered.

Similarly, for the 2D case the second derivatives are approximated by

$$\frac{{\partial }^{2}u}{\partial x^2}\approx \frac{U_{i-1,j}-2U_{i,j}+U_{i+1,j}}{h^2}\mathrm{\text{and}}\frac{{\partial }^{2}u}{\partial y^2}\approx \frac{U_{i,j-1}-2U_{i,j}+U_{i,j+1}}{h^2},$$(5)

and $k^{2}u(x,y)\approx k^2{U}_{i,j}$, which results in the equation for the central node $i,j$ (partial ordering or enumeration)

$$\frac{1}{h^2}(A{U}_{i-1,j}+A{U}_{i,j-1}+B{U}_{i,j}+A{U}_{i+1,j}+A{U}_{i,j+1})=f_{i,j},$$(6)

where $A=1$ and $B=-4+{\left(kh\right)}^2$. In this case, the generating vector of the stencil equation is defined as $G=\frac{1}{h^2}[A,A,B,A,A]$ and the unknown vector as $U=[U_{i-1,j},U_{i,j-1},$ $U_{i,j},U_{i+1,j},U_{i,j+1}]$. This way, the stencil equation (6) for the interior node $i,j$ continues to be the inner product $<G,U>=f_{i,j}$.

Complete centered finite difference method

Consider any generic finite difference scheme with solution $U^h\in S^G$, where $S^G$ denotes the solution space. This space $S^G$ will be called the global approximation space to differentiate it from the local approximation subspace $S^L\subseteq S^G$. The local approximation subspace $S^L$ contains only the components of the vector $U^h$ that are part of the stencil equation, while the global approximation space $S^G$ contains all the components of the vector $U^h$. The dimension of this space is determined by the number of mesh nodes with degrees of freedom, and denoted by $dim\left(S^G\right)=N$. This space contains all vectors with dimension $N$. The interpolant $u^h$, defined as the exact solution of the Helmholtz equation evaluated at the mesh nodes, is an element of this space. Thus, the global error $E^G=U^h-u^h$ can be defined. It is well known that the approximation quality of a finite difference scheme is characterized by three properties: convergence, stability and consistency (Lax & Richtmyer 1956, LeVeque 2007). A finite difference scheme is said to be convergent if $\parallel E^G\parallel \to 0$ when $h\to 0$, where $\parallel \circ \parallel$ denotes some norm in $S^G$. Also, it is well known that the convergence of a scheme is guaranteed if it is stable and consistent.

Let the linear system of algebraic equations of any scheme be denoted by $\mathbf{\text{𝐓}}{U}^h=F$. This finite difference scheme is said to be stable if $\forall h<h_0$ there is the inverse matrix ${\mathbf{\text{𝐓}}}^{-1}$ such that ${\parallel {\mathbf{\text{𝐓}}}^{-1}\parallel }_m<C$, where $C$ is a constant that does not depend on $h$, $h_0>0$ is a specific value of the mesh parameter, and ${\parallel \circ \parallel }_m$ denotes some matrix norm. As mentioned previously, for medium and high $k$ most of the known FDM are not stable for meshes that satisfy the rule of thumb ($k{h}_0\approx 0.6$) (Ihlenburg & Babuška 1995a). Only in meshes that verify $k^{2}h<1$ will there be asymptotic behavior in the convergence of the global error. In the region between $kh<1$ and $k^{2}h<1$, there is pre-asymptotic behavior in the convergence of the global error.

Consider the equation corresponding to an arbitrary interior node $I$ (in global ordering) of the system $\mathbf{\text{𝐓}}{U}^h=F$. This equation is determined by the type of stencil chosen and the finite difference approximations used. Locally at the stencil level, the unknowns of this equation belong to the subspace $S^L\subseteq S^G$ with dimension $dim\left(S^L\right)=L$. That is, by fixing the number of points $L$ on the stencil, the dimension of the local approximation subspace $S^L$ is determined. This space contains all vectors with dimension $L$. Thus, for each central node $I$ of the stencil corresponds to the the component $I$ of the local truncation error vector ${\left[E^L\right]}_I={[\mathbf{\text{𝐓}}u-F]}_I$, where $u$ denotes the exact solution. This finite difference scheme is said to be consistent if the local truncation error $\parallel E^L\parallel \to 0$ when $h\to 0$. On the other hand, a set of $L$ linearly independent vectors forms a basis for this local approximation subspace $S^L$, and any vector in this subspace can be represented as a linear combination of the vectors in this basis. This is the first key idea of the new approach, creating a finite difference method capable of representing any element of the local approximation subspace $S^L$.

In finite difference terms, if the basis vectors of subspace $S^L$ are chosen as the generating vectors $G_l$ of the stencil equation, then this new approach allows representing the equation of any scheme at every interior point $I$ (global ordering or enumeration) as being a linear combination of generating vectors $G_l$ linearly independent, then

$$S_{CCFD}:=⟨\sum _{l=1}^L{\alpha }_l{G}_l,U^h⟩-\left(\sum _{l=1}^L{\alpha }_l\right)f_I=0,$$(7)

where ${\alpha }_l$ are the standard coefficients of the linear combination. These coefficients can be seen as free parameters to be determined considering some criteria for the problem under study, and this is the second key idea of the new approach. In the case of the Helmholtz equation, the criterion used to determine the free parameters will be the minimization of the dispersion relation (Harari & Hughes 1992, Ihlenburg & Babuška 1995b, Babuška et al. 1995).

This new finite difference approach is called the Complete Centered Finite Difference Method (CCFDM) (Nunes 2024). The term ‘complete’ indicates that once a number of stencil points $L$ is fixed, the maximum number of stencil equations necessary to have a complete representation of the local approximation subspace $S^L$ are used, that is, the new approach is capable to generate any element of this subspace $S^L$. The word centered indicates that only centered difference formulas were used in this work, although this approach can also be applied with uncentered difference formulas.

In short, the new formalism is simple and consists of three basic steps. In the first step, a type of stencil with a number of points $L$ is chosen, and with this the dimension of the local approximation subspace $S^L$ is fixed. This subspace determines the number of non-zero diagonals in the matrix $\mathbf{\text{𝐓}}$ of the linear system, and consequently the bandwidth of this matrix. In addition, once $L$ is fixed, the maximum possible order of approximation for the derivatives is also fixed. In the second step, a vector basis is constructed for the local approximation subspace $S^L$. The vectors of this basis are chosen among the generating vectors $G_l$ corresponding to the equation of each scheme for this stencil. The generating vectors $G_l$ are obtained by making different approximations to the terms of the partial differential equation (PDE). Since the dimension of the subspace $S^L$ is $L$, then $L$ linearly independent generating vectors are needed. Just as a reminder, the generating vectors $G_l$ are linearly independent if the system ${\mathbf{\text{𝐆}}}_{L\times L}X=0$ has the unique solution $X=0$, where the columns of the matrix ${\mathbf{\text{𝐆}}}_{L\times L}$ are formed by the generating vectors $G_l$, and the vector $X$ is formed by the coefficients ${\alpha }_l$ of the linear combination. In the third step, the coefficients of the linear combination are determined by requiring or imposing that some specific criterion be checked by the stencil equation. For example, this criterion can be chosen in such a way as to guarantee greater stability and/or accuracy of the method. Once these three steps have been completed, if it is necessary to improve the approximation quality of the method, then it is necessary to increase the dimension of the local approximation subspace $S^L$. In other words, it is necessary to increase the number of points $L$ in the stencil. In the following subsections, the first and second steps of the new formalism will be detailed, and in next section, the third step will be detailed.

CCFDM for 1D case and 3-point stencil

In this case, the local approximation subspace $S^L$ has dimension $dim\left(S^L\right)=L=3$. Consider centered finite difference formulas such that result in three stencil equations with vectors $G_1=[g_{1,1},g_{1,2},g_{1,3}]$, $G_2=[g_{2,1},g_{2,2},g_{2,3}]$ and $G_3=[g_{3,1},g_{3,2},g_{3,3}]$ linearly independent. Thus, the CCFDM equation for the central node $i$ of this stencil is

$$⟨\sum _{l=1}^3{\alpha }_l{G}_l,U^h⟩-\left(\sum _{l=1}^3{\alpha }_l\right)f_i:=A{U}_{i-1}+B{U}_i+C{U}_{i+1}-\left(\sum _{l=1}^3{\alpha }_l\right)f_i=0,$$(8)

where $A={\alpha }_1{g}_{1,1}+{\alpha }_2{g}_{2,1}+{\alpha }_3{g}_{3,1}$, $B={\alpha }_1{g}_{1,2}+{\alpha }_2{g}_{2,2}+{\alpha }_3{g}_{3,2}$ and $C={\alpha }_1{g}_{1,3}+{\alpha }_2{g}_{2,3}+{\alpha }_3{g}_{3,3}$. With this, the first step of formalism is completed. It remains to find three schemes with stencil equations whose generating vectors $G_l$ are linearly independent. These vectors can be found using different approximations for the second derivative and/or for the $k^{2}u$ term of the Helmholtz equation. Here the second derivative is approximated only by the second-order centered finite difference formula (3), and three different approximations are presented for the term $k^{2}u$ of the Helmholtz equation.

The classical centered finite difference scheme (CC) described by the equation (4) generates the first vector $G_1=\frac{1}{h^2}[1,-2+{\left(kh\right)}^2,1]$. A new scheme (NS-1) is generated with the approximation $k^{2}u\approx k^2\frac{[U_{i-1}+U_{i+1}]}{2}$, resulting in the equation for the central node $i$ of the stencil

$$\frac{1}{h^2}[(1+\frac{{\left(kh\right)}^2}{2})U_{i-1}-2U_i+(1+\frac{{\left(kh\right)}^2}{2})U_{i+1}]=f_i.$$(9)

Thus, the generating vector of the equation of this stencil is $G_2=\frac{1}{h^2}[1+\frac{{\left(kh\right)}^2}{2},-2,1+\frac{{\left(kh\right)}^2}{2}]$. This approximation for $k^{2}u$ can be understood as an average of the solution around the central node (Alvarez & Nunes 2024, Nunes 2024). Another new scheme (NS-2) is generated with the approximation $k^{2}u\approx k^2[U_i+\frac{U_{i+1}-U_{i-1}}{2}]$, resulting in the equation

$$\frac{1}{h^2}[(1-\frac{{\left(kh\right)}^2}{2})U_{i-1}+(-2+{\left(kh\right)}^2)U_i+(1+\frac{{\left(kh\right)}^2}{2})U_{i+1}]=f_i,$$(10)

whose generating vector is $G_3=\frac{1}{h^2}[1-\frac{{\left(kh\right)}^2}{2},-2+{\left(kh\right)}^2,1+\frac{{\left(kh\right)}^2}{2}]$. More details and other approximations for $k^{2}u$ can be found in Alvarez & Nunes (2024) and Nunes (2024).

It can be demonstrated that the vectors $G_1$, $G_2$ and $G_3$ are linearly independent (Nunes 2024). With this, the second step of formalism is completed. Therefore, with these vectors $G_1$, $G_2$ and $G_3$, the equation (8) is capable of generating any scheme for this 3-point stencil. Although a complete basis of three vectors was constructed for the subspace $S^{L=3}$, in the case of uniform meshes, the Corollary 1 guarantees that only two vectors will be needed for this basis. This corollary will be seen in the dispersion analysis in the next section. Further on, the dispersion analysis will show that in this local approximation subspace $S^L$ with dimension $L=3$, the CCFDM is capable of eliminating the pollution effect, making it unnecessary to increase the dimension of this subspace $S^L$.

As mentioned in the Introduction, schemes CC (4), NS-1 (9), NS-2 (10) or any other individual scheme are ‘isolated or incomplete’ attempts to improve the numerical solution. This is because it is very difficult for a single scheme with generating vector $G$ to be as close as possible to the ‘best approximation’ in this subspace $S^{L=3}$. However, a basis with three generating vectors $G_1$, $G_2$ and $G_3$ always allows generating the scheme closest to the ‘best approximation’ in this subspace $S^{L=3}$.

CCFDM for 2D case and 5-point stencil

For this case, the local approximation subspace $S^L$ has dimension $dim\left(S^L\right)=5$. Consider centered finite difference formulas such that result in five stencil equations with vectors $G_1$, $G_2$, $G_3$, $G_4$ and $G_5$ linearly independent, where

$$G_l=[g_{l,1},g_{l,2},g_{l,3},g_{l,4},g_{l,5}],l=1,\dots ,5.$$(11)

Thus, the equation for the central node $i,j$ (partial ordering or enumeration) of any 5-point stencil $S_{CCFD}$ can be generated as

$$A_1{U}_{i-1,j}+A_2{U}_{i,j-1}+B{U}_{i,j}+C_1{U}_{i+1,j}+C_2{U}_{i,j+1}=\left(\sum _{l=1}^5{\alpha }_l\right)f_{i,j},$$(12)

where $A_1={\sum }_{l=1}^5{\alpha }_l{g}_{l,1}$, $A_2={\sum }_{l=1}^5{\alpha }_l{g}_{l,2}$, $B={\sum }_{l=1}^5{\alpha }_l{g}_{l,3}$, $C_1={\sum }_{l=1}^5{\alpha }_l{g}_{l,4}$ and $C_2={\sum }_{l=1}^5{\alpha }_l{g}_{l,5}$. With this the first step of formalism is completed. The generating vectors $G_l$ are obtained with five different approximations for the term $k^{2}u$ of the Helmholtz equation and the same approximation (5) for the Laplacian (Nunes 2024).

The first vector $G_1=\frac{1}{h^2}[1,1,-4+{\left(kh\right)}^2,1,1]$ is obtained with the typical centered difference scheme (CC) described by the equation (6). Analogously to the 1D case, a new scheme (NS-1) is generated with the approximation $k^{2}u\approx \frac{k^2}{2}[\frac{U_{i-1,j}+U_{i+1,j}}{2}+\frac{U_{i,j-1}+U_{i,j+1}}{2}]$, resulting in the equation for the central node $i,j$

$$\frac{1}{h^2}[A{U}_{i-1,j}+A{U}_{i,j-1}-4U_{i,j}+A{U}_{i+1,j}+A{U}_{i,j+1}]=f_{i,j},$$(13)

where $A=1+\frac{{\left(kh\right)}^2}{4}$, and the generating vector of the equation of this stencil is $G_2=\frac{1}{h^2}[A,A,-4,A,A]$. Another new scheme (NS-2) is generated with the approximation $k^{2}u\approx k^2[U_{i,j}+\frac{U_{i+1,j}-U_{i-1,j}}{4}+\frac{U_{i,j+1}-U_{i,j-1}}{4}]$, resulting in the equation for the central node $i,j$

$$\frac{1}{h^2}[A{U}_{i-1,j}+A{U}_{i,j-1}+B{U}_{i,j}+C{U}_{i+1,j}+C{U}_{i,j+1}]=f_{i,j},$$(14)

where $A=1-\frac{{\left(kh\right)}^2}{4}$, $B=-4+{\left(kh\right)}^2$ and $C=1+\frac{{\left(kh\right)}^2}{4}$. The generating vector of the equation of this stencil is $G_3=\frac{1}{h^2}[A,A,B,C,C]$. The third new scheme (NS-3) is generated with the approximation $k^{2}u\approx k^2[U_{i,j}+\frac{U_{i-1,j}-U_{i+1,j}}{4}+\frac{U_{i,j+1}-U_{i,j-1}}{4}]$, resulting in the equation for the central node $i,j$

$$\frac{1}{h^2}[A{U}_{i-1,j}+C{U}_{i,j-1}+B{U}_{i,j}+C{U}_{i+1,j}+A{U}_{i,j+1}]=f_{i,j},$$(15)

where $A=1+\frac{{\left(kh\right)}^2}{4}$, $B=-4+{\left(kh\right)}^2$ and $C=1-\frac{{\left(kh\right)}^2}{4}$. The generating vector of the equation of this stencil is $G_4=\frac{1}{h^2}[A,C,B,C,A]$. The last new scheme (NS-4) needed to complete a basis for the subspace $S^L$ is generated with the approximation $k^{2}u\approx k^2[U_{i,j}+\frac{U_{i+1,j}+U_{i-1,j}}{2}-\frac{U_{i,j+1}+U_{i,j-1}}{2}]$, resulting in the equation for the central node $i,j$

$$\frac{1}{h^2}[A{U}_{i-1,j}+C{U}_{i,j-1}+B{U}_{i,j}+A{U}_{i+1,j}+C{U}_{i,j+1}]=f_{i,j},$$(16)

where $A=1+\frac{{\left(kh\right)}^2}{2}$, $B=-4+{\left(kh\right)}^2$ and $C=1-\frac{{\left(kh\right)}^2}{2}$. The generating vector of the equation of this stencil is $G_5=\frac{1}{h^2}[A,C,B,A,C]$. It can be demonstrated that the vectors $G_1$, $G_2$, $G_3$, $G_4$ and $G_5$ are linearly independent (Nunes 2024). With this the second step of formalism is completed.

Although a complete basis of five vectors was constructed for the subspace $S^{L=5}$, in the case of uniform meshes, the Corollary 1 guarantees that only three vectors will be needed for this basis. Further on, the dispersion analysis will show that in this local approximation subspace $S^L$ with dimension $L=5$, the CCFDM is capable of minimizing pollution effect. However, this minimization presents a dispersion relation equivalent to the GLS method (Thompson & Pinsky 1995), making it necessary to increase the dimension of this subspace $S^L$ to obtain a dispersion relation equivalent to the QSFEM (Babuška et al. 1995).

Again, schemes CC (6), NS-1 (52), NS-2 (53), NS-3 (54), NS-4 (55) or any other individual scheme are ‘isolated or incomplete’ attempts to improve the numerical solution. This is because it is very difficult for a single scheme to be as close as possible to the ‘best approximation’ in this subspace $S^{L=5}$. However, a basis with five generating vectors $G_1$, $G_2$, $G_3$, $G_4$ and $G_5$ always allows generating the scheme closest to the ‘best approximation’ in this subspace $S^{L=5}$, and this is done with the optimal choice of linear combination coefficients via dispersion analysis. These statements remain valid if the dimension of the subspace $S^L$ is increased, but as the dimension of this subspace increases it will be more difficult for any individual scheme to be as close as possible to the ‘best approximation’. Furthermore, this is considered impossible for the Helmholtz equation in two and three dimensions, as it is impossible to eliminate the pollution effect in these cases.

CCFDM for 2D case and 9-point stencil

In this case, the local approximation subspace $S^L$ has dimension $dim\left(S^L\right)=9$. Consider centered finite difference formulas such that result in nine stencil equations with vectors $G_1$, $G_2$, $G_3$, $G_4$, $G_5$, $G_6$, $G_7$, $G_8$ and $G_9$ linearly independent. Thus, the equation for the central node $i,j$ (partial ordering or enumeration) of any 9-point stencil $S_{CCFD}$ can be generated as

$$\begin{array}{c}\hfill A_1{U}_{i-1,j}+A_2{U}_{i,j-1}+A_3{U}_{i-1,j-1}+A_4{U}_{i+1,j-1}+B{U}_{i,j}\hfill \\ \hfill +C_1{U}_{i+1,j}+C_2{U}_{i,j+1}+C_3{U}_{i-1,j+1}+C_4{U}_{i+1,j+1}=\left(\sum _{l=1}^9{\alpha }_l\right)f_{i,j}.\hfill \end{array}$$(17)

And so the first step of the new formalism is completed.

The compact 9-point stencil in Figure 1 (red color) can also be decomposed into two 5-point stencils (black and blue) (Vichnevetsky & Bowles 1982). In this stencil, the Laplacian will be approximated in two different ways. The first approximation will use centered formulas of second order (approximate Laplacian with 5-point), and the second approximation will use centered fourth-order formulas (approximate Laplacian with 9 points). In this way, two versions of the CCFDM will be developed. That is, two different bases for the local approximation subspace $S^L$ will be constructed, allowing the development of two versions of the CCFDM. However, it will later be proven that these two versions are equivalent, as they exhibit the same dispersion relation and the same global error behavior. In other words, both versions present a minimum dispersion relation equivalent to QSFEM. In addition, although nine vectors are needed to construct a complete basis for the subspace $S^L$, in the case of uniform meshes, the Corollary 1 guarantees that only five vectors will be needed for this basis. Therefore, in the following development only five vectors will be presented for this basis.

Figure 1
Compact 9-point stencil (red) and two 5-point stencils (black and blue).

CCFDM with 9-point stencil and Laplacian with 5-point

The first version will be denoted by CCFDM-9p-L5p, and consists of decomposing the compact 9-point stencil into two 5-point stencils (Vichnevetsky & Bowles 1982). In Figure 1, these two 5-point stencils are represented as a classic centered 5-point stencil (black color) and another rotated 5-point stencil (blue color). Thus, the Laplacian is approximated by (5) in the 5-point stencil in black, and by

$$\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\approx \frac{U_{i-1,j-1}+U_{i+1,j-1}-4U_{i,j}+U_{i-1,j+1}+U_{i+1,j+1}}{2h^2}$$(18)

in the 5-point stencil in blue (Vichnevetsky & Bowles 1982).

The first three generating vectors are analogous to the classical centered 5-point stencil (black), but now the dimension of these vectors is 9.

$$G_1=\frac{1}{h^2}[1,1,0,0,-4+(kh{)}^2,1,1,0,0],$$(19)

$$G_2=\frac{1}{h^2}[1+\frac{(kh{)}^2}{4},1+\frac{(kh{)}^2}{4},0,0,-4,1+\frac{(kh{)}^2}{4},1+\frac{(kh{)}^2}{4},0,0],$$(20)

$$\begin{array}{cc}\hfill G_3& =\frac{1}{h^2}[1+\frac{{\left(kh\right)}^2}{4},1-\frac{{\left(kh\right)}^2}{4},0,0,-4+{\left(kh\right)}^2,1+\frac{{\left(kh\right)}^2}{4},1-\frac{{\left(kh\right)}^2}{4},0,0].\hfill \end{array}$$(21)

To obtain the fourth and fifth generating vector, the rotated 5-point stencil (blue) is used. The fourth generating vector is obtained by the approximation $k^{2}u\approx \frac{k^2}{2}[\frac{U_{i-1,j+1}+U_{i+1,j-1}}{2}+\frac{U_{i-1,j-1}+U_{i+1,j+1}}{2}]$, where

$$G_4=\frac{1}{h^2}[0,0,\frac{1}{2}+\frac{{\left(kh\right)}^2}{4},\frac{1}{2}+\frac{{\left(kh\right)}^2}{4},-2,0,0,\frac{1}{2}+\frac{{\left(kh\right)}^2}{4},\frac{1}{2}+\frac{{\left(kh\right)}^2}{4}].$$

The fifth generating vector is obtained by the approximation $k^{2}u\approx k^2[U_{i,j}+\frac{U_{i-1,j+1}+U_{i+1,j-1}}{2}-\frac{U_{i-1,j-1}+U_{i+1,j+1}}{2}]$, where

$$G_5=\frac{1}{h^2}[0,0,\frac{1}{2}-\frac{{\left(kh\right)}^2}{2},\frac{1}{2}+\frac{{\left(kh\right)}^2}{2},-2+{\left(kh\right)}^2,0,0,\frac{1}{2}+\frac{{\left(kh\right)}^2}{2},\frac{1}{2}-\frac{{\left(kh\right)}^2}{2}].$$(23)

It can be demonstrated that the vectors $G_1$, $G_2$, $G_3$, $G_4$ and $G_5$ are linearly independent (Nunes 2024). Four more generating vectors ($G_6$, $G_7$, $G_8$ and $G_9$) are missing to complete a basis for the subspace $S^L$. However, as previously stated, the Corrolary 1 guarantees that they will not be necessary in the case of uniform meshes. With this, the second step of the new formalism is completed for this choice of the vector basis for the local approximation subspace.

CCFDM with 9-point stencil and Laplacian with 9-point

The second version will be denoted by CCFDM-9p-L9p, and consists of using a single approximation for the Laplacian operator as in (LeVeque 2007) (compact 9-point stencil red color in Figure 1).

$$\begin{array}{c}\hfill \frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\approx \frac{1}{6h^2}(4U_{i-1,j}+4U_{i+1,j}+4U_{i,j-1}+4U_{i,j+1}\hfill \\ \hfill +U_{i-1,j-1}+U_{i-1,j+1}+U_{i+1,j-1}+U_{i+1,j+1}-20U_{i,j}).\hfill \end{array}$$(24)

The approximations of the term $k^{2}u$ are the same as those used in the previous subsection “CCFDM with 9-point stencil and Laplacian with 5-point“. Thus, the generating vectors $G_1$, $G_2$, $G_3$, $G_4$ and $G_5$ are obtained.

$\begin{array}{cc}\hfill G_1& =\frac{1}{h^2}[\frac{2}{3},\frac{2}{3},\frac{1}{6},\frac{1}{6},-\frac{10}{3}+{\left(kh\right)}^2,\frac{2}{3},\frac{2}{3},\frac{1}{6},\frac{1}{6}],\hfill \end{array}$(25)

$\begin{array}{cc}\hfill G_2& =\frac{1}{h^2}[\frac{2}{3}+\frac{{\left(kh\right)}^2}{4},\frac{2}{3}+\frac{{\left(kh\right)}^2}{4},\frac{1}{6},\frac{1}{6},-\frac{10}{3},\frac{2}{3}+\frac{{\left(kh\right)}^2}{4},\frac{2}{3}+\frac{{\left(kh\right)}^2}{4},\frac{1}{6},\frac{1}{6}],\hfill \end{array}$(26)

$\begin{array}{cc}\hfill G_3& =\frac{1}{h^2}[\frac{2}{3}+\frac{{\left(kh\right)}^2}{4},\frac{2}{3}-\frac{{\left(kh\right)}^2}{4},\frac{1}{6},\frac{1}{6},-\frac{10}{3}+{\left(kh\right)}^2,\frac{2}{3}+\frac{{\left(kh\right)}^2}{4},\frac{2}{3}-\frac{{\left(kh\right)}^2}{4},\frac{1}{6},\frac{1}{6}],\hfill \end{array}$(38)

$\begin{array}{cc}\hfill G_4& =\frac{1}{h^2}[\frac{2}{3},\frac{2}{3},\frac{1}{6}+\frac{{\left(kh\right)}^2}{4},\frac{1}{6}+\frac{{\left(kh\right)}^2}{4},-\frac{10}{3},\frac{2}{3},\frac{2}{3},\frac{1}{6}+\frac{{\left(kh\right)}^2}{4},\frac{1}{6}+\frac{{\left(kh\right)}^2}{4}],\hfill \end{array}$(28)

$$\begin{array}{cc}\hfill G_5& =\frac{1}{h^2}[\frac{2}{3},\frac{2}{3},\frac{1}{6}-\frac{{\left(kh\right)}^2}{2},\frac{1}{6}+\frac{{\left(kh\right)}^2}{2},-\frac{10}{3}+{\left(kh\right)}^2,\frac{2}{3},\frac{2}{3},\frac{1}{6}+\frac{{\left(kh\right)}^2}{2},\frac{1}{6}-\frac{{\left(kh\right)}^2}{2}].\hfill \end{array}$$(29)

It can be demonstrated that the vectors $G_1$, $G_2$, $G_3$, $G_4$ and $G_5$ are linearly independent (Nunes 2024). Analogous to the CCFDM-9p-L5p version, the generating vectors $G_6$, $G_7$, $G_8$ and $G_9$ that complete a basis for the subspace $S^L$ are unnecessary for the case of uniform meshes. With this the second step of the new formalism is completed for the other choice of vector basis.

Dispersion analysis and optimal choise of the free parameters

The dispersion analysis allows estimating the difference between the wavenumber of the exact solution $k$ and the wavenumber of the numerical solution $k^h$. This analysis consists of considering the homogeneous Helmholtz equation ($f=0$) with exact solutions in the form of a plane wave $u(x,y)=e^{\mathbf{\text{𝐢}}\overrightarrow{k}\cdot \overrightarrow{r}}=e^{\mathbf{\text{𝐢}}k(xcos\left(\theta \right)+ysin\left(\theta \right))}$ propagating in the $\theta$-direction, where $\mathbf{\text{𝐢}}=\sqrt{-1}$. It is assumed that the numerical solution should also be a plane wave $U(x,y)=e^{\mathbf{\text{𝐢}}{k}^h(xcos\left(\theta \right)+ysin\left(\theta \right))}$ with discrete wavenumber $k^h$. Thus, substituting this numerical solution into the stencil equation of each scheme, $k^h$ can be obtained as a function of $k$ (Harari & Hughes 1992, Ihlenburg & Babuška 1995b, Babuška et al. 1995, Ihlenburg 1998, Alvarez et al. 2006, Loula et al. 2007, Rochinha et al. 2007, Carmo et al. 2008, Fernandes & Loula 2010).

1D case with 3-point stencil

In this case $\theta =0$ and the numerical solution at point $x_i$ is $U\left(x_i\right):=U_i=e^{\mathbf{\text{𝐢}}{k}^h{x}_i}$. Considering 3-point stencil on uniform meshes, and substituting this numerical solution in the equation (8) the following dispersion relation is obtained,

$$[cos\left(k^{h}h\right)(A+C)+B]+\mathbf{\text{𝐢}}sin\left(k^{h}h\right)[-A+C]=0+\mathbf{\text{𝐢}}0.$$(30)

Therefore, the following two equations are obtained:

$$\mathbf{\text{𝐑𝐞𝐚𝐥 𝐩𝐚𝐫𝐭:}}cos\left(k^{h}h\right)[A+C]+B=0,$$(31)

$$\mathbf{\text{𝐈𝐦𝐚𝐠𝐢𝐧𝐚𝐫𝐲 𝐩𝐚𝐫𝐭:}}sin\left(k^{h}h\right)[C-A]=0.$$(32)

As far as the authors know, this is the first time that two equations appear to be satisfied in the dispersion analysis. In all the works known to us so far, only the equation (20) of the real part has always appeared (Harari & Hughes 1992, Ihlenburg & Babuška 1995b, Babuška et al. 1995, Ihlenburg 1998, Alvarez et al. 2006, Loula et al. 2007, Rochinha et al. 2007, Carmo et al. 2008, Fernandes & Loula 2010). In these works the imaginary part (21) was satisfied automatically, but no mention of this restriction was found. Thus, the equation (21) of the imaginary part is understood as an additional restriction that must satisfy the method to verify the dispersion analysis. In the 1D case, uniform meshes and 3-point stencil, it is necessary that $A=C$. The CC and NS-1 schemes satisfy the restriction (21), but the NS-2 scheme does not satisfy this restriction (21).

In a similar way to all works known to us, the real part equation (20) is used to explain $k^h$ as a function of $k$,

$$k^h=\frac{1}{h}arccos(-\frac{B}{A+C}),$$(33)

where the values of $A$, $B$ and $C$ of each scheme must be replaced. Performing the Taylor series expansion of the function $arccos(\circ )$ around $kh=0$ (Ihlenburg 1998, Fernandes 2009, Nunes 2024) it is possible to obtain a polynomial representation of $k^h-k$ for each scheme

$$\mathbf{\text{𝐂𝐂:}}{k}^h-k=\frac{1}{24}{k}^3{h}^2+\frac{3}{640}{k}^5{h}^4+𝒪\left(k^6{h}^5\right),$$(34)

$$\mathbf{\text{𝐍𝐒-𝟏:}}{k}^h-k=-\frac{5}{24}{k}^3{h}^2+\frac{43}{640}{k}^5{h}^4+𝒪\left(k^6{h}^5\right),$$(35)

$$\mathbf{\text{𝐍𝐒-𝟐:}}{k}^h-k=\frac{1}{24}{k}^3{h}^2+\frac{3}{640}{k}^5{h}^4+𝒪\left(k^6{h}^5\right).$$(36)

For the CCFDM, $A=\frac{1}{h^2}[{\alpha }_1+{\alpha }_2+{\alpha }_3+\frac{{\left(kh\right)}^2}{2}({\alpha }_2-{\alpha }_3)]$, $B=\frac{1}{h^2}[-2({\alpha }_1+{\alpha }_2+{\alpha }_3)+{\left(kh\right)}^2({\alpha }_1+{\alpha }_3)]$ and $C=\frac{1}{h^2}[{\alpha }_1+{\alpha }_2+{\alpha }_3+\frac{{\left(kh\right)}^2}{2}({\alpha }_2+{\alpha }_3)]$. The coefficients ${\alpha }_l$ can be determined to minimize $|k-k^h|$. Thus, it is possible to construct stencils with minimal error in the dispersion relation. In particular, for the 1D case it is possible to reach the absolute minimum ($k^h-k=0$) by making $k^h=k$ in the equation (22), and in this way find the parameters ${\alpha }_l$ that satisfy this dispersion relation.

$${\alpha }_2=(1+{\alpha }_3)\frac{2-2cos\left(kh\right)-k^2{h}^2}{cos\left(kh\right)[2+k^2{h}^2]-2},$$(37)

where ${\alpha }_1=1$ since the coefficients ${\alpha }_l$ are linearly dependent through the dispersion relation, and therefore it is possible to fix one of them (Ihlenburg 1998).

The generating vector $G_3$ does not satisfy the imaginary part of the dispersion relation (21), and must be disregarded by choosing ${\alpha }_3=0$. Thus, of the three generating vectors necessary to complete a basis in subspace $S^L$, the generating vector $G_3$ is unnecessary in the case of uniform meshes. In this way, the coefficients in the equation (8) are determined as $A=C=\frac{1}{h^2}[1+{\alpha }_2(1+\frac{k^2{h}^2}{2})]=\frac{k^4{h}^2}{4-[4+2k^2{h}^2]cos\left(kh\right)}$ and $B=\frac{1}{h^2}[k^2{h}^2-2(1+{\alpha }_2)]=\frac{k^4{h}^{2}cos\left(kh\right)}{-2+[2+k^2{h}^2]cos\left(kh\right)}=-\frac{1}{2}Acos\left(kh\right)$. With this the third step of formalism is completed. This CCFDM eliminates pollution effect ($k^h-k=0$) and there is no need to increase the dimension of the local approximation subspace $S^L$, since its numerical solution coincides with the interpolant. Figure 2 presents $k^h-k$ for the four methods analyzed.

Figure 2
${𝐤}^{𝐡}\mathbf{-}𝐤$ as a function of $𝐤𝐡$ for $𝐤\mathbf{=}\mathrm{𝟏𝟎𝟎}$.

2D case with 5-point stencil

In this case the numerical solution at point $(x_i,y_j)$ is $U_{i,j}=e^{\mathbf{\text{𝐢}}{k}^h(x_{i}cos\left(\theta \right)+y_{j}sin\left(\theta \right))}$. Considering 5-point stencil on uniform meshes, and substituting this numerical solution in the equation (12) the following two equations are obtained:

$$[A_1+C_1]cos(k^{h}hcos\left(\theta \right))+[A_2+C_2]cos(k^{h}hsin\left(\theta \right))+B=0,$$(38)

$$[C_1-A_1]sin(k^{h}hcos\left(\theta \right))+[C_2-A_2]sin(k^{h}hsin\left(\theta \right))=0.$$(39)

Again as in the 1D case, the imaginary part equation (28) introduces a new constraint that can be satisfied by requiring $A_1=C_1$ and $A_2=C_2$. The CC, NS-1 and NS-4 schemes satisfy the restriction (39), but the NS-2 and NS-3 schemes do not satisfy this restriction (39).

The real part (38) is used to explain $k^h$ as a function of $k$. Replacing the values of $A_1$, $A_2$, $B$, $C_1$ and $C_2$ of each scheme in (38) and performing the Taylor series expansion around $kh=0$, after some algebraic manipulations using the software Mathematica¹ (Wolfram 2015), the expressions of $k^h-k$ for each scheme are obtained:

$$\mathbf{\text{𝐂𝐂:}}{k}^h-k=\frac{3+cos\left(4\theta \right)}{96}{k}^3{h}^2+\frac{235+162cos\left(4\theta \right)+35{cos}^2\left(4\theta \right)}{92160}{k}^5{h}^4+𝒪\left(k^7{h}^6\right),$$(40)

$$\mathbf{\text{𝐍𝐒-𝟏:}}{k}^h-k=\frac{-9+cos\left(4\theta \right)}{96}{k}^3{h}^2+\frac{1315-198cos\left(4\theta \right)+35{cos}^2\left(4\theta \right)}{92160}{k}^5{h}^4+𝒪\left(k^7{h}^6\right),$$(41)

$$\mathbf{\text{𝐍𝐒-𝟐:}}{k}^h-k=\frac{3+cos\left(4\theta \right)}{96}{k}^3{h}^2+\frac{235+162cos\left(4\theta \right)+35{cos}^2\left(4\theta \right)}{92160}{k}^5{h}^4+𝒪\left(k^7{h}^6\right),$$(42)

$$\mathbf{\text{𝐍𝐒-𝟑:}}{k}^h-k=\frac{3+cos\left(4\theta \right)}{96}{k}^3{h}^2+\frac{235+162cos\left(4\theta \right)+35{cos}^2\left(4\theta \right)}{92160}{k}^5{h}^4+𝒪\left(k^7{h}^6\right),$$(43)

$$\begin{array}{c}\hfill \mathbf{\text{𝐍𝐒-𝟒:}}{k}^h-k=\frac{3+cos\left(4\theta \right)}{96}{k}^3{h}^2+(\frac{235-1680cos\left(2\theta \right)+8640{cos}^2\left(2\theta \right)}{92160}\hfill \\ \hfill +\frac{162cos\left(4\theta \right)-1200cos\left(2\theta \right)cos\left(4\theta \right)+35{cos}^2\left(4\theta \right)}{92160})k^5{h}^4+𝒪\left(k^7{h}^6\right).\hfill \end{array}$$(44)

For the CCFDM-5p, the dispersion relations (38) and (39) are determined by

$$A_1=\frac{1}{h^2}[{\alpha }_1+(1+\frac{{\left(kh\right)}^2}{4})({\alpha }_2+{\alpha }_4+{\alpha }_5)+(1-\frac{{\left(kh\right)}^2}{4}){\alpha }_3],$$(45)

$$A_2=\frac{1}{h^2}[{\alpha }_1+(1-\frac{{\left(kh\right)}^2}{4})({\alpha }_3+{\alpha }_4+{\alpha }_5)+(1+\frac{{\left(kh\right)}^2}{4}){\alpha }_2],$$(46)

$$B=\frac{1}{h^2}[-4{\alpha }_2+(-4+{\left(kh\right)}^2)({\alpha }_1+{\alpha }_3+{\alpha }_4+{\alpha }_5)],$$(47)

$$C_1=\frac{1}{h^2}[{\alpha }_1+(1+\frac{{\left(kh\right)}^2}{4})({\alpha }_2+{\alpha }_3+{\alpha }_5)+(1-\frac{{\left(kh\right)}^2}{4}){\alpha }_4],$$(48)

$$C_2=\frac{1}{h^2}[{\alpha }_1+(1+\frac{{\left(kh\right)}^2}{4})({\alpha }_2+{\alpha }_3+{\alpha }_4)+(1-\frac{{\left(kh\right)}^2}{4}){\alpha }_5].$$(49)

The coefficients ${\alpha }_l$ are linearly dependent through the dispersion relation, therefore it is possible to fix ${\alpha }_1=1$ (Ihlenburg 1998). The generating vectors $G_3$ and $G_4$ do not satisfy the imaginary part of the dispersion relation (28), and are disregarded making ${\alpha }_3=0$ and ${\alpha }_4=0$. Thus, of the five generating vectors needed to complete a basis in subspace $S^L$, only three are necessary for the case of uniform meshes. These coefficients ${\alpha }_2$ and ${\alpha }_5$ can be determined to minimize $|k-k^h|$. In the 2D case, the absolute minimum $k^h=k$ of the equation (38) cannot be reached for any $\theta$-direction of the plane wave. This absolute minimum will be reached the same number of times as free parameters ${\alpha }_l$ available. Thus, choosing the directions ${\theta }_1$ and ${\theta }_2$ for the plane wave it is possible to obtain

$$\begin{array}{cc}\hfill {\alpha }_2& =\frac{-1}{{\chi }^{5p}}[8(c_2-c_1+s_1-s_2+c_1{s}_2-c_2{s}_1)+2k^2{h}^2(c_1-c_2+s_2-s_1)],\hfill \\ \hfill {\alpha }_5& =\frac{1}{{\chi }^{5p}}\left[k^2{h}^2(c_2-c_1-s_1+s_2)\right],\hfill \end{array}$$(50)

$${\chi }^{5p}=8(c_2-c_1+s_1-s_2+c_1{s}_2-c_2{s}_1)+k^2{h}^2(c_2-c_1+s_2-s_1+2c_1{s}_2-2c_2{s}_1),$$(51)

where $c_1=cos(khcos\left({\theta }_1\right))$, $c_2=cos(khcos\left({\theta }_2\right))$, $s_1=cos(khsin\left({\theta }_1\right))$ and $s_2=cos(khsin\left({\theta }_2\right))$. Therefore, the CCFDM-5p solution matches the interpolant only when $\theta ={\theta }_1$ or $\theta ={\theta }_2$, and

$$\mathbf{\text{𝐂𝐂𝐅𝐃𝐌-𝟓𝐩:}}{k}^h-k=\frac{{\gamma }_1^{5p}{k}^3{h}^2+{\gamma }_2^{5p}{k}^5{h}^4}{82575360{(1+{\alpha }_2+{\alpha }_5)}^2}+𝒪\left(k^7{h}^6\right),$$(53)

$$\begin{array}{cc}\hfill {\gamma }_1^{5p}& =(1+{\alpha }_2+{\alpha }_5)[-2580480(-1+3{\alpha }_2-{\alpha }_5)-20643840{\alpha }_{5}cos\left(2\theta \right)\hfill \\ \hfill & +860160(1+{\alpha }_2+{\alpha }_5)cos\left(4\theta \right)],\hfill \end{array}$$(54)

$$\begin{array}{cc}\hfill {\gamma }_2^{5p}& =2240[101+533{\alpha }_2^2-230{\alpha }_2(1+{\alpha }_5)+{\alpha }_5(202+1829{\alpha }_5)]\hfill \\ \hfill & +107520{\alpha }_5[53{\alpha }_2-19(1+{\alpha }_5)]cos\left(2\theta \right)\hfill \\ \hfill & -16128[-9+11{\alpha }_2^2+2{\alpha }_2(1+{\alpha }_5)-3{\alpha }_5(6+83{\alpha }_5)]cos\left(4\theta \right)\hfill \\ \hfill & -537600{\alpha }_5(1+{\alpha }_2+{\alpha }_5)cos\left(6\theta \right)+15680{(1+{\alpha }_2+{\alpha }_5)}^{2}cos\left(8\theta \right).\hfill \end{array}$$(55)

The dispersion relation of the CCFDM-5p is simplified by replacing the parameters ${\alpha }_2$ equation (50) and ${\alpha }_5$ equation (51) in the equation (53),

$$k^h-k=\frac{1}{48}{k}^3{h}^2[cos\left(2\theta \right)-cos\left(2{\theta }_1\right)][cos\left(2\theta \right)-cos\left(2{\theta }_2\right)]+𝒪\left(k^5{h}^4\right),$$(56)

which is equivalent to the dispersion of the GLS method for ${\theta }_1=\frac{\pi }{8}$ and ${\theta }_2={\theta }_1+\frac{\pi }{4}$ (Thompson & Pinsky 1995).

It is possible to obtain a dispersion equivalent to the equation (53) or equation (56) in another simpler way. For this, it is necessary to require more symmetry in the stencil, which is justified because the Hemlholtz equation in a homogeneous medium has a high degree of symmetry. Restricting the imaginary part of the dispersion relation (28) guarantees stencil symmetry with respect to the $X$ and $Y$ axes. The symmetry of the stencil with respect to the straight lines $y=x$ and $y=-x$ is guaranteed if $A_1=A_2$. The generating vector $G_5$ does not have this symmetry, and can be disregarded by making ${\alpha }_5=0$. Thus, choosing just ${\alpha }_2$ as

$${\alpha }_2=-\frac{2(-4+k^2{h}^2+2c_1+2s_1)}{-8+4c_1+k^2{h}^2{c}_1+4s_1+k^2{h}^2{s}_1},$$(57)

thus, it is possible to obtain

$$\mathbf{\text{𝐂𝐂𝐅𝐃𝐌-𝟓𝐩:}}{k}^h-k=\frac{{\gamma }_1^{5ps}{k}^3{h}^2}{96(1+{\alpha }_2)}+\frac{{\gamma }_2^{5ps}{k}^5{h}^4}{184320{(1+{\alpha }_2)}^2}+𝒪\left(k^7{h}^6\right),$$(58)

$${\gamma }_1^{5ps}=3+cos\left(4\theta \right)+{\alpha }_2[-9+cos\left(4\theta \right)],$$(59)

$$\begin{array}{c}\hfill {\gamma }_2^{5ps}=505+324cos\left(4\theta \right)+35cos\left(8\theta \right)-{\alpha }_2[1150+72cos\left(4\theta \right)-70cos\left(8\theta \right)]\hfill \\ \hfill +{\alpha }_2^2[2665-396cos\left(4\theta \right)+35cos\left(8\theta \right)].\hfill \end{array}$$(60)

This dispersion relation is simplified by replacing the parameters ${\alpha }_2$ (57) in the equation (58),

$$k^h-k=\frac{1}{96}{k}^3{h}^2[cos\left(4\theta \right)-cos\left(4{\theta }_1\right)]+𝒪\left(k^5{h}^4\right),$$(61)

which is equivalent to the dispersion of the GLS method for ${\theta }_1=\frac{\pi }{8}$ (Thompson & Pinsky 1995).

With this the third step of formalism is completed for the 5-point stencil on uniform meshes. It should be highlighted that the equations (53) or (56) and equations (58) or (61) correspond to the same dispersion relation, although they correspond to two versions of the CCFDM-5p. These two versions differ in the degree of symmetry of the stencil. Figures 3 and 4 show the dispersion relation for the analyzed schemes, where ${\theta }_1=\frac{\pi }{8}$ and ${\theta }_2=\frac{3\pi }{8}$. These are the values of ${\theta }_1$ and ${\theta }_2$ that generate the smallest error in the dispersion relation for the CCFDM-5p. With this choice, CCFDM-5p is a scheme with a dispersion relation equivalent to that of the GLS method (Thompson & Pinsky 1995), but with the following difference. The matrix $\mathbf{\text{𝐓}}$ of the linear system of the GLS method has nine non-zero diagonals, while the matrix of the CCFDM has only five non-zero diagonals. However, to have the CCFDM with a dispersion relation equivalent to that of the QSFEM (Babuška et al. 1995) it is necessary to increase the dimension of the local approximation subspace $S^L$ ($L=9$).

Figure 3
Dispersion relations for $$k^h-k{k}^h-k\text{as a function of}khkhfork=100k=100,\theta =\pi /4\theta =\pi /4,{\theta }_1=\pi /8{\theta }_1=\pi /8and{\theta }_2=3\pi /8{\theta }_2=3\pi /8$$.

Figure 4
Dispersion relations for $𝐤\mathbf{=}\mathrm{𝟏𝟎𝟎}$ and $𝐤𝐡\mathbf{=}\mathrm{𝟏}$ as a function of $𝛉$-direction, ${𝛉}_{\mathrm{𝟏}}\mathbf{=}𝛑\mathbf{/}\mathrm{𝟖}$ and ${𝛉}_{\mathrm{𝟐}}\mathbf{=}\mathrm{𝟑}𝛑\mathbf{/}\mathrm{𝟖}$.

2D case with 9-point stencil

Considering 9-point stencil on uniform meshes, and substituting this numerical solution in the equation (17) the two restrictions on the dispersion relation are obtained:

$$\begin{array}{c}\hfill [A_1+C_1]cos(k^{h}hcos\left(\theta \right))+[A_2+C_2]cos(k^{h}hsin\left(\theta \right))\hfill \\ \hfill +[A_3+C_4][cos(k^{h}hcos\left(\theta \right)+k^{h}hsin\left(\theta \right)]\hfill \\ \hfill +[A_4+C_3][cos(k^{h}hcos\left(\theta \right)-k^{h}hsin\left(\theta \right)]+B=0,\hfill \end{array}$$(62)

$$\begin{array}{c}\hfill [C_1-A_1]sin(k^{h}hcos\left(\theta \right))+[C_2-A_2]sin(k^{h}hsin\left(\theta \right))\hfill \\ \hfill +[C_4-A_3][sin(k^{h}hcos\left(\theta \right)+k^{h}hsin\left(\theta \right)]\hfill \\ \hfill +[C_3-A_4][sin(k^{h}hcos\left(\theta \right)-k^{h}hsin\left(\theta \right)]=0.\hfill \end{array}$$(63)

The equation of the imaginary part (63) introduces a restriction that can be satisfied by requiring $A_1=C_1$, $A_2=C_2$, $A_3=C_4$ and $A_4=C_3$. This restriction is not sufficient to guarantee stencil symmetry with respect to the $X$ and $Y$ axes. For this, an additional restriction $A_3=A_4$ is necessary. The symmetry of the stencil with respect to the straight lines $y=x$ and $y=-x$ is guaranteed if $A_1=A_2$. Analogous to the 2D case with 5-point stencil, there are two stencils with different symmetry that generate the same dispersion relation. Here only the case of the stencil with greater symmetry will be presented. That is, the case of the symmetric stencil with respect to the $X$, $Y$ axes and the straight lines $y=x$ and $y=-x$. More details can be found in Nunes (2024).

CCFDM-9p-L5p

For the CCFDM-9p-L5p, the stencil equation (17) is determined by

$$A_1=\frac{1}{h^2}[{\alpha }_1+(1+\frac{{\left(kh\right)}^2}{4}){\alpha }_2+(1+\frac{{\left(kh\right)}^2}{2}){\alpha }_3],$$(64)

$$A_2=\frac{1}{h^2}[{\alpha }_1+(1+\frac{{\left(kh\right)}^2}{4}){\alpha }_2+(1-\frac{{\left(kh\right)}^2}{2}){\alpha }_3],$$(65)

$$A_3=\frac{1}{h^2}[(\frac{1}{2}+\frac{{\left(kh\right)}^2}{4}){\alpha }_4+(\frac{1}{2}-\frac{{\left(kh\right)}^2}{2}){\alpha }_5],$$(66)

$$A_4=\frac{1}{h^2}[(\frac{1}{2}+\frac{{\left(kh\right)}^2}{4}){\alpha }_4+(\frac{1}{2}+\frac{{\left(kh\right)}^2}{2}){\alpha }_5],$$(67)

$$B=\frac{1}{h^2}\{[-4+{\left(kh\right)}^2]({\alpha }_1+{\alpha }_3)-2(2{\alpha }_2+{\alpha }_4)+[-2{\left(kh\right)}^2]{\alpha }_5\},$$(68)

$$C_1=\frac{1}{h^2}[{\alpha }_1+(1+\frac{{\left(kh\right)}^2}{4}){\alpha }_2+(1+\frac{{\left(kh\right)}^2}{2}){\alpha }_3],$$(69)

$$C_2=\frac{1}{h^2}[{\alpha }_1+(1+\frac{{\left(kh\right)}^2}{4}){\alpha }_2+(1-\frac{{\left(kh\right)}^2}{2}){\alpha }_3],$$(70)

$$C_3=\frac{1}{h^2}[(\frac{1}{2}+\frac{{\left(kh\right)}^2}{4}){\alpha }_4+(\frac{1}{2}+\frac{{\left(kh\right)}^2}{2}){\alpha }_5],$$(71)

$$C_4=\frac{1}{h^2}[(\frac{1}{2}+\frac{{\left(kh\right)}^2}{4}){\alpha }_4+(\frac{1}{2}-\frac{{\left(kh\right)}^2}{2}){\alpha }_5].$$(72)

The generating vectors $G_3$ and $G_5$ do not have the required symmetry, and can be disregarded by making ${\alpha }_3={\alpha }_5=0$. Thus, only the free parameters ${\alpha }_2$ and ${\alpha }_4$ remain to be determined, since ${\alpha }_1=1$. These coefficients ${\alpha }_2$ and ${\alpha }_4$ can be determined to minimize $|k-k^h|$. Choosing the directions ${\theta }_1=\pi /16$ and ${\theta }_2=3\pi /16$ for the plane wave it is possible to obtain

$$\begin{array}{cc}\hfill {\alpha }_2& =\frac{8(-c_1+c_2-s_1+s_2+2c_1{s}_1-2c_2{s}_2-c_1{c}_2{s}_1+c_1{c}_2{s}_2-c_1{s}_1{s}_2+c_2{s}_1{s}_2)}{{\chi }_1^{CO}+{\chi }_2^{CO}+{\chi }_3^{CO}}\hfill \\ \hfill & +\frac{4k^2{h}^2(c_1{s}_1-c_2{s}_2-c_1{c}_2{s}_1+c_1{c}_2{s}_2-c_1{s}_1{s}_2+c_2{s}_1{s}_2)+2k^4{h}^4(c_2{s}_2-c_1{s}_1)}{{\chi }_1^{CO}+{\chi }_2^{CO}+{\chi }_3^{CO}},\hfill \end{array}$$(73)

$$\begin{array}{cc}\hfill {\alpha }_4& =\frac{k^4{h}^4(c_1-c_2+s_1-s_2)}{{\chi }_1^{CO}+{\chi }_2^{CO}+{\chi }_3^{CO}},\hfill \end{array}$$(74)

$${\chi }_1^{CO}=8(c_1-c_2+s_1-s_2-2c_1{s}_1+2c_2{s}_2+c_1{c}_2{s}_1-c_1{c}_2{s}_2+c_1{s}_1{s}_2-c_2{s}_1{s}_2),$$(75)

$${\chi }_2^{CO}=2k^2{h}^2(c_1-c_2+s_1-s_2-4c_1{s}_1+4c_2{s}_2+3c_1{c}_2{s}_1-3c_1{c}_2{s}_2+3c_1{s}_1{s}_2-3c_2{s}_1{s}_2),$$(76)

$${\chi }_3^{CO}=k^4{h}^4(c_1{c}_2{s}_1-c_1{c}_2{s}_2+c_1{s}_1{s}_2{c}_2{s}_1{s}_2).$$(77)

The dispersion relation of the CCFDM-9p-L5p is equivalent to the dispersion of the QSFEM for ${\theta }_1=\frac{\pi }{16}$ and ${\theta }_2=\frac{3\pi }{16}$ (Babuška et al. 1995),

$$k^h-k=\frac{k^7{h}^6[cos\left(4\theta \right)-cos\left(4{\theta }_1\right)][cos\left(4\theta \right)-cos\left(4{\theta }_2\right)]}{387072}+𝒪\left(k^9{h}^8\right),$$(78)

where the parameters ${\alpha }_2$ (73) and ${\alpha }_4$ (74) were replaced. Figure 5 shows the dispersion relation for the analyzed schemes.

Figure 5
Dispersion relations for $𝐤\mathbf{=}\mathrm{𝟏𝟎𝟎}$ and $𝐤𝐡\mathbf{=}\mathrm{𝟏}$ as a function of $𝛉$-direction, ${𝛉}_{\mathrm{𝟏}}\mathbf{=}𝛑\mathbf{/}\mathrm{𝟏𝟔}$ and ${𝛉}_{\mathrm{𝟐}}\mathbf{=}\mathrm{𝟑}𝛑\mathbf{/}\mathrm{𝟏𝟔}$.

CCFDM-9p-L9p

For the CCFDM-9p-L9p, the stencil equation (17) is determined by

$$\begin{array}{cc}\hfill A_1& =\frac{1}{h^2}[\frac{2}{3}({\alpha }_1+{\alpha }_4+{\alpha }_5)+(\frac{2}{3}+\frac{{\left(kh\right)}^2}{4}){\alpha }_2+(\frac{2}{3}+\frac{{\left(kh\right)}^2}{2}){\alpha }_3],\hfill \end{array}$$(79)

$$\begin{array}{cc}\hfill A_2& =\frac{1}{h^2}[\frac{2}{3}({\alpha }_1+{\alpha }_4+{\alpha }_5)+(\frac{2}{3}+\frac{{\left(kh\right)}^2}{4}){\alpha }_2+(\frac{2}{3}-\frac{{\left(kh\right)}^2}{2}){\alpha }_3],\hfill \end{array}$$(80)

$$\begin{array}{cc}\hfill A_3& =\frac{1}{h^2}[\left(\frac{{\alpha }_1+{\alpha }_2+{\alpha }_3}{6}\right)+(\frac{1}{6}+\frac{{\left(kh\right)}^2}{4}){\alpha }_4+(\frac{1}{6}-\frac{{\left(kh\right)}^2}{2}){\alpha }_5],\hfill \end{array}$$(81)

$$\begin{array}{cc}\hfill A_4& =\frac{1}{h^2}[\left(\frac{{\alpha }_1+{\alpha }_2+{\alpha }_3}{6}\right)+(\frac{1}{6}+\frac{{\left(kh\right)}^2}{4}){\alpha }_4+(\frac{1}{6}+\frac{{\left(kh\right)}^2}{2}){\alpha }_5],\hfill \end{array}$$(82)

$$\begin{array}{cc}\hfill B& =\frac{1}{h^2}[(\frac{-10}{3}+{\left(kh\right)}^2)({\alpha }_1+{\alpha }_3+{\alpha }_5)-\frac{10}{3}({\alpha }_2+{\alpha }_4)],\hfill \end{array}$$(83)

$$\begin{array}{cc}\hfill C_1& =\frac{1}{h^2}[\frac{2}{3}({\alpha }_1+{\alpha }_4+{\alpha }_5)+(\frac{2}{3}+\frac{{\left(kh\right)}^2}{4}){\alpha }_2+(\frac{2}{3}+\frac{{\left(kh\right)}^2}{2}){\alpha }_3],\hfill \end{array}$$(84)

$$\begin{array}{cc}\hfill C_2& =\frac{1}{h^2}[\frac{2}{3}({\alpha }_1+{\alpha }_4+{\alpha }_5)+(\frac{2}{3}+\frac{{\left(kh\right)}^2}{4}){\alpha }_2+(\frac{2}{3}-\frac{{\left(kh\right)}^2}{2}){\alpha }_3],\hfill \end{array}$$(85)

$$\begin{array}{cc}\hfill C_3& =\frac{1}{h^2}[\left(\frac{{\alpha }_1+{\alpha }_2+{\alpha }_3}{6}\right)+(\frac{1}{6}+\frac{{\left(kh\right)}^2}{4}){\alpha }_4+(\frac{1}{6}+\frac{{\left(kh\right)}^2}{2}){\alpha }_5],\hfill \end{array}$$(86)

$$\begin{array}{cc}\hfill C_4& =\frac{1}{h^2}[\left(\frac{{\alpha }_1+{\alpha }_2+{\alpha }_3}{6}\right)+(\frac{1}{6}+\frac{{\left(kh\right)}^2}{4}){\alpha }_4+(\frac{1}{6}-\frac{{\left(kh\right)}^2}{2}){\alpha }_5].\hfill \end{array}$$(87)

The generating vectors $G_3$ and $G_5$ do not have the required symmetry, and can be disregarded by making ${\alpha }_3={\alpha }_5=0$. Thus, only the free parameters ${\alpha }_2$ and ${\alpha }_4$ remain to be determined, since ${\alpha }_1=1$. These coefficients ${\alpha }_2$ and ${\alpha }_4$ can be determined to minimize $|k-k^h|$. Choosing the directions ${\theta }_1=\pi /16$ and ${\theta }_2=3\pi /16$ for the plane wave it is possible to obtain

$$\begin{array}{cc}\hfill {\alpha }_2& =\frac{8(-c_1+c_2-s_1+s_2+2c_1{s}_1-2c_2{s}_2-c_1{c}_2{s}_1+c_1{c}_2{s}_2-c_1{s}_1{s}_2+c_2{s}_1{s}_2)}{{\chi }_1^{L9p}+{\chi }_2^{L9p}}\hfill \\ \hfill & +\frac{6k^2{h}^2(-c_1{s}_1-c_2{s}_2)}{{\chi }_1^{L9p}+{\chi }_2^{L9p}},\hfill \end{array}$$(88)

$$\begin{array}{cc}\hfill {\alpha }_4& =\frac{2(-c_1+c_2-s_1+s_2+2c_1{s}_1-2c_2{s}_2-c_1{c}_2{s}_1+c_1{c}_2{s}_2-c_1{s}_1{s}_2+c_2{s}_1{s}_2)}{{\chi }_1^{L9p}+{\chi }_2^{L9p}}\hfill \\ \hfill & +\frac{3k^2{h}^2(-c_1-c_2+s_1-s_2)}{{\chi }_1^{L9p}+{\chi }_2^{L9p}},\hfill \end{array}$$(89)

$${\chi }_1^{L9p}=10(c_1-c_2+s_1-s_2-2c_1{s}_1+2c_2{s}_2+c_1{c}_2{s}_1-c_1{c}_2{s}_2+c_1{s}_1{s}_2-c_2{s}_1{s}_2),$$(90)

$${\chi }_2^{L9p}=3k^2{h}^2(c_1{c}_2{s}_1-c_1{c}_2{s}_2+c_1{s}_1{s}_2-c_2{s}_1{s}_2).$$(91)

The dispersion relation of the CCFDM-9p-L9p is equivalent to the dispersion of the QSFEM for ${\theta }_1=\frac{\pi }{16}$ and ${\theta }_2=\frac{3\pi }{16}$ (Babuška et al. 1995),

$$k^h-k=\frac{k^7{h}^6[cos\left(4\theta \right)-cos\left(4{\theta }_1\right)][cos\left(4\theta \right)-cos\left(4{\theta }_2\right)]}{387072}+𝒪\left(k^9{h}^8\right),$$(92)

where the parameters ${\alpha }_2$ (88) and ${\alpha }_4$ (89) were replaced. Figure 6 shows the dispersion relation for the analyzed schemes.

Figure 6
Dispersion relations for $𝐤\mathbf{=}\mathrm{𝟏𝟎𝟎}$ and $𝐤𝐡\mathbf{=}\mathrm{𝟏}$ as a function of $𝛉$-direction, ${𝛉}_{\mathrm{𝟏}}\mathbf{=}𝛑\mathbf{/}\mathrm{𝟏𝟔}$ and ${𝛉}_{\mathrm{𝟐}}\mathbf{=}\mathrm{𝟑}𝛑\mathbf{/}\mathrm{𝟏𝟔}$.

Table I presents the distance $d:=(k^h-k)h$ of the CCFDM and QSFEM for different values of $\theta$ with $kh=1$ fixed. The values presented were calculated using equations (56), (78) and (92). The CCFDM-5p presents zero dispersion for the two angles already mentioned ${\theta }_1=\pi /8$ and ${\theta }_2=3\pi /8$, which is equivalent to the GLS (Thompson Pinsky 1995). Increasing the dimension of the local approximation subspace $S^L$ to $L=9$ two versions CCFDM-9p-L5p and CCFDM-9p-L9p are obtained, which present a dispersion relation equivalent to QSFEM (Babuška et al. 1995).

Table I
Comparison of the dispersion of each method depending on the wave direction.

$\theta$	CCFDM-5p	CCFDM-9p-L5p	CCFDM-9p-L9p	QSFEM
$0$	0.0104166666	0.1291749338$\times {10}^{-5}$	0.1291749338$\times {10}^{-5}$	0.1291749338$\times {10}^{-5}$
$\pi /16$	0.0073656956	0	0	0
$\pi /8$	0	-0.1291749338$\times {10}^{-5}$	-0.1291749338$\times {10}^{-5}$	-0.1291749338$\times {10}^{-5}$
$3\pi /16$	-0.0073656956	0	0	0
$\pi /4$	-0.0104166666	0.1291749338$\times {10}^{-5}$	0.1291749338$\times {10}^{-5}$	0.1291749338$\times {10}^{-5}$
$5\pi /16$	-0.0073656956	0	0	0
$3\pi /8$	0	-0.1291749338$\times {10}^{-5}$	-0.1291749338$\times {10}^{-5}$	-0.1291749338$\times {10}^{-5}$
$7\pi /16$	0.0073656956	0	0	0
$\pi /2$	0.0104166666	0.1291749338$\times {10}^{-5}$	0.1291749338$\times {10}^{-5}$	0.1291749338$\times {10}^{-5}$

Consequence of the uniform mesh in the dispersion analysis and stencil symmetry

The equations (21), (28) and (63) of the imaginary part of the dispersion analysis for uniform meshes imply the following corollary.

Corollary 1. Let it be a uniform mesh with parameter $h=h_x=h_y$. For each central node $I$ (global ordering or enumeration) consider stencil with 3-points in the 1D case, stencil with 5-points and 9-points in the 2D case. Let $G_l$ be the generating vectors of the CCFDM for these 1D and 2D cases. If the equations (21), (28) and (63) of the imaginary part of the dispersion analysis are satisfied, then there are only:

• two linearly independent generating vectors in the 1D case,

• three linearly independent generating vectors in the 2D case and 5-point stencil,

• five linearly independent generating vectors in the 2D case and 9-point stencil.

Proof 1. Item i) In the 1D case and 3-point stencil, the CCFDM is determined by three linearly independent generating vectors $G_1=[A_1,B_1,C_1]$, $G_2=[A_2,B_2,C_2]$ and $G_3=[A_3,B_3,C_3]$. For a uniform mesh, the equation (21) of the imaginary part of the dispersion analysis implies $A=C$, therefore $G_1=[A_1,B_1,A_1]$, $G_2=[A_2,B_2,A_2]$ and $G_3=[A_3,B_3,A_3]$. These vectors are linearly independent, when ${\alpha }_1{G}_1+{\alpha }_2{G}_2+{\alpha }_3{G}_3=0$ only if ${\alpha }_1={\alpha }_2={\alpha }_3=0$. In matrix form this restriction is:

$$\left[\begin{array}{ccc}\hfill A_1\hfill & \hfill A_2\hfill & \hfill A_3\hfill \\ \hfill B_1\hfill & \hfill B_2\hfill & \hfill B_3\hfill \\ \hfill A_1\hfill & \hfill A_2\hfill & \hfill A_3\hfill \\ \hfill \hfill \end{array}\right]\left[\begin{array}{c}\hfill {\alpha }_1\hfill \\ \hfill {\alpha }_2\hfill \\ \hfill {\alpha }_3\hfill \\ \hfill \hfill \end{array}\right]=\left[\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \\ \hfill \hfill \end{array}\right].$$

From which it can be noted that this matrix only has two linearly independent rows. That is, the rank of the matrix is 2 (deficient rank). Items ii) and iii) can be demonstrated in a similar way. More details can be found in Nunes 2024.$▫$

On the other hand, this equation corresponding to the imaginary part is related to the symmetry of the stencil used, thus establishing a connection between the coefficients of the linear system and the symmetry of the stencil. For example, in the 2D case the two stencils analyzed (5-point stencil and 9-point stencil) are shown in the Figure 7.

Figure 7
Two stencils without any constraints: (a) 5-point stencil and (b) 9-point stencil.

The imaginary part of the dispersion equation (28) establishes the constraints $C_1=A_1$ and $C_2=A_2$ for the stencil in the Figure 7a. The imaginary part of the dispersion equation (63) establishes the constraints $C_1=A_1$, $C_2=A_2$, $C_4=A_3$ and $C_3=A_4$ for the stencil in the Figure 7b. From an algebraic point of view, these are restrictions for the coefficients of the linear system of equations $\mathbf{\text{𝐓}}{U}^h=F$. From a geometric point of view, these are restrictions for the symmetry of the stencil. Applying these restrictions, the stencils in the Figure 7 are transformed into the stencils in the Figure 8.

Figure 8
Two stencils with the constraint of the equation of the imaginary part.

In Figure 8a it can be seen that the restriction of the imaginary part of the dispersion relation guarantees symmetry of the stencil with respect to the $X$ and $Y$ axes. However, the symmetry of this stencil with respect to the straight lines $y=x$ and $y=-x$ is only guaranteed if $A_2=A_1$ as shown in Figure 9a. In Figure 8b it can be seen that, even after the restrictions imposed by the imaginary part of the dispersion equation, this stencil still does not have symmetry with respect to the $X$ and $Y$ axes. To guarantee this symmetry it is necessary to impose another restriction, $A_4=A_3$. And the symmetry of this stencil with respect to the straight lines $y=x$ and $y=-x$ is only guaranteed if $A_1=A_2$ as shown in Figure 9b. It should be noted that the symmetry of the stencil with respect to the lines $y=x$ and $y=-x$ will exist if and only if $h_x=h_y$, that is, a uniform mesh. By imposing all these symmetry restrictions, the stencils in Figure 9 are obtained. Thus, this symmetry analysis shows that there is a connection between the symmetry of the problem defined in the continuous domain (equations (1) and (2)) and the discretized problem, and this connection is established through the equation of imaginary part of the dispersion analysis.

Figure 9
Two stencils with all possible symmetry constraints.

Local Truncation Error

The well-known local truncation error $E_I^L$ is determined at each point $I$ (global ordering or enumeration) by replacing the approximate solution $U_I$ with the exact solution $u_I$ in the stencil equation for each scheme. Since $u$ is unknown, it is assumed to be a function smooth enough that it is possible to replace $u$ by its Taylor series expansions around the point $I$ (LeVeque 2007). The errors for each scheme analyzed are presented below. In the case of CCFDM the local truncation error is represented as a function of the free parameters ${\alpha }_l$.

1D case with 3-point stencil

Consider partial ordering or enumeration for the central node $i$. Using the notation $u_i^{\left(n\right)}=\frac{d^{n}u\left(x_i\right)}{d{x}^n}$ and $f_i^{\left(n\right)}=\frac{d^{n}f\left(x_i\right)}{d{x}^n}$ can be obtained:

$$\mathbf{\text{𝐂𝐂:}}{E}_i^L=\frac{f_i^{\left(2\right)}-k^2{u}_i^{\left(2\right)}}{12}{h}^2+𝒪\left(h^3\right),$$(93)

$$\mathbf{\text{𝐍𝐒-𝟏:}}{E}_i^L=\frac{f_i^{\left(2\right)}}{12}{h}^2+\frac{f_i^{\left(2\right)}-k^2{u}_i^{\left(2\right)}}{24}{k}^2{h}^4+𝒪\left(h^5\right),$$(94)

$$\mathbf{\text{𝐍𝐒-𝟐:}}{E}_i^L=k^2{u}_i^{\left(1\right)}h+\frac{f_i^{\left(2\right)}-k^2{u}_i^{\left(2\right)}}{12}{h}^2+\frac{f_i^{\left(1\right)}-k^2{u}_i^{\left(1\right)}}{6}{k}^2{h}^3+𝒪\left(h^4\right),$$(95)

$$\mathbf{\text{𝐂𝐂𝐅𝐃𝐌:}}{E}_i^L=(f_i^{\left(2\right)}-k^2{u}_i^{\left(2\right)})\frac{h^2}{12}+{\alpha }_2[f_i^{\left(2\right)}\frac{h^2}{12}+(f_i^{\left(2\right)}-k^2{u}_i^{\left(2\right)})k^2\frac{h^4}{24}]+𝒪\left(h^5\right).$$(96)

It should be noted that the CC, NS-1 and CCFDM schemes have accuracy order $h^2$. However, the NS-2 scheme has accuracy order $h$. The NS-2 scheme had no impact on the construction of CCFDM because ${\alpha }_3=0$.

2D case with 5-point stencil

Consider partial ordering or enumeration for the central node $i,j$. Using the notation $u_{i,j}^{\left(nx\right)}=\frac{{\partial }^{n}u}{\partial x^n}{|}_{i,j}$ and $u_{i,j}^{\left(ny\right)}=\frac{{\partial }^{n}u}{\partial y^n}{|}_{i,j}$ can be obtained:

$$\mathbf{\text{𝐂𝐂:}}{E}_{i,j}^L=(u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)})\frac{h^2}{12}+𝒪\left(h^3\right),$$(97)

$$\mathbf{\text{𝐍𝐒-𝟏:}}{E}_{i,j}^L=(u_{i,j}^{\left(2x\right)}+u_{i,j}^{\left(2y\right)})\frac{k^2{h}^2}{4}+(u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)})(\frac{h^2}{12}+\frac{k^2{h}^4}{48})+𝒪\left(h^5\right),$$(98)

$$\mathbf{\text{𝐍𝐒-𝟐:}}{E}_{i,j}^L=(u_{i,j}^{\left(1x\right)}+u_{i,j}^{\left(1y\right)})\frac{k^{2}h}{2}+(u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)})\frac{h^2}{12}+(u_{i,j}^{\left(3x\right)}+u_{i,j}^{\left(3y\right)})\frac{k^2{h}^3}{12}+𝒪\left(h^4\right),$$(99)

$$\mathbf{\text{𝐍𝐒-𝟑:}}{E}_{i,j}^L=(u_{i,j}^{\left(1y\right)}-u_{i,j}^{\left(1x\right)})\frac{k^{2}h}{2}+(u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)})\frac{h^2}{12}+(u_{i,j}^{\left(3y\right)}-u_{i,j}^{\left(3x\right)})\frac{k^2{h}^3}{12}+𝒪\left(h^4\right),$$(100)

$$\mathbf{\text{𝐍𝐒-𝟒:}}{E}_{i,j}^L=(u_{i,j}^{\left(2x\right)}-u_{i,j}^{\left(2y\right)})\frac{k^2{h}^2}{2}+(u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)})\frac{h^2}{12}+(u_{i,j}^{\left(4x\right)}-u_{i,j}^{\left(4y\right)})\frac{k^2{h}^4}{12}+𝒪\left(h^5\right),$$(101)

$$\begin{array}{c}\hfill \mathbf{\text{𝐂𝐂𝐅𝐃𝐌:}}{E}_{i,j}^L=(u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)})\frac{h^2}{12}+{\alpha }_2[(u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)})(\frac{h^2}{12}+\frac{k^2{h}^4}{48})\hfill \\ \hfill +(u_{i,j}^{\left(2x\right)}+u_{i,j}^{\left(2y\right)})\frac{k^2{h}^2}{4}]+{\alpha }_5[(u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)})\frac{h^2}{12}+(u_{i,j}^{\left(4x\right)}-u_{i,j}^{\left(4y\right)})\frac{k^2{h}^4}{24}\hfill \\ \hfill +(u_{i,j}^{\left(2x\right)}-u_{i,j}^{\left(2y\right)})\frac{k^2{h}^2}{2}]+𝒪\left(h^5\right).\hfill \end{array}$$(102)

It should be noted that the CC, NS-1, NS-4 and CCFDM schemes have accuracy order $h^2$. However, schemes NS-2 and NS-3 have accuracy order $h$. The NS-2 and NS-3 schemes had no impact on the construction of CCFDM because ${\alpha }_3={\alpha }_4=0$.

2D case with 9-point stencil

CCFDM-9p-L5p

The first three schemes CC, NS-1 and NS-2 generated by the vectors $G_1$, $G_2$ and $G_3$ (black stencil) have the same $E_{i,j}^L$ that the first three schemes of the 2D case with 5-point stencil (equations (97), (98) and (99)). The new schemes NS-3 and NS-4 generated by the vectors $G_4$ and $G_5$ of the rotated 5-point stencil (blue stencil) and the CCFDM-9p-L5p have the following $E_{i,j}^L$:

$$\begin{array}{c}\hfill \mathbf{\text{𝐍𝐒-𝟑:}}{E}_{i,j}^L=[u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)}](\frac{h^2}{12}+\frac{k^2{h}^4}{24})+[u_{i,j}^{\left(2x\right)}+u_{i,j}^{\left(2y\right)}]\frac{k^2{h}^2}{4}+u_{i,j}^{(2x,2y)}(\frac{h^2}{12}+\frac{k^2{h}^4}{4})+𝒪\left(h^5\right),\hfill \end{array}$$(103)

$$\mathbf{\text{𝐍𝐒-𝟒:}}{E}_{i,j}^L=[u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)}](\frac{h^2}{12}-\frac{k^2{h}^4}{3})-2u_{i,j}^{(1x,1y)}{k}^2{h}^2+u_{i,j}^{(2x,2y)}\frac{h^2}{12}+𝒪\left(h^5\right),$$(104)

$$\begin{array}{c}\hfill \mathbf{\text{𝐂𝐂𝐅𝐃𝐌-𝟗𝐩-𝐋𝟓𝐩:}}{E}_{i,j}^L({\alpha }_2,{\alpha }_4)=[u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)}]\frac{h^2}{12}\hfill \\ \hfill +{\alpha }_2\{[u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)}](\frac{h^2}{12}+\frac{k^2{h}^4}{48})+[u_{i,j}^{\left(2x\right)}+u_{i,j}^{\left(2y\right)}]\frac{k^2{h}^4}{4}\}\hfill \\ \hfill +{\alpha }_4\{[u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)}](\frac{h^2}{12}+\frac{k^2{h}^4}{24})+u_{i,j}^{(2x,2y)}(\frac{h^2}{2}+\frac{k^2{h}^2}{4})+[u_{i,j}^{\left(2x\right)}+u_{i,j}^{\left(2y\right)}]\frac{k^2{h}^2}{2}\}.\hfill \end{array}$$(105)

It should be noted that the CC, NS-1, NS-3, NS-4 and CCFDM schemes have accuracy order $h^2$. However, the NS-2 scheme has accuracy order $h$. The NS-2 and NS-4 schemes had no impact on the CCFDM error because ${\alpha }_3={\alpha }_5=0$.

CCFDM-9p-L9p

The five schemes generated by the vectors $G_1$, $G_2$, $G_3$, $G_4$, $G_5$ and the CCFDM-9p-L9p present the following $E_{i,j}^L$:

$$\mathbf{\text{𝐂𝐂:}}{E}_{i,j}^L=[u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)}]\frac{h^2}{12}+\left[u_{i,j}^{(2x,2y)}\right]\frac{h^2}{6}+𝒪\left(h^3\right),$$(106)

$$\mathbf{\text{𝐍𝐒-𝟏:}}{E}_{i,j}^L=[u_{i,j}^{\left(2x\right)}+u_{i,j}^{\left(2y\right)}]\frac{k^2{h}^2}{4}+[u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)}](\frac{h^2}{12}+\frac{k^2{h}^4}{48})+\left[u_{i,j}^{(2x,2y)}\right]\frac{h^2}{6}+𝒪\left(h^5\right),$$(107)

$$\mathbf{\text{𝐍𝐒-𝟐:}}{E}_{i,j}^L=[u_{i,j}^{\left(2x\right)}-u_{i,j}^{\left(2y\right)}]\frac{k^2{h}^2}{2}+[u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)}]\frac{h^2}{12}+[u_{i,j}^{\left(4x\right)}-u_{i,j}^{\left(4y\right)}]\frac{k^2{h}^4}{24}+𝒪\left(h^5\right),$$(108)

$$\begin{array}{c}\hfill \mathbf{\text{𝐍𝐒-𝟑:}}{E}_{i,j}^L=[u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)}](\frac{h^2}{12}+\frac{k^2{h}^4}{24})+u_{i,j}^{(2x,2y)}(\frac{h^2}{6}+\frac{k^2{h}^4}{4})+[u_{i,j}^{\left(2x\right)}+u_{i,j}^{\left(2y\right)}]\frac{k^2{h}^2}{2}+𝒪\left(h^5\right),\hfill \end{array}$$(109)

$$\begin{array}{c}\hfill \mathbf{\text{𝐍𝐒-𝟒:}}{E}_{i,j}^L=[u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)}]\frac{h^2}{12}+u_{i,j}^{(2x,2y)}\frac{h^2}{6}-2u_{i,j}^{(1x,1y)}{k}^2{h}^2-[u_{i,j}^{(3x,1y)}+u_{i,j}^{(1x,3y)}]\frac{k^2{h}^4}{3}+𝒪\left(h^5\right),\hfill \end{array}$$(110)

$$\begin{array}{c}\hfill \mathbf{\text{𝐂𝐂𝐅𝐃𝐌-𝟗𝐩-𝐋𝟗𝐩:}}{E}_{i,j}^L({\alpha }_2,{\alpha }_4)=[u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)}]\frac{h^2}{12}\hfill \\ \hfill +{\alpha }_2\{[u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)}](\frac{h^2}{12}+\frac{k^2{h}^4}{48})+[u_{i,j}^{\left(2x\right)}+u_{i,j}^{\left(2y\right)}]\frac{k^2{h}^2}{4}+u_{i,j}^{(2x,2x)}\frac{h^2}{6}\}\hfill \\ \hfill +{\alpha }_4\{[u_{i,j}^{\left(2x\right)}+u_{i,j}^{\left(2y\right)}]\frac{k^2{h}^2}{2}+[u_{i,j}^{\left(4x\right)}+u_{i,j}^{\left(4y\right)}](\frac{h^2}{12}+\frac{k^2{h}^4}{24})\hfill \\ \hfill +u_{i,j}^{(2x,2y)}(\frac{h^2}{6}+\frac{k^2{h}^4}{4})\}.\hfill \end{array}$$(111)

It should be noted that all schemes have accuracy order $h^2$. The NS-2 and NS-4 schemes had no impact on the CCFDM error because ${\alpha }_3={\alpha }_5=0$.

RESULTS AND DISCUSSION

Two examples in the 1D case and three examples in the 2D case are presented to illustrate the performance of the new approach. As mentioned previously, for the 2D case there are two versions depending on the degree of symmetry of the stencil. The CCDFM-5p has a less symmetrical version (${\alpha }_1=1$, ${\alpha }_3={\alpha }_4=0$, ${\alpha }_2$ equation (50) and ${\alpha }_5$ equation (51)) and a more symmetrical version (${\alpha }_1=1$, ${\alpha }_3={\alpha }_4={\alpha }_5=0$ and ${\alpha }_2$ equation (57)). The CCDFM-9p-L5p and CCDFM-9p-L9p methods also have a less symmetric and a more symmetric version (Nunes 2024). In each case, the two versions have the same dispersion relationship and numerical results. Here the results are presented only for the versions corresponding to the stencil with greater symmetry. Furthermore, the CCDFM-9p-L5p and CCDFM-9p-L9p methods for the 9-point stencil were developed by constructing different bases for the local approximation subspace. These two versions have the same dispersion relation. The numerical results show that the solutions of these two versions are equal to minus the round-off difference.

All finite difference schemes analyzed were implemented in computational codes developed using MATLAB software (Gilat 2014). The linear system of algebraic equations corresponding to each scheme is solved using the ‘backslash’ operator. As the matrix of the linear system is sparse, the MATLAB function ‘spdiags’ was used. This allows to reduce the computational cost in terms of execution time and storage memory. More details about the computational implementation can be found in Nunes 2024, including the computational codes. The numerical results were calculated on a Acer laptop, Linux 64-bit operating system, AMD Ryzen 7 4800H processor 2.9GHz with 8 core and 16 threads up to 4.2GHz, 2TB SSD, 16GB of RAM memory DDR4 - 3200MHz, and NVIDIA GeForce GTX 1650 Ti 4GB graphics card.

Two examples in the 1D case

The first example consists of solving the homogeneous Helmholtz equation (1) on the interval $]0,1[$ subjected to Dirichlet boundary conditions such that the exact solution is $u\left(x\right)=sin\left(kx\right)/sin\left(k\right)$. Figure 10 shows the solutions of the CC and CCFDM methods together with the interpolant for wavenumber $k=4000$. The results are presented only for the $0.55\le x\le 0.60$ interval, as the wavenumber is high. It can be noted that the CCFDM solution coincides with the interpolant. Figure 11 shows the relative error $R{E}^G=\frac{{\parallel U-u\parallel }_2}{{\parallel u\parallel }_2}$ calculated in the $l_2$ vector norm as a function of the number of degrees of freedom $N$. It can be noted that the relative error of the CCFDM method is approximately 13 orders of magnitude smaller than the error of the CC method. In other words, considering the round-off it is possible to state that the CCFDM error is of the order of computer precision.

Figure 10
CC and CCFDM solutions with $𝐟\mathbf{=}\mathrm{𝟎}$, $𝐤\mathbf{=}\mathrm{𝟒𝟎𝟎𝟎}$ and $𝐤𝐡\mathbf{=}\mathbf{0.5}$.

Figure 11
Relative error in the ${𝐥}_{\mathrm{𝟐}}$ norm with $𝐟\mathbf{=}\mathrm{𝟎}$, $𝐤\mathbf{=}\mathrm{𝟏𝟎𝟎𝟎}$ and $\mathrm{𝟓𝟎𝟎}\mathbf{\le }𝐍\mathbf{\le }\mathrm{𝟔𝟎𝟎𝟎}$.

The second example has a source term $f\left(x\right)=2+{\left(kx\right)}^2$ and Dirichlet boundary conditions such that the exact solution is $u\left(x\right)=x^2+cos\left(kx\right)$. Figure 12 shows the solutions of the CC and CCFDM methods together with the interpolant for wavenumber $k=4000$. It can be noted that the CCFDM solution coincides with the interpolant.

Figure 12
CC and CCFDM solutions with $𝐟\mathbf{\left(}𝐱\mathbf{\right)}\mathbf{=}\mathrm{𝟐}\mathbf{+}{\mathbf{\left(}𝐤𝐱\mathbf{\right)}}^{\mathrm{𝟐}}$, $𝐤\mathbf{=}\mathrm{𝟒𝟎𝟎𝟎}$ and $𝐤𝐡\mathbf{=}\mathbf{0.5}$.

Three examples in the 2D case

As a first example, the homogeneous Helmholtz equation is considered over a unitary square domain with Dirichlet boundary conditions such that the exact solution is a plane wave propagating in the $\theta$-direction: $u(x,y)=cos\left(k(xcos\theta +ysin\theta )\right)$. Figure 13 shows the solutions of the two versions of the CCFDM together with the interpolant at section $y=0.5$ for wavenumber $k=1000$, $\theta =0$ and $kh=0.5$. Due to the high wavenumber, three $x$ intervals were zoomed. The $\theta =0$ wave direction was chosen, since in this wave direction the CCDFM-9p does not coincide with the interpolant. As shown by the dispersion analysis, the CCFDM on the 9-point stencil is equivalent to QSFEM. Figures 14 and 15 show the same case only with the difference in $kh$ values, which confirm the good performance of the new approach even on very coarse meshes. The CCFDM-5p results are equivalent to the GLS method and have been omitted, for more details see (Nunes 2024).

Figure 13
CCFDM-9P-L5p and CCFDM-9P-L9p solutions at section $𝐲\mathbf{=}\mathbf{0.5}$ for $𝐟\mathbf{=}\mathrm{𝟎}$, $𝐤\mathbf{=}\mathrm{𝟏𝟎𝟎𝟎}$, $𝛉\mathbf{=}\mathrm{𝟎}$ and $𝐤𝐡\mathbf{=}\mathbf{0.5}$.

Figure 14
CCFDM-9P-L5p and CCFDM-9P-L9p solutions at section $𝐲\mathbf{=}\mathbf{0.5}$ for $𝐟\mathbf{=}\mathrm{𝟎}$, $𝐤\mathbf{=}\mathrm{𝟏𝟎𝟎𝟎}$, $𝛉\mathbf{=}\mathrm{𝟎}$ and $𝐤𝐡\mathbf{=}\mathrm{𝟏}$.

Figure 15
CCFDM-9P-L5p and CCFDM-9P-L9p solutions at section $𝐲\mathbf{=}\mathbf{0.5}$ for $𝐟\mathbf{=}\mathrm{𝟎}$, $𝐤\mathbf{=}\mathrm{𝟏𝟎𝟎𝟎}$, $𝛉\mathbf{=}\mathrm{𝟎}$ and $𝐤𝐡\mathbf{=}\mathrm{𝟐}$.

In the second example, the homogeneous problem with superposition of three plane waves is considered, whose exact solution is $u(x,y)={\sum }_{i=1}^{3}cos\left(k(xcos{\theta }_i+ysin{\theta }_i)\right)$. Figure 16 shows the solutions of the two versions of the CCFDM together with the interpolant at section $y=0.5$ for wavenumber $k=1000$, ${\theta }_1=0$, ${\theta }_2=\frac{\pi }{4}$, ${\theta }_3=\frac{\pi }{8}$ and $kh=0.5$. These three wave directions were chosen because none of them coincide with the wave direction in which the CCDFM-9p shows zero dispersion. Figures 17 and 18 show the same case only with the difference in $kh$ values.

Figure 16
CCFDM-9P-L5p and CCFDM-9P-L9p solutions at section $𝐲\mathbf{=}\mathbf{0.5}$ for $𝐟\mathbf{=}\mathrm{𝟎}$, $𝐤\mathbf{=}\mathrm{𝟏𝟎𝟎𝟎}$, ${𝛉}_{\mathrm{𝟏}}\mathbf{=}\mathrm{𝟎}$, ${𝛉}_{\mathrm{𝟐}}\mathbf{=}\frac{𝛑}{\mathrm{𝟒}}$, ${𝛉}_{\mathrm{𝟑}}\mathbf{=}\frac{𝛑}{\mathrm{𝟖}}$ and $𝐤𝐡\mathbf{=}\mathbf{0.5}$.

Figure 17
CCFDM-9P-L5p and CCFDM-9P-L9p solutions at section $𝐲\mathbf{=}\mathbf{0.5}$ for $𝐟\mathbf{=}\mathrm{𝟎}$, $𝐤\mathbf{=}\mathrm{𝟏𝟎𝟎𝟎}$, ${𝛉}_{\mathrm{𝟏}}\mathbf{=}\mathrm{𝟎}$, ${𝛉}_{\mathrm{𝟐}}\mathbf{=}\frac{𝛑}{\mathrm{𝟒}}$, ${𝛉}_{\mathrm{𝟑}}\mathbf{=}\frac{𝛑}{\mathrm{𝟖}}$ and $𝐤𝐡\mathbf{=}\mathrm{𝟏}$.

Figure 18
CCFDM-9P-L5p and CCFDM-9P-L9p solutions at section $𝐲\mathbf{=}\mathbf{0.5}$ for $𝐟\mathbf{=}\mathrm{𝟎}$, $𝐤\mathbf{=}\mathrm{𝟏𝟎𝟎𝟎}$, ${𝛉}_{\mathrm{𝟏}}\mathbf{=}\mathrm{𝟎}$, ${𝛉}_{\mathrm{𝟐}}\mathbf{=}\frac{𝛑}{\mathrm{𝟒}}$, ${𝛉}_{\mathrm{𝟑}}\mathbf{=}\frac{𝛑}{\mathrm{𝟖}}$ and $𝐤𝐡\mathbf{=}\mathrm{𝟐}$.

In the third example, the source term $f(x,y)=4+k^2(x^2+y^2)$ and Dirichlet boundary conditions are considered such that the exact solution is $u(x,y)=x^2+y^2+sin\left[k(xcos\theta +ysin\theta )\right]$. Figure 19 shows the solutions of the two versions of the CCFDM together with the interpolant at sections $y=0.5$ for wavenumber $k=1000$, $\theta =\frac{\pi }{4}$ and $kh=0.5$. Figures 20 and 21 show the same case only with the difference in $kh$ values.

Figure 19
CCFDM-9P-L5p and CCFDM-9P-L9p solutions at section $𝐲\mathbf{=}\mathbf{0.5}$ for $𝐟\mathbf{(}𝐱\mathbf{,}𝐲\mathbf{)}\mathbf{=}\mathrm{𝟒}\mathbf{+}{𝐤}^{\mathrm{𝟐}}\mathbf{(}{𝐱}^{\mathrm{𝟐}}\mathbf{+}{𝐲}^{\mathrm{𝟐}}\mathbf{)}$, $𝐤\mathbf{=}\mathrm{𝟏𝟎𝟎𝟎}$, $𝛉\mathbf{=}\frac{𝛑}{\mathrm{𝟒}}$ and $𝐤𝐡\mathbf{=}\mathbf{0.5}$.

Figure 20
CCFDM-9P-L5p and CCFDM-9P-L9p solutions at section $𝐲\mathbf{=}\mathbf{0.5}$ for $𝐟\mathbf{(}𝐱\mathbf{,}𝐲\mathbf{)}\mathbf{=}\mathrm{𝟒}\mathbf{+}{𝐤}^{\mathrm{𝟐}}\mathbf{(}{𝐱}^{\mathrm{𝟐}}\mathbf{+}{𝐲}^{\mathrm{𝟐}}\mathbf{)}$, $𝐤\mathbf{=}\mathrm{𝟏𝟎𝟎𝟎}$, $𝛉\mathbf{=}\frac{𝛑}{\mathrm{𝟒}}$ and $𝐤𝐡\mathbf{=}\mathrm{𝟏}$.

Figure 21
CCFDM-9P-L5p and CCFDM-9P-L9p solutions at section $𝐲\mathbf{=}\mathbf{0.5}$ for $𝐟\mathbf{(}𝐱\mathbf{,}𝐲\mathbf{)}\mathbf{=}\mathrm{𝟒}\mathbf{+}{𝐤}^{\mathrm{𝟐}}\mathbf{(}{𝐱}^{\mathrm{𝟐}}\mathbf{+}{𝐲}^{\mathrm{𝟐}}\mathbf{)}$, $𝐤\mathbf{=}\mathrm{𝟏𝟎𝟎𝟎}$, $𝛉\mathbf{=}\frac{𝛑}{\mathrm{𝟒}}$ and $𝐤𝐡\mathbf{=}\mathrm{𝟐}$.

Figure 22 shows the relative error $R{E}^G=\frac{{\parallel U-u\parallel }_2}{{\parallel u\parallel }_2}$ calculated in the $l_2$ vector norm as a function of the mesh parameter (lower axis) or $kh$ (upper axis). The CCFDM-5p shows an error similar to the CC scheme, where the main highlight is that the pre-asymptotic behavior appears only from $kh0.1$. This corresponds to a mesh refinement much lower than rule of thumb ($kh\approx 0.6$). This was to be expected because the CCFDM-5p has a dispersion relation equivalent to the GLS method. On the other hand, the CCFDM-9p-L5p and CCFDM-9p-L9p have a dispersion relationship equivalent to the QSFEM method. The main highlight is that they present pre-asymptotic convergence in the region $0.79kh1$ (red straight line with slope $m_1$) and asymptotic in the region $0.2kh0.79$ (black straight line with slope $m_2$). For meshes with $kh0.2$, round-off errors impact the solution, and it would be necessary to use computers with quadruple precision (Strouboulis et al. 2006).

Figure 22
Relative error in the ${𝐥}_{\mathrm{𝟐}}$ norm with $𝐟\mathbf{=}\mathrm{𝟎}$, $𝐤\mathbf{=}\mathrm{𝟑𝟎𝟎}$, $𝛉\mathbf{=}\mathrm{𝟎}$ and $\mathrm{𝟏𝟎𝟎}\mathbf{\times }\mathrm{𝟏𝟎𝟎}\mathbf{\le }𝐍\mathbf{\le }\mathrm{𝟐𝟎𝟎𝟎}\mathbf{\times }\mathrm{𝟐𝟎𝟎𝟎}$.

Figure 23 shows the relative error $R{E}^G$ keeping the $kh=0.5$ value constant for $f=0$, $\theta =0$ and $50\le k\le 1000$. The main highlight is that CCFDM-9p-L5p and CCFDM-9p-L9p are equivalent and have a much smaller error than CCFDM-5p and CC scheme.

Figure 23
Relative error $𝐑{𝐄}^{𝐆}$ with $𝐤𝐡\mathbf{=}\mathbf{0.5}$ constant for $𝛉\mathbf{=}\mathrm{𝟎}$ and $\mathrm{𝟓𝟎}\mathbf{\le }𝐤\mathbf{\le }\mathrm{𝟏𝟎𝟎𝟎}$.

Figure 24 shows the relative error $R{E}^G$ as a function of wave direction for the case $f=0$, $k=200$ and $kh=0.5$. Again, the main highlight is that CCFDM-9p-L5p and CCFDM-9p-L9p are equivalent and have a much smaller error than CCFDM-5p and CC scheme. As expected, the CCFDM-5p shows zero error for the angles ${\theta }_1=\pi /8$ and ${\theta }_1=3\pi /8$ similarly to the GLS method. The CCFDM-9p-L5p and CCFDM-9p-L9p methods show zero error for angles ${\theta }_1=\pi /16$, ${\theta }_2=3\pi /16$, ${\theta }_3=5\pi /16$ and ${\theta }_4=7\pi /16$ similarly to the QSFEM method. The numerical value $R{E}^G$ is not exactly zero due to round-off (Strouboulis et al. 2006).

Figure 24
Relative error as a function of wave direction for the cases $𝐟\mathbf{=}\mathrm{𝟎}$, $𝐤\mathbf{=}\mathrm{𝟐𝟎𝟎}$, $𝐤𝐡\mathbf{=}\mathbf{0.5}$ and $\mathrm{𝟎}\mathbf{\le }𝛉\mathbf{\le }𝛑\mathbf{/}\mathrm{𝟐}$.

All the theoretical analysis and numerical experiments carried out show that the CCFDM-9pL-5p and CCFDM-9pL-9p have the same dispersion relation and the same global error behavior, although they differ in the order of approximation of the Laplacian. This allows us to state that for the new approach CCFDM the order of approximation of the derivatives is not the most important thing, if the number of stencil points is kept the same (dimension of the subspace $S^L$), and if a basis for this subspace is constructed. In other words, the most important issue in terms of stability and accuracy is the ability to completely represent the local approximation subspace. Hence the term ‘Complete’ used in the name of the new approach, since it was possible to build a basis for this subspace in the finite difference framework, that is, any element of this subspace can be generated.

CONCLUSIONS

A new approach in the finite difference framework was developed. This approach fully considers the local approximation subspace determined by the type of stencil chosen, since a vector basis is constructed for this subspace using different schemes generated by linearly independent vectors. Thus, the new method can represent any element of this subspace. That is, any finite difference scheme that belongs to this subspace. Furthermore, some of the vectors used to form this basis are new schemes developed by modifying only the $k^{2}u$ term of the PDE.

In essence, the new approach can be formalized in three simple steps: choosing the dimension of the local approximation subspace, constructing a vector basis for this subspace, and determining the coefficients of the linear combination. When applied to the Helmholtz equation on uniform meshes, CCFDM can always minimize the dispersion relation for all stencils in all dimensions. In the 1D case and 3-point stencil, the absolute minimum is achieved with zero dispersion. In the 2D case and 5-point stencil, the absolute minimum is reached in the directions $\theta ={\theta }_1=\frac{\pi }{8}$ and $\theta ={\theta }_1+\frac{\pi }{4}$ with a dispersion relation equivalent to GLS, where the linear system matrix has only five non-zero diagonals. In the 2D case and 9-point stencil, the CCFDM-9p-L5p and CCFDM-9p-L9p versions were developed using two different bases for the local approximation subspace. Both versions reach the absolute minimum in the directions $\theta ={\theta }_1=\frac{\pi }{16}$ and $\theta ={\theta }_1+\frac{\pi }{8}$ with a dispersion relation equivalent to QSFEM. That is, both versions are equivalent, and this shows that regardless of the basis chosen for the local approximation subspace, the new approach is capable of generating any element of this subspace.

For the first time in dispersion analysis, two equations emerge that must be satisfied. The equation corresponding to the imaginary part is related to the symmetry of the stencil used, thereby establishing a link between the algebraic question (coefficients of the linear system) and the geometric question of the mesh (stencil symmetry). For each scheme, dispersion analysis, calculation of the local truncation error and experimental analysis of the global error were carried out, showing that the CCFDM is consistent. The numerical results demonstrate that the new approach, in addition to presenting good performance, is very flexible and promising for application to other PDEs.

Based on the analysis presented, it is clear that this new approach introduces a paradigm shift in the framework of finite differences. Future research should investigate this approach on non-uniform 1D, 2D, and 3D meshes, and its application to the diffusion-convection equation with dominant convection. Furthermore, it would be beneficial to extend these ideas to other numerical methods involving sparse matrices, such as the finite element method and the finite volume method.

ACKNOWLEDGMENTS

This work was supported by the Universidade Federal Fluminense and the Coordena~{a}o de Aperfeioamento de Pessoal de N'{i}vel Superior (CAPES, Brasil). The authors are grateful for the constructive comments made by the two reviewers. Dedicated to the Memory of Professor Ivo M. Babuka.

Dedicated to the Memory of Professor Ivo M. Babuška.

• 1
  https://www.wolfram.com/mathematica

References

• ALVAREZ GB, LOULA AFD, CARMO EGD ROCHINHA FA. 2006. A discontinuous finite element formulation for Helmholtz equation. Comput Methods Appl Mech Engrg 195(33-36): 4018-4035.
• ALVAREZ GB NUNES HF. 2024. Novos Esquemas de Diferenças Finitas para a Equação de Helmholtz. REMAT: Revista Eletrônica da Matemática 10: e4001.
• BABUŠKA I, IHLENBURG F, PAIK ET SAUTER SA. 1995. A Generalized Finite Element Method for solving the Helmholtz equation in two dimensions with minimal pollution. Comput Methods Appl Mech Engrg 128: 325-359.
• BABUŠKA IM SAUTER SA. 1997. Is the Pollution Effect of the FEM Avoidable for the Helmholtz Equation Considering High Wave Numbers? SIAM J Numer Anal 34(6): 2392-2423.
• CARMO EGD, ALVAREZ GB, LOULA AFD ROCHINHA FA. 2008. A nearly optimal Galerkin projected residual finite element method for Helmholtz problem. Comput Methods Appl Mech Engrg 197(13-16): 1362-1375.
• FERNANDES DT LOULA AFD. 2010. Quasi optimal finite difference method for Helmholtz problem on unstructured grids. Int J Numer Methods Eng 82(10): 1244-1281.
• GILAT A. 2014. MATLAB: An Introduction with Applications. Wiley, 416 p.
• HARARI I HUGHES TJR. 1991. Finite element methods for the Helmholtz equation in an exterior domain: Model problems. Comput Methods Appl Mech Engrg 87(1): 59-96.
• HARARI I HUGHES TJR. 1992. Galerkin/least-squares finite element methods for the reduced wave equation with non-reflecting boundary conditions in unbounded domains. Comput Methods Appl Mech Engrg 98(3): 411-454.
• IHLENBURG F. 1998. Finite element analysis of acoustic scattering. New York: Springer-Verlag, 226 p.
• IHLENBURG F BABUŠKA I. 1995a. Finite element solution of the Helmholtz equation with high wave number Part I: The h-version of the FEM. Comput Math Appl 30(9): 9-37.
• IHLENBURG F BABUŠKA I. 1995b. Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation. Int J Numer Methods Eng 38(22): 3745-3774.
• LAX PD RICHTMYER RD. 1956. Survey of the Stability of Linear Finite Difference Equations. Commun Pure Appl Math 9(2): 267-293.
• LEVEQUE RJ. 2007. Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems. Philadelphia: SIAM - Society for Industrial and Applied Mathematics, 341 p.
• LOULA AFD, ALVAREZ GB, CARMO EGD ROCHINHA FA. 2007. A discontinuous finite element method at element level for Helmholtz equation. Comput Methods Appl Mech Engrg 196(4-6): 867-878.
• NABAVI M, SIDDIQUI MHK DARGAHI J. 2007. A new 9-point sixth-order accurate compact finite-difference method for the Helmholtz equation. J Sound Vib 307(3): 972-982.
• NUNES HF. 2024. Método Completo de Diferenças Finitas Centradas para a Equação de Helmholtz. Volta Redonda: Master’s Thesis, Universidade Federal Fluminense, 133 p. (Unpublished).
• FERNANDES DT. 2009. Métodos de Diferenças Finitas e Elementos Finitos para o Problema de Helmholtz. Petrópolis: Ph.D. Thesis, Laboratório Nacional de Computação Científica, 126 p.
• ROCHINHA FA, ALVAREZ GB, CARMO EGD LOULA AFD. 2007. A locally discontinuous enriched finite element formulation for acoustics. Commun Numer Meth Engng 23: 623-637.
• SINGER I TURKEL E. 1998. High-order finite difference methods for the Helmholtz equation. Comput Methods Appl Mech Engrg 163: 343-358.
• STROUBOULIS T, BABUŠKA I HIDAJAT R. 2006. The generalized finite element method for Helmholtz equation: Theory, computation, and open problems. Comput Methods Appl Mech Engrg 195(37-40): 4711-4731.
• SUTMANN G. 2007. Compact finite difference schemes of sixth order for the Helmholtz equation. J Comput Appl Math 311(1): 15-31.
• THOMPSON L PINSKY P. 1995. A galerkin least squares finite element method for the two-dimensional Helmholtz equation, Int J Numer Methods Eng 38(3): 371-397.
• VICHNEVETSKY R BOWLES JB. 1982. Fourier Analysis of Numerical Approximations of Hyperbolic Equations. Philadelphia: SIAM - Society for Industrial and Applied Mathematics, 152 p.
• WOLFRAM S. 2015. An Elementary Introduction to the Wolfram Language. Champaign: Wolfram Media Inc, 324 p.
• WU T. 2017. A dispersion minimizing compact finite difference scheme for the 2D Helmholtz equation. J Comput Appl Math 203: 497-512.
• WU T XU R. 2018. An optimal compact sixth-order finite difference scheme for the Helmholtz equation. Comput Math Appl 75(7): 2520-2537.

Publication Dates

• Publication in this collection
  25 Nov 2024
• Date of issue
  2024

History

• Received
  18 May 2024
• Accepted
  11 Aug 2024