Simulations and analysis of underwater acoustic wave propagation model based on Helmholtz equation

Abstract

In this study, we explore the effectiveness of elliptic partial differential equations (PDEs) in two and three dimensional space based on Helmholtz equation for simulations of acoustic sound in a very complex environment of propagation. For this purpose, we use an advance and robust numerical technique by utilizing the properties of shifted Chebyshev spectral collocation method. This technique is an extension of the traditional Chebyshev polynomials, incorporating a shift in their argument to enhance flexibility across a wider domain, while retaining an extraordinary numerical characteristic such as orthogonality and spectral convergence making them exceptionally effective in finding the approximate solutions. The exponential order of convergence of the proposed approach is shown both through theoretical and numerical approaches. We provide a number of numerical experiments to verify the theoretical results. The spectral convergence has been substantially enhanced by these numerical examples. The exponential order is further validated by numerical error behaviour in both $L_2$ and $L_{\mathrm{\infty }}$ norms.

Keywords:

• Acoustic wave propagation model
• Helmholtz equation
• shifted Chebyshev polynomials
• error analysis
• numerical simulations

Mathematics Subject Classifications:

• 41A10
• 39A14
• 65M70

1. Introduction

Understanding the complex dynamics of underwater acoustic propagation is crucial for various applications, including but not limited to sonar systems and environmental monitoring. It also plays a significant role in understanding and predicting the behaviour of sound waves in the ocean environment. Further ocean acoustic underwater propagation models have been extensively used in many applications over several decades, including but not limited to anti-submarine warfare, global-scale underwater sound propagation, characterization of acoustic communication channels, ocean acoustic tomography and various other oceanographic applications. In early 1950, a considerable efforts were made in comprehending and forecasting sound propagation in the ocean, resulting in the accumulation of substantial knowledge concerning underwater sound propagation and its dependencies on the environment. Mathematical models and their numerical solutions play a crucial role in the study and simulation of acoustic wave propagation. Acoustic wave propagation refers to the transmission of sound waves through various media, such as air, water, or solid materials. Understanding and predicting how sound waves behave in different environments is essential in various fields, including engineering, physics, environmental science. Recent advancements across various scientific disciplines have enhanced our confidence in predicting the acoustic field for a specific source. Much of this progress is due to the developments in computational acoustics, which offered solutions to the model equation, enabling the computation of the acoustic field for a given environmental context [1].

Exploring wave propagation in underwater acoustics has consistently been a prominent subject. Numerical simulation serves as an effective tool for examining the propagation of sound waves within complex media. Understanding the propagation characteristics of sound waves in diverse ocean environments constitutes a primary focus in the analysis of ocean sound fields. The complex distribution of sound waves in the ocean is attributed to the influence of several factors including the water surface, the aquatic environment, and the seafloor [2]. The principles of ocean acoustic propagation dictate that sound waves at a fixed frequency conform to the elliptic partial differential Helmholtz equation. Over the past few decades, the computation methods for analysing sound fields in acoustic media have seen significant advancements. Particularly in underwater acoustics, a range of techniques has been effectively used, as extensively detailed in a notable reference [3]. Each of these methods, however, comes with its own set of constraints. Well-established approaches like acoustic ray theory, normal mode analysis, simplified parabolic equations, and Green's functions for layered media stand out for their widespread application. The evolution of more precise wave theories, combined with the advent of high-speed computing technologies and advances in numerical physics, has led to the emergence of more accurate and varied solutions for specific underwater acoustic challenges [4–6].

Pioneering work by Dawson and Fawcett [7] involved using the Boundary Element Method for studying acoustic scattering in the open ocean, particularly over flat waveguide surfaces with localized areas of deformation. Grilli et al. subsequently introduced a hybrid model in Boundary Element Method (BEM) applications, combining the conventional method in an inner region with variable bathymetry and an eigenfunction expansion in an outer region of constant depth [8]. Santiago and Wrobel applied the sub-region technique within a boundary element framework to analyse two-dimensional acoustic wave propagation in shallow waters with irregular seabed topography [8,9]. Implementing boundary conditions is crucial in confining the Helmholtz equation to accurately determine sound pressure across the entire spatial domain [10]. In recent developments, Wang et al. effectively used a one-dimensional spectral method for solving the parabolic equation model [11].

The spectral method, known for its high accuracy and rapid convergence, has significantly advanced the field of ocean acoustic calculations. Nevertheless, challenges persist in the computational process. One such challenge arises from the initial simplification of the Helmholtz equation into a parabolic model, facilitating numerical solutions but severely limiting its applicability and introducing approximation errors. The parabolic equation model is suited primarily for scenarios with slow horizontal variations in physical parameters and an adoption of the far-field assumption [12–14]. This approach ignores the contribution of reverse acoustic waves, limiting the solution to one-way acoustic propagation [15]. The normal mode model, on the other hand, is suitable only when physical parameters remain constant and cannot handle scenarios with variable parameters over distance, as seen in actual ocean settings [16]. The ray model applies to high-frequency situations, where sound wave wavelengths are considerably smaller than any relevant physical scale [17]. The mirror image method relies on superimposing real and mirrored sound source points to obtain the actual sound field, but it demands stringent requirements on sound speed, limiting its use to constant sound speed profiles in ocean environments. In [18], a wave number integration method is only appropriate for horizontally layered media and employs numerical integration to determine the sound field [19]. A higher order compact difference scheme utilizing uniform mesh sizes across various coordinate directions, is applied for discretizing both two- and three-dimensional Helmholtz equations is used in [20]. Apart from this in science and engineering, partial differential equations in general and fractional differential equations in particular are essential mathematical tools used in describing a wide range of phenomena, especially in mathematical modelling of infectious diseases in biological sciences [21–26].

Spectral methods used a global approach for approximating unknown functions, where computations at any point depend on the entire domain and neighbouring points. This is in contrast to local approaches, where computations rely solely on neighbouring points. These methods are commonly used for discretizing special coordinates, and their formulation is heavily reliant on approximation and test functions, often referred to as trial and weight functions, respectively. In this approach, trial functions are combined linearly to approximate an unknown function, thereby providing an estimated solution for differential equations. The trial functions' truncated series expansion is utilized to ensure the differential equations meet specific conditions, such as initial or boundary conditions. The goal is to minimize the error between the exact solution and the solution derived from the truncated series. This is achieved by minimizing a suitable norm, which involves implicit or explicit considerations of the residuals. Spectral methods, often viewed as a specialized form of the method of residuals, are distinguished by the residuals satisfying orthogonality conditions with certain test functions. Early spectral methods are set apart from finite element and finite difference methods by their choice of trial functions. These methods, classified into Galerkin, Tau, and collocation methods, differ based on the test functions used. Fourier spectral methods are preferred for periodic functions, while orthogonal polynomials like Legendre or Chebyshev polynomials are ideal for non-periodic domains.

The primary objective of this paper is to develop an effective numerical scheme for two-dimensional elliptic partial differential equations (PDEs) based on the Helmholtz equation, by employing a shifted Chebyshev spectral method with appropriate boundary conditions. The choice of using shifted Chebyshev polynomials in the spectral method is driven by their high accuracy and minimal phase error in solving the specified problem. A significant advantage of spectral methods is their exponentially decaying error, implying an infinite convergence rate in space even with a minimal number of collocation points. These methods offer flexibility in choosing basis functions, which can be selected based on the problem at hand. However, their implementation can be challenging, mainly because they use global functions as basis functions. This characteristic makes them less suitable for handling local features and sharp gradients, as seen in phenomena like Gibbs one [27–41].

The remaining structure of the paper includes method and mathematical formulation in Section 2. Some preliminaries are introduced in Section 3 needed for the analysis for the proposed scheme which is presented in Section 4. Numerical examples are presented in Section 5 to confirm the exponential convergence, while a conclusion is drawn in Section 6.

2. Preliminaries of shifted Chebyshev polynomials

Chebyshev polynomials and their orthogonal properties over the interval $[-1,1]$ make them a most suitable candidate to solve differential equations arising in the mathematical modelling of many real-world problems. They are often used as basis functions in spectral methods while solving PDEs, which is the most powerful numerical method in terms of their accuracy and efficiency as compared to other classical methods like finite difference and final element methods with a very fever grids. Depending on the specific problem, the range $[0,1]$ is often more convenient as compare to the range $[-1,1].$ The modified version of Chebyshev polynomials is more useful while dealing with the functions that have discontinuity or singularities at x = 0, or x = 1, or at both points. It is also a natural choice if the problem has Dirichlet boundary conditions which are given at $x=0,$ and x = 1. In a broader context, we can establish Chebyshev polynomials which are appropriate for any specified finite interval $[u,v]$ for the variable ζ. This is achieved by aligning the chosen interval with the range $[-1,1]$ through a linear transformation involving a new variable ζ [42]. (1) $\zeta =\frac{2x-(u+v)}{v-u}.$(1) We are more interested in the special case when $[u,v]=[0,1]$. Then the linear transformation given in Equation (1(1) $\zeta =\frac{2x-(u+v)}{v-u}.$(1) ) becomes: $\zeta =2x-1.$Thus we obtain the modified Chebyshev polynomial which is denoted by ${{\mathcal{T}}^{\mathcal{\ast }}}_n(\zeta )$ and stratifying the recursive relation of the form: ${{\mathcal{T}}^{\mathcal{\ast }}}_n(\zeta )=2\zeta \cdot {{\mathcal{T}}^{\mathcal{\ast }}}_{n-1}(\zeta )-{{\mathcal{T}}^{\mathcal{\ast }}}_{n-2}(\zeta ),$where ${{\mathcal{T}}^{\mathcal{\ast }}}_0(\zeta )=1$ and ${{\mathcal{T}}^{\mathcal{\ast }}}_1(\zeta )=2\zeta -1$. The power form representation of ${{\mathcal{T}}^{\mathcal{\ast }}}_n(\zeta )$ is one of the most essential formulas and is given by: ${{\mathcal{T}}^{\mathcal{\ast }}}_n(\zeta )=n\sum _{k=0}^n(-1{)}^{n-k}\frac{(n+k-1)!2^{2k}}{(n-k)!(2k)!}{\zeta }^k.$The modified shifted polynomials retain the orthogonality condition: ${\int }_0^1{T}_m^{\ast }(\zeta )T_n^{\ast }(\zeta )\frac{1}{\sqrt{\zeta -{\zeta }^2}}\mathrm{d}\zeta ={\delta }_{\mathrm{mn}}{C}_n,$where $C_n=h_n\frac{\pi }{2}$ with $h_0=2,h_1=1,n\ge 1$ and ${\delta }_{\mathrm{nm}}$ is the Kronecker delta, which is defined as: ${\delta }_{\mathrm{mn}}=\{\begin{array}{ll}1& \mathrm{if}m=n,\\ 0& \mathrm{if}n\ne m.\end{array},$The first few modified Chebyshev polynomials are: $\begin{array}{rl}{\mathcal{T}}_0^{\ast }(\zeta )& =1\\ {\mathcal{T}}_1^{\ast }(\zeta )& =2\zeta -1\\ {\mathcal{T}}_2^{\ast }(\zeta )& =8{\zeta }^2-8\zeta +1\\ {\mathcal{T}}_3^{\ast }(\zeta )& =32{\zeta }^3-48{\zeta }^2+18\zeta -1\\ {\mathcal{T}}_4^{\ast }(\zeta )& =128{\zeta }^4-256{\zeta }^3+160{\zeta }^2-32\zeta +1\\ {\mathcal{T}}_5^{\ast }(\zeta )& =512{\zeta }^5-1280{\zeta }^4+1120{\zeta }^3-400{\zeta }^2+50\zeta -1\\ {\mathcal{T}}_6^{\ast }(\zeta )& =2048{\zeta }^6-6144{\zeta }^5+6912{\zeta }^4-3584{\zeta }^3+840{\zeta }^2-72\zeta +1\\ {\mathcal{T}}_7^{\ast }(\zeta )& =8192{\zeta }^7-28672{\zeta }^6+36864{\zeta }^5-23296{\zeta }^4+7648{\zeta }^3-1312{\zeta }^2+98\zeta -1.\end{array}$

3. Mathematical model description and solution methodology

In various scientific and engineering applications, second-order PDEs, whether parabolic, elliptic, or hyperbolic, and whether linear or nonlinear, play a crucial role. Helmholtz equation, a fundamental PDEs, serve as a key mathematical model for understanding the dynamics of sound wave propagation which is a fascinating and critical area of study in various scientific and engineering disciplines, offering insights into the behaviour of sound waves as they traverse through different media, including scattering and radiation problems, muffler systems, car cabin acoustics etc. The time-independent Helmholtz equation is also used to model a harmonic sound pressure field at a specific angular frequency, where the dependent variable is the sound pressure and the sound pressure wave is propagating in a medium with density at the speed of sound. It can be can also be derived from the heat conduction equation, Schrödinger equation, telegraph and other wave-type, or evolutionary, equations. We consider the two-dimensional in-homogeneous acoustic underwater wave propagation mathematical model based on elliptic partial differential equations of Helmholtz type of the form: (2) $\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}+{\lambda }^2\mathrm{\Psi }(\eta ,\xi )=f(\eta ,\xi ).(\eta ,\xi )\in \mathrm{\Gamma }$(2) subject to general Dirichlet boundary conditions: (3) $\{\begin{array}{ll}\mathrm{\Psi }(\eta ,\alpha )=f_1(\eta ),& \mathrm{\Psi }(\eta ,\beta )=f_2(\eta ),\\ \mathrm{\Psi }(\gamma ,\xi )=f_3(\xi ),& \mathrm{\Psi }(\mu ,\xi )=f_4(\xi ).\end{array}$(3) or with general Neumann boundary conditions: (4) $\{\begin{array}{ll}{\displaystyle \frac{\mathrm{\partial \Psi }}{\mathrm{\partial }\eta }(\eta ,\alpha )=h_1(\eta ),}& {\displaystyle \frac{\mathrm{\partial \Psi }}{\mathrm{\partial }\eta }=h_2(\eta ),}\\ {\displaystyle \frac{\mathrm{\partial \Psi }}{\mathrm{\partial }\xi }=h_3(\xi ),}& {\displaystyle \frac{\mathrm{\partial \Psi }}{\mathrm{\partial }\xi }=h_4(\xi ).}\end{array}$(4) For ease of understanding, we assume a section of the ocean characterized by a flat surface and a flat bottom, resulting in the region Γ being shaped as a rectangle as shown in Figure 1.

Figure 1. Model problem geometry region.

In Equation (2(2) $\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}+{\lambda }^2\mathrm{\Psi }(\eta ,\xi )=f(\eta ,\xi ).(\eta ,\xi )\in \mathrm{\Gamma }$(2) ), $\mathrm{\Psi }(\eta ,\xi )$ is a function need to be workout. λ represents the acoustics wave number, $\mathrm{\Gamma }=[0,1]\times [0,1]$ denotes the convex symmetric domain in two dimension and $f(\eta ,\xi )$ is a source term, which is assumed to continues on Γ and $\alpha ,\beta ,\gamma ,\mu \in \mathrm{\Gamma }.$

Due to the oscillatory nature of Equation (2(2) $\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}+{\lambda }^2\mathrm{\Psi }(\eta ,\xi )=f(\eta ,\xi ).(\eta ,\xi )\in \mathrm{\Gamma }$(2) ), it is very hard to find the approximate solution especially when the wave number λ is very large. Also the variable wave number in the Helmholtz equation adds significant complexity, posing considerable challenges for numerical computation. The study of numerical solutions for the Helmholtz equation is of vital theoretical and practical importance. Various numerical methods have been employed to tackle this equation, including finite element, boundary element, finite difference, and meshless methods, where the finite element method is notably the most prevalent, yet its accuracy tends to diminish significantly as the wave number in the equation increases. Conversely, the spectral methods based on shifted Chebyshev polynomials are gaining popularity among researchers due to several benefits. These include straightforward discretization, ease in deriving formulas, higher computational accuracy, and strong convergence (exponential), when the unknown in the given equation is a smooth function.

To discretize Equation (2(2) $\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}+{\lambda }^2\mathrm{\Psi }(\eta ,\xi )=f(\eta ,\xi ).(\eta ,\xi )\in \mathrm{\Gamma }$(2) ), subject to the boundary condition given by Equation (3(3) $\{\begin{array}{ll}\mathrm{\Psi }(\eta ,\alpha )=f_1(\eta ),& \mathrm{\Psi }(\eta ,\beta )=f_2(\eta ),\\ \mathrm{\Psi }(\gamma ,\xi )=f_3(\xi ),& \mathrm{\Psi }(\mu ,\xi )=f_4(\xi ).\end{array}$(3) ) using shifted Chebyshev spectral method, first the square domain $[-1,1]\times [-1,1]$ is mapped into to the domain of interest $\mathrm{\Gamma }=[0,1]\times [0,1]$ and a mesh of $(N+1)\times (M+1)$ points is formulated, where the collocation points are given by: $\begin{array}{rl}{\eta }_i& =\frac{1}{2}(1-\cos \left(\frac{\mathrm{\pi i}}{N}\right)),i=0,1,\dots ,N.\\ {\xi }_j& =\frac{1}{2}(1-\cos \left(\frac{\mathrm{\pi j}}{M}\right)),j=0,1,\dots ,M.\end{array}$A global interpolant using Chebyshev polynomials, with order N along x-axis and order M along the y-axis is applied to interpolate the values at the mesh points. (5) $\mathrm{\Psi }(\eta ,\xi )=\sum _{n=0}^N\sum _{m=0}^M{a}_{\mathrm{nm}}{T}_n^{\ast }(\eta )T_m^{\ast }(\xi ).$(5) Derivatives at the collocation points are obtained by first differentiating the interpolant and then evaluating this differentiated polynomial at those points. There are two primary methods to compute these derivatives [43]. In this context, we express the differentiation process as the multiplication of a differentiation matrix with a vector, specifically for each line parallel to the y-axis, the derivative in the x-direction is calculated in the same manner. Kronecker product, which is a special case of tensor product for vectors is used to compute partial derivatives.

Given two matrices, $\mathcal{A}$ and $\mathcal{B}$, where $\mathcal{A}$ is of size $m\times n$ and $\mathcal{B}$ is of size $p\times q$, the Kronecker product $\mathcal{A}\otimes \mathcal{B}$ results in a matrix of size $\mathrm{mp}\times \mathrm{nq}$. The operation is performed by multiplying each element of A by the entire matrix B.

For matrices: $\begin{array}{rl}\mathcal{A}& =\left(\begin{array}{cccc}{\alpha }_{11}& {\alpha }_{12}& \cdots & {\alpha }_{1n}\\ {\alpha }_{21}& {\alpha }_{22}& \cdots & {\alpha }_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ {\alpha }_{m1}& {\alpha }_{m2}& \cdots & {\alpha }_{\mathrm{mn}}\end{array}\right)\\ \mathcal{B}& =\left(\begin{array}{cccc}{\beta }_{11}& {\beta }_{12}& \cdots & {\beta }_{1q}\\ {\beta }_{21}& {\beta }_{22}& \cdots & {\beta }_{2q}\\ ⋮& ⋮& \ddots & ⋮\\ {\beta }_{p1}& {\beta }_{p2}& \cdots & {\beta }_{\mathrm{pq}}\end{array}\right)\end{array}$The Kronecker product is: $\mathcal{A}\otimes \mathcal{B}=\left(\begin{array}{cccc}{\alpha }_{11}\mathcal{B}& {\alpha }_{12}\mathcal{B}& \cdots & {\alpha }_{1n}\mathcal{B}\\ {\alpha }_{21}\mathcal{B}& {\alpha }_{22}\mathcal{B}& \cdots & {\alpha }_{2n}\mathcal{B}\\ ⋮& ⋮& \ddots & ⋮\\ {\alpha }_{m1}\mathcal{B}& {\alpha }_{m2}\mathcal{B}& \cdots & {\alpha }_{\mathrm{mn}}\mathcal{B}\end{array}\right)$In this expanded form, each ${\alpha }_{\mathrm{ij}}\mathcal{B}$ is $\mathrm{\alpha p}\times q$ matrix obtained by multiplying the scalar $a_{\mathrm{ij}}$ from matrix $\mathcal{A}$ with every element of matrix $\mathcal{B}$. The partial derivative appearing in Equation  (2(2) $\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}+{\lambda }^2\mathrm{\Psi }(\eta ,\xi )=f(\eta ,\xi ).(\eta ,\xi )\in \mathrm{\Gamma }$(2) ) takes the form: $\begin{array}{rl}& \frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}=D^{\ast 2}\otimes I_n,\\ & \frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}=I_m\otimes D^{\ast 2},\end{array}$where $I_n$, $I_m$ are the identity matrix of order n and m $D^{\ast 2}$ is the square of the shifted Chebyshev differentiation matrix. Its elements are given by: (6) $D_{\mathrm{ij}}^{\ast }=\{\begin{array}{ll}{\displaystyle \frac{2N^2+1}{6},}& \mathrm{if}i=j=0\mathrm{or}i=j=N,\\ {\displaystyle -\frac{x_i}{2(1-x_i^2)},}& \mathrm{if}i=j\mathrm{and}i\ne 0,N,\\ {\displaystyle \frac{c}{x_i-x_j},}& \mathrm{if}i\ne j,\end{array}$(6) where $c=(-1{)}^{i+j}$.

The Helmholtz Equation (2(2) $\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}+{\lambda }^2\mathrm{\Psi }(\eta ,\xi )=f(\eta ,\xi ).(\eta ,\xi )\in \mathrm{\Gamma }$(2) ) after discretize at the collocation points, leads to a system of linear equations: (7) $\sum _{k=0}^n\sum _{l=0}^m(D^{\ast 2}\otimes I_n+I_m\otimes D^{\ast 2}+{\lambda }^2{I}_m\otimes I_n){\mathrm{\Psi }}_{\mathrm{kl}}=f_{\mathrm{nm}},$(7) The Dirichlet boundary conditions are applied by modifying the system of equations at the boundary points: $\begin{array}{rl}{\mathrm{\Psi }}_{i,0}& =f_1({\eta }_i),i=0,1,\dots ,N,\\ {\mathrm{\Psi }}_{i,N}& =f_2({\eta }_i),i=0,1,\dots ,N,\\ {\mathrm{\Psi }}_{0,j}& =f_3({\xi }_j),j=0,1,\dots ,M,\\ {\mathrm{\Psi }}_{M,j}& =f_4({\xi }_j),j=0,1,\dots ,M.\end{array}$Let $\begin{array}{rl}{\lambda }^2& =\mathrm{diag}({\lambda }_{11}^2,\dots ,{\lambda }_{1n}^2,{\lambda }_{21}^2,\dots {\lambda }_{2n}^2,\dots ,{\lambda }_{m1}^2,\dots ,{\lambda }_{\mathrm{mn}}^2),\\ \mathrm{\Psi }& =({\psi }_{11},\dots ,{\psi }_{1n},{\psi }_{21},\dots {\psi }_{2n},\dots ,{\psi }_{m1},\dots ,{\psi }_{\mathrm{mn}}),\\ F& =(f_{11},\dots ,f_{1n},f_{21},\dots ,f_{2n},\dots ,f_{m1},\dots ,f_{\mathrm{mn}}).\end{array}$Equation (7(7) $\sum _{k=0}^n\sum _{l=0}^m(D^{\ast 2}\otimes I_n+I_m\otimes D^{\ast 2}+{\lambda }^2{I}_m\otimes I_n){\mathrm{\Psi }}_{\mathrm{kl}}=f_{\mathrm{nm}},$(7) ) can then be written as: (8) $\mathcal{L}\mathrm{\Psi }=F,$(8) where $\mathcal{L}=(D^{\ast 2}\otimes I_n+I_m\otimes D^{\ast 2}+{\lambda }^2{I}_m\otimes I_n)$The linear system Equation (8(8) $\mathcal{L}\mathrm{\Psi }=F,$(8) ) is then solved to find the values of Ψ at the collocation points.

4. Error analysis

In this section, we first present some basic lemmas, which is essentials for the theoretical investigation of the model Equation (2(2) $\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}+{\lambda }^2\mathrm{\Psi }(\eta ,\xi )=f(\eta ,\xi ).(\eta ,\xi )\in \mathrm{\Gamma }$(2) ). These initial results are not only crucial for understanding the error analysis of the numerical scheme being considered but also play a role in the numerical discretization process [44–47]. It's crucial to understand that in assessing the convergence of a spectral method across three spatial dimensions, all spatial discretizations need to converge for the overall method to be considered convergent. This requirement stems from our approach of discretizing the three spatial dimensions independently. Therefore, our analysis of convergence will focus on the case of a single spatial dimension, bearing in mind that the convergence criteria must be met for both x-axis, y-axis and z-axis discretizations, provided that the unknown function is smooth. The degree of convergence is significantly influenced by the precision of our modified Chebyshev spectral differentiation.

Lemma 4.1

[45], Interpolation error estimation for Chebyshev polynomial

Let $\mathrm{\Psi }\in H^m[a,b]$ and denote ${\mathcal{I}}_N\mathrm{\Psi }$ the polynomial interpolation polynomial with respect to $(N+1)$-Gauss type points $\{x_k{\}}_{k=0}^N$. Then $\begin{array}{rl}& \Vert \mathrm{\Psi }-{\mathcal{I}}_N\mathrm{\Psi }{\Vert }_{L^2(I)}\le C{N}^{-m}\Vert \mathrm{\Psi }{\Vert }_{{\tilde{H}}_{m,N}(I)},\\ & \Vert \mathrm{\Psi }-{\mathcal{I}}_N\mathrm{\Psi }{\Vert }_{L^{\mathrm{\infty }}(I)}\le C{N}^{1/2-m}\Vert \mathrm{\Psi }{\Vert }_{{\tilde{H}}_{m,N}(I)}.\end{array}$

Lemma 4.2

[45], Truncation error estimation for Chebyshev polynomial

Assume that ${\mathcal{P}}_N\mathrm{\Psi }$ is $L_w^2(I)$ orthogonal projection onto the space of Chebyshev polynomials with degree at most $(N+1)$ and C then for positive constant C and for any function $\mathrm{\Psi }\in C(I)\cap H_m^w(I)$, the following error estimates hold $\begin{array}{rl}& \Vert \mathrm{\Psi }-{\mathcal{P}}_N\mathrm{\Psi }{\Vert }_{L^2(I)}\le C{N}^{-m}\Vert \mathrm{\Psi }{\Vert }_{{\tilde{H}}_{m,N}(I)},\\ & \Vert \mathrm{\Psi }-{\mathcal{P}}_N\mathrm{\Psi }{\Vert }_{L^{\mathrm{\infty }}(I)}\le C{N}^{-m}(1+\ln N)\Vert \mathrm{\Psi }{\Vert }_{{\tilde{H}}_{m,N}(I)}.\end{array}$

Lemma 4.3

[45], Exponential convergence of Chebyshev spectral method

Let w is the mth spectral derivative of Chebyshev polynomials with respect to Ψ $(m\ge 1)$ then we have the following error estimate of exponential order of convergence since the solution Ψ is analytic. $\Vert w_j-{\mathrm{\Psi }}^{(m)}(x_j)\Vert =O(K^{-N}),$

Lemma 4.4

[48]

Assume that ${\mathcal{F}}_i(t)$ be the ith Lagrange interpolation polynomial with the $(N+1)$-point Gauss Chebyshev, Gauss-Radau Chebyshev or Gauss-Lobatto Chebyshev points $\{x_k{\}}_{k=0}^N$. Then $\underset{t\in I}{max}\sum _{i=0}^N|{\mathcal{F}}_i(t)|\le C{N}^{\frac{1}{2}}.$

Using the above Lemmas 4.1–4.4, one can easily state and prove the following main result.

Theorem 4.5

Assume that $\mathrm{\Psi }\in H_w^p[0,1]$ with $p>1/2$ and 0<p<m is the exact solution of the model Equation (2(2) $\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}+{\lambda }^2\mathrm{\Psi }(\eta ,\xi )=f(\eta ,\xi ).(\eta ,\xi )\in \mathrm{\Gamma }$(2) ) and ${\mathrm{\Psi }}_N$ is the approximate solution which is obtained using Equation (7(7) $\sum _{k=0}^n\sum _{l=0}^m(D^{\ast 2}\otimes I_n+I_m\otimes D^{\ast 2}+{\lambda }^2{I}_m\otimes I_n){\mathrm{\Psi }}_{\mathrm{kl}}=f_{\mathrm{nm}},$(7) ), then there exist a positive constant ${\mathcal{C}}^{\prime }$ and ${\mathcal{C}}^{″}$ independent of N, such that $\Vert \mathrm{\Psi }-{\mathrm{\Psi }}_N{\Vert }_{H_w^p[0,1]}\le {\mathcal{C}}^{\prime }{N}^{2p-m}\Vert \psi {\Vert }_{H_w^p[0,1]}+{\mathcal{C}}^{″}{N}^{-m}\Vert \psi {\Vert }_{H_w^p[0,1]}.$

5. Numerical illustrations

In this section, we illustrate three numerical examples, two of them for 2-dimensional case, while one for 3-dimensional case for the confirmation of spectral rate of convergence for the proposed numerical method.

5.1. Example 1

(9) $\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}+{\lambda }^2\mathrm{\Psi }(\eta ,\xi )=({\lambda }^2-2{\pi }^2)\sin (\mathrm{\pi \eta })\sin (\mathrm{\pi \xi })$(9) subject to the boundary conditions of pure Dirichlet form: (10) $\{\begin{array}{l}\mathrm{\Psi }(0,\xi )=0,\mathrm{\Psi }(1,\xi )=0,\\ \mathrm{\Psi }(\eta ,0)=0,\mathrm{\Psi }(\eta ,1)=0.\end{array}$(10) The exact solution of Equation (9(9) $\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}+{\lambda }^2\mathrm{\Psi }(\eta ,\xi )=({\lambda }^2-2{\pi }^2)\sin (\mathrm{\pi \eta })\sin (\mathrm{\pi \xi })$(9) ) is given by $\mathrm{\Psi }(\eta ,\xi )=\sin (\mathrm{\pi \eta })\sin (\mathrm{\pi \xi })$. The error behaviour for $L_2$ and $L_{\mathrm{\infty }}$ norm for N = 20 and N = 30 collocation points is shown in Figure 2. One can see that the error decays exponentially in both norms, which further authenticates the result of Theorem 4.5. Comparison between the exact solution and approximate solution for different collocation points and wavelength λ are given in Figures 3–6.

Figure 2. Example 1: Error behaviour in $L_2$ and $L_{\mathrm{\infty }}$ norm for N = 20 (left) and N = 30 (right).

Figure 3. Example 1: Numerical solution (left) vs exact solution (right) at N = 15 $\lambda =50$.

Figure 4. Example 1: Numerical solution (left) vs exact solution (right) at N = 30 $\lambda ={10}^2$.

Figure 5. Example 1: Numerical solution (left) vs exact solution (right) at N = 50 $\lambda ={10}^3$.

Figure 6. Example 1: Numerical solution (left) vs exact solution (right) at N = 100 $\lambda ={10}^4$.

5.2. Example 2

Consider the following example of Helmholtz equation: (11) $\begin{array}{rl}\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}+{\lambda }^2\mathrm{\Psi }(\eta ,\xi )& =-2[(1-6{\eta }^2)({\xi }^2-{\xi }^4)+(1-6{\xi }^2)({\eta }^2-{\eta }^4)\\ & +{\lambda }^2[({\eta }^2-{\eta }^4)({\xi }^4-{\xi }^2)]\end{array}$(11) subject to the same boundary condition given in Equation (10(10) $\{\begin{array}{l}\mathrm{\Psi }(0,\xi )=0,\mathrm{\Psi }(1,\xi )=0,\\ \mathrm{\Psi }(\eta ,0)=0,\mathrm{\Psi }(\eta ,1)=0.\end{array}$(10) ) on $[0,1]\times [0,1]$:

The exact solution of Equation (11(11) $\begin{array}{rl}\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}+{\lambda }^2\mathrm{\Psi }(\eta ,\xi )& =-2[(1-6{\eta }^2)({\xi }^2-{\xi }^4)+(1-6{\xi }^2)({\eta }^2-{\eta }^4)\\ & +{\lambda }^2[({\eta }^2-{\eta }^4)({\xi }^4-{\xi }^2)]\end{array}$(11) ) is given by $\mathrm{\Psi }(\eta ,\xi )=({\eta }^2-{\eta }^4)({\xi }^4-{\xi }^2)$. Numerical solution is plotted against the exact solution using different values of collocation points and wavelength λ is given in Figures 7–10. It is clear from the figures that exact and numerical solutions are getting closer and closer by increasing the collocation points and wavelength λ. The error between exact and numerical solutions in different norms is shown in Table 1, which further confirms our theoretical justification given in Theorem 4.5.

Figure 7. Example 2: Numerical solution (left) vs exact solution (right) at N = 15 $\lambda =80$.

Figure 8. Example 2: Numerical solution (left) vs exact solution (right) at N = 30 $\lambda ={10}^2$.

Figure 9. Example 2: Numerical solution (left) vs exact solution (right) at N = 50 $\lambda ={10}^3$.

Figure 10. Example 2: Numerical solution (left) vs exact solution (right) at N = 100 $\lambda ={10}^4$.

Table 1. Example 2: Error behaviour between exact and numerical solution.

N	$L_2$ error	$L_{\mathrm{\infty }}$ error	N	$L_2$ error	$L_{\mathrm{\infty }}$ error
6	$2.283e^{-02}$	$6.137e^{-02}$	20	$1.489e^{-14}$	$3.977e^{-14}$
8	$8.384e^{-04}$	$2.543e^{-03}$	22	$2.138e^{-16}$	$6.031e^{-16}$
10	$9.218e^{-06}$	$2.806e^{-05}$	24	$2.898e^{-16}$	$8.368e^{-16}$
14	$5.723e^{-09}$	$1.512e^{-08}$	26	$2.924e^{-16}$	$8.507e^{-16}$
16	$1.312e^{-12}$	$2.676e^{-10}$	28	$2.836e^{-16}$	$7.533e^{-16}$
18	$1.300e^{-12}$	$1.516e^{-12}$	30	$4.364e^{-16}$	$1.606e^{-16}$

5.3. Example 3

Consider the Helmholtz in 3-dimension of the form: Consider the following example of Helmholtz equation: (12) $\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\sigma }^2}+{\lambda }^2\mathrm{\Psi }(\eta ,\xi ,\sigma )=(1-{\lambda }^2)(\sin (\kappa {\eta }^2\pi )\sin (\mathrm{\kappa \pi }{\xi }^2)\sin (\kappa {\sigma }^2\pi ))$(12) subject to the boundary conditions: (13) $\{\begin{array}{l}\frac{\mathrm{\partial \Psi }}{\mathrm{\partial }\eta }(0,\xi ,\sigma )=0,\mathrm{\Psi }(1,\xi ,\sigma )=-\sin (\mathrm{\pi \xi })\sin (\mathrm{\pi \sigma }),\mathrm{\Psi }(\eta ,0,\sigma )=0,\\ \mathrm{\Psi }(\eta ,1,\sigma )=0\mathrm{\Psi }(\eta ,\xi ,0)=0,\mathrm{\Psi }(\eta ,\xi ,1)=0.\end{array}$(13) The exact solution of Equation (12(12) $\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\eta }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\xi }^2}+\frac{{\mathrm{\partial }}^2\mathrm{\Psi }}{\mathrm{\partial }{\sigma }^2}+{\lambda }^2\mathrm{\Psi }(\eta ,\xi ,\sigma )=(1-{\lambda }^2)(\sin (\kappa {\eta }^2\pi )\sin (\mathrm{\kappa \pi }{\xi }^2)\sin (\kappa {\sigma }^2\pi ))$(12) ) is given by $\mathrm{\Psi }(\eta ,\xi ,\sigma )=\cos (\mathrm{\pi \eta })\sin (\mathrm{\pi \xi })\sin (\mathrm{\pi \sigma })$. Figures 11–14 are plotted using the fixed value of $\kappa =1$, while increasing the values of collocation points and wavelength as shown in figures captions. The approximate and exact solution confirm the high accuracy of our numerical scheme. The error behaviour in both $L_2$ and $L_{\mathrm{\infty }}$ norm between exact and approximate solution is shown in Table 2, which conforms to the exponential convergence of our proposed scheme.

Figure 11. Example 3: Numerical solution (left) vs exact solution (right) at N = 10 $\lambda =15$.

Figure 12. Example 3: Numerical solution (left) vs exact solution (right) at N = 15 $\lambda =30$.

Figure 13. Example 3: Numerical solution (left) vs exact solution (right) at N = 20 $\lambda =40$.

Figure 14. Example 3: Numerical solution (left) vs exact solution (right) at N = 25 $\lambda =50$.

Table 2. Example 3: Error behaviour between exact and numerical solution.

N	$L_2$ error	$L_{\mathrm{\infty }}$ error	N	$L_2$ error	$L_{\mathrm{\infty }}$ error
6	$2.165e^{-01}$	$7.255e^{-01}$	20	$2.837e^{-08}$	$2.134e^{-07}$
8	$4.907e^{-02}$	$1.971e^{-02}$	22	$1.294e^{-09}$	$1.877e^{-08}$
10	$4.538e^{-03}$	$2.469e^{-03}$	24	$6.308e^{-11}$	$8.475e^{-09}$
14	$5.745e^{-05}$	$2.854e^{-04}$	26	$1.593e^{-11}$	$4.481e^{-10}$
16	$5.549e^{-06}$	$3.316e^{-05}$	28	$1.539e^{-12}$	$1.107e^{-11}$
18	$3.790e^{-07}$	$2.174e^{-06}$	30	$9.675e^{-12}$	$1.516e^{-12}$

6. Conclusion

In this study spectral methods based on shifted Chebyshev polynomials and their properties are used to successfully address the complexities involved in modelling acoustic wave propagation in underwater environments, governed by the elliptic partial differential equations of Helmholtz type. The proposed numerical scheme significantly enhances the efficiency of simulating underwater acoustic wave propagation, which not only reduces computational time but also maintains high accuracy. Spectral methods exhibit exponential convergence, which means that the error in the numerical solution decreases exponentially as the number of basis functions increases. As a consequences a very high accuracy with relatively few degrees of freedom can be achieved as compared to other numerical methods like finite difference or finite element methods, particularly for problems with smooth solutions. A number of numerical experiments have been performed to confirm the theoretical justification of exponential convergence. Our future goal is to apply the spectral method to the Helmholtz equation in cylindrical coordinates and its fractional form.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Deanship of Scientific Research, Vice Presidency for Graduate Studies and Scientific Research, King Faisal University, Saudi Arabia [grant number 6154].