{'fraction_non_alphanumeric': 0.043658220558594876, 'fraction_numerical': 0.027569824359343506, 'mean_word_length': 4.565905448717949, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An optimized dispersion-relation-preserving combined compact difference scheme to solve advection equations

In this study, we first present an improved version of the classical sixth-order combined compact difference (CCD6) scheme to enhance the convective stability of advection equations through an increased dispersion accuracy. This improved fifth-order dispersion-relation-preserving combined compact difference scheme (DRPCCD5) has been rigorously analyzed through the dispersion, phase speed anisotropy and stability analyses. We then couple the DRPCCD5 scheme with the previous fifth-order compact-reconstruction weighted essentially non-oscillatory (CRWENO5) scheme using a novel hybrid strategy based on the monotonicity-maintenance criteria. To verify the resulting"optimized"hybrid scheme (ODRPCCD5), several benchmark problems with available exact solution are investigated. The comparison to the previous fifth-order WENO (WENO5) scheme shows that the ODRPCCD5 avoids numerical oscillation around discontinuities, handles large gradients well, and is much faster at the same accuracy because a coarser mesh can be used.


Introduction
Numerical simulations of advection equations are commonly found in many applications of practical importance, such as shock waves, shallow water flow, magnetohydrodynamics, and two-phase flow models. When numerically solving such convection-dominated partial differential equations (PDEs), it is desirable to minimize the indispensable dispersion error, which is defined as the discrepancy between the numerical and actual wavenumbers, because this enhances convective stability and allows for accurate capturing of small length scales in the wave phase [1].
For this purpose, dispersion-relation-preserving (DRP) approaches have been developed to enhance convective stability by rigorously preserving the dispersion relation [1][2][3][4][5][6]. Furthermore, compact difference schemes offer spectral accuracy with fewer grid points to improve convective stability [7][8][9][10]. These compact difference schemes have been extended to combined compact difference schemes (CCD) [11], in which first and second derivative terms are simultaneously evaluated in an implicit manner, making the scheme more compact and accurate. CCD schemes suffer from stability issues of boundary conditions when solving the PDE. In fact, these schemes need special treatment at the boundary nodes, in particular when simulating thin boundary layer problems. Hence, the boundary closures have been improved [12] to obtain better numerical properties, and the corresponding dissipation and de-aliasing properties have been discussed [13].
High spectral resolution schemes, such as the compact difference and CCD schemes, inevitably produce numerical oscillations near discontinuities and lead to failure of the flow simulation. In order to avoid numerical oscillations, high resolution schemes often use flux/slope limiters to bound the solution gradient around shocks or discontinuities [14,15]. Some representative schemes belonging to this class of methods include the essentially non-oscillatory (ENO) scheme [16,17] and their weighted variants, known as the weighted ENO (WENO) [18,19]. It's well known that the ENO and WENO schemes may be too dissipative for compressible turbulence simulations and aero-acoustics problems. Hence, the compact-reconstruction weighted essentially non-oscillatory (CRWENO) scheme [20] has been presented, in which compact sub-stencils are identified at each interface and combined using the WENO weights. WENO schemes have been intensively used for problems containing both shocks and complicated smooth solution structures [21,22].
Algorithms with high accuracy are required to capture small wavelengths and nonoscillatory behaviors across discontinuities like shock waves. For this purpose, special finite difference schemes have been introduced [23]. Also, the hybrid finite difference scheme based on the minimized dispersion and controllable dissipation (MDCD) technique has been developed to solve advection equations. This MDCD technique has been coupled with an optimized WENO scheme to make discontinuity capturing possible [24]. Many researchers have also proposed various alternative ways to improve the numerical schemes [25][26][27][28]. However, accuracy still remains a challenge because, to our knowledge, most if not all existing numerical schemes suffer from the drawback that they switch to a noncompact scheme at and near discontinuities, resulting in a loss of resolution.
In this study, a fifth-order dispersion-relation-preserving combined compact difference (DRPCCD5) scheme which has better DRP properties than previously reported compact difference schemes over a considerable range of wavenumbers is proposed. This scheme ensures that resolved energy components propagate closer to the correct physical speed, and that complex phenomena, involving interactions among different wavelength scales, can be captured. Furthermore, the DRPCCD5 scheme is coupled with the CRWENO5 scheme using a novel hybrid strategy based on the monotonicity-maintenance criteria. The numerical properties of the resulting "optimized" hybrid scheme (ODRPCCD5) are then 2 rigorously analyzed using several benchmark problems. This paper is organized as follows. Section 2 describes discretization of a standard advection equation and the time marching method, which is used in the present study. The schemes construction is carried out in Section 3. Section 4 includes the fundamental analysis of dispersion, dissipation, phase speed anisotropy, numerical group velocity, and numerical phase velocity for the proposed DRPCCD5 scheme. Several benchmark tests are performed in Section 5 to validate the ODRPCCD5 scheme. Section 6 draws concluding remarks based on the results presented in Section 5.

Time marching method
The one-dimensional linear wave equation can be expressed as ∂u ∂t where t is time, x the spatial coordinate, u the field variable, and f = cu with c the constant propagation speed of the wave. A conservative finite difference discretization of Eq. (1) results in an ordinary differential equation, which can be expressed as where h is the grid spacing andf i+ 1 2 is the numerical approximation of flux between points x i and x i+1 . In the present study, we apply the fourth-order Runge-Kutta (RK4) scheme and the sixth-order symplectic Runge-Kutta (SRK6) scheme [29] for time evolution. The explicit RK4 scheme reads (1) ) + 2F(u (2) )), For the SRK6 scheme, given the solution u n at t = n∆t, the solution u n+1 is obtained from the following iteration. We start with computing u ( j) and F ( j) = F(u ( j) ), where j=1 to 3, by numerically solving the following equations iteratively: where c = 1 2 3 5 . These updated values correspond to the times t = n + ( 1 2 + c)∆t, t = n + 1 2 ∆t, and t = n + ( 1 2 − c)∆t, respectively. Upon reaching the user's specified tolerance (10 −8 ), the solution at t = (n + 1)∆t is obtained as 3 The RK4 scheme is mainly used to run the numerical tests in this study because the implicit SRK6 scheme provides nearly the same results, but is very time-consuming (see results of linear advection problem #1 in Table 1 and Fig. 7).

Fifth-order non-compact difference scheme
The numerical flux can be reconstructed using a left or right biased interpolation [20]. The appropriate interpolation is chosen based on the sign of the wave speed, which in the case of a scalar PDE is given bŷ where the superscripts + and − denote left and right biased interpolations respectively. Note that the approximation of the left biased numericalû i+1/2 is described in this section. The first derivative term ∂u ∂x can be approximated to the desired order (r), reading ∂u ∂x whereû i+1/2 for odd r is computed using the linear reconstruction on a stencil Here b k is the coefficient and u i+k = u x i +kh . For r = 5, it readŝ This scheme has fifth-order spatial accuracy according to the derived modified equation given below The drawback of this scheme is that the magnitude of the leading error of the resulting scheme is too large. In addition, this scheme must adopt a very fine mesh to correctly capture the important advection flow structures. A space-time accurate numerical simulation of advection problems requires higher spatial resolution and dispersion-relationpreserving (DRP) properties. Such schemes act as an important numerical tool to solve complex physical problems displaying a large bandwidth of spatio-temporal scales. For this reason, we develop a new fifth-order dispersion-relation-preserving combined compact difference (DRPCCD5) scheme in the Section 3.2.

Fifth-order dispersion-relation-preserving combined compact difference (DRPCCD5) scheme
In this section, we present an improved upwind combined compact difference scheme. The first and the second derivative terms ( ∂u ∂x and ∂ 2 u ∂x 2 ) in a four-point grid stencil are ap-proximated as The coefficients shown in Eq. (14) are derived through Taylor series expansion. Elimination of the leading truncation error terms in the modified equation analysis enables us to get the formal accuracy order of six [11].
Derivation of the coefficients in Eq. (13) is started from performing Taylor series expansion on the terms with respect to u i to get the modified equation. The six leading truncation error terms derived in the modified equation analysis are then eliminated to get a set of six algebraic equations Derivation of two further algebraic equations are needed to determine all eight coefficients in Eq. (13). One way of deriving the two equations so as to get a better approximation of ∂u ∂x is to reduce numerical error of the accumulative type. We can then expect to retain the theoretical dispersive property of ∂u ∂x [2]. Our strategy of achieving the goal of reducing numerical dispersion error is to match the exact and numerical wavenumbers. Use of this underlying approach amounts to equating the effective wavenumbers α ′ and α ′′ to those shown on the right-hand sides of Eqs. (21) and (22) [2]. Following this line of derivation, we are led to get the two equations for α ′ h and α ′′ h as follows Equations (21) and (22) are solved to get the expression for α ′ h which has been used subsequently to minimize the dispersion error. The real and imaginary parts of α ′ h provide information regarding the dispersion error (phase error) and dissipation error (amplitude error), respectively. To improve the dispersive accuracy for α ′ , the exact value αh should be very close to To achieve the goal of improving solution accuracy, the positive-value error function E(α) defined below should be very small over the following integration interval for the modified wavenumber αh In Eq. (23) the weighting function W is chosen to be the denominator of This choice facilitates us to integrate E(α) exactly. To make the error function defined in 0 ≤ αh ≤ 7π 8 to be positive and minimal, two extreme conditions given by are enforced. These two constraint equations enforced for maximizing the dispersion accuracy are used together with the other six algebraic equations derived from the modified equation analysis to get not only a smaller dissipation error but also an improved dispersion accuracy. Note that several integration ranges have been numerically determined so as to find the best one that renders the smallest value of E. The resulting eight introduced unknown coefficients can be determined as a 1 = 0.8873686, 9692791 and c 3 = 1.9528828 from the above reduction of dispersion and dissipation errors. The upwinding scheme developed theoretically in four stencil points i− 2, i− 1, i and i + 1 for ∂u ∂x has the spatial accuracy of order fifth according to the derived modified equation given below It is noted that, unlike our strategy, Zhou et al. [28] chose the coefficient c 3 as free parameter so that the other seven coefficients are expressed as the linear functions of c 3 by Taylor's expansion. Then these eight coefficients were numerically determined by the standard sequential quadratic programming (SQP) method [30] to minimize the error function shown in Eq. (23). However, this optimization result is highly sensitive to the initial guess of c 3 , as pointed by Zhou et al. [28]. Define first the values of u at the half nodal points i ± 1 2 as follows: and One can then substitute them into Eq. (9) to get the algebraic equation for ∂u ∂x at the node i. Derivation of α i , β i and γ i is then followed by comparing the coefficients derived in Eq. (13) for ∂u ∂x | i . After a term-by-term comparison of Eq. (9), we are led to get the coefficients as follows: α 1 = 0.8873686, α 2 = 0.0491178, β 1 = 0.1495320, β 2 = −0.2507682, β 3 = −0.0123598, γ 1 = −0.0163964, γ 2 = 1.9528828. In brief,û i+1/2 of DRPCCD5 scheme for c i+1/2 ≥ 0 is given bŷ Thus, the magnitude of the leading error term in the compact interpolation is less than the corresponding non-compact interpolation on the same order (see Eqs. (12) and (26)).û i+1/2 of DRPCCD5 scheme for c i+1/2 < 0 can be similarly derived:

Weighted essentially non-oscillatory (WENO) scheme
Advection equations admit discontinuous solutions. Weighted essentially non-oscillatory schemes are designed to achieve the high order of accuracy at smooth regions and switch to lower order interpolation to avoid oscillations near discontinuities.

Fifth-order WENO (WENO5) scheme
The form of the interface flux reconstructed by the WENO5 scheme [19] readŝ In the above equation, we write The optimal weights are c 1 = 1 10 , c 2 = 6 10 and c 3 = 3 10 . A very small number (ε = 10 −6 ) is used to prevent division by zero. The smoothness indicators β k are given to detect large discontinuities and automatically switch to the stencil that generates the least oscillatory reconstruction by The WENO5 scheme gives fifth-order accurate results in smooth regions of the solution and is non-oscillatory near discontinuities.

Fifth-order compact-reconstruction WENO (CRWENO5) scheme
The drawback of higher order WENO schemes is the increasingly wide stencil when increasing the order of accuracy. Therefore, the CRWENO5 has been constructed using three third-order compact interpolations as candidates [20]. The CRWENO5 scheme can be expressed as Note thatf i+1/2 in Eq. (34) is the approximation of the left biased numerical fluxf CRW ENO+ Since the weights ω k in Section 3.3.1 are overly dissipative, they are determined, as suggested in the literature [31,32], using α k as Here, the optimal weights are c 1 = 1 5 , c 2 = 1 2 and c 3 = 3 10 . τ is simply defined as the absolute difference between β 0 and β 2 .

Fifth-order optimized dispersion-relation-preserving combined compact difference scheme (ODRPCCD5)
In this section, we briefly present the hybrid strategy to couple CCD with WENO schemes proposed by [24,27,28] and our novel hybrid strategy based on the monotonicitymaintenance criteria. Both strategies are compared with each other in Section 5.1.2.
Step 2: and perform Step 3 to Step 7 according to Fig. 3(b).
Coupling DRPCCD5 with CRWENO5 using this hybrid strategy leads to the ODR-PCCD5 scheme.

Dispersion and dissipation errors
The solution for the model equation is given by where i ≡ √ −1 andû α is the Fourier mode of the wave number α. Differentiation of the above equation leads to The approximated derivative term ∂u ∂x can be similarly written as Here, K r and K i , denoting the real and imaginary parts of α ′ h (cf. Eq. (21)), account for the dispersion and dissipation errors, respectively. Fig. 4 shows the dispersion and dissipation characteristics of fifth-order non-compact finite difference scheme (FD5) [20], fifth-order compact difference (CD5) scheme [20], eighth-order optimized compact difference (OCD8) scheme [33] and our proposed DR-PCCD5 scheme. It can be seen that the DRPCCD5 scheme has a better spectral resolution than the OCD8 scheme. The dispersion property of the DRPCCD5 scheme is better than those of the other schemes because of the improved dispersive accuracy. Furthermore, at frequencies with low dispersion error, the DRPCCD5 scheme has less dissipation than the other schemes.

Assessment of the phase speed anisotropy
In anisotropic two-dimensional problems, first-order differencing schemes tend to produce phase space errors [7,11]. To evaluate the phase space error of our DRPCCD5 scheme, we take the following two-dimensional advection equation into consideration u t + c x u x + c y u y = 0.
Here, c x = c cos(θ) and c y = c sin(θ) denote the velocity components along the x and y directions, respectively. For a wave propagating at the angle θ (≡ tan −1 ( c y c x )) with respect to the x-axis, the numerical phase speed anisotropy can be derived as follows [11] ℜ( c * c ) =
One can find from Fig. 5 that our proposed scheme reproduces phase speed anisotropies much better than the sixth-order combined compact difference (CCD6) scheme [11] at all scaled wavenumbers.

Amplification factor, numerical group velocity and numerical phase velocity
The properties, such as amplification factor, numerical group velocity and numerical phase velocity [9], of the present DRPCCD5 scheme are analyzed by solving the onedimensional wave equation, where the fourth-order accuracy Runge-Kutta (RK4) scheme is applied in time evolution. The present scheme is compared with the previous sixth-order combined compact difference (CCD6) scheme [11]. The general numerical solution of Eq. (47) is identified as such that the initial solution is given by Note that the u(x m ,t n ) can be obtained by substituting the above initial condition as [9] u(x m ,t n ) = A 0 (α)(G 2 In Eq.(55), the numerical amplification factor G(α) is defined as G(α) = G r +iG i = U(α,t n+1 ) U(α,t n ) . The term β is obtained as tanβ = − G i G r . The numerical group speed and numerical phase velocity are obtained as where Cr = c∆t h = ω∆t αh denotes the Courant number.

One-dimensional problems
The ODRPCCD5 scheme is tested to solve three linear advection equations and one inviscid Burgers' equation. The L 2 -errors and their corresponding spatial rates of convergence are tested for the linear advection problem#1. The computational costs are compared using different spatial discretization schemes for the linear advection problems#2. Two hybrid strategies descried in section 3.4 are used to solve linear advection problem#2. Finally, we extend ODRPCCD5 scheme to solve the one-dimensional Euler equations of the polytropic gas dynamics.

Linear advection problem #1
The problem with the smooth initial condition u(x, 0) = sin(2πx) for Eq. (1) with c = 1 is solved. Periodic boundary conditions are applied at two boundaries of the region 0 ≤ x ≤ 1. To compare the computational efficiency of time evolution, we solve this problem by using the sixth-order implicit symplectic Runge-Kutta scheme (SRK6) and the fourth-order explicit Runge-Kutta scheme (RK4). The twin-tridiagonal coefficient matrix for the DR-PCCD5 scheme is solved by the computationally effective solver, including twin-forward elimination and twin-backward substitution techniques, which is described in [11]. All the computational times are obtained using a Core i7, 3.40 GHz computer with 64.0 GB of RAM. Table 1 shows that the SRK6 scheme costs more CPU time than the RK4 scheme when the same spatial scheme and grid are used. Fig. 7 shows that the computational errors mainly come from the spatial discretization, by comparing RK4/WENO5 and RK4/DRPCCD5. Therefore, we employ the RK4 scheme for time evolution in the following numerical cases. The L 2 -errors and their corresponding spatial rates of convergence, by using DRPCCD5, WENO5, and CRWENO5 schemes, are given in Table 2 with time step ∆t = 1 × 10 −5 . It can be seen that all schemes can approximately achieve their theoretical order of accuracy.

Linear advection problem #2
The one-dimensional linear equation u t + u x = 0 is solved considering the following initial condition [34]: The constants are taken as a = 0.5, z = −0.7, δ = 0.005, α = 10, and β=(log 2)/36δ 2 . This initial condition consists of a discontinuous square wave, an exponential wave, a triangular wave, and a parabolic wave. Periodic boundary conditions are imposed here. The time step is chosen as ∆t = 0.05h. Fig. 8 shows the exact waveform and the waveform obtained by the WENO5 and ODRPCCD5 scheme on a grid with 200 points at t = 2 and t = 4. Figs. 9 and 10 show the magnified solution for the exponential and square waves at t = 4. In Fig. 9, one can see that the ODRPCCD5 show less clipping at the extreme than the WENO5 in the case of the exponential wave. In Fig. 10, the ODRPCCD5 scheme is less dissipative than the WENO5 scheme across the discontinuities. Since the DRPCCD5 scheme is not classified to be a non-oscillatory scheme, the predicted kinks near the root of square wave is computationally inevitable. Comparing the magnitude of errors produced by WENO5, DRPCCD5 and ODRPCCD5 for this test problem shows that ODRPCCD5 performs better. The computational costs using WENO5, DRPCCD5 and ODRPCCD5 schemes are compared based on different grids, as shown in Table 3. The ODRPCCD5 scheme needs more CPU time than the other two schemes if the same grid is used because this scheme is hybrid. However, the spectral properties of the ODRPCCD5 scheme imply that it may apply a coarser grid to achieve the same resolution as the WENO5 scheme at the same order of convergence. As shown in Fig. 11, the ODRPCCD5 scheme with 600 grids reaches a better resolution than WENO5 scheme with 1600 grids. Meanwhile, it only needs 2.37s in comparison with 3.9s by WENO5 scheme.
The two hybrid strategies introduced in section 3.4 are used to solve the advection equations, and the numerical results are plotted in Fig. 12 with 400 grids and ∆t = 0.05h at t = 2.0. In Fig. 12, we can see that the solution is not damped when using the previous hybrid strategy by [24,27,28] when r c = 0.1. Therefore, this hybrid strategy needs an appropriate trial parameter (r c ) to damp the oscillation. In contrast, our hybrid strategy is based on the monotonicity-maintenance criteria, which automatically limits oscillations and captures discontinuities, as shown in Fig. 12(b).

Linear advection problem #3
We solve the linear equation u t + u x = 0, −1 ≤ x ≤ 1, with periodic boundary condition [35]. The initial condition reads The predicted results in the domain with 200 grid points are plotted in Fig. 13 at t = 20. It can be seen that ODRPCCD5 scheme performs better than the WENO5 scheme.

Non-linear advection problem
We solve the Burgers' equation u t + (0.5u 2 ) x = 0, −1 ≤ x ≤ 1, with periodic boundary condition. The initial condition is u(x, 0) = 2 + sin(π(x + 1)). The solution to Burgers' equation is smooth for t < 1 π and it develops shocks for t = 1 π . The results obtained at t = 0.3 (before shock) and t = 0.35 (after shock) are plotted in Fig. 14 in the domain with 200 grid points. The time step is chosen as ∆t = 0.1h in this computation. We observe that ODRPCCD5 gives better results than the DRPCCD5 scheme at t = 0.35.

The Shu-Osher problem
In this case, we solve the one-dimensional Euler equations of gas dynamics [24] ∂ ∂t 13 where ρ, q, p and E are the density, velocity, pressure and total energy of the conserved fluid, respectively. The initial conditions are (ρ, u, p) = (3.857143, 2.629369, 10.3333) , i f x ≤ 1 (1 + 0.2 sin(5x), 0, 1) , otherwise This test case leads to very strong shock waves and is employed to validate the shockcapturing capability of the proposed ODRPCCD5 scheme. Reflective boundary conditions are applied at both x = 0 and x = 10. Since the exact solution for this problem is not available, the solution computed in 10000 grids is considered as the exact solution. Fig. 15 shows waveforms at t = 0.45, t = 0.9, t = 1.35 and t = 1.8 (grid spacing h = 1 40 , time step ∆t = 0.05h). It can be seen that the shock-waves are well reproduced by our proposed ODRPCCD5 scheme.

Two-dimensional problems
In this subsection, we illustrate the capacity of the ODRPCCD5 scheme through twodimensional numerical simulations.

Vortex flow problem
The equation φ t + (uφ) x + (vφ) y = 0 is solved using an initial circle shape in a square of unit length, within which the vortex flow field (u, v) is given by [36] u = −sin 2 (πx)sin(2πy), The radius of the circle is 0.15 located at the center (0.5, 0.75). At t = T the flow field was reversed, so that the exact solution at t = 2T should coincide with the initial condition. This problem has been known to be computationally challenging since its solution is stretched and torn by the vortex flow where a very thin filament having a scale of single mesh size can be generated. Computations were performed for T = 2.5 and ∆t = 1 1000 . The predicted results of WENO5 and ODRPCCD5 are compared for the calculation of φ = 0. The results obtained in 100 × 100 grids at t = 1.5, 2.5, 4, 5 are plotted in Fig. 16. It is clear that the solution computed using the ODRPCCD5 scheme is maintained within a thin and elongated filament on the scale of one grid spacing. On the contrary, the WENO5 scheme results in a considerable reduction of the area at the head and tail of the filament. In Fig. 16(d), one can see that the solution computed using our proposed scheme returns to its initial state. In Fig. 17, the ODRPCCD5 scheme using a 100 × 100 mesh can reach the same resolution at t = 2.5 as the WENO5 scheme using a 200 × 200 mesh. Hence, the ODRPCCD5 needs less CPU time (8.30s) than the WENO5 scheme (17.32s).

Zalesak's problem
The Zalesak's problem [37,38] is one of the best known benchmark cases for testing the developed advection scheme. The slotted disk has a radius of 15 and a slot width of 5. It is initially located at (50,75) in the domain of size (100,100). The prescribed velocity field is given as 14 The results predicted for 100 × 100 grid points at t = 50π, t = 100π, t = 150π and t = 200π are plotted in Fig. 18(a). The results are also plotted in Fig. 18(b) in the domain with 200 × 200 grid points. The solution computed with the proposed scheme is in good agreement with the exact (or initial) solution as shown in Fig. 18(b).

Concluding remarks
In this paper, a fifth-order dispersion-relation-preserving combined compact difference (DRPCCD5) scheme has been proposed, which shows increased dispersion accuracy and improved dispersion-relation-preserving properties compared to the CCD6 [11] scheme. To make discontinuity capturing possible and handle large gradients, an optimized DRPCCD5 scheme (ODRPCCD5), which couples the DRPCCD5 and CRWENO5 schemes, is constructed using a novel hybrid strategy based on the monotonicity-maintenance criteria. The numerical solutions of linear problems show that our ODPRCCD5 scheme performs very well and is much faster than the previous WENO5 scheme at the same accuracy. In addition, the ODPRCCD5 scheme produces non-oscillatory solutions of the Euler equations in domains with discontinuities, and it can handle sharp resolutions when solving the twodimensional vortex flow and Zalesak's problems. We plan to apply our algorithm to solve the three-dimensional Navier-Stokes equations for the simulation of two-phase flows in future studies.     (b) Figure 4: Comparison of K i (αh) and K r (αh) amongst the proposed DRPCCD5 scheme, WENO5 scheme [19], CD5 scheme [20], and OCD8 scheme [33]. (a) K i ; (b) K r .  Figure 5: Comparison of the predicted phase speed anisotropy, which is plotted against θ, for the proposed DRPCCD5 scheme and the CCD6 scheme of Chu and Fan [11]. (a) DRPCCD5 scheme; (b) CCD6 scheme [11].      (b) Figure 11: The predicted results for linear advection problem#2 are plotted at t = 4.