Construction of Compact Finite Difference Schemes by Classic Differential Quadrature

Using classic differential quadrature formulae and uniform grids, this paper systematically constructs a variety of high-order finite difference schemes, and some of these schemes are consistent with the so-called boundary value methods. The derived difference schemes enjoy the same stability and accuracy properties with correspondent differential quadrature methods but have a simpler form of calculation; thus, they can be seen as a compact format of classic differential quadrature methods. Through systematic Fourier stability analysis, the characteristics such as the dissipation, dispersion and resolution of the different schemes were studied and compared.


Introduction
Many physical phenomena possess a range of space and time scales, and direct numerical simulations of these processes require all the relevant scales to be properly represented in the numerical model.These requirements have led to the development of the so-called compact finite difference scheme [1].Usually, the compact schemes have much smaller numerical dispersion and dissipation errors than traditional finite difference schemes of the same order of accuracy on the same mesh.So far, researchers have developed a variety of compact difference schemes [2][3][4][5][6].
This paper presents a simple way for constructing the compact finite difference schemes by the use of the differential quadrature method (DQM).The DQM was firstly proposed by Bellman and Casti in the early 1970s [7,8].The basic idea of DQM is to extend analogously the numerical integration rule to the numerical differentiation, where the derivative of a function on a coordinate direction is expressed as the weighted linear combination of function values for all discrete points along this direction.
At present, the differential quadrature method has been widely developed.According to the used trial or test function, the differential quadrature method generally includes the classical differential quadrature method, harmonic differential quadrature method, Fourier differential quadrature method, and so on.The classic differential quadrature method mainly uses the general polynomial or Lagrangian interpolation function as a test function, the harmonic differential quadrature method uses the trigonometric function as the trial function, while the Fourier differential quadrature method uses Fourier series expansion as the trial function [9,10].Theoretically, the above differential quadrature rules can be used either for spatial discretization or for time discretization.In other words, the differential quadrature method can be used for the solution of both the partial differential equation and the ordinary differential initial problem [11][12][13][14][15].In summary, differential quadrature methods have the advantages of high precision, good numerical stability and convenient programming, and have been successfully applied in many fields [16,17].
Using classic differential quadrature formulae and uniform grids, we systematically constructed a variety of finite difference schemes, and some of these schemes are consistent with the boundary value methods (BVM) proposed by Brugnano et al. [18,19].Considering that the derived difference schemes enjoy the same stability and accuracy properties with correspondent DQM but have a simpler form of calculation, the derived difference schemes can be regarded as the compact expression of classic DQM.Naturally, these compact difference methods could be used instead of classic differential quadrature methods for engineering applications [20].

The Classic Differential Quadrature Methods
Suppose function f (x) is sufficiently smooth in the considered interval, and the first-order derivative of f (x) at each grid point c p , p ∈ (1, s) is approximated by a linear sum of all the function values in the whole domain, that is The above Equation ( 1) is simply the basic expression of DQM, within which g pq represents the weighting coefficient of DQM; s is the grid number of the partition of each interval; and c p , p ∈ (0, s) represent the grid points.
The weighting coefficients of DQM only depend on the trial function and the discrete sample points, namely, the distribution of grid points, and have nothing to do with a specific problem.Usually, we would like to choose the general polynomial or Lagrange interpolation polynomial as the trial function, and use these functions to approximate the value of f (1) (x) at each grid points, which is defined as the classic DQM in this paper.The researchers have proved that the weighting coefficients derived by these two trial functions are consistent [17].The concrete calculation expression of weighting coefficients for classic DQM can be described as follows: where, Obviously, if the distribution of grid points is centro-symmetric in the normalized interval [0, 1], it can be deduced that it can be known from Equation (5) that the weighting coefficient matrix G is centro-antisymmetric.
Regarding the selection of grid points, there are some commonly used approaches, including the uniform grids, Chebyshev grids, Chebyshev-Gauss-Lobatto grids, and Legendre grids.These types of grids can be concretely described as follows.
The Legendre grids are the zero points of d s dx s x s (x − 1) s = 0; and c s+1 = 1.
Take the following ordinary differentia initial value problem as an example Supposing x n and x n+1 represent respectively the starting point and ending point of an arbitrary interval, h = x n+1 − x n is the interval length or the so-called stencil length.Let c = (x − x n )/h and make the interval x ∈ [x n , x n+1 ] normalized as c ∈ [0, 1], then the Equation ( 10) can be converted into On this basis, we can apply the s-stage differential quadrature method (1) to Equation (11), and deduce that [15,16] . . .
in which y n+z = y(x n + c z h), z ∈ (0, s), If the initial value of y n is given, the value of state variable y at different grid points, that is y n+z ∈ (1, s), can be obtained by solving the aforementioned nonlinear Equation (12).
It should be noted that, although the ordinary differential Equation ( 10) is used as a trial case, the variable x in this equation can either represent the space or the time, and naturally the stencil length h can be either a spatial discrete step size or a time discrete step.For ordinary differential initial value problems, it has been proven in reference [15] that, when using uniform grids, Chebyshev grids and Chebyshev-Gauss-Lobatto grids, the classical DQM-i.e., Equation ( 12)-are all s-stage s-order methods with A-stability, while the classic DQM combined with Legendre grids can achieve the accuracy of 2s-2 order.

The Generalized Backward Differential Formulae
As shown in Figure 1, we divide the interval x ∈ [x n , x n+1 ] into s segments, namely the grids, in which h = x n+1 − x n , x(x n + c 0 h) ≡ x n , x(x n + c s h) ≡ x n+1 .

The Generalized Backward Differential Formulae
As shown in Figure 1, we divide the interval [ ] into s segments, namely the grids, in which

Defining
( 2) 2, for even ( 1) 2, for odd and taking the υ -th equation from Equation ( 12), we can obtain 0 ( ) in which, , (0, ) It can be shown that the above numerical method (15) has the accuracy of s-order.
The numerical method (15) has a similar form to the traditional backward difference formulae (BDF), but we will later show that this method has the better numerical stability then traditional BDF.Thus, for the sake of simplicity, we define this numerical method as the generalized backward difference formulae (GBDF).
When selecting the uniform grids, the coefficient values ( z α ) of GBDF are shown in the following Tables 1 and 2.
Table 1.The coefficients of generalized backward difference formulae (GBDF) using uniform grids (for even s).Defining υ = (s + 2)/2, for even s (s + 1)/2, for odd s (14) and taking the υ-th equation from Equation ( 12), we can obtain in which, It can be shown that the above numerical method (15) has the accuracy of s-order.
The numerical method ( 15) has a similar form to the traditional backward difference formulae (BDF), but we will later show that this method has the better numerical stability then traditional BDF.Thus, for the sake of simplicity, we define this numerical method as the generalized backward difference formulae (GBDF).
When selecting the uniform grids, the coefficient values (α z ) of GBDF are shown in the following Tables 1 and 2.
Table 1.The coefficients of generalized backward difference formulae (GBDF) using uniform grids (for even s).Table 2.The coefficients of GBDF using uniform grids (for odd s).
−5 For the uniform grids, we can define that therefore, the Equation ( 15) can be rewritten as the following formula It can be verified from Tables 1 and 2 that, the numerical method ( 18) is consistent with the GBDF proposed by Brugnano et al. [19].In reference [21], it has been shown that the methods (18) are both A-stable numerical methods for time-domain numerical integration.

The Central Compact Difference Schemes
As is shown in Figure 1, let s be an even number.By defining and selecting the υ-th equation in Equation ( 12), it follows that in which, It can be shown that the above numerical method (20) has the accuracy of s-order.Furthermore, it can be verified that the coefficients α z , z ∈ (0, s) are centro-antisymmetrical about α υ .If the distribution of grid points is centro-symmetrical in the normalized interval [0, 1].So, the numerical method (20) can be called the central difference scheme.
When selecting the uniform grids, the coefficient values (α z ) of the method (20) are shown in the following Table 3.
Table 3.The coefficients of central difference schemes using uniform grids.

The Generalized Adams Methods
As is shown in Figure 2, let x n−1 < x < x n , and divide the interval then we can obtain the following formula according to Equation ( 12) in which ( 1) ( 1) ) , ( 0, ) n z z y y x c h z s In additional, define Since it has been proved that [12,15] where, e is the unit column vector, therefore, the Equation ( 23) can be rewritten as Subtracting the υ -th equation from the It can be shown that the above numerical method (27) has the accuracy of ( 1) + s -order.In additional, define Since it has been proved that [12,15] where, e is the unit column vector, therefore, the Equation ( 23) can be rewritten as Subtracting the υ-th equation from the (υ + 1)-th equation in the matrix Equation (26) yields the following equation where, It can be shown that the above numerical method (27) has the accuracy of (s + 1)-order.The numerical method (27) has a similar form to the traditional Adams methods, but we will later show that this method has the better numerical stability then traditional Adams methods.We call this numerical method the generalized Adams methods (GAMs).
When selecting the uniform grids, the coefficient values (β z ) of GAMs are shown in the following Tables 4 and 5 respectively.
For the uniform grids, we can define that therefore, the method (27) could be rewritten as According to Tables 4 and 5, we can verify that the proposed method (30) is consistent with the GAMs proposed by Brugnano et al. [19].In reference [21], it has been shown that, for even s, the methods (30) are both A-stable, and for odd s the methods (30) are also A-stable while the stability region is just the complete left half plane.
The same as Equation ( 12), Equation (32) is also the basic calculation scheme of DQM, however, the scheme of Equation ( 12) is more widely used in engineering application.In general, the matrix G in Equation (32) is not a full rank matrix, so its inverse matrix does not exist.
If the distribution of grid points is centro-symmetrical in the normalized interval [0, 1], the matrix G in Equation (32) will be a centro-antisymmetrical matrix.Since the inverse matrix of G does not exist, it is natural to consider its generalized inverse matrix.As is well known, the generalized inverse matrix-namely, the pseudo-inverse matrix G T of centro-antisymmetrical matrix with even order-will also be a centro-antisymmetrical matrix. Define then the Equation (32) can be written as Appl.Sci.2017, 7, 284 8 of 13 Subtracting the υ-th equation from the (υ + 1)-th equation in the matrix Equation (34), we can obtain the following equation It is obvious that the aforementioned Equation ( 35) is the GAMs with (s + 1)-order, in which Considering that s = 2υ − 1 and A defined in Equation ( 33) is a centro-antisymmetrical matrix with even dimension, it can be easily proved that ρ z in Equation ( 36) is centro-symmetrical.Therefore, Equation ( 35) can be further written as When selecting the uniform grids, the values ρ z are listed in the following Table 6.Obviously, when s = 1, the aforementioned Equation (37) is just the classic implicit trapezoidal rule, so, we call this numerical method the extended trapezoidal rules of the first kind (ETR 1 ).
Although ETR 1 are simply the GAMs, the deduction method of these two methods is slightly different.Therefore, comparing Table 6 with Table 4, it can be clearly seen that the concrete results-i.e., the coefficients-are different.
It has been shown in reference [21] that the stability regions of ETR 1 methods are the same as that of the classic implicit trapezoidal integration rule, and they are all A-stable methods.

Extended Trapezoidal Rules of the Second Kind
Let s be an odd number, and the expression of υ is also the same with Equation (31).According to matrix Equation ( 12), it can be drawn that . . .
In the above matrix equation, add (υ − 1)-th equation to the υ-th equation, we can get where, Considering that s = 2υ − 1, as well as the weighting coefficient g pq is centro-antisymmetrical and α z is centro-antisymmetrical, the Equation (39) could lastly be written as where, When selecting the uniform grids, the coefficient values (ω z ) in Equation ( 41) are shown in the following Table 7.
Table 7. Coefficients of extended trapezoidal rules of the second kind (ETR 2 ) using uniform grids.Obviously, when s = 1, the aforementioned method is just the implicit trapezoidal method; thus, this numerical method can be named the extended trapezoidal rules of the second kind (ETR 2 ).
According to Table 7, it can be verified that the method (41) is consistent with ETR 2 methods that proposed by Brugnano et al. [19].The stability regions of ETR 2 methods are the same with that of classic implicit trapezoidal integration rule, and they are all A-stable methods.

Fourier Stability Analysis for the Derived Compact Difference Schemes
Stability analysis is an important part of the numerical integration approach.In general, we use linear stability analysis by the use of a standard scalar test equation ( .x = λx) [21][22][23].However, the linear stability analysis does not provide information about the performance of schemes at different length scales.Unlike this, the Fourier stability analysis provides an effective way to quantify the resolution characteristics such as dispersion or phase error of the differencing approximations.In this section, we analyze the dispersion and dissipation characteristics of the derived difference schemes using Fourier stability analysis [1,24].
By the Fourier stability analysis, different numerical schemes have different estimates of the modified scaled wavenumber, and it is in general a complex quantity [25].It can easily be shown that the real part of the modified scaled wavenumber tends to amplify/attenuate the corresponding Fourier mode while the imaginary part introduces a phase error.Thus the deviations of the real and imaginary parts of the modified scaled wavenumber from the exact scaled wavenumber can be respectively taken to be the indicators of numerical dissipation and dispersion.Furthermore, for numerical stability of any scheme, one must look at the imaginary part of the modified scaled wavenumber.The imaginary part of the modified scaled wavenumber represents numerical dissipation only when it is negative.Any scheme that produces a positive imaginary part of the modified scaled wavenumber is numerically unstable because a positive imaginary part is equivalent to adding anti-diffusion [25].
Figure 3 shows the modified scaled wavenumber of GBDF based on uniform grids.
respectively taken to be the indicators of numerical dissipation and dispersion.Furthermore, for numerical stability of any scheme, one must look at the imaginary part of the modified scaled wavenumber.The imaginary part of the modified scaled wavenumber represents numerical dissipation only when it is negative.Any scheme that produces a positive imaginary part of the modified scaled wavenumber is numerically unstable because a positive imaginary part is equivalent to adding anti-diffusion [25].
Figure 3 shows the modified scaled wavenumber of GBDF based on uniform grids.Based on the above results, we can see that GBDF are all stable numerical methods, and this kind of method has a strong numerical dissipation and dispersion.
Figure 4 shows the modified scaled wavenumber of central compact difference schemes (CCDS) based on uniform grids.Based on the above results, we can see that GBDF are all stable numerical methods, and this kind of method has a strong numerical dissipation and dispersion.
Figure 4 shows the modified scaled wavenumber of central compact difference schemes (CCDS) based on uniform grids.Based on the above results, we can see that CCDS have a strong numerical dispersion but have almost no numerical dissipation.
Figure 5a,b show the modified scaled wavenumber of GAMs based on uniform grids and respectively for even s and odd s.Based on the above results, we can see that CCDS have a strong numerical dispersion but have almost no numerical dissipation.
Figure 5a,b show the modified scaled wavenumber of GAMs based on uniform grids and respectively for even s and odd s.Based on the above results, we can see that GBDF are all stable numerical methods, and this kind of method has a strong numerical dissipation and dispersion.
Figure 4 shows the modified scaled wavenumber of central compact difference schemes (CCDS) based on uniform grids.Based on the above results, we can see that CCDS have a strong numerical dispersion but have almost no numerical dissipation.
Figure 5a,b show the modified scaled wavenumber of GAMs based on uniform grids and respectively for even s and odd s.Based on the above results, we can see that the ETR1 are roughly equivalent to the GAMs with odd s in terms of numerical dissipation and dispersion.In fact, as mentioned earlier, these two methods are essentially the same.
Figure 7 shows the modified scaled wavenumber of ETR2 based on uniform grids.Based on the above results, we can see that the ETR2 are roughly equivalent to the ETR1 and GAMs with odd s in terms of numerical dissipation and dispersion.
To sum up, from the perspective of numerical dissipation, the CCDS, GAMs with odd s, ETR1 and ETR2 methods perform better than the GAMs with even s and GBDF methods, but from the Based on the above results, we can see that the ETR 1 are roughly equivalent to the GAMs with odd s in terms of numerical dissipation and dispersion.In fact, as mentioned earlier, these two methods are essentially the same.
Figure 7 shows the modified scaled wavenumber of ETR 2 based on uniform grids.Based on the above results, we can see that the ETR1 are roughly equivalent to the GAMs with odd s in terms of numerical dissipation and dispersion.In fact, as mentioned earlier, these two methods are essentially the same.
Figure 7 shows the modified scaled wavenumber of ETR2 based on uniform grids.Based on the above results, we can see that the ETR2 are roughly equivalent to the ETR1 and GAMs with odd s in terms of numerical dissipation and dispersion.
To sum up, from the perspective of numerical dissipation, the CCDS, GAMs with odd s, ETR1 and ETR2 methods perform better than the GAMs with even s and GBDF methods, but from the perspective of phase dispersion, the high-order GBDF methods behave slightly better than other methods.In the wavenumber range 1 h ω ≤ , the proposed methods both have a high Based on the above results, we can see that the ETR 2 are roughly equivalent to the ETR 1 and GAMs with odd s in terms of numerical dissipation and dispersion.
To sum up, from the perspective of numerical dissipation, the CCDS, GAMs with odd s, ETR1 and ETR2 methods perform better than the GAMs with even s and GBDF methods, but from the perspective of phase dispersion, the high-order GBDF methods behave slightly better than other methods.In the wavenumber range ωh ≤ 1, the proposed methods both have a high numerical resolution.
From the practical application of the point of view, the above different schemes have different practical application scenarios.For example, the GBDF method may be a better choice for electromagnetic transient simulation of electric power systems, because the dissipation of this method could avoid or effectively dampen fictitious numerical oscillations caused by abrupt phenomena.However, for electromechanical transient simulation of electric power systems, the GBDF method is likely to yield the wrong result due to the strong numerical dissipation: a stable system of its own, the corresponding simulation results may be stable.Based on the above analysis results combined with our experience, the high-order scheme of GAMs with even s and CCDS could be a better choice for electromechanical transient simulation of electric power systems.

Conclusions
Using the classic differential quadrature formulae while adopting the uniform grid points, this paper has symmetrically constructed a variety of compact difference schemes.Although some of the constructed methods in this paper are consistent with the boundary value method proposed by Brugnano et al., they are completely different in the means of construction.In fact, based on the methods described in this paper, it is easy to construct the compact difference schemes based on second-order derivatives.
Using Fourier stability analysis method, the characteristics such as the dissipation, dispersion and resolution of the different schemes were studied and compared.Based on these results, it is found that the CCDS, GAMs with odd s, ETR 1 and ETR 2 have almost no numerical dissipation, and GAMs, ETR 1 and ETR 2 have a high resolution when the wave number is not very large.In summary, the derived compact difference schemes in this paper could be used instead of classic differential quadrature methods for engineering applications.

Figure 1 .
Figure 1.Schematic diagram of the grid points for generalized backward difference formulae (GBDF).

Figure 1 .
Figure 1.Schematic diagram of the grid points for generalized backward difference formulae (GBDF).

Figure 2 .
Figure 2. Schematic diagram of the grid points for generalized Adams methods (GAMs).

Figure 2 .
Figure 2. Schematic diagram of the grid points for generalized Adams methods (GAMs).

3. 4 .
The Extended Trapezoidal Rules3.4.1.Extended Trapezoidal Rules of the First KindAs is shown in Figure1, and let s be an odd number.Define υ = (s + 1

Figure 3 .
Figure 3. Plot of real and imaginary parts of the modified wavenumber for GBDF.

Figure 3 .
Figure 3. Plot of real and imaginary parts of the modified wavenumber for GBDF.

Figure 4 .
Figure 4. Plot of real and imaginary parts of the modified wavenumber for central compact difference schemes (CCDS).

Figure 4 .
Figure 4. Plot of real and imaginary parts of the modified wavenumber for central compact difference schemes (CCDS).

Figure 4 .
Figure 4. Plot of real and imaginary parts of the modified wavenumber for central compact difference schemes (CCDS).

Figure 5 .
Figure 5. (a) Plot of real and imaginary parts of the modified wavenumber for GAMs with even s; (b) plot of real and imaginary parts of the modified wavenumber for GAMs with odd s.

Figure 5 .
Figure 5. (a) Plot of real and imaginary parts of the modified wavenumber for GAMs with even s; (b) plot of real and imaginary parts of the modified wavenumber for GAMs with odd s.

Figure 7 .
Figure 7. Plot of real and imaginary parts of the modified wavenumber for ETR2.

Figure 6 .
Figure 6.Plot of real and imaginary parts of the modified wavenumber for ETR 1 .

Figure 6 Figure 6 .
Figure6shows the modified scaled wavenumber of ETR1 based on uniform grids.

Figure 7 .
Figure 7. Plot of real and imaginary parts of the modified wavenumber for ETR2.

Figure 7 .
Figure 7. Plot of real and imaginary parts of the modified wavenumber for ETR 2 .

Table 4 .
Coefficients of generalized Adams methods (GAMs) using uniform grids (for even s).

Table 5 .
Coefficients of GAMs using uniform grids (for odd s).

Table 6 .
Coefficients of extended trapezoidal rules of the first kind (ETR 1 ) using uniform grids.