# Study of the stability of a meshless generalized finite difference scheme applied to the wave equation

^{1}Civil Engineering Faculty, Universidad Michoacana de San Nicolás de Hidalgo, Morelia, Michoacán, Mexico^{2}Faculty of Physics and Mathematics, UMSNH, Morelia, Michoacán, Mexico

When designing and implementing numerical schemes, it is imperative to consider the stability of the applied methods. Prior research has presented different results for the stability of generalized finite-difference methods applied to advection and diffusion equations. In recent years, research has explored a generalized finite-difference approach to the advection-diffusion equation solved on non-rectangular and highly irregular regions using convex, logically rectangular grids. This paper presents a study on the stability of generalized finite difference schemes applied to the numerical solution of the wave equation, solved on clouds of points for highly irregular domains. The stability analysis presented in this work provides significant insights into the proper discretizations needed to obtain stable and satisfactory results. The proposed explicit scheme is *conditionally stable*, while the implicit scheme is *unconditionally stable*. Notably, the stability analyses presented in this paper apply to any scheme which is at least second order in space, not just the proposed approach. The proposed scheme offers effective means of numerically solving the wave equation, particularly for highly irregular domains. By demonstrating the stability of the scheme, this study provides a foundation for further research in this area.

## 1. Introduction

The wave equation is a well-known partial differential equation that is commonly used to model the behavior of a scalar function $u(\overrightarrow{r},t)$, where $\overrightarrow{r}$ denotes the spatial coordinates and *t* represents time. This scalar function can describe various physical phenomena, such as pressure waves on media or the displacement of particles from their equilibrium positions. The wave equation can be expressed as:

*c* is a real non-zero coefficient representing wave propagation velocity, and ∇^{2} denotes the spatial Laplacian operator.

Numerically solving the wave equation is a classical problem in partial differential equations, and many methods have been developed for this purpose. However, existing approaches often require rectangular or symmetrical regions, or involve a high computational cost.

Some finite-difference-based schemes have been developed recently to obtain numerical solutions to different partial differential equations. The main advantage of these schemes is their ability to achieve satisfactory results with low computational cost, making them relatively easy to implement. Some examples of such schemes can be found in [1–4].

After constructing a finite difference method, it is crucial to determine whether the scheme provides a stable approximation. Stability is essential for convergence, and several challenges arise in stability analysis when the scheme becomes more complex. Some authors have proposed bounds for specific problems. For instance, Alcrudo offered a practical selection of the temporal discretization Δ*t* for the 2D scalar advection equation to ensure stability in [5]. This can be achieved by considering the inequality

where Δ*x* and Δ*y* are the spatial discretizations, and *a* and *b* represent the advective velocities in the *x* and *y* directions, respectively.

For its part, other advances have been made in designing and applying finite difference schemes, including the corresponding stability analysis of these schemes. Appadu presents in [6, 7] different schemes to numerically solve the 1D advection-diffusion equation with very satisfactory results. In addition, the work presented in [8] shows a complete stability analysis, obtaining stability regions of each of the proposed methods.

Similar studies can be found for the wave equation. For instance, a proper selection for Δ*t* is presented in [9] as

where *h* is the grid size, *a*_{1} is the sum of the absolute values of weights for the finite difference operators for ∂^{2}*u*/∂*t*^{2}, and *a*_{2} is the sum of the absolute values of weights for the finite difference operators for ∇^{2}*u*. This bound ensures the production of stable results for the wave equation; however, as can be seen, it is limited to regular discretizations over a rectangle. Another remarkable example of the application of finite differences to obtain the numerical solution of the wave equation can be found in [10], where some CTCS schemes are presented for solving the 1+1*D* case, with great performance.

For the case of interest in this work, irregular 2D spatial domains (see Figure 1) are considered. Due to this, the classical bounds cannot be applied since the spatial discretizations are no longer a regular set of points on a grid.

In this context, a generalized finite difference scheme is proposed in this paper, building upon the ideas presented in [11]. However, despite the utility of generalized finite differences, only some advances in stability analysis exist for these schemes. In previous works such as [12–15], the authors present stability analyses for generalized finite difference schemes applied to a modified Lax-Wendroff scheme for advection equation, and the pure advection, diffusion, and advection-diffusion equations. Given the importance of establishing stability conditions for the wave equation, this paper presents a Von-Neumann stability analysis for the proposed scheme.

## 2. Generalized finite differences applied to wave equation

The stability analysis presented in the following section arises from the problem of obtaining a finite-difference scheme for the wave equation problem over a simply connected planar domain Ω with a positively oriented Jordan polygon as its boundary. The problem is defined as:

In order to address this problem, first, the equation is partially discretized in time, followed by the discretization of the spatial derivatives. Both discretizations are presented in the following subsections.

### 2.1. Temporal discretization

The stability of the scheme for the wave equation relies heavily on the discretization of the second-order partial time derivative. To achieve this, it is possible to follow the steps provided below.

1. First, the time derivative can be expressed as its finite difference approximation. This is,

where ${u}_{i}^{k}$ stands for the approximation of *u* at the grid node *p*_{i}, at a time level *k*.

2. For an arbitrary node on the cloud, *p*_{0} = (*x*_{0}, *y*_{0}), the equation can be written as

or, solving for ${u}_{0}^{k+1}$,

As seen, this approximation requires two different time steps to be computed. Due to this, it is essential to properly approximate the second time step.

3. It is possible to consider the condition for the first partial derivative in time

and apply central finite differences to get

now, solving for ${u}_{0}^{-1}$, it can be rewritten as

4. Let us assume that (4) holds at *k* = 0, this is,

therefore, it is possible to replace (5) into (6) to obtain,

With this, the discretizations (4) and (7) can compute all the time steps required for the scheme.

### 2.2. Spatial discretization

Next, the spatial discretization is carried out for the partially discretized equation in time, which is accomplished using a generalized finite differences scheme.

To begin, a second-order linear operator

where *A*, *B*, *C*, *D*, *E*, and *F* are spatial functions specified by the operator, is considered. At a central node *p*_{0} = (*x*_{0}, *y*_{0}), this operator can be approximated for a cloud points distribution, as illustrated in Figure 2, by using the values of *u* at nearby nodes *p*_{i} = (*x*_{i}, *y*_{i}), *i* = 1, 2, …, *q*, that are sufficiently close, which will be further discussed later in the manuscript.

A finite-difference scheme at the node *p*_{0} can then be expressed as a linear combination

where Γ_{0}, Γ_{1}, … , Γ_{q} are appropriate weights.

In [16, 17], it is stated that a finite difference scheme *L*_{0} is consistent with the linear operator *L* if the local truncation error τ satisfies

when *p*_{1}, *p*_{2}, … , *p*_{q}→*p*_{0}.

The Taylor's series expansion of the consistency condition, up to the second order, yields

in this case, Δ*x*_{i} = *x*_{i}−*x*_{0} and Δ*y*_{i} = *y*_{i}−*y*_{0}. Now, according to the consistency condition (8), it is required that

These conditions define the linear system

An important remark must be made, system (9) must be full-row-ranked to produce suitable approximations of the waves, and therefore, this fact must be taken into account for choosing the neighbors in the clouds.

Now, the scheme defined by (9) can be used to approximate the linear operator,

at each time level, taking *A* = *C* = (*cΔt*)^{2}, and *B* = *D* = *E* = *F* = 0. The resulting Γ_{i} coefficients, with the time discretizations (4) and (7) define the Generalized Finite-Difference Schemes

and

for the wave equation, where *k* represent the time level and ${u}_{i}^{k}$ is the approximation to the solution in the point *p*_{i} = (*x*_{i}, *y*_{i}) at time *kΔt*.

Equation (11) is the algebraic representation of an explicit scheme, following the idea presented on [18]; it is possible to extend this scheme, adding a parameter, λ, to involve two different time levels:

where λ can take values in [0, 1]. Particularly, for λ = 0.5, the scheme is a Crank-Nicolson-like scheme.

As mentioned earlier, a critical issue for this method is the number of neighbor nodes *q* used in the scheme. Extensive tests have shown that considering all the nodes in a space lower or equal to δ, where δ is the average node distance along the boundary, an average of 6 neighbors are chosen for each approximation. This is a significant reduction in the number of nodes compared with other schemes in the literature [19, 20]. Figure 3 shows an example of this selection.

## 3. Stability analysis

Once the scheme is defined, a question that could naturally arise is whether the scheme produces a stable approximation since this is crucial to achieve convergence. Sousa illustrates in [21] the difficulties in stability analysis when the schemes become more complex, as is the case in clouds of points; taking this into account, in this section, analyses for the proposed schemes are presented to get theoretical bounds for the stability of the methods.

For the case of the explicit scheme (11),

following the Von-Neumann's stability analysis [22, 23], a generic component of the error at an arbitrary node *p*_{0} = (*x*_{0}, *y*_{0}) at the *k*-th time-step could be measured as,

where φ represents the amplification factor of the error. Since ${\Phi}_{0}^{k}$ satisfies the finite difference scheme, then,

so, in this case, φ must satisfy

Then, the amplification factor satisfies

By applying the consistency conditions up to the second order, it is possible to write

or

with α = 2−Δ*t*^{2}*c*^{2}(*r*^{2}+*s*^{2}). Let us notice in this expression that α ≤ 2, and let us choose −2 ≤ α ≤ 2.

The solution of the quadratic equation can be expressed as

in this case, since α^{2} ≤ 4, then,

and then, the scheme is * conditionally stable*.

Let us consider the case α^{2} = 4; for this value, we get the bound

Following the same idea, for the implicit scheme (12),

where the error can be written as

Then, the amplification factor satisfies

Once again, by applying the consistency conditions up to the second order,

now, considering

it is possible to write

this is,

and, solving for φ,

now, considering λ^{2}β ≤ 4, and using ∥·∥^{2},

With this, in general, the implicit scheme turns out to be * stable* and, for the case λ = 0, it turns out to be

*.*

**unconditionally stable**## 4. Numerical tests

To show the performance of the explicit and the implicit proposed schemes, four different regions were selected for the numerical tests: The unitary square [0, 1] × [0, 1], denoted as CUA, and three irregular, non-rectangular, non-symmetrical geometries designated as ENG, HAB, and PAT. Arbitrary clouds of points were generated for each region using an algorithm designed by Persson and Strang and presented in [24]. In Figure 4, it is possible to see one of the generated clouds of points for each region.

Two examples were solved for each region to show the method's accuracy. In both examples, the time interval [0, 1] was subdivided into 1, 000 subintervals, i.e., Δ*t* = 0.001, to assure stability in most tests with the explicit scheme. Then, the norm of the quadratic error was computed as

where ${U}_{i}^{k}$ and ${u}_{i}^{k}$ are the approximated and theoretical solutions, respectively, at the *i*-th node and at a time *k*; and *n* is the total number of cloud nodes. Numerical solutions were computed for the explicit scheme (λ = 1) and implicit scheme with λ = 0.75, 0.50, 0.25, 0.00.

### 4.1. Example 1

Considering, as presented in [25], the equation

with initial conditions given by

and the following boundary condition

### 4.2. Example 2

For this second test, following the idea on [11], the equation

is considered, with initial conditions given by

and the boundary condition

Tables 1, 2 depict the error norms for the tests. Two different values of the error norm are considered in each case: in the first set of results, the mean of the error norm, over all time steps, is presented in order to show the overall performance of the schemes; in the second set of results, the error norm for the last time step is presented.

For the case of the implicit scheme, with λ = 0.75, Figures 5, 6 present a comparison between the numerical solution (on the left) and the theoretical solution (on the right) of the tests, in three different time levels for PAT region. It is possible to notice that both the approximated and the theoretical solutions present the same behavior.

The empirical convergence order of the method was estimated both for time and space. For the spatial case, the means of the quadratic error were calculated for clouds with a different number of nodes, for λ = 0. Denoting *m* as the total amount of nodes, Δ_{m} = 1/*m*, and *e*_{Δm} as the mean of the error computed then, given two pairs (*e*_{Δi}, Δ_{i}) and (*e*_{Δj}, Δ_{j}), the spatial empirical convergence order *p*_{s} was computed as

similarly, for the temporal convergence order, *p*_{t} , the computations were performed, denoting *e*_{Δt} as the error computed for the cloud using a Δ_{t} value. The results for these computations are reported in Table 3.

The method implementation was developed in Python, and executed with Anaconda in Visual Studio Code. All tests were performed on an iMac (21.5-inch, Late 2013) with a 2.7 GHz quad-core Intel Core i5 processor and 8 GB of 1,600 MHz DDR3 RAM. The execution times, for all cases, were measured using the *time* library in Python, measuring the actual execution time instead of processor time. In the most time-consuming case, the execution took 3.015 s, while in the fastest case, it took 0.241 s.

Due to its stability, the implicit scheme can be applied to compute results for problems where there is no theoretical solution, even in non-simply-connected domains. In Figures 7, 8, it is possible to see the propagation of a wave, over CUA and HAB regions, with reflexive boundary conditions and a hole inside the region.

## 5. Discussion and future work

In this work, two generalized finite differences schemes were presented to numerically solve the wave equation on highly irregular regions. Both schemes arise from the Taylor series expansion of the consistency condition introduced in [16, 17] for all finite difference schemes. The classical Lax-Richtmyer Equivalence Theorem [26] states that,

since the proposed schemes arise from a consistency condition, stability is an important issue to study for a convergent scheme.

The presented Von-Neumann's stability analyses show that it is possible to obtain stability for both schemes; this agrees with the numerical results. For the case of the explicit scheme that turned out to be * conditionally stable*, it is possible to use the presented bound for Δ

*t*to get stable results. On the other hand, the implicit scheme proved to be

*for the case of λ = 0, and no bound for Δ*

**unconditionally stable***t*is required to compute stable results; even more, extensive tests have shown that with relatively big time-steps, it is possible to obtain stable results. Furthermore, the stability analyses were performed without considering any particular data structures. Due to this, they can be applied to all and every up-to-second-order finite difference scheme, not only to the ones presented in this work.

On the other hand, the numerical results show that the proposed schemes produce stable and accurate approximations of the exact solution of the wave equation, even in highly-irregular ones, such as HAB and PAT. The computed norms of the quadratic errors presented in this work show that the proposed schemes could be effective numerical methods to model different real-life scenarios involving the wave equation. And even more, the performed tests show that the results obtained are accurate without considering any weight functions, as in other methods.

A crucial issue to be remarked on is that the irregularity of the region needs to be considered since the number of nodes on the boundary needs to be enough to reproduce the area that is being studied correctly; due to this, a proper selection of the number of neighbor points and a correct choice of these might affect the accuracy when the region boundaries are highly irregular, and this has to be studied in the future.

Finally, it is worth highlighting the necessity of exploring the adaptability of the proposed methods to different physical problems governed by partial differential equations, such as shallow-water problems, where Robin boundary conditions are required to describe real-life phenomena.

## Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

## Author contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

## Funding

Universidad Michoacana de San Nicolás de Hidalgo and CONACyT provide funds for open-access publication fees, and the International Center for Numerical Methods in Engineering provides technical and academic support.

## Acknowledgments

The authors thank AULA CIMNE-Morelia, CIC UMSNH, and CONACyT for supporting this research.

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## References

1. Barrera-Sánchez P, Castellanos-Noda L, Domínguez-Mota FJ, González-Flores G, Pérez-Domínguez A. Adaptive discrete harmonic grid generation. *Math Comput Simul.* (2009) 79:1792–809. doi: 10.1016/j.matcom.2007.04.015

2. Domínguez-Mota FJ, Mendoza-Armenta S, Tinoco-Guerrero G, Tinoco-Ruiz JG. Numerical solution of poisson-Like equations with Robin boundary conditions using a finite difference scheme defined by an optimality condition. In: *IMACS Series in Computational and Applied Mathematics: MASCOT 2011*. Rome. (2012).

3. Tinoco-Ruiz JG, Domínguez-Mota FJ, Mendoza-Armenta S, Tinoco-Guerrero G. Solving stokes equation in plane irregular region using an optimal consistent finite difference scheme. In: *20th International Meshing Roundtable*. Paris. (2011).

4. Tinoco-Ruiz JG, Domínguez-Mota FJ, Tinoco-Guerrero G, Fernández-Valdés PM, Mendoza-Armenta S. Numerical solution of differential equations in irregular plane regions using quality structured convex grids. *Int J Model Simul Sci Comput.* (2013) 4:1350004. doi: 10.1142/S1793962313500049

5. Alcrudo F, García-Navarro P. A high-resolution Godunov-type scheme in finite volumes for the 2D shallow-water equations solution Godunov-type scheme in finite volumes for the 2D shallow-water equations. *Int J Num Methods Fluids.* (1993) 16:489–505.

6. Appadu AR. Numerical solution of the 1D advection-diffusion equation using standard and nonstandard finite difference schemes. *J Appl Math.* (2013) 2013:734374. doi: 10.1155/2013/734374

7. Appadu AR. Optimized composite finite difference schemes for atmospheric flow modeling. *Num Methods Partial Diff Equat.* (2019) 35:2171–92. doi: 10.1002/num.22407

8. Appadu AR, Djoko JK, Gidey HH. A computational study of three numerical methods for some advection-diffusion problems. *Appl Math Comput.* (2016) 272:629–47. doi: 10.1016/j.amc.2015.03.101

9. Lines LR, Slawinski R, Bording RP. A recipe for stability of finite-difference wave-equation computations. *Geophysics.* (1999) 64:967–9.

10. Appadu AR, Waal GN. CTCS schemes for second order wave equation: numerical results and spectral analysis. *Int J Eng Res Afr.* (2021) 55:47–65. doi: 10.4028/www.scientific.net/JERA.55.47

11. Tinoco-Guerrero G, Arias-Rojas H, Guzmán-Torres JA, Roman-Gutiérrez R, Tinoco-Ruiz JG. A meshless finite difference scheme applied to the numerical solution of wave equation in highly irregular space regions. *Comput Math Appl.* (2023) 126:25–33. doi: 10.1016/j.camwa.2023.01.035

12. Tinoco-Guerrero G, Domínguez-Mota FJ, Gaona-Arias A, Ruiz-Zavala ML, Tinoco-Ruiz JG. A stability analysis for a generalized finite-difference scheme applied to the pure advection equation. *Math Comput Simul*. (2018) 147:293–300. doi: 10.1016/j.matcom.2017.06.001

13. Tinoco-Guerrero G, Domínguez-Mota FJ, Tinoco-Ruiz JG. Stability aspects of a modified Lax-Wendroff scheme for irregular 2D regions. In: *IMACS Series in Computational and Applied Mathematics: MASCOT 2015*. Rome. (2015). p. 161–70.

14. Tinoco-Guerrero G, Domínguez-Mota FJ, Tinoco-Ruiz JG. A study of the stability for a generalized finite-difference scheme applied to the advection-diffusion equation. *Math Comput Simul.* (2020) 176:301–11. doi: 10.1016/j.matcom.2020.01.020

15. Tinoco-Ruiz JG, Domínguez-Mota FJ, Lucas-Martínez JS, Tinoco-Guerrero G. Study of the stability of a generalized finite difference scheme applied to the diffusion equation in irregular 2D space regions using convex grids. In: *IMACS Series in Computational and Applied Mathematics: MASCOT 2018*. Rome. (2019).

16. Strikwerda JC.*Finite Difference Schemes and Partial Differential Equations*. Philadelphia, PA: Society for Industrial and Applied Mathematics (2004).

17. Thomas JW. *Numerical Partial Differential Equations: Finite Difference Methods. vol. 2 of Texts in Applied Mathematics*. New York, NY: Springer (1995).

18. Tinoco-Guerrero G, Domínguez-Mota FJ, Tinoco-Ruiz JG, Ruiz-Dí|az E. An implicit modified Lax-Wendroff scheme for irregular space 2D regions. In: *IMACS Series in Computational and Applied Mathematics: MASCOT 2013*. El Escorial. (2013). p. 18.

19. Huang J, Fan CM, Chen JH, Yan J. Meshless generalized finite difference method for the propagation of nonlinear water waves under complex wave conditions. *Mathematics.* (2022) 10. doi: 10.3390/math10061007

20. Salete E, Vargas AM, García A, Benito JJ, Ureña F, Ureña M. An effective numeric method for different formulations of the elastic wave propagation problem in isotropic medium. *Appl Math Modell.* (2021) 96:480–96. doi: 10.1016/j.apm.2021.03.015

21. Sousa E. The controversial stability analysis. *Appl Math Comput.* (2003) 145:777–94. doi: 10.1016/S0096-3003(03)00274-1

22. Grebenkov DS, Nguyen BT. Geometrical structure of Laplacian eigenfunctions. *SIAM Rev.* (2013) 55:601–67. doi: 10.1137/120880173

23. Strauss WA. Chapter 8: Computations of Solutions. In:Rosatone L, Corliss S, , editors. *Partial Differential Equations: An Introduction.* 2nd ed. John Wiley & Sons (2009). p. 199–227.

24. Persson PO, Strang G. A simple mesh generator in Matlab. *SIAM Rev.* (2004) 46:329–45. doi: 10.1137/S0036144503429121

25. Flores J, García A, Negreanu M, Salete E, Ureña F, Vargas AM. A spatio-temporal fully meshless method for hyperbolic PDEs. *J Comput Appl Math.* (2023) 2023:115194. doi: 10.1016/j.cam.2023.115194

Keywords: stability analysis, meshless method, generalized finite difference, wave equation, numerical solution of PDE

Citation: Tinoco-Guerrero G, Domínguez-Mota FJ, Guzmán-Torres JA, Román-Gutiérrez R and Tinoco-Ruiz JG (2023) Study of the stability of a meshless generalized finite difference scheme applied to the wave equation. *Front. Appl. Math. Stat.* 9:1214022. doi: 10.3389/fams.2023.1214022

Received: 28 April 2023; Accepted: 26 June 2023;

Published: 11 July 2023.

Edited by:

Edgar O. Reséndiz-Flores, Tecnológico Nacional de México/IT Saltillo, MexicoReviewed by:

Bishnu Lamichhane, The University of Newcastle, AustraliaAppanah Rao Appadu, Nelson Mandela University, South Africa

Jorge Eliecer Ospino Portillo, Universidad del Norte, Colombia, Colombia

Copyright © 2023 Tinoco-Guerrero, Domínguez-Mota, Guzmán-Torres, Román-Gutiérrez and Tinoco-Ruiz. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Gerardo Tinoco-Guerrero, gerardo.tinoco@umich.mx

^{†}These authors have contributed equally to this work and share last authorship