Radial basis function-generated finite differences with Bessel weights for the 2D Helmholtz equation

TL;DR

This paper introduces a novel RBF-FD method using Bessel function weights to solve the 2D Helmholtz equation, effectively reducing dispersion and pollution effects in numerical solutions.

Contribution

The paper presents a new RBF-FD scheme with Bessel weights and a regularization technique to improve stability and accuracy for the 2D Helmholtz equation.

Findings

Mitigates dispersion effects in numerical solutions

Controls condition number via regularization

Demonstrates convergence and efficiency through numerical tests

Abstract

In this paper we obtain approximated numerical solutions for the 2D Helmholtz equation using a radial basis function-generated finite difference scheme (RBF-FD), where weights are calculated by applying an oscillatory radial basis function given in terms of Bessel functions of the first kind. The problem of obtaining weights by local interpolation is ill-conditioned; we overcome this difficulty by means of regularization of the interpolation matrix by perturbing its diagonal. The condition number of this perturbed matrix is controlled according to a prescribed value of a regularization parameter. Different numerical tests are performed in order to study convergence and algorithmic complexity. As a result, we verify that dispersion and pollution effects are mitigated.

Tables

Table 1. Table 1 : This table shows values of the condition number of the matrices 𝐉 k subscript 𝐉 𝑘 \mathbf{J}_{k} and 𝐉 ~ k subscript ~ 𝐉 𝑘 \widetilde{\mathbf{J}}_{k} for several stencil sizes n 𝑛 n , and relative local truncation errors of the approximation ( ∂ S , k ) ⁡ u 2 ​ ( 𝐱 ) ≈ ∂ x u 2 ​ ( 𝐱 ) subscript 𝑆 𝑘 subscript 𝑢 2 𝐱 subscript 𝑥 subscript 𝑢 2 𝐱 (\partial_{S,k})u_{2}(\mathbf{x})\approx\partial_{x}u_{2}(\mathbf{x}) . We used stencils as in Fig. 1 . We compared relative errors | ( ∂ S , k ) ⁡ u 2 ​ ( 𝐱 ) − ∂ x u 2 ​ ( 𝐱 ) | / | ∂ x u 2 ​ ( 𝐱 ) | subscript 𝑆 𝑘 subscript 𝑢 2 𝐱 subscript 𝑥 subscript 𝑢 2 𝐱 subscript 𝑥 subscript 𝑢 2 𝐱 |(\partial_{S,k})u_{2}(\mathbf{x})-\partial_{x}u_{2}(\mathbf{x})|/|\partial_{x}u_{2}(\mathbf{x})| produced by using LDL T , MDI and ITMDI with 15 15 15 iterations. The function u 2 subscript 𝑢 2 u_{2} as in ( 30 ).

Stencil size ($n$)	$\kappa \left({𝐉}_k\right)$	$\kappa \left({\widetilde{𝐉}}_k\right)$	Error LDLT	Error MDI	Error ITMDI
9	7.8e+14	7.8e+08	0.0036	0.0036	0.0036
13	1.34e+17	1.34e+11	0.000138	0.000139	0.000138
25	2.38e+18	2.38e+12	1.46e-07	7.93e-07	7.72e-07
49	1.61e+18	1.61e+12	6.44e-07	4.28e-08	5.79e-09
69	8.93e+17	8.93e+11	2.57e-06	1.14e-07	7.18e-09
81	1.44e+18	1.44e+12	3.98e-07	7.06e-08	1.45e-08

Table 2. Table 2 : This table shows values of the condition number of the matrices 𝐉 k subscript 𝐉 𝑘 \mathbf{J}_{k} and 𝐉 ~ k subscript ~ 𝐉 𝑘 \widetilde{\mathbf{J}}_{k} , respect to the stencil size n 𝑛 n , and relative local truncation errors of the approximation ( ∂ S , k ) ⁡ u 2 ​ ( 𝐱 ) ≈ ∂ x u 2 ​ ( 𝐱 ) subscript 𝑆 𝑘 subscript 𝑢 2 𝐱 subscript 𝑥 subscript 𝑢 2 𝐱 (\partial_{S,k})u_{2}(\mathbf{x})\approx\partial_{x}u_{2}(\mathbf{x}) . We used stencils as in Fig. 2 . We compared relative errors | ( ∂ S , k ) ⁡ u 2 ​ ( 𝐱 ) − ∂ x u 2 ​ ( 𝐱 ) | / | ∂ x u 2 ​ ( 𝐱 ) | subscript 𝑆 𝑘 subscript 𝑢 2 𝐱 subscript 𝑥 subscript 𝑢 2 𝐱 subscript 𝑥 subscript 𝑢 2 𝐱 |(\partial_{S,k})u_{2}(\mathbf{x})-\partial_{x}u_{2}(\mathbf{x})|/|\partial_{x}u_{2}(\mathbf{x})| produced by using LDL T , MDI and ITMDI with 15 15 15 iterations. We used the u 2 subscript 𝑢 2 u_{2} in ( 30 ).

Stencil size ($n$)	$\kappa \left({𝐉}_k\right)$	$\kappa \left({\widetilde{𝐉}}_k\right)$	Error LDLT	Error MDI	Error ITMDI
7	9.79e+10	9.79e+04	0.00311	0.00437	0.00311
13	1.9e+16	1.9e+10	3.42e-05	3.7e-05	3.42e-05
19	2.45e+16	2.45e+10	1.48e-08	8.36e-05	4.84e-05
31	1.61e+17	1.62e+11	1.15e-06	4.33e-07	1.18e-07
37	1.71e+17	1.71e+11	9.29e-08	3.23e-08	2.69e-08
61	9.75e+18	9.78e+12	1.07e-07	9.56e-09	5.28e-09

Table 3. Table 3 : Results for approximated solutions of ( 31 ). We used a square uniform grid in Ω ∩ ∂ Ω Ω Ω \Omega\cap\partial\Omega . For inner nodes the stencil size is n = 13 𝑛 13 n=13 , at boundary nodes n b = 15 subscript 𝑛 𝑏 15 n_{b}=15 , the number of nodes per wavelength is kept constant with N g = 6 subscript 𝑁 𝑔 6 N_{g}=6 .

$\frac{k}{2\pi }$	$\frac{1}{h}$	Nodes ($N$)	$\kappa \left(𝐇\right)$	${\Vert u_1-{\tilde{u}}_1\Vert }_{\infty }$	${\Vert u_2-{\tilde{u}}_2\Vert }_{\infty }$
10	60	3721	1.40e+04	1.97e-04	1.79e-04
20	120	14641	7.54e+04	1.95e-04	1.56e-04
40	240	58081	4.17e+05	1.96e-04	1.10e-04
80	480	231361	2.38e+06	1.98e-04	1.72e-04
120	720	519841	6.55e+06	1.96e-04	1.21e-04

Table 4. Table 4 : Results for approximated solutions of ( 31 ). We used an hexagonal uniform grid in Ω ∩ ∂ Ω Ω Ω \Omega\cap\partial\Omega . For inner nodes the stencil size is n = 13 𝑛 13 n=13 , at boundary nodes n b = 25 subscript 𝑛 𝑏 25 n_{b}=25 , the number of nodes per wavelength is kept constant with N g = 6 subscript 𝑁 𝑔 6 N_{g}=6 .

$\frac{k}{2\pi }$	$\frac{1}{h}$	Nodes ($N$)	$\kappa \left(𝐇\right)$	${\Vert u_1-{\tilde{u}}_1\Vert }_{\infty }$	${\Vert u_2-{\tilde{u}}_2\Vert }_{\infty }$
10	60	4237	2.52e+04	1.32e-04	3.30e-05
20	120	16752	1.09e+05	4.49e-05	3.35e-05
40	240	66861	5.25e+05	9.30e-05	3.59e-05
80	480	266680	3.46e+06	6.69e-05	5.62e-05
120	720	599458	8.21e+06	9.05e-05	4.25e-05

Table 5. Table 5 : Results for approximated solutions of ( 31 ). We used a square uniform grid in Ω ∩ ∂ Ω Ω Ω \Omega\cap\partial\Omega . For inner nodes the stencil size is n = 19 𝑛 19 n=19 , at boundary nodes n b = 25 subscript 𝑛 𝑏 25 n_{b}=25 , the number of nodes per wavelength is kept constant with N g = 6 subscript 𝑁 𝑔 6 N_{g}=6 .

$\frac{k}{2\pi }$	$\frac{1}{h}$	Nodes ($N$)	$\kappa \left(𝐇\right)$	${\Vert u_1-{\tilde{u}}_1\Vert }_{\infty }$	${\Vert u_2-{\tilde{u}}_2\Vert }_{\infty }$
10	60	4237	8.18e+04	1.97e-05	1.23e-05
20	120	16752	8.58e+05	2.22e-05	1.77e-05
40	240	66861	6.12e+05	1.94e-05	1.47e-05
80	480	266680	5.22e+06	1.78e-05	1.28e-05
120	720	599458	9.53e+06	1.85e-05	9.86e-06

Table 6. Table 6 : Here k 2 ​ π = 20 𝑘 2 𝜋 20 \frac{k}{2\pi}=20 . In a square uniform grid, we took stencils of size: n = 9 𝑛 9 n=9 for inner nodes, and n b = 15 subscript 𝑛 𝑏 15 n_{b}=15 for boundary nodes. By applying a linear regression we have that log 10 ⁡ ( ‖ u 1 − u ~ 1 ‖ ∞ ) ≈ 5.76 ​ log 10 ⁡ ( h ) + 8.3 subscript 10 subscript norm subscript 𝑢 1 subscript ~ 𝑢 1 5.76 subscript 10 ℎ 8.3 \log_{10}(\|u_{1}-\widetilde{u}_{1}\|_{\infty})\approx 5.76\log_{10}(h)+8.3 and log 10 ⁡ ( ‖ u 2 − u ~ 2 ‖ ∞ ) ≈ 6.132 ​ log 10 ⁡ ( h ) + 8.932 subscript 10 subscript norm subscript 𝑢 2 subscript ~ 𝑢 2 6.132 subscript 10 ℎ 8.932 \log_{10}(\|u_{2}-\widetilde{u}_{2}\|_{\infty})\approx 6.132\log_{10}(h)+8.932 .

NPW=$N_g$	$\frac{1}{h}$	Nodes ($N$)	$\kappa \left(𝐇\right)$	${\Vert u_1-{\tilde{u}}_1\Vert }_{\infty }$	${\Vert u_2-{\tilde{u}}_2\Vert }_{\infty }$
6.0	120.0	14641	7.38e+04	2.54e-04	1.71e-04
8.6	171.4	29584	1.03e+05	2.67e-05	1.56e-05
12.2	244.9	60025	1.47e+05	3.07e-06	1.79e-06
17.5	350.0	122500	2.06e+05	3.55e-07	2.06e-07
25.0	500.0	251001	2.75e+09	7.65e-08	2.64e-08

Table 7. Table 7 : Here k 2 ​ π = 20 𝑘 2 𝜋 20 \frac{k}{2\pi}=20 . In a hexagonal grid, we took stencils of size: n = 13 𝑛 13 n=13 for inner nodes, and n b = 25 subscript 𝑛 𝑏 25 n_{b}=25 for boundary nodes. By applying a linear regression we have that log 10 ⁡ ( ‖ u 1 − u ~ 1 ‖ ∞ ) ≈ 6.01 ​ log 10 ⁡ ( h ) + 10.01 subscript 10 subscript norm subscript 𝑢 1 subscript ~ 𝑢 1 6.01 subscript 10 ℎ 10.01 \log_{10}(\|u_{1}-\widetilde{u}_{1}\|_{\infty})\approx 6.01\log_{10}(h)+10.01 and log 10 ⁡ ( ‖ u 2 − u ~ 2 ‖ ∞ ) ≈ 5.82 ​ log 10 ⁡ ( h ) + 9.33 subscript 10 subscript norm subscript 𝑢 2 subscript ~ 𝑢 2 5.82 subscript 10 ℎ 9.33 \log_{10}(\|u_{2}-\widetilde{u}_{2}\|_{\infty})\approx 5.82\log_{10}(h)+9.33 .

NPW=$N_g$	$\frac{1}{h}$	Nodes ($N$)	$\kappa \left(𝐇\right)$	${\Vert u_1-{\tilde{u}}_1\Vert }_{\infty }$	${\Vert u_2-{\tilde{u}}_2\Vert }_{\infty }$
6.0	120.0	16752	1.34e+24	9.18e-03	3.31e-03
8.6	171.4	33960	9.84e+05	1.22e-04	1.05e-04
12.2	244.9	69338	1.69e+16	7.37e-05	1.35e-05
17.5	350.0	141753	4.04e+32	6.89e-06	7.51e-06
25.0	500.0	289291	2.14e+20	8.50e-07	3.82e-07

Table 8. Table 8 : Results in computing solutions corresponding to the smooth velocity model in ( 36 ).

$\omega /2\pi$	Nodes ($N$)	$\kappa \left(𝐇\right)$	Time (s) for ${𝐇}_{\Omega }$	Time (s) for ${𝐇}_{\partial \Omega }$	Time (s) for LU
2.5	1072	1.09e+04	0.37	0.24	0.02
5	4404	1.46e+05	1.18	0.45	0.10
10	17563	6.25e+05	4.72	1.01	0.45
20	70585	1.93e+06	18.98	2.58	3.00
40	283458	1.27e+07	75.10	9.05	17.14

Table 9. Table 9 : Results in computing solutions corresponding to the 2004 BP velocity model.

$\omega /2\pi$	Nodes ($N$)	$\kappa \left(𝐇\right)$	Time for ${𝐇}_{\Omega }$ (s)	Time for ${𝐇}_{\partial \Omega }$ (s)	Time for LU (s)
2	38488	9.01e+04	20.59	4.72	0.99
4	153223	3.43e+05	74.27	10.52	5.42
6	345389	9.93e+05	173.78	19.75	16.04
8	613529	7.22e+06	192.71	33.50	35.76

Equations

$$\left\{\begin{array}[]{rcll}-\Delta u(\mathbf{x})-\omega^{2}c(\mathbf{x})^{-2}u(\mathbf{x})&=&f(\mathbf{x}),&\mbox{ in }\Omega\\ b\frac{\partial}{\partial\mathbf{n}}u(\mathbf{x})+\mathrm{i}\omega c(\mathbf{x})^{-1}\mathcal{B}u(\mathbf{x})&=&g(\mathbf{x}),&\mbox{ on }\Gamma=\partial\Omega\end{array}\right.$$

$${\phi }_k^{(d)}(r)=\frac{J_{d/2-1}(kr)}{(kr{)}^{d/2-1}},d=1,2,\dots ,$$

$$S(\mathbf{x})=\sum _{j=1}^n{\alpha }_j{\phi }_k^{(d)}(\Vert \mathbf{x}-{\mathbf{x}}_j\Vert )$$

$${\varphi }_{C,k}(r)=J_0(k\sqrt{r^2+C^2}),$$

$${\varphi }_k(r)=J_0(kr),$$

$$\underset{\delta \to \infty }{\lim }{2}^{\delta }\delta !\frac{J_{\delta }(2\sqrt{\delta }r)}{2\sqrt{\delta }r}=e^{-r^2}$$

$$\mathbf{b}=(\mathbf{A}+\beta \mathbf{I})\tilde{\mathbf{y}},$$

$$\{\begin{array}{ll}\mathcal{L}u(\mathbf{x})=f(\mathbf{x}),& \text{mbox}if\mathbf{x}\in \Omega \\ \mathcal{B}u(\mathbf{x})=g(\mathbf{x}),& \text{mbox}if\mathbf{x}\in \partial \Omega ,\end{array}$$

$${\mathcal{P}}_{X,\epsilon }u(\mathbf{x})=\sum _{j=1}^N{\alpha }_j{\Phi }_{\epsilon }(\mathbf{x},{\mathbf{x}}_j),$$

$$u({\mathbf{x}}_k)=\sum _{j=1}^N{\alpha }_j{\Phi }_{\epsilon }({\mathbf{x}}_k,{\mathbf{x}}_j),\text{mbox}withk=1,\dots ,N.$$

$$\left(\begin{array}[]{c}\alpha_{1}\\ \vdots\\ \alpha_{N}\\ \end{array}\right)=(\Phi_{X,\varepsilon})^{-1}u|_{X},$$

$$u(\mathbf{x})\approx \tilde{u}(\mathbf{x})={\mathcal{P}}_{S_i,{\epsilon }_i}u(\mathbf{x})=\sum _{j=1}^{n_i}{\alpha }_j^i{\Phi }_{{\epsilon }_i}(\mathbf{x},{\mathbf{x}}_j^i).$$

$${\Phi }_{S_i,{\epsilon }_i}{\mathbf{\alpha }}^i=U_i$$

$${\mathbf{\alpha }}^i={\Phi }_{S_i,{\epsilon }_i}^{-1}{U}_i.$$

$$\mathcal{L}{\mathcal{P}}_{S_i,{\epsilon }_i}u({\mathbf{x}}_i)$$

$$=\mathcal{L}{\Phi }_{S_i,{\epsilon }_i}^1{\mathbf{\alpha }}^i$$

$$=\mathcal{L}{\Phi }_{S_i,{\epsilon }_i}^1{\Phi }_{S_i,{\epsilon }_i}^{-1}{U}_i$$

$$\mathcal{L}\Phi_{S_{i},\varepsilon_{i}}^{1}=\left(\begin{array}[]{cccc}\mathcal{L}\Phi_{\varepsilon_{i},}(\mathbf{x},\mathbf{x}_{1}^{i})|_{\mathbf{x}=\mathbf{x}_{i}}&\mathcal{L}\Phi_{\varepsilon_{i},}(\mathbf{x},\mathbf{x}_{2}^{i})|_{\mathbf{x}=\mathbf{x}_{i}}&\cdots&\mathcal{L}\Phi_{\varepsilon_{i},}(\mathbf{x},\mathbf{x}_{n_{i}}^{i})|_{\mathbf{x}=\mathbf{x}_{i}}\end{array}\right)$$

$$\mathcal{B}{\mathcal{P}}_{S_i,{\epsilon }_i}u({\mathbf{x}}_i)=\mathcal{B}{\Phi }_{S_i,{\epsilon }_i}^1{\Phi }_{{}_{S_i,{\epsilon }_i}}^{-1}{U}_i.$$

$$W{\mathcal{L}}_{S_i,{\epsilon }_i}{U}_i$$

$$W{\mathcal{B}}_{S_i,{\epsilon }_i}{U}_i$$

$$\tilde{u}(\mathbf{x})=\sum _{j=1}^{n_i}{\alpha }_j^i{J}_0(k_i\Vert \mathbf{x}-{\mathbf{x}}_j^i\Vert ).$$

$$U_i={\mathbf{J}}_{k_i}{\mathbf{\alpha }}_i,$$

$$\Delta \tilde{u}(\mathbf{x}){|}_{\mathbf{x}={\mathbf{x}}_i}$$

$${\Delta }_{S_i,k_i}u({\mathbf{x}}_i)$$

$$(\mathbf{I}-\beta {\tilde{\mathbf{J}}}_k^{-1}{)}^{-1}=\sum _{m=0}^{\infty }(\beta {\tilde{\mathbf{J}}}_k^{-1}{)}^m$$

$$\kappa ({\tilde{\mathbf{J}}}_k)<\kappa ({\mathbf{J}}_k).$$

$$\begin{array}[]{ccc}{\mathbf{J}}_{k}^{-1}-\widetilde{\mathbf{J}}_{k}^{-1}&=&\widetilde{\mathbf{J}}_{k}^{-1}(\mathbf{I}-\beta\widetilde{\mathbf{J}}_{k}^{-1})^{-1}-\widetilde{\mathbf{J}}_{k}^{-1}\\ &=&\widetilde{\mathbf{J}}_{k}^{-1}\left((\mathbf{I}-\beta\widetilde{\mathbf{J}}_{k}^{-1})^{-1}-\mathbf{I}\right)\end{array}$$

$$\mathbf{\alpha }-\tilde{\mathbf{\alpha }}=\frac{1}{\beta }(\beta {\tilde{\mathbf{J}}}_k^{-1})\left((\mathbf{I}-\beta {\tilde{\mathbf{J}}}_k^{-1}{)}^{-1}-\mathbf{I}\right)U.$$

$$\Vert \mathbf{\alpha }-\tilde{\mathbf{\alpha }}{\Vert }_2\le \frac{1}{{\lambda }_n}\left(\frac{\beta }{{\lambda }_n+\beta }\right)\Vert U{\Vert }_2.$$

$${\mathbf{J}}_k^{-1}$$

$$\mathbf{\alpha }=\frac{1}{\beta }\sum _{m=1}^{\infty }(\beta {\tilde{\mathbf{J}}}_k^{-1}{)}^{m}U$$

Full text

††journal: Elsevier

1 Introduction

The Helmholtz equation is an elliptic Partial Differential Equation (PDE) which represents time-independent solutions of the wave equation. This equation models a wide variety of physical phenomena. These include among others, acoustic wave scattering, time harmonic acoustic, electromagnetic fields, water wave propagation, membrane vibration and radar scattering. One of the objectives of numerical solutions of Helmholtz equation is to build a solver dealing with (i) a wide range of wave numbers and (ii) decrease the accumulation of spurious dispersion in computation due to the pollution effect. Given that an increase in wave number requires an appropriate increase in the mesh resolution for maintaining the level of accuracy, most numerical methods face difficulties for tackling the pollution effect.

In this work we consider to find numerical solutions for the 2$D$ Helmholtz equation given by

[TABLE]

where $\omega$ is the angular frequency, $c(\mathbf{x})>0$ is the sound speed of the continuous media, $f(\mathbf{x})$ is the source term, $\mathbf{n}$ is unitary normal vector to the boundary $\Gamma$, $b$ takes values zero or one, $\mathcal{B}$ is a certain linear operator and $g(\mathbf{x})$ is certain exact data on $\Gamma$. Here $\mathrm{i}=\sqrt{-1}$.

We apply a basic idea inspired in a combination of Trefftz method and Radial basis functions (RBF). Trefftz method consider linear combinations of solutions of the equation itself to solve a PDE. In this case the solutions we apply are given by oscillatory radial basis functions in terms of Bessel functions of the first kind. The weights of the linear combination are found by considering local stencils in the way of finite difference method(FD). This joint formulation of RBF and FD is well known as RBF-FD method. RBF-FD methods have been widely applied in the solution of partial differential equations. For a better idea the reader may consult [7, 10] and references therein; in particular, Fornberg [12] studied the family of oscillatory Bessel RBF’s we use here.

When solving a differential equation is very important to apply methods adaptable to the geometry or node distribution in the solution domain. RBFs are an appropriate meshless tool, given that they only depend on the distance between points, do not require any prescribed structure on them. This method consider solutions in the form $S(\mathbf{x})=\sum_{j=1}^{n}\varphi(\|\mathbf{x}-\mathbf{x}_{j}\|)+p(\mathbf{x})$; where $\varphi(r)$ is a radial function such as the Gaussian family of radial basis functions (GRBF) $\varphi(r)=e^{-(\varepsilon r)^{2}}$, $\|\cdot\|$ is the Euclidean norm and $X=\{\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{n}\}$ is a set of scattered points in the domain $\Omega$. Depending of the chosen RBF the low degree polynomial term $p(\mathbf{x})$ could be or not included.

From the first publications of Kansa [16, 17] until now, there has been an increasing interest and success in these methods with a wide range of applications [7, 13, 28]. In particular, several RBFs methods have been developed and applied to obtain numerical solutions of PDEs. There exist a large number of these functions, in fact a theorem from Bochner (1932) [15] shows that (under certain conditions) if the Fourier transform $\widehat{\Phi}$ of $\Phi$ is positive on ${\mathbb{R}^{d}}$ with $\widehat{\Phi}>0$ then $\Phi$ is positive definite in ${\mathbb{R}^{d}}$. Nevertheless, only some of them are usually chosen, depending on every particular application.

Usually, RBFs as multiquadrics $\phi(r)=\sqrt{1+(\varepsilon r)^{2}}$ contain a shape parameter $\varepsilon$ which decides the flatness of the function and by consequence the condition of the interpolation matrix. As $\varepsilon\rightarrow 0$ the shape of $\phi(r)$ goes from very peaked($\varepsilon$ large)to nearly flat ($\varepsilon$ small), in this last case the interpolation has shown to be remarkably accurate [10]. Until recently, the literature only have shown non-oscillatory RBFs, nevertheless the family of Bessel oscillatory RBFs (2) applied here provides existence and uniqueness of the interpolation problem, they do not diverge in the limit of flat basis functions for any node geometry and have exact polynomial reproduction of arbitrary order [8]. In this paper the role of the shape parameter is taken by the wave number $k$. The main argument for applying these functions in Helmholtz problem is that they are in themselves solutions of the equation.

The rest of the paper is organized as follows. In the next section we define the oscillatory Bessel function we are going to work with. In section 3 we give the general setting of the RBF-FD method. In Section 4 we describe the method of diagonal increments which deals with the ill-conditioning of the interpolation matrix. Sections 5,6 show results obtained in testing our methods with some well-known problems and benchmarks of current literature on Helmholtz equation.

2 Oscillatory RBF

There exist a wide number of Trefftz methods for the Helmholtz problems that have been surveyed in [14]. These are schemes of type finite elements where test and trial functions are local solutions of the differential equation to solve. Inspired by Trefftz methods, in this paper we work with a family of oscillatory RBF’s whose members are solutions of the homogeneous Helmholtz equation. Besides, given the oscillatory behavior of solutions of Helmholtz equation, it makes sense to consider such a family, whose members are given in terms of Bessel functions of the first kind.

The oscillatory RBF class $\varphi_{k}^{(d)}(r)$ (BRBF) is the family of radial basis functions given by

[TABLE]

which are detailed studied in [12]. Here $J_{\alpha}(r)$ is denoting the Bessel function of the first kind and order $\alpha$. Two remarkable properties of these functions are:

the non-singularity of the interpolation matrix for arbitrarily scattered data in up to $d$ dimensions, when $d>1$,

• 2.

and that the Laplace eigenvalue problem $\Delta\varphi+k^{2}\varphi=0$ has as bounded solutions, at the origin, the functions given in (2), thus any interpolant of the form

[TABLE]

will also satisfy $\Delta S+k^{2}S=0$.

In the case $d=2$ the current literature shows very few applications of the oscillatory RBF (2); this has been because the function (3) implies that $\Delta S=-k^{2}S$ so, by the weak maximum principle, the function $S$ cannot have local maximum at points where for some neighborhood it is negative; a fact that put restrictions to be used for general 2D interpolation. But in this work such a feature becomes a strength, since we are just interpolating solutions of Helmholtz problems, which locally can be seen as plane waves satisfying the homogeneous Helmholtz equation. In early works, as in [18], oscillatory RBF’s based on Bessel functions have been employed to solve the 2D Helmholtz equation with constant wavenumber within the approach of global collocation method and using the RBF

[TABLE]

which has two shape parameters with $k$ corresponding to the wavenumber and $C$ is empirically chosen. The ill-conditioning of the interpolation matrix that arises from (4) is overcome by way of a regularized singular value decomposition method.

For our interest, the 2D Helmholtz problem with large wavenumber, we take the special case $d=2$. So we work with the oscillatory RBF

[TABLE]

such that in the approach RBF-FD the shape parameter $k$ will be evaluated at the wave number $k(\mathbf{x})=\omega/c(\mathbf{x})$ corresponding to the center of the stencil.

Among the strengths of the oscillatory RBF family $\varphi_{k}^{(d)}(r)$ over other radial functions we must remark that the Gaussian family is contained in the Bessel RBF class in the limiting case

[TABLE]

In fact, all other RBFs could suffer divergence when $\varepsilon\to 0$. It was shown in [9] that such divergence can never occur when using GRBF, independently of the node distribution.

It is well known that for assembling the sparse matrix, which discretizes the Helmholtz problem is necessary to solve a small linear equation system at each node. As it will be seen, interpolation matrices are ill-conditioned and we deal with this issue by the Method of Diagonal Increments (MDI) [20], [22] adding to the diagonal entries a small regularization parameter $\beta>0$, thus we solve, instead of the linear system $\mathbf{b}=\mathbf{A}\mathbf{y}$, the equation

[TABLE]

where $\mathbf{I}$ is the identity matrix. An explanation of MDI will be given in the section 4, where it is shown that the matrix $\widetilde{\mathbf{A}}=\mathbf{A}+\beta\mathbf{I}$ is better conditioned than $\mathbf{A}$ and $\widetilde{\mathbf{y}}\approx\mathbf{y}$. Now are described the properties of discretizing Helmholtz problems with Bessel RFB.

3 Discretization by RBF-FD method

Under the RBF interpolation framework, we want to approximate the solution of boundary value problems in the form111Helmholtz problems we are dealing with in this paper, can be seen as particular cases of (7).

[TABLE]

where $\mathcal{L}$ and $\mathcal{B}$ are linear partial differential operators whose coefficients have a good enough regularity, and $\Omega$ is a bounded, open and connected set in $\mathds{R}^{d}$. It is assumed that (7) it is a well-posed problem.

In interpolation with Radial Basis Functions the goal is to reconstruct a real-valued or complex-valued function $u$ defined on a bounded domain $\Omega\subset\mathds{R}^{d}$ from the values $u(\mathbf{x}_{k})$ of $u$ on a finite set of $N$ scattered nodes $X=\{\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{N}\}\subset\Omega\subset\mathds{R}^{d}$, where $d$ is a positive integer. A radial basis function with shape parameter $\varepsilon$ is defined as a function $\Phi_{\varepsilon}:\mathds{R}^{d}\times\mathds{R}^{d}\rightarrow\mathds{R}$ such that $\Phi_{\varepsilon}(\mathbf{x},\mathbf{y})=\phi(\varepsilon\|\mathbf{x}-\mathbf{y}\|)$, where $\phi:[0,\infty)\rightarrow\mathds{R}$ is a single variable function [25, 7, 23]. A sufficiently smooth function $u:\Omega\subset\mathds{R}^{d}\rightarrow\mathds{R}$, with $\Omega$ an open set whose boundary is regular enough, can be approximated by the interpolant

[TABLE]

by forcing the condition $\mathcal{P}_{X,\varepsilon}u(\mathbf{x}_{k})=u(\mathbf{x}_{k})$ for $k=1,\ldots,N$, where weights $\alpha_{j}$ can be determined by solving the linear system

[TABLE]

Provided that the interpolation matrix $\Phi_{X,\varepsilon}=(\Phi_{\varepsilon}(\mathbf{x}_{k},\mathbf{x}_{j}))_{1\leq k,j\leq N}$ is non-singular, we have the solution

[TABLE]

where u|_{X}=\left(\begin{array}[]{cccc}u(\mathbf{x}_{1})&u(\mathbf{x}_{2})&\cdots&u(\mathbf{x}_{N})\end{array}\right)^{T}, hence the interpolant $\mathcal{P}_{X,\varepsilon}u$ in (8) is known.

When is necessary to compute solutions of (7) on large set of nodes, the resultant matrix for collocation method is dense, huge and ill-conditioned, carrying a prohibited computational cost. A variant of collocation method that allows to deal with large domains is the local version [24, 26]. Making a local interpolation is possible to obtain a sparse matrix which discretizes the linear partial differential operator. We now give a description of the RBF-FD method.

For an open and bounded set $\Omega\subset\mathds{R}^{d}$, let $X=\{\mathbf{x}_{i}\}_{i=1}^{N}\subset\Omega\cup\partial\Omega$ be a set of interpolation nodes. For every $\mathbf{x}_{i}\in X$, we create an influence domain $S_{i}\subset X$, which is formed by the $n_{i}$ nearest neighbor interpolation nodes, where $n_{i}$ is a positive integer, for $i=1,2,\ldots,N$. That is, we consider an $n_{i}$-stencil $S_{i}=\{\mathbf{x}_{j}^{i}\}_{j=1}^{n_{i}}\subset X$, where $\mathbf{x}_{1}^{i}\equiv\mathbf{x}_{i}$ and we denote the convex hull of the stencil $S_{i}$ by $\Omega_{i}$, i.e., $\Omega_{i}=\mbox{ConvexHull}(S_{i})$. Thus we have a collection of subsets $\{S_{i}\}_{i=1}^{N}$ formed by nodes of $X$. By RBF interpolation, for any $\mathbf{x}\in\Omega$ we choose an $\Omega_{i}$ such that $\mathbf{x}\in\Omega_{i}$, so we can approximate $u(\mathbf{x})$ as

[TABLE]

Note we have taken the shape parameter $\varepsilon_{i}$ depending of the location $\mathbf{x}_{i}$, thus the shape parameter of the RBF is conveniently manipulated, according to local known information related with the PDE.

By collocating the $n_{i}$ nodes of the stencil $S_{i}$, we obtain a small linear system

[TABLE]

where {U}_{i}=\left(\begin{array}[]{cccc}\widetilde{u}(\mathbf{x}_{1}^{i})&\widetilde{u}(\mathbf{x}_{2}^{i})&\cdots&\widetilde{u}(\mathbf{x}_{n_{i}}^{i})\end{array}\right)^{T}, $\Phi_{S_{i},\varepsilon_{i}}=\left(\Phi_{\varepsilon_{i}}(\mathbf{x}_{j}^{i},\mathbf{x}_{k}^{i})\right)_{1\leq j,k\leq n_{i}}$ is the local interpolation matrix and $\bm{\alpha}^{i}=\left(\alpha_{1}^{i}\ \ \alpha_{2}^{i}\ \ \cdots\ \ \alpha_{n_{i}}^{i}\right)^{T}$.

The unknown coefficients $\bm{\alpha}^{i}$ in (12) can be expressed in terms of the function values at the local interpolation nodes as

[TABLE]

The inverse matrix $\Phi_{S_{i},\varepsilon_{i}}^{-1}$ exists provided that $\Phi_{S_{i},\varepsilon_{i}}$ is positive definite, which is true for Bessel RBF (2).

Now, with the aim to get a local discretized version for (7) we consider $\mathbf{x}_{i}\in\Omega\cap X$ (or $\mathbf{x}_{i}\in\partial\Omega\cap X$). In both cases it must be applied a linear partial differential operator, either $\mathcal{L}$ or $\mathcal{B}$, to equation (11). For $\mathbf{x}_{i}\in\Omega$, with (13), we have

[TABLE]

where

[TABLE]

is a row matrix. Similarly, for $\mathbf{x}_{i}\in\partial\Omega\cap X$ we have

[TABLE]

We denote $W\mathcal{L}_{S_{i},\varepsilon_{i}}=\mathcal{L}\Phi_{S_{i},\varepsilon_{i}}^{1}\Phi_{S_{i},\varepsilon_{i}}^{-1}$ and $W\mathcal{B}_{S_{i},\varepsilon_{i}}=\mathcal{B}\Phi_{S_{i},\varepsilon_{i}}^{1}\Phi_{S_{i},\varepsilon_{i}}^{-1}$. From (14) and (15) it follows a discretized local version of (7)

[TABLE]

for $1\leq i\leq N$.

The above system of linear equations can be assembled forming a sparse matrix $\mathbf{H}$ of size $N\times N$ where the $i-$th row, associated to $\mathbf{x}_{i}\in X$, has at most $n_{i}$ nonzero entries, and the unknown column matrix is given by U=\left(\begin{array}[]{cccc}\widetilde{u}(\mathbf{x}_{1})&\widetilde{u}(\mathbf{x}_{2})&\cdots&\widetilde{u}(\mathbf{x}_{N})\end{array}\right)^{T} (we recall that $\mathbf{x}_{i}\equiv\mathbf{x}_{1}^{i}$). We can then obtain an approximated solution $\tilde{u}(\mathbf{x})$ at all interpolation nodes by solving $\mathbf{H}U=F$. The $\mathbf{H}$ can be thought as a discretized version of the problem (7).

3.1 Bessel RBF-FD

Suppose that $u$ is a solution of the Helmholtz equation $\Delta u(\mathbf{x})+k(\mathbf{x})^{2}u(\mathbf{x})=0$, for $\mathbf{x}\in\Omega$, and $u(\mathbf{x})$, with $\mathbf{x}\in\partial\Omega$, it satisfies certain boundary condition. If $X=\{\mathbf{x}_{i}\}_{i=1}^{N}\subset\Omega\cup\partial\Omega$ is a set of nodes, for $\mathbf{x}_{i}\in X\cap\Omega$ we take a stencil $S_{i}=\{\mathbf{x}_{j}^{i}\}_{j=1}^{n_{i}}\subset X$ based on $\mathbf{x}_{i}$, with $\mathbf{x}_{1}^{i}=\mathbf{x}_{i}$. For $\mathbf{x}\in\Omega_{i}=\mbox{ConvexHull}(S_{i})$ we define, with $k_{i}=k(\mathbf{x}_{i})$, the interpolant

[TABLE]

With the local interpolation matrix, $\mathbf{J}_{k_{i}}=(J_{0}(k_{i}\|\mathbf{x}_{l}^{i}-\mathbf{x}_{j}^{i}\|))_{1\leq l,j\leq n_{i}}$, which is positive definite [12], and forcing the condition $\widetilde{u}(\mathbf{x}_{l}^{i})=u(\mathbf{x}_{l}^{i})$, then from (16) we have the linear equation

[TABLE]

where U_{i}=\left(\begin{array}[]{cccc}u(\mathbf{x}_{1}^{i})&u(\mathbf{x}_{2}^{i})&\cdots&u(\mathbf{x}_{n_{i}}^{i})\end{array}\right)^{T} and \bm{\alpha}_{i}=\left(\begin{array}[]{cccc}\alpha_{1}^{i}&\alpha_{2}^{i}&\cdots&\alpha_{n_{i}}^{i}\end{array}\right)^{T}. In view that $\phi_{k}$, defined in (5), satisfies the homogeneous Helmholtz equation, then

[TABLE]

Hence the interpolant (16) satisfies the homogeneous Helmholtz equation. On the other hand, applying $\Delta_{S_{i},k_{i}}$ to the solution $u$, we have

[TABLE]

Note that $\Delta_{S_{i},k_{i}}u(\mathbf{x}_{i})-\Delta\widetilde{u}(\mathbf{x}_{i})=-k_{i}^{2}(u(\mathbf{x}_{i})-\widetilde{u}(\mathbf{x}_{i}))$, thus, for solutions of the homogeneous Helmholtz problem the local truncation error for the Laplace operator has a theoretical error depending of wavenumber at $\mathbf{x}_{i}$ and of the error of the local interpolant. The error of the approximated solutions is produced by the interpolant (16) and by the ill-conditioning of the matrix $\mathbf{J}_{k_{i}}$, in solving the linear system (17). Next we will deal with solutions of these systems.

4 Method of Diagonal Increments (MDI)

The interpolation matrix $\mathbf{J}_{k}$ is ill-conditioned, especially for certain node distributions. In literature there are several methods for dealing with the ill-conditioning when the shape parameter is small [11], but in our case we are taking the shape parameter as the wavenumber $k$, which can be large. So we have chosen the MDI. For our case $\mathbf{J}_{k}$ will be considered ill-conditioned when the condition number in the spectral norm 222The spectral norm matches with the matrix norm induced by Euclidean norm for vectors, i.e.,

$\|\mathbf{A}\|_{2}=\sup\{\|\mathbf{A}\bm{\alpha}\|_{2}\ :\ \bm{\alpha}\in\mathds{R}^{d}\mbox{ with }\|\bm{\alpha}\|_{2}=1\}$. , $\kappa(\mathbf{J}_{k})=\|\mathbf{J}_{k}\|_{2}\|\mathbf{J}_{k}^{-1}\|_{2}$, satisfies $\kappa(\mathbf{J}_{k})>10^{15}$, which hinders that the solution $\bm{\alpha}$ of $U=\mathbf{J}_{k}\bm{\alpha}$ be accurately calculated, in double precision, through Cholesky factorization. An alternative to compute $\bm{\alpha}$ with better tolerance respect to large condition numbers, allowing roughly 2 orders of magnitude more, i.e., up to $\kappa(\mathbf{J}_{k})\sim 10^{17}$, is the Block-$LDL^{T}$-decomposition (LDLT) [1, 6]. When $\mathbf{J}_{k}$ is ill-conditioned we solve the better conditioned problem $U=(\mathbf{J}_{k}+\beta\mathbf{I})\widetilde{\bm{\alpha}}$ instead, where $\mathbf{I}$ is the identity matrix and $\beta$ a small positive real number. Next, we will give some important aspects about the spectrum of $\mathbf{J}_{k}+\beta\mathbf{I}$. The following development is based on the Riley’s method [21].

Remark 1**.**

$\mathbf{J}_{k}$* is positive definite [12], thus its spectrum is real and positive. If $\displaystyle\{\lambda_{m}\}_{m=1}^{n}$ is its spectrum, with $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{n}$, then $\{\lambda_{m}+\beta\}_{m=1}^{n}$ is the spectrum of $\widetilde{\mathbf{J}}_{k}=\mathbf{J}_{k}+\beta\mathbf{I}$ and $\{\frac{\beta}{\lambda_{m}+\beta}\}_{m=1}^{n}$ is the spectrum of $\beta\widetilde{\mathbf{J}}_{k}^{-1}$, hence we have the spectral norms, $\|\beta\widetilde{\mathbf{J}}_{k}^{-1}\|_{2}=\frac{\beta}{\lambda_{n}+\beta}$ and $\|(\mathbf{I}-\beta\widetilde{\mathbf{J}}_{k}^{-1})^{-1}\|_{2}=\frac{\lambda_{n}+\beta}{\lambda_{n}}$. The above implies that the Neumann series $\sum_{m=0}^{\infty}(\beta\widetilde{\mathbf{J}}_{k}^{-1})^{m}$ converges and the equality*

[TABLE]

holds.

Remark 2**.**

If $\{\lambda_{m}\}_{m=1}^{n}$ is the spectrum of $\mathbf{J}_{k}$, as in remark 1, then the condition number of $\mathbf{J}_{k}$ is given by $\kappa(\mathbf{J}_{k})=\frac{\lambda_{1}}{\lambda_{n}}$ and $\kappa(\widetilde{\mathbf{J}}_{k})=\frac{\lambda_{1}+\beta}{\lambda_{n}+\beta}$, which implies that

[TABLE]

With this, the matrix $\widetilde{\mathbf{J}}_{k}$ is better conditioned than $\mathbf{J}_{k}$.

Remark 3**.**

Given that $\widetilde{\mathbf{J}}_{k}=\mathbf{J}_{k}+\beta\mathbf{I}$, then $\mathbf{J}_{k}^{-1}=\widetilde{\mathbf{J}}_{k}^{-1}(\mathbf{I}-\beta\widetilde{\mathbf{J}}_{k}^{-1})^{-1}$. If $\bm{\alpha}$ is the true solution of the equation $\mathbf{J}_{k}\bm{\alpha}=U$ and $\widetilde{\bm{\alpha}}$ is the solution for the perturbed system $\widetilde{\mathbf{J}}_{k}\widetilde{\bm{\alpha}}=U$, we can compare $\bm{\alpha}$ and $\widetilde{\bm{\alpha}}$. Note the following:

[TABLE]

thus,

[TABLE]

Since $\|(\mathbf{I}-\beta\widetilde{\mathbf{J}}_{k}^{-1})^{-1}-\mathbf{I}\|_{2}=\frac{\beta}{\lambda_{n}}$, finally we have,

[TABLE]

Now, with the purpose of obtaining closer solutions to the true one $\bm{\alpha}$, and improve the error bound (21), we consider the following. By using the Neumann series (18) we have

[TABLE]

If $\mathbf{J}_{k}\bm{\alpha}=U$ and $\widetilde{\mathbf{J}}_{k}\widetilde{\bm{\alpha}}=U$ then, from the remark 3,

[TABLE]

and

[TABLE]

If we truncate the series in (22) up to order $M$, we obtain an approximation of the true solution $\bm{\alpha}$, we denote it by

[TABLE]

From (22) and (24), the error of the approximation $\widetilde{\bm{\alpha}}_{M}$ can be bounded by using the formula

[TABLE]

Finally, by taking the Euclidean norm, we have in terms of the spectral norm,

[TABLE]

From the remark 1 we can conclude that

[TABLE]

An iterative procedure to compute (24), with the better conditioned matrix $\beta\widetilde{\mathbf{J}}_{k}$, can be obtained just by noting that, with $\widetilde{\bm{\alpha}}=\widetilde{\mathbf{J}}_{k}^{-1}U$,

[TABLE]

Hence we can compute $\widetilde{\bm{\alpha}}_{M}$ as:

[TABLE]

We call this algorithm the Iterative MDI (ITMDI) and its origin goes back to [21].

Since $0<\frac{\beta}{\lambda_{n}+\beta}<1$, it is important to note that the parameter $\beta$ needs to be selected to be large enough to improve conditioning but small enough so that the convergence of the method is faster [22, 21]. An undesirable situation is when $\lambda_{n}$ is near to the machine epsilon,333In double precision the machine epsilon is approximately 2.22e-16. in theses cases the ratio $\frac{\beta}{\lambda_{n}+\beta}$ is very close to 1 and in the convergence may be too slow. However, we have seen that with a small $M$ (for example $M=15$) is enough to improve the error (21). All codes were typed in Matlab R2016a and run in a laptop with Core i7 processor at 2.8 Ghz with 12 GiB of RAM.

4.1 Local truncation error and condition number

As can be seen in Figures 1 and 2, the condition number of $\mathbf{J}_{k}$ may becomes very large, then we adopt to handle values into a computationally acceptable range and to use this fact for obtaining the regularization parameter $\beta$. For small stencils we take $10^{7}\leq\kappa_{0}\leq 10^{14}$, with $\kappa_{0}=10^{7+\sqrt{n}}$ where $n$ is the size of the stencil, and we take

[TABLE]

as regularization parameter, ensuring, from remark 2, that $\kappa(\widetilde{\mathbf{J}}_{k})\approx\kappa_{0}$, which is an adequate condition number to work in double precision. Now, $\kappa_{0}$ must be taken such that $\kappa_{\min}\leq\kappa_{0}\leq\kappa_{\max}$, thus

[TABLE]

It is important to remark that $\beta$ is recalculated for every node, so its values may have a wide range; in our numerical tests and according to (28), has been about $0\leq\beta\leq 10^{-6}$.

We have noted empirically that the matrix $\mathbf{J}_{k}$ is worse conditioned for stencils with nodes collocated symmetrically on a regular grid, e.g. with square and hexagonal grids. See figures 1 and 2, where we can observe that severe ill-conditioning begins with symmetric stencils of $13$ nodes. However, in this case, with a small perturbation in the position of the nodes, its associated interpolation matrix $\mathbf{J}_{k}$ has a better condition number.

For numerical tests we considered the solutions for Helmholtz equation, $u_{1}$ and $u_{2}$, given by

[TABLE]

and

[TABLE]

where $H_{0}^{(1)}$ is the Hankel function of the first kind and $k$ is a constant wavenumber. The solution $u_{1}$ corresponds to the solution of the single source problem located at $\mathbf{x}_{s}=(2,2)$, whereas the solution $u_{2}$ corresponds to a solution of the problem with four single sources located at $\mathbf{x}_{s_{1}}=(-20,-20)$, $\mathbf{x}_{s_{2}}=(20,20)$, $\mathbf{x}_{s_{3}}=(-20,20)$, and $\mathbf{x}_{s_{4}}=(20,-20)$. In the first test, we considered the function $u_{2}$ to see the behavior of the local truncation error of the approximation $(\partial_{x})_{S,k}$. In addition we validate the conditioning of the matrix $\mathbf{J}_{k}$ and the better conditioning of $\widetilde{\mathbf{J}}_{k}$. With stencils $S$ as in Fig. 1 and Fig. 2 we approximate the solutions of the system $\mathbf{J}_{k}\bm{\alpha}=U$, with $U=u_{2}|_{S}$, by using LDLT, MDI, and ITMDI with 15 iterations ($M=15$ in (26)). We used these results to compute the approximation $(\partial_{x})_{S,k}u_{2}(\mathbf{x})$. In tables 1 and 2 we can see a comparison of the relative local truncation error, given by $|(\partial_{x})_{S,k}u_{2}(\mathbf{x})-\partial_{x}u_{2}(\mathbf{x})|/|\partial_{x}u_{2}(\mathbf{x})|$.

4.2 Pollution-effect and convergence

To see the impact of the pollution effect in numerical solutions, the standard procedure is: to compare the errors in several solutions obtained by increasing the wavenumber and keeping constant the Number of nodes Per Wavelength (NPW), i.e., the product $hk=\frac{2\pi}{Ng}$ should be constant [3].

4.2.1 Test 1

In this test we calculate the approximated solution for the problem

[TABLE]

with the known data $g(\mathbf{x})$, $\Omega=(-0.5,0.5)\times(-0.5,0.5)$ and $k=\omega c^{-1}$, with constant wave speed $c\equiv 1$. Results are verified with solutions $u_{1}$ and $u_{2}$ in (29) and (30), respectively. Tables 3, 4 and 5 show errors when the resolution is kept constant at $N_{g}=6$ NPW. In calculations we have taken the perturbed matrices, $\widetilde{\mathbf{J}}_{k}$, such that $\kappa_{0}=\kappa(\widetilde{\mathbf{J}}_{k})=10^{-6}\kappa(\mathbf{J}_{k})$. In the three cases, for uniform square and hexagonal grids, we see that the order of the error remains constant, i.e., $\|u_{1}-\widetilde{u}_{1}\|_{\infty}\sim\mathcal{O}(1)$ and $\|u_{2}-\widetilde{u}_{2}\|_{\infty}\sim\mathcal{O}(1)$, as $h\rightarrow 0$ with $hk=\frac{2\pi}{6}$. Hence, in these examples, pollution effects are mitigated.

Results of convergence tests are summarized in tables 6 and 7 and Fig 3. For these tests we choose $\widetilde{\mathbf{J}}_{k}=\beta\mathbf{I}+\mathbf{J}_{k}$, with $\beta$ according to (27), such that the condition number $\kappa_{0}=\kappa(\widetilde{\mathbf{J}}_{k})=10^{-4}\kappa(\mathbf{J}_{k})$.

4.2.2 Test 2

In this example we consider the problem

[TABLE]

with $\Omega=(0,1)\times(0,1)$, whose solution is given by the plane wave $u(x,y;k,\theta)=e^{\mathrm{i}k(x\cos\theta+y\sin\theta)}$ when the data $g$ on the impedance boundary condition is given by

[TABLE]

with $k_{1}=k\cos\theta$ and $k_{2}=k\sin\theta$. This example is standard for testing numerical dispersion of solvers for Helmholtz equation. The approximated solutions for this problem ware calculated using Gaussian RBF-FD on hexagonal grids with 7-stencil (GRBF-FD-7p), Bessel RBF-FD with 9-stencils (BRBF-FD-9p) and 13-stencils (BRBF-FD-13p) on uniform Cartesian grids and we compare with results reported in [3]. In all methods it has been fixed NPW$=2\pi$to keep constant resolution when the wavenumber k is increasing. For GRBF-FD-7p we have used the approximations in [19] with its respective shape parameter $\varepsilon_{op}$. We use it to approximate the Laplace operator and all partial derivative operators involved in the boundary condition. To solve the local interpolations in BRBF-FD9p and BRBF-FD13p we used condition numbers $\kappa_{0}=10^{-4}\kappa(\mathbf{J}_{k})$ and $\kappa_{0}=10^{-6}\kappa(\mathbf{J}_{k})$, respectively.

Results and comparisons can be seen in Fig. 4. On the left frame we can see that errors of BRBF-FD9p and BRBF-FD13p are smaller than GRBF-FD7p, ROT-FD9p and OP-FD9p, at least in two orders of magnitude. Besides, we see that there is less anisotropy in the error of BRBF-FD13p, where over all propagation angles we have improved the error at least three orders of magnitude. We can see that the behavior of GRBF-FD7p is similar to OP-FD9p. On the right frame we can see that the dispersion and pollution effects are mitigated with BRBF-FD13p because, when the wavenumber increases while we keep a fixed resolution with $kh=\frac{2\pi}{Ng}=1$, the error remains almost constant.

5 Examples with some Helmholtz problems

In this section we test our BRBF-FD scheme by computing numerical solutions of some Helmholtz problems with second and third order Absorbing Boundary Conditions (ABC) [5].

5.1 Approximated fundamental solutions

It is known that the problem $-\Delta u-k^{2}u=\delta(\mathbf{x})$ in the free-space $\mathds{R}^{2}$ has a unique solution when it is imposed the Sommerfeld radiation condition

[TABLE]

Particularly, the associated Green’s function, which is solution of $-\Delta u(\mathbf{x})-k^{2}u(\mathbf{x})=\delta(\mathbf{x}-\mathbf{x}_{0})$, [4] is given by $u(\mathbf{x})=G(\mathbf{x},\mathbf{x}_{0})$, with the fundamental solution

[TABLE]

We compute approximated Green’s functions in the free space truncated to a bounded domain $\Omega\subset\mathds{R}^{2}$ for $\mathbf{x}\neq\mathbf{x}_{0}$, where $\mathbf{x}_{0}\in\Omega$, through the boundary value problem

[TABLE]

where $\mathcal{B}=1+\frac{3}{4k^{2}}\frac{\partial^{2}}{\partial\bm{\tau}^{2}}-\frac{\mathrm{i}}{4k^{3}}\frac{\partial^{3}}{\partial\mathbf{n}\partial\bm{\tau}^{2}}$. The boundary condition corresponds to the ABC in the Padé approximation [5]. The single source is given by the Gaussian function

[TABLE]

where $\sigma$ is a value such that $\int_{\Omega}\widetilde{\delta}(\mathbf{x}-\mathbf{x}_{0})d\mathbf{x}\approx 1$. To solve (33) we have used the BRFB-FD9p scheme in $\Omega=(0,1)\times(0,1)$, by using square $9-$stencils at inner nodes, and $19$-stencils at boundary nodes. The results shown in Fig. 5 and Fig. 6 were calculated with $k=500$ and $N_{g}=6$  NPW, i.e., $h=\frac{2\pi}{6k}$. We point out that results show a good accuracy at $\mathbf{x}\neq\mathbf{x}_{0}$. Besides, for the source located at the center of the square domain, the wavelength of the numerical solution on the boundary, matches very close to the exact one. This is a good indication that dispersion errors are not significant. However in Fig. 6 (bottom) we see that amplitude has a considerable discrepancy with respect to the exact one, this is due to the approximation of ABC.

5.2 Heterogeneous medium

5.2.1 Smooth medium

For this qualitative test, we have calculated approximated solutions of the problem

[TABLE]

where $\mathcal{B}$ is the operator $\mathcal{B}=1+\frac{c(\mathbf{x})^{2}}{2\omega^{2}}\frac{\partial^{2}}{\partial\bm{\tau}^{2}}$ corresponding to the ABC of second order and $\Omega=(-0.5,0.5)\times(-0.5,0.5)$. Here we perform two examples with the velocity functions

[TABLE]

and

[TABLE]

For these velocity models, nodes distributions are sketched on the left column of Fig. 7. On center and right columns it can be seen the real part of the approximated solution for two different single sources. Table 8 shows results for required times to assembly sparse matrices $\mathbf{H}=\mathbf{H}_{\Omega}+\mathbf{H}_{\Gamma}$ and for solution of the system $-\mathbf{H}U=\mathbf{F}$ by LU factorization.

5.2.2 Non smooth medium (Test in the 2004 BP velocity model)

We consider the 2004 BP velocity benchmark, which is a popular model in research for velocity estimation methods in seismic imaging, which is presented as a challenge due to its complexity and large scale [2]. The velocity model $c(\mathbf{x})$ can be seen in the density plot on the middle-top in Fig. 8. Roughness of the velocity model such as hard interfaces and sharp transitions generates strong reflections that hinder the efficiency of known iterative methods due to the ill-conditioning of the matrix, when the number of iterations increases [27]. In addition, for large $\omega$, the interaction of high frequency waves with short wavelength structures such as discontinuities, increases the reflections, further deteriorating the convergence rate. In BRBF-FD schemes the local interpolation matrices $\mathbf{J}_{k}$ becomes dramatically ill-conditioned. However, in our tests, we have found empirically that in keeping the condition number of $\widetilde{\mathbf{J}}_{k}$ to the value $\kappa_{0}\approx 10^{(1+\sqrt{n_{s}})}$ remains the condition number of $\mathbf{H}$ in the range: $10^{4}\leq\kappa(\mathbf{H})\leq 10^{9}$. Hence, in this situation, it is feasible to perform LU factorization. We see in Fig. 8 that wavelengths of the wavefield have the expected behavior according to wave speed. In this test we solved (35) with ABC of second order. Table 9 resumes the computational complexity of the method by showing the execution times for key routines depending of frequency values and number of nodes.

6 Conclusions and future research

In this paper we perform local interpolation with a shape parameter-free oscillatory RBF based on Bessel functions of the first kind to obtain a higher order RBF-FD scheme for solving Helmholtz equation. In this approach the shape parameter is substituted for the local wavenumber. However due to the local interpolation the resulting matrices are extremely ill-conditioned even for stencils with a low number of nodes. We overcome this issue with a regularization method that introduce small perturbations in the diagonal o the matrix. In some tests we have achieved convergence rates of third and sixth order, this is a performance in accordance with the state of the art methods for Helmholtz problems.

Among the pending problems and future research that may improve our scheme it is very important:

To explore, in an analytic way, the behavior of the numerical solutions of Helmholtz equation with discontinuous and piecewise constant coefficients.

• 2.

To improve the choice of regularization parameters.

• 3.

To exploit the meshless features of the BRBF-FD method in 3D Helmholtz problems.

Acknowledgments

This work is supported by Colombian Oil Company ECOPETROL and COLCIENCIAS as a part of the research project grant No. 0266-2013.