An efficient method for 3D Helmholtz equation with complex solution

1.   Introduction

This paper is devoted to establishing a highly accurate method for solving the Helmholtz equation that has many applications in time-harmonic wave propagation phenomena, such as the acoustic cavity and radiation wave ^[1]. Some of the numerical methods that have been used to solve the Helmholtz equation in recent years are: Legendre wavelet method ^[2], least-squares finite element method ^[3], meshless Chebyshev collocation method ^[4], implicit finite difference method ^[5], staggered discontinuous Galerkin method ^[6], Fourier-Bessel method ^[7], radial basis function-generated finite difference scheme ^[8] and hybrid approach proposed in ^[9].

Nowadays, researchers use orthogonal polynomials instead of the Fourier functions to solve various classes of equations numerically, because these types of basis functions have spectral accuracy property and their computational volume is less than other basis functions. For instance, see ^{[10,11,12,13,14]}. Orthogonal polynomials are defined into two categories based on the definition of their internal product ^[15], discrete and continuous. In numerical algorithms, the expansion coefficients for orthogonal continuous polynomials are achieved by integrating, but the expansion coefficients for orthogonal discrete polynomials are got using a finite summation. In recent years, researchers have been successfully used discrete polynomials to provide numerical methods for solving ordinary and fractional differential equations. Some differential equations that are solved using discrete polynomials are: variable-order fractional extended Fisher-Kolmogorov equation ^[16], time-delay fractional optimal control problems ^[17], fractional Volterra partial integro-differential equations ^[18], fractional viscoelastic model ^[19], fractional reaction-advection-diffusion equations ^[20], etc..

In this paper, we utilize the orthonormal shifted discrete Chebyshev polynomials (OSDCPs) to obtain the numerical solution of the 3D Helmholtz equation.

The main objectives of this study are briefly expressed in the following:

● Investigating the 1D and 2D OSDCPs and their derivative matrices.

● Establishing a computational method based on the 1D and 2D OSDCPs to obtain the numerical solution of the 3D Helmholtz equation.

To solve the intended equation, we first approximate the unknown solution in the main problem by the 1D and 2D OSDCPs. Next, using these polynomials, we expand the second-order derivatives. So, the derivative matrices of the second order are obtained. In the following, we replace approximations obtained in the previous steps into the main problem. By this replacing and the collocation technique, the primary problem is turned into an algebraic equations system. Some of the most important advantages of the established method are listed in the following:

● The proposed method has spectral accuracy and low computational volume.

● A small number of the OSDCPs are required to obtain highly accurate results.

● The developed method transforms the solution of the main problem into the solution of a linear system of algebraic equations, which is easily solvable.

In the sequel, we provide five sections as follows: the 1D and 2D OSDCPs and some of their properties are provided in Section 2. A numerical method is adopted by applying the 2D OSDCPs in Section 3. Section 4 provides some numerical examples for investigating the accuracy of the proposed method. Section 5 explains the conclusion of the article.

2.   The Helmholtz equation and discrete Chebyshev polynomials

Here, we review the 3D Helmholtz equation and introduce the 1D and 2D OSDCPs together with some of their attributes.

2.1.   The 3D Helmholtz equation

In this work, we focus on the 3D Helmholtz equation

with the boundary conditions

or

where $\eta$ is the wave number, $f$ and $h_{\varsigma}$ for $\varsigma = 1, 2, \ldots, 6$ are given complex functions, $L_{1}$–$L_{3}$ are real constants. Also, $\vartheta$ is an unknown complex function.

2.2.   The 1D basis polynomials

The 1D OSDCPs are defined over $[0, z_{f}]$ in the following form ^[21]:

where $S^{(r)}_{\kappa}$ is the first kind Stirling numbers and

The set $\left\lbrace \mathcal{CP}_{z_{f}, \ell}(z, \widetilde{M})\right\rbrace ^{\widetilde{M}}_{\ell = 0}$ produces an orthogonal set over $[0, z_{f}]$ based on the following internal product definition:

So, we have

Any function $\vartheta(z)\in [0, z_{f}]$ can be expanded by the 1D OSDCPs as follows:

where

in which

and

Theorem 2.1. ^[21] Let $\mathbb{CP}_{z_{f}, \widetilde{M}}(z)$ is the vector introduced in (2.10). Then, the $\rho$th derivative of this vector can be gotten as

where components of the square matrix ${D}^{(\rho, z_{f})}_{z} = [{d}^{(\rho, z_{f})}_{ij}]$ are calculated by the following formula:

The 2D OSDCPs are constructed over $[0, 1]\times[0, 1]$ by the 1D OSDCPs as

Utilizing the 2D OSDCPs, any continuous function $\vartheta\left(\xi, \gamma \right)$ can be expanded over $[0, 1]\times[0, 1]$ in the following form:

where

with

and

Theorem 2.2. Suppose that $\mathbb{CP}_{{N}_{1}\, {N}_{2}}\left(\xi, \gamma\right)$ is the 2D OSDCPs vector demonstrated in (2.17). Then, the $\rho$th derivative of this vector can be calculated as follow:

where

in which ${I}_{\left({N}_{2}+1\right)\times \left({N}_{2}+1\right) } = diag(1, 1, \cdots, 1)$ and

Proof. According to (2.13), the $\rho$th derivative of the vector $\mathcal{CP}_{ij}\left(\xi, \gamma; {\rm{N}}_{1}, {\rm{N}}_{2}\right)$ $(i = 1, 2, \ldots, {\rm{N}}_{1}$, $j$ = $1, 2, \ldots, {\rm{N}}_{2})$ can be written in the following form:

By computing $\dfrac{d^{\rho}}{d\, \xi^{\rho}}\mathcal{CP}_{i}(\xi; {\rm{N}}_{1})\, (i = 1, 2, \ldots, {\rm{N}}_{1})$ and substituting in the above relation, we get

Then, by applying Theorem 2.1, we have

where

or equivalently

□

As a numerical example, we want to compute matrix $\text{D}^{(1)}_{\xi}$ expressed in the Theorem 2.2. By assuming that $\text{N}_{1} = \text{N}_{2} = 3$, we have

and

Then, we obtain

Similarly, the matrix $\text{D}^{(2)}_{\xi}$ is computed as

Theorem 2.3. Suppose that $\mathbb{CP}_{{N}_{1}\, {N}_{2}}\left(\xi, \gamma\right)$ is the 2D OSDCPs vector demonstrated in (2.17). Then, the $\rho$th derivative of this vector can be calculated as follow:

where

such that ${O}$ demonstrates the zero matrix of order $({N}_{2}+1)\times({N}_{2}+1)$ and

Proof. According to (2.13), the $\rho$th derivative of the vector $\mathcal{CP}_{ij}\left(\xi, \gamma; \text{N}_{1}, \text{N}_{2}\right)\, (\, i = 1, 2, \ldots, \text{N}_{1},$ $j$ = $1, 2, \ldots, \text{N}_{2})$ can be written in the following form:

By computing $\dfrac{d^{\rho}}{d\, \gamma^{\rho}}\mathcal{CP}_{j}(\gamma; \text{N}_{2})\, (j = 1, 2, \ldots, \text{N}_{2})$ and substituting in the above relation, we get

Then, by applying Theorem 2.1, we have

where

or equivalently

□

As a numerical example, for $\left(\text{N}_{1} = \text{N}_{2} = 3 \right)$, we get

3.   Numerical scheme

In this section, a numerical scheme is provided for the 3D Helmholtz Eq (2.1) by utilizing the OSDCPs. To this end, we represent the complex function $\vartheta(\xi, \gamma, z)$ as follow:

where $\omega$ and $\nu$ are real functions. Moreover, we consider the complex functions $f(\xi, \gamma, z)$, $g(\xi, \gamma)$ and $h_{\ell}(., .)\, (\ell = 1, 2, 3, 4)$ in the following forms:

where $f_{k}(., ., .)$, and $h_{lk}(., .)$ for $l = 1, 2, 3, 4, 5, 6, \, \, k = 1, 2$ are given real functions. By putting (3.1) into (2.1), we obtain

Also, by replacing (3.2) into (2.2) or (2.3), the boundary conditions are obtained as follows:

or

In the following, using (2.8) and (2.14), we expand the unknown functions $\omega(\xi, \gamma, z)$ and $\nu(\xi, \gamma, z)$ as follows:

where $\Upsilon$ and $\Lambda$ are matrices with $\left(\text{N}_{1}+1 \right)\left(\text{N}_{2}+1 \right)\times \left(\widetilde{M}+1 \right)$ unknown elements. Regarding Theorems 2.1–2.3, we have

and

The functions expressed in (3.4) or (3.5) can be expanded by the OSDCPs as

where $\mathcal{G}_{kj}$ and $\mathcal{H}_{ij}(.)$ for $k = 1, 2, \, i = 1, 2, 3, 4, \, j = 1, 2$ are given vectors. Using (2.3), (3.6) and (3.9), we obtain the following relations:

The residual functions can be defined via (2.1) and (3.6)–(3.8) as follows:

Finally, we get a system containing $2 (\text{N}_{1}+1)(\text{N}_{2}+1)\times (\mathcal{M}+1)$ equations by (3.10) and (3.11) as

where $\left(\xi_{\kappa}, \gamma_{\ell}, z_{m} \right)$ are defined as follows:

By solving (3.12), we determine the elements of the matrices $\Upsilon$ and $\Lambda$, and subsequently using (3.6), we derive a solution for primary problem (2.1). In this work, we have used the "linsolve" command of Matlab R2020b to solve this system.

The step-by-step algorithm of the proposed method is given in the following:

4.   Test problem

This section presents four test problems to demonstrate the validity of the suggested method.

We evaluate the precision of our algorithm using $L_{\infty}$ norm as follows:

where $\tilde{\omega}$ and $\tilde{\nu}$ are the approximations of $\omega$ and $\nu$, respectively. In addition, we compute the convergence order (CO) as

where $N_{1}$ and $N_{2}$ are the number of the 2D OSDCPs applied in the first and second implementations, respectively. It should be noted that in all examples we assume that $\text{N}_{1} = \text{N}_{2}$.

Example 4.1. Consider the 3D Helmholtz equation

with the boundary conditions

where $\left(L_{1}, L_{2}, L_{3}\right) = \left(\eta\cos\mu\cos\varpi, \eta\sin\mu\cos\varpi, \eta\sin\varpi \right)$. The true solution of this example is

We apply the method proposed in Section 3 for solving this example with some values $(\mu, \varpi)$, $(N, \widetilde{M})$ and $\eta$. The obtained results are presented in Tables 1 and 2 and Figures 1–3. These results indicate that the proposed method is sufficiently accurate.

Example 4.2. Consider the 3D Helmholtz equation

with the boundary conditions

where $\left(L_{1}, L_{2}, L_{3}\right) = \left(\eta\cos\mu\cos\varpi, \eta\sin\mu\cos\varpi, \eta\sin\varpi \right)$. The analytic solution of this example is

that is a plane wave. Using the above analytic solution, we can determine functions $h_{\ell}(., .), \, \ell = 1, 3, 5$. We apply the method proposed in Section 3 for solving this example with $(\mu, \varpi) = \left(\dfrac{\pi}{4}, 0\right)$. Results obtained for $\eta = \sqrt{2}, \, 5\sqrt{2}, \, 10\sqrt{2}, \, 15\sqrt{2}$ are presented in Tables 3 and 4, Figures 4–7. The results listed in Tables 3 and 4 indicate that the proposed method is sufficiently accurate even for large wave numbers. Figure 8 shows the logarithm of errors for various values of wave number that decrease with increasing $N$.

Example 4.3. Consider the 3D Helmholtz equation

with the boundary conditions

where $\left(L_{1}, L_{2}, L_{3}\right) = \left(\eta\cos\mu\cos\varpi, \eta\sin\mu\cos\varpi, \eta\sin\varpi \right)$. The analytic solution of this example is $\vartheta(\xi, \gamma, z) = e^{i\left(L_{1}\xi+L_{2}\gamma+L_{3}z\right) }$. We apply the method proposed in Section 3 for solving this example with $\left(\mu = \dfrac{\pi}{8}, \varpi = 0\right)$ and some values of $(N, \widetilde{M})$ and $\eta$. The obtained results are presented in Tables 5, 6 and Figures 9–12. These results indicate that the proposed method is sufficiently accurate.

Example 4.4. Consider the 3D Helmholtz equation

with

where the exact solution is $\vartheta(\xi, \gamma, z) = \dfrac{\sin(\eta \xi)\sin(\eta \gamma)\sin(\eta z)}{\eta^{2}}$. Other required information can be derived from the expressed exact solution. We apply the method proposed in Section 3 for solving this example with some values of $(N, \widetilde{M})$ and $\eta$. The obtained results are presented in Table 7, Figures 13 and 14. These results indicate that the proposed method is sufficiently accurate.

5.   Conclusions

In this paper, we provided a numerical scheme using the orthonormal shifted discrete Chebyshev polynomials. Using the proposed method, we convert the 3D Helmholtz equation into a system of algebraic equations which can be easily solved. The results obtained of solving some numerical examples confirmed the high accuracy of the presented algorithm. Note that the proposed scheme can be developed for other types of partial differential equations, such as the Kdv-Burgers-Kuramoto equation, the Schrödinger equation, and the Benjamin-Bona-Mahony equation.

Conflict of interest

The authors declare that they have no competing interests.