Numerical Simulation of Acoustic Equation Using Radial Point Interpolation Method with Discontinuous Galerkin Time Integration

The Numerical methods are research and industrial strategies commonly used by the finite difference method (FDM), finite element method (FEM) and finite volume method (FVM). The technique is a mesh based or formation of the domain. Owing to the complicated and time-consuming nature of the mesh method in the complex domain, it encounters numerous inconsistencies. One way of eluding this is by the use of a meshless method. This technique eliminates the use of but rather makes use of nodes in the distribution of its domain. This paper introduces the use of the radial point interpolation method (RPIM) to approximate the acoustic equations using the discontinuous Galerkin method (DGM) time integration. In order to determine the numerical behaviour, its results were simulated with the exact solution. The DGM time integration and order of accuracy is also compared with some commonly used procedures, such as the backward Euler and trapezoid methods. The size of support domain responsible for the numerical accuracy is also examined. Finally a comparison of numerical simulations of the exact results obtained during a specific time snapshot is displayed.


Introduction
In recent times, the meshless method has attracted lots of attention from researches seeking for ways to develop and utilize it for numerical modelling in engineering. The advantage of this method compared to the conventional predecessors such as finite difference method (FDM), finite element method (FEM) and finite volume method (FVM) is that it does not make use of mesh. The use of mesh in solving complex space domain shape is complicated and time consuming [1,2].
One of the earliest commonly used in astrophysics simulations is smoothed particle hydrodynamics (SPH), [3]. However, Diffuse-element method (DEM) is the first meshless procedure that makes use of moving least square on the Galerkin method [4]. EFG uses the principle of moving least square (MLS) on Galerkin [1,5]. There is also a reproducing kernel particle method (RKPM) which uses kernel approximation [6].
The above mentioned procedures make use of the moving-least-square principle, which lacks the Kronecker delta properties. This causes difficulty in applying boundary condition (BC). To remediate the difficulty, several methods are proposed, some of which are the Lagrange multipliers [1], penalty method [7], coupling with FEM [8], collocation [9] and point interpolation method (PIM) [10]. PIM possesses the property of Discontinuous Galerkin Time Integration approximation function which passes through each node in the support domain, and shape function. It also possesses the Kronecker delta properties, which eliminates any form of difficulties inherent the application of boundary condition to PIM [9]. In the beginning, PIM used polynomial as the basis function, but the appearance of singularity inspired the usage of radial function as its basis, thereby, converting it into the radial point interpolation method (RPIM) [10,11].
While numerical modelling is used for time-dependent cases, the meshless procedures are used for space-discretization. As for time-discretization, either explicit or implicit method should be used. One of the difficulties associated with the forward Euler and Runge-Kutta techniques, it is conditional stable. To maintain stability, the width of time step (t) should be smaller, thereby, rendering the method efficient for numerical calculations. Furthermore, there is the implicit method such as the backward Euler, and the trapezoidal methods (Crank-Nicholson method), which are unconditionally stable and as such requires greater time step [12]. The use of meshless discretization and Crank-Nicholson method was carried out by [13]. The procedure was used to calculate dissipation of consolidation, and the Biot process [14]. Nevertheless, the choice of a large time step could lower the accuracy of the calculation. Unfortunately, the accuracy of backward Euler method is merely order-1, while Crank-Nicholson is order-2, for this reason, the accuracy order of the implicit methods needs to be developed, to improve its accuracy. One of the high accuracy implicit techniques used for time integration is the discontinuous Galerkin method (DGM) [15]. This paper will discuss the development of a time and space discretization method between radial point interpolation method (RPIM) and discontinuous Galerkin method (DGM). The model equation used is the one-dimensional acoustic equation. Derivation of weak form and time integration was being done. Simulation was carried out by varying the number of nodes, time step, and the size of support domain. The result of discretization calculation using discontinuous Galerkin method (DGM) is compared to several other standard implicit techniques such as the Euler backward method, the trapezoidal method and exact solution, to create a resulting numerical behaviour.

Acoustic Equation
Acoustic equation, also known as acoustic wave or scalar wave for one dimension, is presented on Eq. 1. This equation is used to describe sounds in one dimension or 1D. It is also similar to the SH wave problem [2]. This equation possesses an exact solution, which could be used to test the numerical solutions being used.
Where is the dependent variable, which in this case is pressure, is the constant of acoustic velocity, is the space axis and represents time. Eq. 1 possesses the exact solution which is expressed on Eq. 2 [2].
Simulation was carried out on initial condition (IC) in the form of Gaussian function, which exact solution would appear as Eq. 3:

Radial Point Interpolation Method (RPIM)
RPIM is a method of approach where the interpolation function passes through nodes spread throughout the support domain. If the value point (commonly Gauss integration points) is being approached by RPIM, its accuracy is being influenced by the number of nodes in the support domain. To calculate the dimension of the support domain on point , Eq. 4 should be utilized.

= (4)
Where the dimension of the support domain is, is the size of support domain expressed in dimensionless quantity, is the average of distance between nodes. RPIM interpolates ( ) using nodal value on nodes inside support domain on point according to Eq. 5 [16].
Where ( ) is constant, is the number of nodes on support domain from point (Gauss integration point) and basis function is expressed in the form of Eq. 6.
Basis function for one dimension on axis uses multiquadratic (MQ) type as expressed in Eq. 7.
Constant ( ) is the result of interpolations through all nodes on support domain, which is expressed as matrix form in Eq. (9).
Where moment matrix is expressed as a matrix in Eq. (10).
The value of is expressed in Eq. 11.
The value of ( ) was obtained by reversing Eq. 9 which result to Eq. 12.
Where the shape function is expressed in Eq. 14.
Weak Form and Spatial Discretization To anticipate the time integration method using the backward Euler method, trapezoidal method and DGM, the used weak form is adjusted to produce first order ordinary differential equation (ODE) [17]. Numerical solution using RPIM on acoustic equation was carried out by changing Eq. 1 into two first-order equations, which are respectively being displayed in Eq. 15a and Eq. 15b respectively, and by introducing a new variable, ( ).
The RPIM solution requires the approximation of the RPIM on the involved variables according to Eq. 16 and Eq. 17.
Where is the number of nodes on the support domain. By approximating the unknown and in Eq. 16 and Eq. 17 using ℎ and ℎ and by mutliplying and integrating it with test function and mesh Ω (as linear integral), we obtained Eq. 18a and Eq. 18b.
To obtain the weak form, Eq. 18a and Eq. 18b was undertaken to reduce the order of the equations, which produces Eq. 19a and Eq. 19b.
Eq. 20a and Eq. 20b could be expressed by matrix notations as displayed on Eq. 21a and Eq. 21b.

Time Integration
Time integration is carried out in standard form using backward Euler method, trapezoid method and discontinuous Galerkin method (DGM). It was carried out by multiplying each side Eq. 21a and Eq. 21b with − , thus we obtained Eq. 28a and Eq. 28b.
Using matrix notations, Eq. 29a and Eq. 29b could be expressed as in Eq.30, where is the identity matrix.
The RPIM simulation with backward Euler is carried out by solving Eq. 30. The process is done by determining the value of and at time , after which Eq. 30 is inverted, thereby, generating and values.
Step time + 1 using an unknown variable is repeated to obtain the final time.
Time integration using other implicit methods other than backward Euler could be conducted using the trapezoidal method. Here time integration for Eq. 28a and Eq. 28b is being carried out using an average gradient time step and + 1, as displayed in Eq. 31a and Eq. 31b.
Eq. 31a and Eq. 31b could be expressed in matrix as seen on Eq. 32.
The value of and are unknown in Eq. 33a and Eq. 33b with an independent variable time of . The value of and were approximated by Eq. 34 and Eq. 35, giving rise to an independent variable with coordinate .
Thus Eq. 37a and Eq. 37b could be solved for 1 = 1 and 2 = 2 , and expressed in matrix as seen on Eq. 38.

Numerical Example
This chapter will show the solutions of 1-dimensional acoustic equation using RPIM-DGM and first-order time integration. Simulation was carried out using RPIM and time integration with backward Euler process, trapezoidal method, and DGM. Simulation on axis was carried out with ≤ ≤ where = −8.0 (left) and  In Fig. 1 to Fig. 3 the simulation result displays the order of accuracy from three different time integration. The order of accuracy indicates the rate of convergence from numerical method into an exact solution. The simulation results which indicates the order of accuracy on varying number of nodes relative to log 10 ‖error‖ 2 where error = exact -numeric is being displayed from Fig. 1 to Fig. 3.
According to these figures, the rate of convergence of RPIM-DGM is higher than that of RPIM-trapezoid and RPIM-backward Euler for time increments Δ = 0.5; 0.2; and 0.1. The results of the rate of convergence for ∆t = 0.1, RPIM-Euler backward is 0.95, that of the RPIM-trapezoid is 2.58, while that of the RPIM-DGM is 3.01 as shown in Table 1. For ∆t = 0.2, RPIM-Euler backward is 0.22, RPIM-trapezoid is 1.81 and RPIM-DGM is 3.05. At ∆t = 0.5, the backward RPIM-Euler is 0.01, the RPIM-trapezoid is 0.48, and the RPIM-DGM has a value of 1.75. Table 1 shows that for a larger time increment, order of accuracy will diminish. On the other hand, smaller time increment will result in higher order of accuracy.
The result of error accumulation relative to time for RPIM-backward Euler, RPIM-Euler trapezoid and RPIM-DGM is displayed from Fig. 4 to Fig. 6. From these figures it could be observed that simulations with RPIM-backward Euler, results to huge errors, while that of RPIM-trapezoid gives a smaller accumulation error, with RPIM-DGM having the smallest accumulation error rate. Smaller accumulation error rate will result in a gentle graphic line slope.  The accuracy of numerical simulation using RPIM-DGM approach is also influenced by the size of the support domain s as displayed in Fig. 7, 8, and 9. A minimum and maximum support domain size will result in larger error.
For Δ = 0.5 with 20 nodes, the maximum accuracy for s ranges from 7 to 10, while for 30 and 50 nodes, the maximum accuracy for s ranges from 2 to 12, as shown in Fig. 7. For Δ = 0.2 with 20, 30 and 50 nodes, the maximum accuracy for s ranges from 5 to 11, as shown in Fig. 8. While for Δ = 0.1 with 20, 30 and 50 nodes, the maximum accuracy for s ranges from 7 to 11, as shown in Fig. 9. In general, the more the number of nodes used, the higher the accuracy in the wider s range. The calculation result of RPIM-DGM which takes the form of acoustic wave is being compared directly to the exact solution and displayed for snapshot on time =0; 1; 2; 3; and 4 as displayed in Figure 10 to 14. The pressure which is the calculated result of RPIM-DGM is displayed as dotted lines, while the exact solution is illustrated using solid lines.
It could be observed from Fig.10 to 14 that the lines representing both RPIM-DGM with coincident solution. From Fig.10 to 14, it is observed that stability problems do not occur with the integration of time with DGM.

Conclusions
A numerical approach with RPIM and time integration using DGM or RPIM-DGM has been carried out to approximate a one-dimensional acoustic equation in order to produce good results. The RPIM-DGM numerical simulation is also compared to other commonly used time integration approaches such as Euler backward and trapezoid methods. In addition, the numerical method is compared with the exact method so that its behaviour can be identified. The results of accumulated errors for each method show different values. The RPIM-DGM results produces the smallest accumulated number of errors, while the biggest accumulated number of errors is produced using the RPIM-Euler backward technique. The size of the domain support affects the accuracy of the RPIM-DGM simulation. The optimal size of the support domain is 7 to 10 from the average distance between nodes. The numerical simulation results show that the simulation does not encounter stability inconsistencies.