# An Improved Collocation Meshless Method Based on the Variable Shaped Radial Basis Function for the Solution of the Interior Acoustic Problems

Authors: Shuang Wang; Shande Li; Qibai Huang; Kun Li
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

As an efficient tool, radial basis function (RBF) has been widely used for the multivariate approximation, interpolating continuous, and the solution of the particle differential equations. However, ill-conditioned interpolation matrix may be encountered when the interpolation points are very dense or irregularly arranged. To avert this problem, RBFs with variable shape parameters are introduced, and several new variation strategies are proposed. Comparison with the RBF with constant shape parameters are made, and the results show that the condition number of the interpolation matrix grows much slower with our strategies. As an application, an improved collocation meshless method is formulated by employing the new RBF. In addition, the Hermite-type interpolation is implemented to handle the Neumann boundary conditions and an additional sine/cosine basis is introduced for the Helmlholtz equation. Then, two interior acoustic problems are solved with the presented method; the results demonstrate the robustness and effectiveness of the method.


---

## Body

## 1. Introduction

As compared to traditional design criteria, such as the function, structure, and manufacturing technology, the comfort of a mechanical system has taken up a more and more important place. In order to improve the properties of the vibration and acoustic of machines, especially for those who have cavities, such as the aircrafts, vehicles, and rail trucks, it is very necessary to determine the interior acoustic field. Until now, the FEM and boundary element boundary (BEM) are still the most widely used numerical methods for solving the acoustic problems [1, 2].In the FEM, the computational domain of the problem is divided into lots of small elements, and the value of the field points can be obtained by a superposition of polynomial shape functions (usually simple and low-order). Thus, it is very suitable for solving the problems in bounded domains, such as the interior acoustic problems. On the other hand, BEM is more applicable to deal with the problems in unbounded domain, but it can hardly compare with the FEM for solving the bounded problems. The reason for this is that the system matrices may become fully polluted when the number of the elements increased, and if the collocation scheme is employed, the matrices may even be not symmetric [3, 4]. To improve the ability of BEM in solving these problems, several treatment methods have been proposed, such as the fast multipole BEM [2, 5, 6] and the wave boundary elements [7], which may make the BEM more efficient for solving interior acoustic problems. And during the last decades, considering the robust of the FEM, most researchers preferred to use it to solve the interior acoustic problems.It is well known that interior acoustic problems are generally governed by the Helmholtz equation or the wave equation. Nevertheless, the solution of these equations solved by standard FEM suffers from the so-called “pollution effect,” which is directly related to the dispersion [8]. In other words, the wave number of the numerical solution is different from the one of the exact solution, especially for the high wave numbers. In order to get reliable solutions for high wave numbers, a refined mesh or a higher order of the polynomial approximation can be applied. Through these schemes, the solutions can be improved, but this will result in large systems of equation; thus, the computation time will greatly increase. Thereby, the minimization or even the elimination of the pollution effect is strongly desirable. In the past years, a number of efforts have been spent on the development of FEM to achieve this aim. Thompson and Pinsky presented a Galerkin least squares FEM to reduce the dispersion [9], and Oberai and Pinsky developed the residual-based FEM to reach this goal [10]. Partition of unity finite element method (PUFEM) is another method, in which a priori knowledge about the solution of the problem is introduced into the FEM [11]. Unfortunately, this method will lead to ill-conditioned system matrices; hence, the use of the classical iterative solvers becomes complicated.In parallel with the FEM and BEM, meshless methods are relatively new numerical methods. As compared with FEM, no mesh is needed in meshless method for the discretization of the problems. On the contrary, the governing equations are discreted by arbitrarily distributed points. In [12–14], meshless methods have been presented to solve practical engineering problems. However, using the meshless methods to solve the acoustic problems, especially the interior acoustic problem is relatively few. Element free Galerkin method (EFG) was used to reduce the dispersion effect of the Helmholtz equation by  Bouillard et al.  [15] In addition,  Wenterodt  [16] applied the radial point interpolation method (RPIM) to solve the interior acoustic problem governed by Helmholtz equation. And the results showed that the RPIM has a better performance than EFG in reducing the dispersion effect; thus, better solutions can be obtained.The interpolation function in RPIM is RBF, which is widely used in numerical simulation for its efficiency and simplicity. However, the RBF interpolation may suffer from the contradiction between the accuracy and the stability, which can be described as uncertainly principle [17]. The condition number of the RBF interpolation matrix becomes very large when the interpolation points are irregularly or density arranged. And the ill-conditioned matrix will limit the application of the RBF, especially for large-scale problems. In order to guarantee the robustness and stability of interpolation with RBFs, many numerical treatments have been proposed, such as compactly supported RBF [12], precondition method [18], domain decomposition method [19], RBF with variable shape parameter [20], and node adaptively method [21]. In this paper, we focus on the RBF with variable shape parameter.The main concept of variable shaped RBFs is to apply different shape parameters in RBFs according to the local density of its corresponding interpolation point. Therefore, the columns of the interpolation matrix elements are more distinct and the condition number becomes relatively smaller. This numerical scheme has been investigated by some researchers. The accuracy of multiquadrics (MQ) interpolation was improved by using variable shape parameters according to the work of Kansa and Carlson [22, 23].  Sarra and Sturgill  proposed a random variation scheme for shape parameter selection [24]. And most researchers put their attention on the improvement of the accuracy rather than stability by applying variable shaped RBFs in the previous work. On the contrary, Zhou et al. focused on the great performance of variable shaped RBFs on the stability [20]. In addition, they proposed a quadric scheme for the variable shaped MQ and a linear scheme for the variable shaped inverse MQ (IMQ) to improve the interpolation matrix condition number.Since the variation strategy is the key factor for the variable shaped RBF. In this paper, we proposed several new variation strategies for the MQ and IMQ to improve their performance. Numerical results show that the variable shaped MQ with our strategies can improve both accuracy and stability, and the variation strategies for IMQ can improve the stability of the IMQ under an acceptable accuracy. Based on the variable shaped MQ, an improved collocation meshless method is presented, and instead of the traditional polynomials terms added for RBF, we proposed an additional sine/cosine basis for the Helmholtz equations. Then, two interior acoustic problems are solved with the proposed method. Since there are no analytical solutions for the two problems, refined FEM solutions are presented as the reference solutions. Comparison with the reference solutions shows the efficiency of the proposed method in solving the interior acoustic problems.The paper is organized as follows: in Section2, the interior acoustic problems are briefly introduced. Section 3 demonstrates the outline of the RBF interpolation and its error estimation. The RBF with variable shape parameters and their variation strategies are presented in Section 4. The improved collocation meshless method is described in Section 5. In Section 6, two examples are solved with the proposed method. The paper ends with conclusion in Section 7.


## 2. Interior Acoustic Problem

As shown in Figure1, we consider the interior acoustic problem in domain Ω with boundaries ΓD, ΓN, and ΓR. The propagation of the acoustic in the domain is governed by the following wave equation:
(2.1)ΔP-1c2∂2P∂t2=0,
where P is the acoustic pressure (small perturbations around a steady state) and c denotes the speed of sound in the medium. Considering a steady harmonic wave of angular frequency ω, the pressure can be described as
(2.2)P=p·e-jωt,
where p denotes the complex amplitude of acoustic pressure. By introducing the wave number k=ω/c and substituting (2.2) into (2.1), the following Helmholtz equation can be obtained:
(2.3)Δp+k2p=0.
In addition, the particle velocity v is another important quantity in the interior acoustic analysis, which can be obtained through the equation
(2.4)jρckv+∇p=0,
where ρ denotes the density of the medium. Associated with the governing equation, there are three typical boundary conditions encountered in the interior acoustic problem.Figure 1
The interior acoustic problem;Ω: computational domain; ΓD: sound soft boundary; ΓN: sound hard boundary; ΓR: impedance boundary.Sound soft boundary condition:(2.5)p=0.Sound hard boundary condition:(2.6)∂p∂n=0.
Impedance boundary condition:
(2.7)∂p∂n=-jωpZ,
where Z is the acoustic input impedance of the external domain. Note that in the two opposite limits Z→∞ and Z→0, the sound-hard and sound-soft boundary conditions are obtained.


## 3. Outline of the RBF Approximation

### 3.1. Basis Formulation

Consideringn distinct points x1,…,xn in thereal number set ℝd (d is the dimension) and a continuous function u(x), the approximation of the u(x) at the point x with RBF can be expressed as
(3.1)u^(x)=∑i=1nR(ri)ai+∑k=1mpk(x)ck,
where
(3.2)ri=∥x-xi∥2.
And R is the radial basis function, ai, ck are coefficients to be determined, and pk is the polynomial of an order not exceeding m. To ensure the u^(x) approximates the u(x), the following conditions should be satisfied:
(3.3)∑i=1nR(rij)ai+∑k=1mpk(xj)ck=u(xj).
To ensure the unique solution of the coefficients, additional equations should be considered
(3.4)∑i=1npk(xi)ai+∑j=1nDBpk(xj)ck=0,k=1,2,…,m.
Combining (3.3) and (3.4), the following matrix can be obtained:
(3.5)[ABTB0]{ac}={u0},
where the matrix element for A is R(rij), for B is pk(xi), and a,c are the coefficient vectors. Solving (3.5), the coefficients ai and ck can be gotten, and the approximation function is then given by (3.1). It should be noticed that the additional polynomials are not necessary if the RBF is strictly positive definite (SPD).


### 3.2. Typical RBFs

There are several advantages of the radial basis function (RBF) to be used for interpolation. First, there are generally few restrictions on the way the data prescribed and it can easily be applied in almost any dimension. Second, it has the naturally mesh-free property. And high-order accuracy is its other merit. Until now, many RBFs have been developed, and the most widely used global supported RBFs are listed in Table1.Table 1
Typical radial basis functions with dimensionless shape parameters.
Name
Expression
Shape parameter
Multiquadrics (MQ)
R
i
(
r
)
=
(
r
2
+
(
ε
l
)
2
)
q
ε
≥
0
,
q
Inverse multiquadrics (IMQ)
R
i
(
r
)
=
(
r
2
+
(
ε
l
)
2
)
-
q
ε
≥
0
,
q
Gaussian (GA)
R
i
(
r
)
=
exp
[
-
ε
(
r
/
l
)
2
]
ε
Thin plate spline (TPS)
R
i
(
r
)
=
r
η
η
Note:l is a characteristic length that relates to the nodal spacing in the computational domain.In these RBFs, the TPS is conditional positive definite, and others are SPD RBFs. Thus, when the TPS is used for interpolation, the polynomial terms should be employed to avoid the singularity.


### 3.3. Error Estimation and Stability of the RBF Interpolation

For the advantage of the RBF interpolation, it has been applied to a wide range of engineering problems. However, works concerning the error estimation and stability of this interpolation with mathematical theories are relatively less. The most noticeable work was completed by Schaback [17]. For the error estimation, the formulation is as follows:
(3.6)|u(x)-ua(x)|≤|u|f·P(x),
where P(x) is called the power function, which is the norm of the error function on F evaluated at x. And F is a native space, which can be described as
(3.7)F={f∈L2(Rn)∣∫Rd|u^(ω)|ψ(ω)dω<∞}.
In addition, u^ is in L1(Rd) and
(3.8)u|f2=∫Rd|u^(ω)|ψ(ω)dω<∞,
where ψ(ω) denotes the d-variate generalized Fourier transform of φ. Furthermore, Schaback proposed that the error and the sensitivity of the RBF interpolation can be described by P(x) and the eigenvalue λ of the interpolation matrix, which seems to be intimately related. P(x) and λ are bounded by F(h(x)) and G(h(x)) separately, the exact forms of which can be found in [17].And the shape parameters are directly contained inG(h(x)) but not in F(h(x)). Here, we list the G(h(x)) as follows:
(3.9)MQ:G(h)=he(-12.76εd/h),IMQ:G(h)=h-1e(-12.76εd/h),Guassian:G(h)=h-de(-40.71ε2d3/h2),
where h=maxx∈Ωminxi∈X∥x-xi∥, Ω is the interpolation domain, and X denotes the set of the interpolation points, d represents the dimensions, and ε is the shape parameter of the corresponding RBF.On the other hand, the error bound that contains the shape parameter was described in Madych’s theoretical analysis. His conclusions are proposed for two different classes of functions,Bσ and Eσ,
(3.10)Bσ={f∈L2(Rn):f~(ξ)=0if∥ξ∥>σ},Eσ={f∈L2(Rn):f~(ξ)<∞},
where σ is a positive constant,
(3.11)f~(ξ)=∫Rnf(x)e-i(x,ξ)dx
is the multivariate Fourier transform of function f(x), and
(3.12)∥f~∥Eσ2=∫Rn|f~(ξ)|2e|ξ|2/σdξ.
Considering the different functions, Madych presented the following error estimation for a class of RBF interpolations, which include Gaussian, MQ, and IMQ:
(3.13)forf∈Bσ:εr~O(eaελε/h)for0<λ<1anda>0,forf∈Eσ:εr~O(eaε2λε/h)for0<λ<1anda>0,
where εr is the error. It is worth noting that there is an intimate relation between the accuracy and the stability of the RBF interpolation. In the following section, we will present some strategies to balance the accuracy and the stability.


## 3.1. Basis Formulation

Consideringn distinct points x1,…,xn in thereal number set ℝd (d is the dimension) and a continuous function u(x), the approximation of the u(x) at the point x with RBF can be expressed as
(3.1)u^(x)=∑i=1nR(ri)ai+∑k=1mpk(x)ck,
where
(3.2)ri=∥x-xi∥2.
And R is the radial basis function, ai, ck are coefficients to be determined, and pk is the polynomial of an order not exceeding m. To ensure the u^(x) approximates the u(x), the following conditions should be satisfied:
(3.3)∑i=1nR(rij)ai+∑k=1mpk(xj)ck=u(xj).
To ensure the unique solution of the coefficients, additional equations should be considered
(3.4)∑i=1npk(xi)ai+∑j=1nDBpk(xj)ck=0,k=1,2,…,m.
Combining (3.3) and (3.4), the following matrix can be obtained:
(3.5)[ABTB0]{ac}={u0},
where the matrix element for A is R(rij), for B is pk(xi), and a,c are the coefficient vectors. Solving (3.5), the coefficients ai and ck can be gotten, and the approximation function is then given by (3.1). It should be noticed that the additional polynomials are not necessary if the RBF is strictly positive definite (SPD).


## 3.2. Typical RBFs

There are several advantages of the radial basis function (RBF) to be used for interpolation. First, there are generally few restrictions on the way the data prescribed and it can easily be applied in almost any dimension. Second, it has the naturally mesh-free property. And high-order accuracy is its other merit. Until now, many RBFs have been developed, and the most widely used global supported RBFs are listed in Table1.Table 1
Typical radial basis functions with dimensionless shape parameters.
Name
Expression
Shape parameter
Multiquadrics (MQ)
R
i
(
r
)
=
(
r
2
+
(
ε
l
)
2
)
q
ε
≥
0
,
q
Inverse multiquadrics (IMQ)
R
i
(
r
)
=
(
r
2
+
(
ε
l
)
2
)
-
q
ε
≥
0
,
q
Gaussian (GA)
R
i
(
r
)
=
exp
[
-
ε
(
r
/
l
)
2
]
ε
Thin plate spline (TPS)
R
i
(
r
)
=
r
η
η
Note:l is a characteristic length that relates to the nodal spacing in the computational domain.In these RBFs, the TPS is conditional positive definite, and others are SPD RBFs. Thus, when the TPS is used for interpolation, the polynomial terms should be employed to avoid the singularity.


## 3.3. Error Estimation and Stability of the RBF Interpolation

For the advantage of the RBF interpolation, it has been applied to a wide range of engineering problems. However, works concerning the error estimation and stability of this interpolation with mathematical theories are relatively less. The most noticeable work was completed by Schaback [17]. For the error estimation, the formulation is as follows:
(3.6)|u(x)-ua(x)|≤|u|f·P(x),
where P(x) is called the power function, which is the norm of the error function on F evaluated at x. And F is a native space, which can be described as
(3.7)F={f∈L2(Rn)∣∫Rd|u^(ω)|ψ(ω)dω<∞}.
In addition, u^ is in L1(Rd) and
(3.8)u|f2=∫Rd|u^(ω)|ψ(ω)dω<∞,
where ψ(ω) denotes the d-variate generalized Fourier transform of φ. Furthermore, Schaback proposed that the error and the sensitivity of the RBF interpolation can be described by P(x) and the eigenvalue λ of the interpolation matrix, which seems to be intimately related. P(x) and λ are bounded by F(h(x)) and G(h(x)) separately, the exact forms of which can be found in [17].And the shape parameters are directly contained inG(h(x)) but not in F(h(x)). Here, we list the G(h(x)) as follows:
(3.9)MQ:G(h)=he(-12.76εd/h),IMQ:G(h)=h-1e(-12.76εd/h),Guassian:G(h)=h-de(-40.71ε2d3/h2),
where h=maxx∈Ωminxi∈X∥x-xi∥, Ω is the interpolation domain, and X denotes the set of the interpolation points, d represents the dimensions, and ε is the shape parameter of the corresponding RBF.On the other hand, the error bound that contains the shape parameter was described in Madych’s theoretical analysis. His conclusions are proposed for two different classes of functions,Bσ and Eσ,
(3.10)Bσ={f∈L2(Rn):f~(ξ)=0if∥ξ∥>σ},Eσ={f∈L2(Rn):f~(ξ)<∞},
where σ is a positive constant,
(3.11)f~(ξ)=∫Rnf(x)e-i(x,ξ)dx
is the multivariate Fourier transform of function f(x), and
(3.12)∥f~∥Eσ2=∫Rn|f~(ξ)|2e|ξ|2/σdξ.
Considering the different functions, Madych presented the following error estimation for a class of RBF interpolations, which include Gaussian, MQ, and IMQ:
(3.13)forf∈Bσ:εr~O(eaελε/h)for0<λ<1anda>0,forf∈Eσ:εr~O(eaε2λε/h)for0<λ<1anda>0,
where εr is the error. It is worth noting that there is an intimate relation between the accuracy and the stability of the RBF interpolation. In the following section, we will present some strategies to balance the accuracy and the stability.


## 4. Variable Shaped RBF

### 4.1. The Role of the Shape Parameter in RBF Interpolation

Since the shape parameter in the RBF greatly affects the interpolation error and the stability, it has drawn many attentions of the researchers. For the MQ case, Hardy’s original work showed that the larger shape parameter caused the unstable results [25]; therefore, he advised the use of small shape parameter. However, contrary to Hardy’s belief, both  Tarwate  [26] and Kansa [27] observed that the accuracy of the MQ interpolation actually increases with increasing shape parameter. Furthermore, Huang et al. [28] established that the theoretically larger shape parameter often results in high accuracy. Nevertheless, as the shape parameter is continuously increasing, the interpolation function becomes flat. As the function flattened, it becomes insensitive to the distance r; thus, the elements of the interpolation matrix become nearly identical. The interpolation matrix may become ill-conditioned, which leads to unstable results.Once the role of the shape parameter is realized, many efforts followed to choose the optimal shape parameter for the RBF interpolation. Even though many algorithms have been developed, there is no widely accepted theory for choosing the proper shape parameter for various applications; therefore, it is still an open problem.In this paper, we try to use variable shape parameter strategies to improve the performance of RBF for interpolation.


### 4.2. Why RBF with Variable Shape Parameter

There are at least three advantages of using the RBF with variable shape parameter.(1)
It can improve the accuracy of the RBF interpolation. Kansa’s work showed that the accuracy of the RBF interpolation can be greatly improved if spatially variable shape parameters are employed. In addition, Sturgill proposed a random variable shape parameter strategy, which can greatly improve the numerical accuracy [24].(2)
The condition number of the interpolation matrix can be greatly reduced and the stability of the numerical methods based on RBF can be improved. In [20], the author presented several variation schemes for the RBF with variable shape parameters, and then introduced the variable shaped RBF into the dual reciprocity boundary face method. Numerical results showed the robustness and better stability of the new method.(3)
Runge errors encountered in the interpolation can be overcome by the spatially variable shape parameter. Fornberg detailed and described the Runge phenomenon in the RBF interpolation and demonstrated the effectiveness of the spatially variable shape parameters in dealing with the Runge errors both theoretically and numerically [29].


### 4.3. Variation Strategies

Heuristically, it has been demonstrated that using a variable shape parameter is an excellent idea. However, as mentioned in Section3, there is a contradiction between the accuracy and the sensitivity generally, which can also be described by Schaback as “uncertainly principle.” Thus, the variation strategies are crucially important to the performance of the variable shaped RBF. Until now, there are several variable shape parameter strategies that have been proposed. In [27], Kansa presented the following variation scheme:
(4.1)εj=[εmin2(εmax2εmin2)(j-1)/(N-1)]1/2j=1,…,N
to improve the accuracy of the MQ, where εj is the jth shape parameter in MQ, εmin and εmax are the minimum and maximum values of the shape parameter separately, and N is the number of the points used for interpolation. Furthermore, Kansa pointed out that the very accurate numerical result could be gotten if εmin and εmax varied by several orders of magnitude [22].In addtion, a linearly varying formula and a random variable shaped parameter strategy are also investigated.(
1
) The linearly varying formula:
(4.2)εj=εmin+(εmax-εminN-1)j,j=1,…,N.(
2
) The random variation scheme:
(4.3)εj=εmin+(εmax-εmin)×rand(1,N),j=1,…,N.
What is more, the strategy that one constant shape parameter applied in the domain and another constant shape parameter implemented on the boundaries were also suggested to guarantee the stability of the numerical approach. It should be noticed that for different RBFs, the variation strategies should be different. In this paper, we focus on the widely used MQ and IMQ.For MQ, we consider the upper bound ofλ-1, which can represent the sensitivity of the RBF interpolation, and the formulation is as follows:
(4.4)λ-1≤h-1e25.25ε/hin2D.
We propose the following two variation strategies for the MQ:
(4.5)εi=hi3,(4.6)εi=2hi2,
where hi=mini≠j∥xi-xj∥ and xi and xj denote the different center points for interpolation. Compared with the constant shape parameter, we substitute the expression (4.5) into (4.4); the new bound can be obtained
(4.7)λ-1≤h-1e25.25h2.
It should be noticed that h-1e25.25h2=O(h-1) when h goes to zero. In addition, the new error estimations as follows are obtained:
(4.8)forf∈Bσ:εr~O(eahλh),forf∈Bσ:εr~O(eah6λh).
In the same way, we can get the error estimations and the bound of the λ-1 for the variation strategy εi=2hi3. Furthermore, the balance between the accuracy and the stability may be found through the variation strategies (4.5) and (4.6).For the IMQ, the bound ofλ-1 is given by
(4.9)λ-1≤he25.25ε/hin2D.
We propose the following variation strategy for this RBF:
(4.10)εi=hi/5.
Substituting (4.10) into (4.9), we can get the new bound for the λ-1:
(4.11)λ-1≤he5.05.
Furthermore, the error estimations can be formulated by
(4.12)forf∈Bσ:εr~O(eah/5λ1/5),forf∈Eσ:εr~O(eah2/25λ1/5).
Through the above expressions, the balance between the accuracy and stability for the IMQ may be found. In order to verify the effectiveness of our variation strategy, we consider the following interpolation example:
(4.13)f=sinxsinysinzΩ:(x,y,z)∈[-1,1]×[-1,1]×[-1,1].
There are total N points that are used for interpolation. Both the MQ and IMQ with variation strategies are tested, and comparisons with the corresponding RBFs with constant shape parameters are shown in Figures 2–5. In addition, the parameters of constant-shaped RBFs are evaluated with the average value of all the variable parameters.Figure 2
Relative error versus the number of points for MQ; variable shapedMQ(1): εj=hi3; (2) variable shaped MQ(2): εj=2hi2.Figure 3
Condition number versus the number of points for MQ; variable shapedMQ(1): εj=hi3; (2) variable shaped MQ(2): εj=2hi2.Figure 4
Relative error versus the number of points for IMQ.Figure 5
Condition number versus the number of points for IMQ.The relative error of the numerical solution is defined as(4.14)e=∑i=1N(uiexact-uinum)2∑i=1N(uiexact)2,
where uiexact and uinum refer to the exact solution and the numerical solution at the point xi, respectively, and N is the number of the points used for interpolation. In order to measure the stability of the variable shaped RBF, condition number is employed, which can be defined as
(4.15)Cond(A)=∥A∥∥A-1∥,
where A is the interpolation matrix, A-1 and ∥A∥ are its inversion and norm separately.As shown in Figure2, when the number of the interpolation points is too small (N<7), the accuracy is not too good, which is similar to the FEM. However, when the number of the points increase, and N≥7, the accuracy can be greatly improved due to the properties of the MQ. Furthermore, the more points used for interpolation, the more accurate results can be obtained both for the constant shaped MQ and the variable shaped MQ, and the MQ with εi=2hi2 can get the best accuracy as the number of points increases. Figure 3 shows the condition numbers of the interpolation matrices with different RBFs. It can be observed that the condition number can be well controlled with the variable shaped MQ presented. Thus, we can figure out that the variable shaped MQ with the proposed variation strategies outperform the constant shaped MQ both in the accuracy and the stability.Figures3 and 4 show the comparisons of the constant shaped IMQ and the variable shaped IMQ in the accuracy and the stability. Even though the convergence of the constant shaped IMQ is better than the variable shaped IMQ, the condition number with the variable shaped IMQ grows much slower than the constant shaped IMQ.All the numerical results show that both the variable shaped MQ and IMQ with proper variation strategies can greatly reduce the condition number of the interpolation matrix and improve the performance of RBFs. At the same time, the accuracy of the MQ can also be improved with our strategies. Thus, in the following section, an improved collocation meshless method will be formulated by employing the variable shaped MQ.


## 4.1. The Role of the Shape Parameter in RBF Interpolation

Since the shape parameter in the RBF greatly affects the interpolation error and the stability, it has drawn many attentions of the researchers. For the MQ case, Hardy’s original work showed that the larger shape parameter caused the unstable results [25]; therefore, he advised the use of small shape parameter. However, contrary to Hardy’s belief, both  Tarwate  [26] and Kansa [27] observed that the accuracy of the MQ interpolation actually increases with increasing shape parameter. Furthermore, Huang et al. [28] established that the theoretically larger shape parameter often results in high accuracy. Nevertheless, as the shape parameter is continuously increasing, the interpolation function becomes flat. As the function flattened, it becomes insensitive to the distance r; thus, the elements of the interpolation matrix become nearly identical. The interpolation matrix may become ill-conditioned, which leads to unstable results.Once the role of the shape parameter is realized, many efforts followed to choose the optimal shape parameter for the RBF interpolation. Even though many algorithms have been developed, there is no widely accepted theory for choosing the proper shape parameter for various applications; therefore, it is still an open problem.In this paper, we try to use variable shape parameter strategies to improve the performance of RBF for interpolation.


## 4.2. Why RBF with Variable Shape Parameter

There are at least three advantages of using the RBF with variable shape parameter.(1)
It can improve the accuracy of the RBF interpolation. Kansa’s work showed that the accuracy of the RBF interpolation can be greatly improved if spatially variable shape parameters are employed. In addition, Sturgill proposed a random variable shape parameter strategy, which can greatly improve the numerical accuracy [24].(2)
The condition number of the interpolation matrix can be greatly reduced and the stability of the numerical methods based on RBF can be improved. In [20], the author presented several variation schemes for the RBF with variable shape parameters, and then introduced the variable shaped RBF into the dual reciprocity boundary face method. Numerical results showed the robustness and better stability of the new method.(3)
Runge errors encountered in the interpolation can be overcome by the spatially variable shape parameter. Fornberg detailed and described the Runge phenomenon in the RBF interpolation and demonstrated the effectiveness of the spatially variable shape parameters in dealing with the Runge errors both theoretically and numerically [29].


## 4.3. Variation Strategies

Heuristically, it has been demonstrated that using a variable shape parameter is an excellent idea. However, as mentioned in Section3, there is a contradiction between the accuracy and the sensitivity generally, which can also be described by Schaback as “uncertainly principle.” Thus, the variation strategies are crucially important to the performance of the variable shaped RBF. Until now, there are several variable shape parameter strategies that have been proposed. In [27], Kansa presented the following variation scheme:
(4.1)εj=[εmin2(εmax2εmin2)(j-1)/(N-1)]1/2j=1,…,N
to improve the accuracy of the MQ, where εj is the jth shape parameter in MQ, εmin and εmax are the minimum and maximum values of the shape parameter separately, and N is the number of the points used for interpolation. Furthermore, Kansa pointed out that the very accurate numerical result could be gotten if εmin and εmax varied by several orders of magnitude [22].In addtion, a linearly varying formula and a random variable shaped parameter strategy are also investigated.(
1
) The linearly varying formula:
(4.2)εj=εmin+(εmax-εminN-1)j,j=1,…,N.(
2
) The random variation scheme:
(4.3)εj=εmin+(εmax-εmin)×rand(1,N),j=1,…,N.
What is more, the strategy that one constant shape parameter applied in the domain and another constant shape parameter implemented on the boundaries were also suggested to guarantee the stability of the numerical approach. It should be noticed that for different RBFs, the variation strategies should be different. In this paper, we focus on the widely used MQ and IMQ.For MQ, we consider the upper bound ofλ-1, which can represent the sensitivity of the RBF interpolation, and the formulation is as follows:
(4.4)λ-1≤h-1e25.25ε/hin2D.
We propose the following two variation strategies for the MQ:
(4.5)εi=hi3,(4.6)εi=2hi2,
where hi=mini≠j∥xi-xj∥ and xi and xj denote the different center points for interpolation. Compared with the constant shape parameter, we substitute the expression (4.5) into (4.4); the new bound can be obtained
(4.7)λ-1≤h-1e25.25h2.
It should be noticed that h-1e25.25h2=O(h-1) when h goes to zero. In addition, the new error estimations as follows are obtained:
(4.8)forf∈Bσ:εr~O(eahλh),forf∈Bσ:εr~O(eah6λh).
In the same way, we can get the error estimations and the bound of the λ-1 for the variation strategy εi=2hi3. Furthermore, the balance between the accuracy and the stability may be found through the variation strategies (4.5) and (4.6).For the IMQ, the bound ofλ-1 is given by
(4.9)λ-1≤he25.25ε/hin2D.
We propose the following variation strategy for this RBF:
(4.10)εi=hi/5.
Substituting (4.10) into (4.9), we can get the new bound for the λ-1:
(4.11)λ-1≤he5.05.
Furthermore, the error estimations can be formulated by
(4.12)forf∈Bσ:εr~O(eah/5λ1/5),forf∈Eσ:εr~O(eah2/25λ1/5).
Through the above expressions, the balance between the accuracy and stability for the IMQ may be found. In order to verify the effectiveness of our variation strategy, we consider the following interpolation example:
(4.13)f=sinxsinysinzΩ:(x,y,z)∈[-1,1]×[-1,1]×[-1,1].
There are total N points that are used for interpolation. Both the MQ and IMQ with variation strategies are tested, and comparisons with the corresponding RBFs with constant shape parameters are shown in Figures 2–5. In addition, the parameters of constant-shaped RBFs are evaluated with the average value of all the variable parameters.Figure 2
Relative error versus the number of points for MQ; variable shapedMQ(1): εj=hi3; (2) variable shaped MQ(2): εj=2hi2.Figure 3
Condition number versus the number of points for MQ; variable shapedMQ(1): εj=hi3; (2) variable shaped MQ(2): εj=2hi2.Figure 4
Relative error versus the number of points for IMQ.Figure 5
Condition number versus the number of points for IMQ.The relative error of the numerical solution is defined as(4.14)e=∑i=1N(uiexact-uinum)2∑i=1N(uiexact)2,
where uiexact and uinum refer to the exact solution and the numerical solution at the point xi, respectively, and N is the number of the points used for interpolation. In order to measure the stability of the variable shaped RBF, condition number is employed, which can be defined as
(4.15)Cond(A)=∥A∥∥A-1∥,
where A is the interpolation matrix, A-1 and ∥A∥ are its inversion and norm separately.As shown in Figure2, when the number of the interpolation points is too small (N<7), the accuracy is not too good, which is similar to the FEM. However, when the number of the points increase, and N≥7, the accuracy can be greatly improved due to the properties of the MQ. Furthermore, the more points used for interpolation, the more accurate results can be obtained both for the constant shaped MQ and the variable shaped MQ, and the MQ with εi=2hi2 can get the best accuracy as the number of points increases. Figure 3 shows the condition numbers of the interpolation matrices with different RBFs. It can be observed that the condition number can be well controlled with the variable shaped MQ presented. Thus, we can figure out that the variable shaped MQ with the proposed variation strategies outperform the constant shaped MQ both in the accuracy and the stability.Figures3 and 4 show the comparisons of the constant shaped IMQ and the variable shaped IMQ in the accuracy and the stability. Even though the convergence of the constant shaped IMQ is better than the variable shaped IMQ, the condition number with the variable shaped IMQ grows much slower than the constant shaped IMQ.All the numerical results show that both the variable shaped MQ and IMQ with proper variation strategies can greatly reduce the condition number of the interpolation matrix and improve the performance of RBFs. At the same time, the accuracy of the MQ can also be improved with our strategies. Thus, in the following section, an improved collocation meshless method will be formulated by employing the variable shaped MQ.


## 5. The Improved Collocation Meshless Method

### 5.1. Collocation Scheme

According to the formulation procedures of the governing equations, various meshless methods can be divided into three types: meshless methods based on weak form, meshless methods based on strong form (collocation form), and meshless methods based on the combination of the weak and strong form. Considering the straightforward and truly meshless  method  of the collocation form, it is introduced into the approach of the interior acoustic problems. In order to demonstrate its application in the solution of the boundary value problems, the following equation is considered:(5.1)A∂2u∂x2+B∂2u∂y2+C∂2u∂x∂y+D∂u∂x+E∂u∂y+Fu=G,Ω,
with the boundary conditions:Dirichlet boundary(5.2)u=HonΓD,
Neumann boundary
(5.3)∂u∂n=IonΓN.Here,Ω is the computational domain. A, B, C, D, E, F, G, H, and I can be constant or depend on x and y. Consider that there are nΩ points that are distributed in the domain, nD points on the Dirichlet boundary, and nN points on the Neumann boundary.By substitutingu in (5.1) with the approximation u(xi) expressed in (3.1), the following nΩ linear equations can be obtained:
(5.4)Ai∂2u(xi)∂x2+Bi∂2u(xi)∂y2+Ci∂2u(xi)∂x∂y+Di∂u(xi)∂x+Ei∂u(xi)∂y+Fiu(xi)=Gi,
wherei=1,2,…,nΩ, xi denotes the ith points in the domain Ω and note that the xi=(xi,yi). In the same way, the following nD+nN equations are satisfied on the Dirichlet boundary ΓD and the Neumann boundary:
(5.5)u(xi)=Hii=1,2,…,nD,∂u(xi)∂n=Iii=1,2,…,nN.


### 5.2. Treatment for Neumann Boundary Conditions

It is well known that meshless methods based on collocation scheme are weak in dealing with Neumann boundary conditions; in order to overcome this issue, a Hermite-type interpolation is employed [30]. The formula is as follows:
(5.6)u(x)=∑i=1nRi(x)ai+∑j=1nN∂RjN(x)∂nbj+∑k=1mpk(x)ck.Here, theRi(x) and RjN(x) denote the RBF, and pk(x) is the polynomial terms, with which a linear field can be reproduced exactly and the accuracy of the interpolation can be improved in most cases. n is the total number of the points both in the domain and on the boundaries for interpolation. nN represents the number of the interpolation points on the Neumann boundaries and m is the number of the polynomial. In addition, ai, bj, and ck are the constant coefficients that need to be determined. n is the vector of the unit outward normal on the Neumann boundaries, and in (5.6), the derivative terms can be obtained through the following expression:
(5.7)∂RjN(x)∂n=lxj∂RjN(x)∂x+lyj∂RjN(x)∂y,
where lxj=cos(n,x) and lyj=cos(n,y) denote the direction cosines in x and y directions separately.In (5.6), there are totally n+nN+m unknown qualities that need to be determined; however, by enforcing the interpolation functions pass through n+nN points both in the domain and on the boundaries, only n+nN equations can be gotten. Therefore, to ensure the unique solution of the coefficients, the following constraint condition is added:
(5.8)∑i=1npk(xi)ai+∑j=1nNpk(xj)bj=0,k=1,2,…,m.Here,xi and xj represent the points in the domain and on the Neumann boundary separately.


### 5.3. The Interpolation Function with Sine/Cosine Basis for the Helmholtz Equations

As mentioned in Section3, the additional polynomial terms are not always necessary for the interpolation if the RBF is strictly positive definite (SPD), such as MQ and IMQ. Nevertheless, considering the property of the Helmholtz equation, we add sine and cosine expressions into the RBF, which are better suitable to reproduce the wave characteristics of the Helmholtz equation. Therefore, the polynomial terms are substituted by the following expression:
(5.9)p(x)=[1kxcosθ+kysinθkxsinθ+kycosθ],
where k is the wave number and θ denotes the angle of propagation.


## 5.1. Collocation Scheme

According to the formulation procedures of the governing equations, various meshless methods can be divided into three types: meshless methods based on weak form, meshless methods based on strong form (collocation form), and meshless methods based on the combination of the weak and strong form. Considering the straightforward and truly meshless  method  of the collocation form, it is introduced into the approach of the interior acoustic problems. In order to demonstrate its application in the solution of the boundary value problems, the following equation is considered:(5.1)A∂2u∂x2+B∂2u∂y2+C∂2u∂x∂y+D∂u∂x+E∂u∂y+Fu=G,Ω,
with the boundary conditions:Dirichlet boundary(5.2)u=HonΓD,
Neumann boundary
(5.3)∂u∂n=IonΓN.Here,Ω is the computational domain. A, B, C, D, E, F, G, H, and I can be constant or depend on x and y. Consider that there are nΩ points that are distributed in the domain, nD points on the Dirichlet boundary, and nN points on the Neumann boundary.By substitutingu in (5.1) with the approximation u(xi) expressed in (3.1), the following nΩ linear equations can be obtained:
(5.4)Ai∂2u(xi)∂x2+Bi∂2u(xi)∂y2+Ci∂2u(xi)∂x∂y+Di∂u(xi)∂x+Ei∂u(xi)∂y+Fiu(xi)=Gi,
wherei=1,2,…,nΩ, xi denotes the ith points in the domain Ω and note that the xi=(xi,yi). In the same way, the following nD+nN equations are satisfied on the Dirichlet boundary ΓD and the Neumann boundary:
(5.5)u(xi)=Hii=1,2,…,nD,∂u(xi)∂n=Iii=1,2,…,nN.


## 5.2. Treatment for Neumann Boundary Conditions

It is well known that meshless methods based on collocation scheme are weak in dealing with Neumann boundary conditions; in order to overcome this issue, a Hermite-type interpolation is employed [30]. The formula is as follows:
(5.6)u(x)=∑i=1nRi(x)ai+∑j=1nN∂RjN(x)∂nbj+∑k=1mpk(x)ck.Here, theRi(x) and RjN(x) denote the RBF, and pk(x) is the polynomial terms, with which a linear field can be reproduced exactly and the accuracy of the interpolation can be improved in most cases. n is the total number of the points both in the domain and on the boundaries for interpolation. nN represents the number of the interpolation points on the Neumann boundaries and m is the number of the polynomial. In addition, ai, bj, and ck are the constant coefficients that need to be determined. n is the vector of the unit outward normal on the Neumann boundaries, and in (5.6), the derivative terms can be obtained through the following expression:
(5.7)∂RjN(x)∂n=lxj∂RjN(x)∂x+lyj∂RjN(x)∂y,
where lxj=cos(n,x) and lyj=cos(n,y) denote the direction cosines in x and y directions separately.In (5.6), there are totally n+nN+m unknown qualities that need to be determined; however, by enforcing the interpolation functions pass through n+nN points both in the domain and on the boundaries, only n+nN equations can be gotten. Therefore, to ensure the unique solution of the coefficients, the following constraint condition is added:
(5.8)∑i=1npk(xi)ai+∑j=1nNpk(xj)bj=0,k=1,2,…,m.Here,xi and xj represent the points in the domain and on the Neumann boundary separately.


## 5.3. The Interpolation Function with Sine/Cosine Basis for the Helmholtz Equations

As mentioned in Section3, the additional polynomial terms are not always necessary for the interpolation if the RBF is strictly positive definite (SPD), such as MQ and IMQ. Nevertheless, considering the property of the Helmholtz equation, we add sine and cosine expressions into the RBF, which are better suitable to reproduce the wave characteristics of the Helmholtz equation. Therefore, the polynomial terms are substituted by the following expression:
(5.9)p(x)=[1kxcosθ+kysinθkxsinθ+kycosθ],
where k is the wave number and θ denotes the angle of propagation.


## 6. Solving the Interior Acoustic Problems by the Improved Collocation Meshless Method

To verify the efficiency of the proposed method in solving the interior acoustic problems, two numerical examples are carried out. The method is first applied to a problem with anL-shape domain. Then an irregular domain with complex boundary conditions is further investigated. Since both of the two examples have no analytical solutions, comparisons with the refined FEM solutions are made. Considering the better performance of the variable shaped MQ with the εi=2hi2 variation strategy, it will be employed in the proposed method.Example 6.1 (L-shape interior acoustic problem).
As shown in Figure6, an L-shape acoustic domain is considered, a harmonic plane wave (p=p0e-jwt) is applied on the boundary A, and on all other boundaries impedance boundary condition (∂p/∂n=jkp) is considered. L=1 and p0=1 are considered in this example. To approach the problem with the proposed method, total 100 points are employed, including 40 of them on the boundaries. Since there is no analytical solution for the problem, refined FEM solutions with 2500 elements are presented as a reference. It should be noticed that the number of the FEM elements should satisfy the rule of thumb of at least 10 grid points per wave length.
Figure7 shows a pressure variation of a fixed field point (0.4, 0.6) with respect to the wave number. The good agreement with the reference solution can be seen from the figure. We can also figure out that when k≥17, the solution becomes a little inaccurate. The reason for this is that there is a pollution effect caused by dispersion in solving the Helmholtz equations with numerical methods for the high wave number. The relatively small number of the interpolation points is another reason for the error. As described, only 100 points are used for interpolation and 2500 elements are employed in the reference solution. Furthermore, our numerical tests showed that very accurate results (for k=17~20) can be obtained by increasing very few additional interpolation points. In addition, we believe that the accuracy can also be improved by using better-shape parameter strategies, which is an ongoing work.
In addition, another two tests with different wave numbers (k=10, 20) are conducted. Figures 8 and 9 show the numerical solution along the diagonal of the computational domain (x=y). Comparison with the reference solutions shows the efficiency of the proposed method in dealing with the interior acoustic problems. Even though relatively few points are used for interpolation, very accurate results are obtained.Figure 6
TheL-shape interior acoustic problem.Figure 7
Comparison of the numerical solutions at the field point(0.3,0.5) with respect to wave number.Figure 8
Comparison of numerical solutions along the diagonalx=y for k=10 (545.9 Hz).Figure 9
Comparison of numerical solutions along the diagonalx=y for k=20 (1091.8 Hz).Example 6.2 (irregular shape domain with complex boundary conditions).
As well known, the boundary conditions play an important role in the boundary value problems. Just as described in Section2, there are generally three kinds of boundary conditions for the interior acoustic problems. In order to verify the effectiveness of the proposed method in dealing with these boundaries, an irregular shape domain with complex boundary conditions is considered.
The geometry is shown in Figure10, AB=1.8, BC=0.4, DE=1.0, and EA=0.8, the unit is meter. The impedance boundary conditions are imposed on the AB and DE, where Z=1.25×343(Pa·s/m). On the CD and AE, sound-hard boundary and sound-soft boundary are considered separately. In addition, the sound source is on the BC (p=p0e-jwt), where p0=1 in this example. Including 60 points on the boundaries, there are totally 400 points arranged for simulation. For comparison, the FEM solutions with 5000 elements are also given. Figure 11 shows a pressure variation of a fixed field point (0.6, 0.4) with respect to the wave number. The same phenomenon can be seen that when k≥17, the solution becomes a little inaccurate. The reasons for this are the same as the ones described in Example 6.1.
The comparisons of the numerical solutions fork=10 and k=20 are also illustrated in Figures 12–13. All the numerical results show that the proposed method can well handle the boundary conditions and very accurate solutions can be obtained.Figure 10
The computational domain for the Example6.2.Figure 11
Comparison of the numerical solutions at the field point(0.6and0.4) with respect to wave number.Figure 12
Comparison of numerical solutions along the oblique liney=-4x/9+4/5 for k=10 (545.9 Hz).Figure 13
Comparison of numerical solutions along the linex(y=0.4) for k=20 (1091.8 Hz).


## 7. Conclusion

In this paper, several new variation strategies are proposed for the RBF with variable shape parameters in order to improve the performance of the RBF interpolation. Numerical results show that the variable shaped RBF can usually improve the stability of the RBF interpolation. What is more, the accuracy can also be improved with some proper variation strategies.Then, an improved collocation meshless method is formulated by employing the new variable shaped RBF. Considering the property of the Hemholtz equation, a sine/cosine basis is presented for interpolation. Then interior acoustic problems with different geometries and different boundary conditions are solved with the present method. Since there are no analytical solutions for the examples, the refined FEM solutions are considered as reference solutions. Compared with the FEM solutions, even though relatively few points are arranged, very accurate results can be obtained with the proposed method. In addition, points for interpolation can be chosen freely; thus, it can easily be extended to more complex practical problems.



---

*Source: 632072-2012-10-17.xml*