Next Article in Journal
New Integral Inequalities via Generalized Preinvex Functions
Next Article in Special Issue
Numerical Scheme Based on the Implicit Runge-Kutta Method and Spectral Method for Calculating Nonlinear Hyperbolic Evolution Equations
Previous Article in Journal
The Darboux Transformation and N-Soliton Solutions of Gerdjikov–Ivanov Equation on a Time–Space Scale
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Reconstructed Interpolating Differential Operator Method with Arbitrary Order of Accuracy for the Hyperbolic Equation

1
College of Meteorology and Oceanography, National University of Defense Technology, Changsha 410000, China
2
Department of Mathematics, National University of Defense Technology, Changsha 410000, China
*
Author to whom correspondence should be addressed.
Axioms 2021, 10(4), 295; https://doi.org/10.3390/axioms10040295
Submission received: 15 September 2021 / Revised: 27 October 2021 / Accepted: 2 November 2021 / Published: 6 November 2021

Abstract

:
We propose a family of multi-moment methods with arbitrary orders of accuracy for the hyperbolic equation via the reconstructed interpolating differential operator (RDO) approach. Reconstruction up to arbitrary order can be achieved on a single cell from properly allocated model variables including spatial derivatives of varying orders. Then we calculate the temporal derivatives of coefficients of the reconstructed polynomial and transform them into the temporal derivatives of the model variables. Unlike the conventional multi-moment methods which evolve different types of moments by deriving different equations, RDO can update all derivatives uniformly via a simple linear transform more efficiently. Based on difference in introducing interaction from adjacent cells, the central RDO and the upwind RDO are proposed. Both schemes enjoy high-order accuracy which is verified by Fourier analysis and numerical experiments.

1. Introduction

In the last few decades, high-order numerical methods are becoming increasingly popular in the research community due to the advantages over the low-order methods in achieving higher resolution more efficiently with lower computational cost [1]. The conventional high-order approaches based on high-order finite difference (FD) and the finite volume (FV) methods only use one degree of freedom (DOF) on each node or cell. For these methods, a high-order approximation requires a wide stencil which introduces huge communication costs between computational nodes on supercomputers.
To narrow the stencil, a natural and feasible approach is to include two or more DOFs per cell so as to construct a higher-order polynomial. The constrained interpolation profile (CIP) scheme [2,3] evolves two moments, i.e., the physical variable and the spatial derivative simultaneously and independently according to the governing equations in different forms. Specifically, CIP describes the spatial profile using the third-order Hermite interpolation and updates the profile according to the local analytic solution in via the semi-Lagrangian approach. Following this pioneering work, the so-called CIP-conservative semi-Lagrangian (CIP-CSL) schemes [4,5] add the cell-averaged value as an extra moment to ensure the conservativeness for the scalar conservative advection transport. Successive studies have resulted in a more general framework, the so-called CIP/multi-moment finite volume method (CIP/MM FVM) which has been applied to various fluid dynamic simulations [6,7,8]. The interpolated differential operator (IDO) method [9] is similar to CIP but IDO updates the model variables using the Eulerian method instead of the semi-Lagragian method. A conservative variant of IDO scheme also predicts the cell-averaged moment separately [10]. Both CIP and IDO require an additional equation for the spatial derivative moment which can be difficult to obtain when treating nonlinear and high-dimensional problems.
A compact CIP/MM-FVM formulation that uses high-order derivative moments can reach arbitrary accuracy order has been presented in [11]. This approach predicts the high-order derivative moments in the Eulerian representation via solving a linearized high-order derivative Riemann derivative problem [12]. As a more flexible variant of CIP/MM-FVM, the multi-moment constrained finite volume (MCV) [13] treats the cell-averaged moments as constraints instead of explicit model variables. Following this principle, the MCV method has been generalized to the unstructured meshes [14,15] recently.
High-order schemes based on local flux reconstruction are another representative class of high-order schemes which also use compact stencil for spatial discretization and are increasingly attractive. Examples includes the discontinuous Galerkin(DG) [16], the spectral volume [17] and the flux reconstruction [18] method.
This article presents a new efficient compact multi-moment method termed reconstructed interpolating differential operator (RDO) method which does not involve the cross-cell reconstruction for the hyperbolic equations. This work can be viewed as an generalization of the original IDO scheme [9] to arbitrary orders, non-uniform meshes and higher dimensions. Moreover, this work proposes a novel updating rules for model variables via a direct linear transform, which does not require introducing additional equations for different types of moments. The core idea is simple and direct. In each cell, the solution is approximated by a polynomial to interpolate the moments including the point values and high-order derivatives at the cell boundaries. The next step computes the temporal derivatives of the polynomial coefficients with the aid of solution points. Then the temporal derivatives are retrieved from the temporal derivatives of polynomial coefficient via a linear transform. To deal with the non-uniqueness of temporal derivatives computed form different cells at cell boundaries, the central and the upwind RDO are proposed. Fourier analysis and numerical tests indicate that our scheme can achieve arbitrary orders of accuracy.
The remaining part of this paper is organized as follows. Section 2 gives the basic formulation of RDO in one dimension. Section 3 describes the generalization of RDO to higher dimensions. Then Section 4 investigates the accuracy order and stability by classical Fourier analysis. Numerical tests in Section 5 validate the Fourier analysis results and compare the computational efficiency and numerical performance of RDO of various orders. Section 6 ends this paper with a few conclusions and discussions.

2. Basic Formulation for the One-Dimensional Problem

Consider the 1D scalar hyperbolic equation
u t + f ( u ) x = 0 ,
where x is the spatial coordinate, t denotes time, u is the conserved variable and f ( u ) is the flux.
To begin with, we introduce some notations here. The bold lower-case letter such as a , u stands for a column vector. The bold capital letter like A denotes a matrix. The column vector a 1 , , a m T is denoted by a k k = 1 m for simplicity. A T denotes the transpose of A . The p norm of the vector a is denoted by a p and a is a simplified form of a 2 throughout this paper. The matrix norm of H is defined by H = max z = 1 H z / z .

2.1. Spatial Discretization

We divide the computational domain x min , x max into N non-overlapping cells or elements, among which the j-th cell is I j : = [ x j 1 / 2 , x j + 1 / 2 ] . For convenience, we denote the middle point of I j by x j = ( x j 1 / 2 + x j + 1 / 2 ) / 2 and the length of I j by Δ x j = x j + 1 / 2 x j 1 / 2 .
Model variables including point-based value and derivatives are independently evolved in our scheme. For simplicity, the RDO scheme that have M model variables in each cell is denoted by RDO-M. The allocation of model variables has a slight difference between M = 2 K and M = 2 K + 1 . For RDO- 2 K , 2 K model variables in each cell are all located at the cell boundaries. These model variables are { u j ± 1 / 2 , ( u x ) j ± 1 / 2 , , ( u x K 1 ) j ± 1 / 2 } , where ( u x k ) j 1 / 2 is the approximation of u x k ( x j 1 / 2 ) at x = x j 1 / 2 .
For RDO- ( 2 K + 1 ) , we add the function value at the middle of the cell, i.e., u j u ( x j ) as an extra model variable and use a 2 K -thorder polynomial to represent the solution polynomial. Specifically, the ( 2 K + 1 ) model variables used in RDO- ( 2 K + 1 ) are { u j ± 1 / 2 , ( u x ) j ± 1 / 2 , , ( u x K 1 ) j ± 1 / 2 } u j .

2.2. Updating Model Variables

This section only presents the detailed formulation of RDO- 2 K ( K 2 ). The procedure for RDO- ( 2 K + 1 ) is almost the same and hence omitted here for simplicity. The procedure to compute time derivatives of model variables can be divided into three stages as follows.

2.2.1. First Stage: Reconstructing Polynomial from Model Variables

This stage finds a solution polynomial h j ( x ) on I j = [ x j 1 / 2 , x j + 1 / 2 ] to interpolate model variables defined at cell boundaries. The interpolation polynomial is of ( 2 K 1 ) -th order since there are totally 2 K model variables. We first introduce a local coordinate to map I j to the standard cell [ 1 , 1 ] ,
s ( x ) : = 2 ( x x j ) Δ x j [ 1 , 1 ] , for x [ x j 1 / 2 , x j + 1 / 2 ] .
The usage of s ( x ) unifies the interpolation templates and hence reduces the computational cost in polynomial reconstruction. The solution polynomial h j ( x ) under the local coordinate s is then h ˜ j ( s ) = h j ( x ( s ) ) . d h ˜ j d s and d h j d x are connected by the chain rule
d h ˜ j ( s ) d s = d x d s d h d x = Δ x j 2 d h d x .
Suppose h ˜ j ( s ) = i = 0 2 K 1 a i s i . Then the polynomial coefficients a i i = 0 2 K 1 can be determined by the following linear system:
d k d s k h ˜ j ( s ) | s = ± 1 = i = k 2 K 1 i ! ( i k ) ! a i ( ± 1 ) i = Δ x j 2 k · ( u x k ) j ± 1 / 2 ,
for k = 0 , , K . This is a linear system about 2 K unknowns a 0 , , a 2 K 1 and we can write it into the matrix form
Ha = d j ,
where H is the coefficient matrix and d j denotes the rescaled model variables
d j = [ u j 1 / 2 , Δ x j 2 ( u x ) j 1 / 2 , , Δ x j 2 K 1 ( u x K 1 ) j 1 / 2 , u j + 1 / 2 , Δ x j 2 ( u x ) j + 1 / 2 , , Δ x j 2 K 1 ( u x K 1 ) j + 1 / 2 ] T .
The coefficients of the solution polynomial is then
a = H 1 d j .
Through the high-order Hermite interpolation (4), the solution polynomial is ( K 1 ) -th smooth globally. Introducing the local coordinate makes the inverse of the coefficient matrix H the same for all cells so as to reduce the computational costs in solving the linear system in each cell. When the condition number of this linear system is too high as K increases, one can use the symbolic tool in Maple or Matlab to obtain a sufficiently accurate inverse of the matrix to avoid the loss of accuracy.

2.2.2. Second Stage: Computing Temporal Derivatives on Solution Points

This part computes the temporal derivative of the obtained solution polynomial h j ( x ) . For this, we choose 2 K solution points x j , k = x j + ξ k Δ x j , k = 1 , , 2 K on I j with ξ k [ 1 , 1 ] . { ξ k } k = 1 2 K are the same for each element. Through numerical and theoretical analysis later, we find that a simple choice of uniformly distributed solution points is sufficient for good numerical performance, i.e.,
ξ k = 1 + 2 ( k 1 ) 2 K 1 , k = 1 , , 2 K .
We can now easily update the solution polynomials via updating solution values located on selected solution points. Current solution values and approximated spatial derivatives are, respectively,
u j , k = h ˜ j ( s ( x j , k ) ) , ( u x ) j , k = Δ x j 2 1 d h ˜ j ( s ) d s | s = s ( x j , k ) , k = 1 , , 2 K .
where d h ˜ j d s can be obtained by
d h ˜ j d s | s = s ( x j , k ) k = 1 2 K = D a , D = 0 ξ 1 ( 2 K 1 ) ξ 1 2 K 2 0 ξ 2 ( 2 K 1 ) ξ 2 2 K 2 0 ξ M ( 2 K 1 ) ξ M 2 K 2 .
According to (1), the temporal derivatives of solution values become
d u j , k d t = f u ( u j , k ) · ( u x ) j , k , k = 1 , , 2 K .
Then we use d u j , k d t to retrieval the temporal derivatives of the coefficients of solution polynomials h j ( s ) . This can be done by solving the following linear systems
d h ˜ j s ( x j , k ) d t = i = 0 2 K 1 d a i d t s ( x j , k ) i = d u j , k d t , k = 1 , , 2 K .
Note that s ( x j , k ) = ξ k . This is a 2 K × 2 K linear system about d a i d t and the coefficient matrix is the Vandermonde matrix
V = 1 ξ 1 ξ 1 2 K 1 1 ξ 2 ξ 2 2 K 1 1 ξ 2 K ξ 2 K 2 K 1 .
This coefficient matrix remains the same for all cells. Therefore, we also can store the inverse in advance to reduce computational costs.

2.2.3. Third Stage: Retrieving the Temporal Derivatives

Using d a i d t , the temporal derivatives of model variables can be directly retrieved by the following linear transform,
d d t ( u x k ) j ± 1 / 2 j = Δ x j 2 k i = k 2 K 1 i ! ( i k ) ! d a i d t ( ± 1 ) i , k = 0 , , K 1 .
The subscript j in the left-hand side of (14) indicates that this temporal derivative is computed from the cell I j alone. At the cell boundary x = x j + 1 / 2 , we obtain the temporal derivatives from the I j and I j 1 at the same time, which can take different values. This leaves us the problem about how to define the proper temporal derivatives of ( u x k ) j + 1 / 2 k = 0 K properly.
A natural thought is choosing the one from the upwind direction so as to obey the physical property of the hyperbolic equation. We denote d ϕ j + 1 / 2 d t computed in the cell I j by d ϕ j + 1 / 2 d t j for the model variable ϕ . Then the upwind temporal derivative of ϕ is
d d t ϕ j + 1 / 2 = d ϕ j + 1 / 2 d t j up ,
where j up denotes the index in the upwind direction which is determined by the sign of f u ( u j + 1 / 2 ) . To be more specific,
j up = j , if f u ( u j + 1 / 2 ) 0 ; j + 1 , otherwise .
This scheme is termed upwind RDO scheme due to the upwind property.
Another idea to determine d d t ϕ j 1 / 2 is simply computing the algebraic average of the temporal derivatives from two cells just like the classical central difference scheme. Then the central temporal derivative of some model variables ϕ is
d d t ϕ j 1 / 2 = 1 2 d ϕ j + 1 / 2 d t j + d ϕ j + 1 / 2 d t j + 1
One may consider a sum weighted by cell length in (17) when encountering the non-uniform mesh like the spectral element method [19]. However, (17) works sufficiently well, as shown in the numerical tests on non-uniform meshes.

2.3. Time Integration

We have obtained the spatial discretization and temporal derivatives of model variables, which can be described in the semi-discrete form
d Φ d t = L ( Φ ) ,
where Φ denotes the collection of all model variables. Suppose the value of Φ at n-th time step is Φ n , then Φ n + 1 at the next time step is computed by the fourth-order Runge–Kutta method [20] to ensure the numerical accuracy in time. We have
Φ n + 1 = Φ n + Δ t 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) ,
where Δ t is the time step length and
k 1 = L ( Φ n ) , k 2 = L ( Φ n + 1 2 k 1 Δ t ) , k 3 = L ( Φ n + 1 2 k 2 Δ t ) , k 4 = L ( Φ n + k 3 Δ t ) .

2.4. Discussion

The involved polynomial interpolation procedures (4) and (14) are implemented within the cell. In other words, no cross-cell polynomial interpolation is needed in the proposed schemes, which means the scheme is compact and makes it possible to solve the problems on the non-uniform mesh more flexibly than the conventional IDO methods.
We can also construct RDO-2 and RDO-3 by assigning merely function values at cell boundaries without derivative moments. Then central and the upwind RDO-2 are, respectively, the classical central difference and the upwind difference method. Besides, the spectral element method [21] with three element base function in one dimension is also the special case of central RDO-3.

3. Generalization to Higher Dimensions

Now we extend the algorithm to rectangular meshes to solve the two dimensional conservative equations and briefly explain how to extend this problem into three dimension. The core issue is to define proper model variables on cell boundaries and find a benign polynomial space to interpolate them. Furthermore, the reconstructed piecewise polynomials should also be at least differentiable at the cell boundaries. Therefore we also use Hermite interpolation to represent the solution polynomials.

3.1. Spatial Discretization

Suppose the computational domain is a rectangle [ a , b ] × [ c , d ] and is divided into N 1 × N 2 non-overlapping cells. We denote the cell at the i-th row and j-th column by I i , j = [ x i 1 / 2 , x i + 1 / 2 ] × [ y j 1 / 2 , y j + 1 / 2 ] . Denote Δ x i = x i + 1 / 2 x i 1 / 2 and Δ y j = y j + 1 / 2 y j 1 / 2 .
To explain the intuition of spatial discretization, we first consider the bicubic Hermite interpolation on the cell I i , j . In this case, 16 DOFs are needed as 4 DOFs are assigned on each vertex. These 4 DOFs are u , u x , u y and the second-order mixed derivative u x y . A bit different from the 1D situation, u x y are requited here to guarantee that the reconstructed piecewise Hermite polynomials possesses continuous normal derivatives along cell edges. A similar allocation of model variables can be seen in the IDO method for solving the two-dimensional Poission Equations [22].
For general RDO- 2 K in 2D, the needed model variables on each vertex are i + j u x i y j 0 i , j K 1 , which means that each element has totally 4 K 2 model variables. However, since the adjacent cells share the same model variables on common vertices, the number of effective DOFs on each element is only K 2 . We use the following polynomial spaces for representing the solution polynomial is the following full polynomial spaces [23],
Q 2 K = { q ( x , y ) | q ( x , y ) = 0 l , m 2 K 1 a l , m x l y m } .
Its space dimension is dim ( Q 2 K ) = ( 2 K ) 2 . Q 2 K is the tensor product of the 1D Hermite polynomials.
For RDO- ( 2 K + 1 ) , the reconstructed polynomial space becomes Q 2 K + 1 which has ( 2 K + 1 ) 2 dimensions. Hence, each element needs additionally ( 2 K + 1 ) 2 ( 2 K ) 2 = 4 K + 1 model variables, which include u i , j that denotes the function value on the barycenter, and K model variables, u , u n , u n 2 , , u n K 1 on four edge centers. Here u n k denotes the kth-order normal derivative on the edge center.
Figure 1a,b illustrate the choices of model variables for RDO-4 and RDO-5, which are representative examples of RDO- ( 2 K ) and RDO- ( 2 K + 1 ) , respectively. Dots represent the point-valued variables. Single arrows represent the first derivatives. Paired arrows denote the mixed derivatives u x y .

3.2. First Stage: Reconstructing Solution Polynomial

The first step interpolates the solution polynomials h i , j ( x , y ) on I i , j . To begin with, we also introduce a local coordinate to map the cell I i , j to the standard cell [ 1 , 1 ] × 1 , 1 ,
s ( x , y ) = 2 ( x x i ) / Δ x i [ 1 , 1 ] ; r ( x , y ) = 2 ( y y j ) / Δ y j [ 1 , 1 ] .
Similar to the one-dimensional case, h i , j ( x , y ) is also transformed into h ˜ i , j ( r , s ) under the local coordinate. Correspondingly, 4 K 2 model variables serves as constraints. The solution polynomial
h ˜ i , j ( r , s ) = α = 0 2 K 1 β = 0 2 K 1 a α , β r α s β
is then determined by
l + m r l s m h ˜ ( r , s ) | r = ± 1 , s = ± 1 = α = l 2 K 1 β = m 2 K 1 a α , β l ! ( l α ) ! m ! ( m β ) ! ( ± 1 ) α ( ± 1 ) β = Δ x i 2 l Δ y j 2 m ( u x l y m ) i ± 1 / 2 , j ± 1 / 2 , 0 l , m 2 K 1 .
The inverse of the coefficient matrix of the linear system (23) about a α , β can also be solved in advance so as to effectively save computational costs in solving (23).

3.3. Second Stage: Computing Temporal Derivatives on Solution Points

This stage represents the solution polynomial on solution points and update it. The chosen solution points in 2D are tensor products of 1D points. Specifically, these 4 K 2 points in the local coordinate are in the form of
( s k , r l ) = ( ξ k , ξ l ) , k , l 1 , , 2 K .
Figure 1c,d illustrate the distribution of solution points for RDO-4 and RDO-5. The solution values and spatial derivatives can be computed as
u i , j ; k , l = h ˜ i , j ( s , r ) ,
( u x ) i , j ; k , l = ( Δ x j 2 ) 1 s h ˜ j ( s , r ) ,
( u y ) i , j ; k , l = ( Δ y j 2 ) 1 r h ˜ j ( s , r ) .
Then the temporal derivatives of solution values become
d u i , j ; k , l d t = f u ( u i , j ; k , l ) ( u x ) i , j ; k , l g u ( u i , j ; k , l ) ( u y ) i , j ; k , l .
Then the temporal derivatives of coefficients of the solution polynomial, d a α , β d t can be obtained by solving the following linear system,
d h ˜ j ( s ( x i , j ; k , l ) ) d t = l = 0 2 K 1 m = 0 2 K 1 d a α , β d t l ! ( l α ) ! m ! ( m β ) ! ( ± 1 ) α ( ± 1 ) β .

3.4. Third Stage: Retrieving the Model Variables

From d a i , j d t , we can obtain the temporal derivatives of model variables by the following relationships for RDO- 2 K
d d t u x l y m i ± 1 / 2 , j ± 1 / 2 i , j = ( Δ x i 2 ) l ( Δ y i 2 ) m α = l 2 K 1 β = m 2 K 1 l ! ( l α ) ! m ! ( m β ) ! d a i , j d t ( 1 ) α + β
for l , m 0 , , 2 K 1 . The subscript i , j in RHS of (30) specifies that the temporal derivatives are computed from I i , j and do not involve any external information. To properly introduce the external interaction, both central and upwind schemes can be reconstructed similar to the one dimensional case.

3.5. Three Dimensional Case

In three dimensions, RDO- 2 K on the cuboid mesh can also be constructed similarly via the polynomial space
C 2 K = { p ( x , y , z ) | 0 i , j , k K 1 a i , j , k x i y j z k } ,
and chooses i + j + k u x i y j z k 0 i , j , z K 1 on each vertex as model variables. There are totally ( 2 K ) 3 model variables in each cell. The remaining procedure is exactly the same as in lower dimensional space and we do not repeat the formulation here.

4. Fourier Analysis

This part analyses the accuracy and the stability of presenting scheme in solving the linear advection equation. For concreteness, we only give a detailed analysis for RDO-4 since the generalized RDO-M can be analyzed following the same routine.

4.1. Formulation of the Amplification Matrix

Consider the linear scalar advection equation
u t + u x = 0 .
Suppose the spatial domain is 0 , L and is discretized by a uniform mesh spacing Δ x . The initial condition is periodic, u ini = e i w x / Δ x where i denotes 1 and the scaled wavenumber is w = 2 π k Δ x / L [ 0 , π ) . Then the initial model variables for the upwind RDO-4 on I j = x j 1 / 2 , x j + 1 / 2 are
u j ± 1 / 2 = e i w ( x j ± 1 / 2 Δ x ) / Δ x ,
( u x ) j ± 1 / 2 = i w Δ x e i w ( x j ± 1 / 2 Δ x ) / Δ x .
The core in Fourier analysis is reformulating the algorithm into the matrix form
d d t d = 1 Δ x S d
using the periodicity of u ini and S denotes the amplification matrix. We do this step by step as follows.
Stage 1 described in Section 2 reconstructs the Hermite polynomial coefficients a j from the rescaled model variables d j in the cell I j ,
a j = H 1 d j ,
where H is defined by (5) and d j is the composition of the model variables after the coordinate transform,
d j = [ u j 1 / 2 , ( u x ) j 1 / 2 Δ x / 2 , u j + 1 / 2 , ( u x ) j + 1 / 2 Δ x / 2 ] T = [ 1 , i w 2 , 1 , i w 2 ] T e i w x j / Δ x .
In Stage 2, the solution values and derivatives on solution points are, respectively,
u j , k k = 1 4 = VH 1 d j ,
and
( u x ) j , k k = 1 4 = 2 Δ x DH 1 d j ,
where V and D can be found in (13) and (10). According to (31), the temporal derivatives on the solution points are
d d t u j , k k = 1 4 = 2 Δ x DH 1 d j .
Then the temporal derivatives on solution points are transformed back to the temporal derivatives of model variables, which is the inverse problem of (37):
d d t d j j = H V 1 d d t u j , k k = 1 4 = 2 Δ x H V 1 DH 1 d j = 1 Δ x S ˜ d j ,
where the subscript j in the left-hand side indicates that the temporal derivatives are computed locally form j-th cell, and S ˜ = 2 H V 1 D H 1 .
In the linear case, V 1 D is an invariant about the choice of solution points. Furthermore, the matrix V 1 D has a simple structure
V 1 D = 0 1 0 0 0 0 2 0 0 0 0 3 0 0 0 0 .
Therefore, it can be easily concluded that RDO schemes using different sets of solution points are identical in solving the linear advection equation as the following theorem tells.
Theorem 1.
For RDO-M schemes ( M 4 ), the numerical result is independent of the solution points chosen in solving the linear advection equation on uniform meshes.
Proof. 
This can be proved by showing that V 1 D is invariant about the choice of solution points ξ 1 , , ξ M . In fact, one can easily check that
V 1 D = 1 ξ 1 ξ 1 M 1 1 ξ 2 ξ 2 M 1 1 ξ M ξ M M 1 1 0 ξ 1 ( M 1 ) ξ 1 M 2 0 ξ 2 ( M 1 ) ξ 2 M 2 0 ξ M ( M 1 ) ξ M M 2 = 0 1 0 2 0 3 M 1 0 .
Moreover, it is noticed that the in-cell temporal derivatives (39) only relies on H , V 1 D and d j . H and d j are independent with the choice of solution points. Therefore, we can conclude that the semi-discrete scheme is independent of solution points chosen. □
Equation (39) gives the in-cell numerical approximation of derivatives of model derivatives. Denote the temporal derivatives p j = d d t d j j . From the periodicity of the initial profile, we know that
d j ± 1 = e ± i w d j ,
p j ± 1 = e ± i w p j
Then for the upwind discretization, the periodicity of the initial profile (42) tells that
d d t d j = p j 1 ( 3 ) , p j 1 ( 4 ) , p j ( 3 ) , p j ( 4 ) T = e i w p j ( 3 ) , e i w p j ( 4 ) , p j ( 3 ) , p j ( 4 ) T .
Hence, we can obtain the amplification matrix S up for the upwind scheme as follows,
d d t d j = 1 Δ x T up S ˜ d j = 1 Δ x S up d j , S up T up S ˜
with
T up = 0 0 e i w 0 0 0 0 e i w 0 0 1 0 0 0 0 1 .
Similarly, the temporal derivatives of d j for the central RDO-4 are
d d t d j = 1 2 p j 1 ( 3 ) + p j ( 1 ) , p j 1 ( 4 ) + p j ( 2 ) , p j ( 3 ) + p j + 1 ( 1 ) , p j ( 4 ) + p j + 1 ( 2 ) T = 1 2 e i w p j ( 3 ) + p j ( 1 ) , e i w p j ( 4 ) + p j ( 2 ) , e i w p j ( 1 ) + p j ( 3 ) , e i w p j ( 2 ) + p j ( 4 ) T .
Combining (39) and (47) we have
d d t d j = 1 Δ x T central S ˜ d j = 1 Δ x S central d j , S central : = T central S ˜
where
T central = 1 / 2 0 e i w / 2 0 0 1 / 2 0 e i w / 2 e i w / 2 0 1 / 2 0 0 e i w / 2 0 1 / 2 .

4.2. Accuracy Analysis

We can see that both the central and the upwind RDO-4 have fourth spatial local discretization errors. However, this does not suffice to ensure the third order accuracy in numerical experiments. In fact, the order accuracy is determined by the principal eigenvalues σ ( w ) according to [18], with σ ( w ) being the principle eigenvalue of the amplification matrix S . Other eigenvalues of S are spurious ones. Specifically, if σ ( w ) can be expanded as a Taylor series in w,
σ ( w ) = i w + O ( w m + 1 ) ,
then the scheme is m-th order accurate. Since it is difficult to calculate an expression manually in practice, we use the symbolic computation tools in Matlab ® .
For the upwind RDO-4, we can obtain the principal eigenvalue of the amplification matrix S upwind is
σ ( w ) = i w 1 72 w 4 + O ( w 5 ) .
This indicates that the upwind RDO-4 is third order accurate in space. Besides, the term 1 72 w 4 introduces a small amount of dissipation which tends to smear the solution, which is certificated in the numerical experiments.
The principal eigenvalue of the amplification matrix S central for central RDO-4 is
σ ( w ) = i w i 16 w 3 + O ( w 5 ) ,
and we can see that only second order spatial accuracy is maintained for the central scheme. It is also worthy noting that i 16 w 3 introduces slight dispersion, which means that phase speed varies as w changes. Therefore, during the long-term simulation of a complicated initial profile that contains waves of different frequency, the profile will be distorted finally. This is also verified in numerical experiments in Section 5.
The accuracy analysis for higher-order RDO schemes is exactly the same as RDO-4. However, a concise Taylor expansion of principal eigenvalue like (51) and (52) is very difficult (if not impossible) to obtain even using the mathematical software. Here we try to calculate the truncation order of σ ( w ) of a scheme by the following approximation. First, choose a proper w, say w = π / 4 . Then the error is
δ ( w ) = σ ( w ) + i w .
Next, we halve the wave number w and obtain the error corresponding to w / 2
δ ( w / 2 ) = σ ( w / 2 ) + i w / 2 .
Noticing that
O w / 2 1 2 m + 1 O ( w m + 1 ) ,
for the m-th order accurate scheme, one has
δ ( w / 2 ) 1 2 m + 1 δ ( w ) .
Hence the order of accuracy can be approximated by a simple logarithmic operation
m log ( | δ ( w ) | / | δ ( w / 2 ) | ) log ( 2 ) 1 .
The modulus operation is required here since the error result can be complex valued. Errors and estimated accuracy orders for the upwind and the central schemes can be seen in Table 1 and Table 2. The errors for the upwind schemes agree with the local discretization error well. The upwind scheme using M model variables per cell achieves ( M 1 ) -th accuracy order strictly. However, the accuracy orders for the central schemes varies quite irregularly. Specifically, when M = 5 and 9, the accuracy orders for central schemes remain consistent with that of upwind schemes. When M = 6 and 7, central schemes achieve higher accuracy orders than upwind schemes. Nevertheless, there also some choices of M like M = 4 and 8 such that the accuracy order of the central scheme is lower than the corresponding upwind scheme. It is remarkable that when M = 7 , the central scheme enjoys two orders higher accuracy than that of the local discretization error.
Table 1 and Table 2 and also contain abundant information about the dissipation and dispersion of the schemes. It should be noted that a positive real part of δ ( w ) will cause the scheme eventually blow up if no additional dissipation is added in time stepping. In Table 2, the real parts of δ ( w ) are negative for all upwind schemes which introduces a certain amount of spurious dissipation so as to stabilize the scheme. From Table 1 we can see that δ ( w ) are pure imaginary numbers with all choices of M, which means that the dispersion is dominant in the numerical error for the central schemes. However, it can be easily observed that the dispersion and dissipation error decays rapidly as M increases for the proposed schemes.
In (39), one can see that the numerical error in computing H 1 dominates the error of temporal derivatives, since H and V 1 D have explicit expression. However, the error of computing H 1 is inevitably increasing as M increases, which indicates that there is a limit of accuracy for the presented method. Denote the numerical inverse of H by H ^ , which is obtained by Matlab. It is well known that the error H 1 H ^ hinges on the condition number of H . Moreover, a more accurate estimate H 1 H ^ can be bounded by the following inequality,
H 1 H ^ = 1 H H H 1 H ^ 1 H I H H ^ : = E inv .
Then E inv determines the smallest error that RDO-M can reach through the grid refinement. The results of cond ( H ) and E inv with M varying form 4 to 11 are shown in Table 3. We can see that cond ( H ) and E inv grow quickly as M increases. Hence, there is a tradeoff between the higher convergence rate and E inv . In this work, we suggest that M 9 is ideal for the practical numerical experiments.

4.3. Stability Analysis

Now we investigate the maximum tolerant time step length when using the Runge–Kutta time integration scheme. Let the time step length be Δ t = c Δ x and c is the Courant number for (31). The amplification matrix for the proposed schemes using the fourth order Runge Kutta method for time integration (19) can be simplified as
R = I + c S + c 2 S 2 2 + c 3 S 3 6 + c 4 S 4 24 .
To ensure the stability, the Courant number c should be chosen such that the spectral radius of the amplification matrix R , i.e., the largest absolute value of its eigenvalues, is not greater than 1 for all wavenumbers w. The spectral radius of R for the upwind schemes and the central schemes are presented in Figure 2 and Figure 3, respectively, which can also tell the Courant number limitation related to the stability condition.
Both types of schemes have a fairly large permissible range of c. As more moments are used, the stability condition becomes more strict for both types of RDO schemes. Generally, the Courant number limitation of the central scheme is larger than that of the upwind schemes. For example, the Courant number limit for upwind RDO-5 is about 0.25 and this number is approximately 0.39 for the central RDO-5, although both schemes are fourth-order accurate. The Courant number limits for RDO schemes are presented in Table 4.

5. Numerical Results

This section provides several numerical tests to verify the formal accuracy and other numerical properties like diffusion and dispersion of the proposed schemes.

5.1. Accuracy Test

5.1.1. One-Dimensional Linear Case

We first test the accuracy when solving the scalar hyperbolic Equation (31) analysed in Section 4. The initial condition is u ( x , 0 ) = sin ( x ) on [ 0 , 2 π ) with the periodical boundary condition. The exact solution is u ( x , t ) = sin ( x t ) .
Both the central and the upwind RDO-M with the number of moments M varying from 4 to 9 are tested on uniform meshes. To measure the accuracy, two typical norms for errors are used here,
E 1 = j = 1 N | u j u j true | ,
E = max 1 j N | u j u j true | .
We study the convergence rate by recording the errors at t = 2 π and gradually refining the grids. The time step length is set as a very small number such that the temporal error does not influence the accuracy results. The error results are presented in Table 5. In general, the orders of accuracy computed from the grid refinement agree with the Fourier analysis results in Table 1 and Table 2 well, with same trends in both orders of accuracy and error magnitudes. It can be observed that the upwind RDO-M achieves ( M 1 ) -th order of accuracy for all M, which agrees with the local discretization error. On the other hand, the loss of accuracy order when M = 4 , 8 and the supravergence phenomenon when M = 6 , 7 indicated by Fourier analysis are also observed for the central RDO schemes in this test.
Besides the uniform meshes, we also consider the non-uniform meshes that are randomly generated by setting cell boundaries as
x j 1 2 = x j 1 2 + η ζ j Δ x , ζ j Unif ( 1 , 1 ) , j = 1 , , N ,
where N denotes the number of cells, Δ x = L / N , ζ j j = 1 N are independently random variables, η represents the magnitude of perturbation and Unif ( 1 , 1 ) denotes the uniform distribution over 1 , 1 . We set η = 0.1 in the present work. In this setting, the maximal and the minimal possible values of Δ x j are, respectively, 0.8 Δ x and 1.2 Δ x .
The accuracy results of RDO schemes over non-uniform meshes are listed in Table 6. It can be observed that the lack of uniformity in the mesh causes slightly larger errors and does not influence the convergence rate essentially. Therefore, the proposed schemes works consistently well on both uniform and non-conform meshes.

5.1.2. Two-Dimensional Linear Case

Consider the following two-dimensional scalar hyperbolic equation
u t + u x + u y = 0 .
The initial condition is u ( x , y ; 0 ) = sin ( x + y ) on 0 , 2 π × 0 , 2 π with the periodical boundary condition. The exact solution is u ( x , y ; t ) = sin ( x + y 2 t ) . The errors and the orders of accuracy for the central and the upwind RDO-M with M changing from 4 to 8 are shown in Table 7.
The errors and accuracy orders of the two-dimensional schemes are similar to that of the one-dimensional schemes, which illustrates that both the upwind and the central RDO schemes work well in two dimensions. Specifically, central RDO schemes achieves significantly lower errors compared with upwind RDO schemes when M = 6 , 7 . However, for any other presented values of M, upwind RDO performs better.

5.1.3. Nonlinear Case

This case considers solving the one-dimensional Burger’s equation
u t + u u x = 0 .
The initial condition is u ( x , 0 ) = 0.5 + sin x with the periodical boundary condition for x 0 , 2 π . We compute the numerical solution at t = 0.2 π when the profile is still smooth and record the errors under grid refinement. The results can be seen in Table 8. Both schemes lose accuracy due to nonlinearity compared with results of the linear cases shown in Table 5. The accuracy orders of central RDO schemes coincide the accuracy orders of local polynomial reconstruction better when M 5 . However, upwind schemes generally enjoy smaller errors in this test.

5.2. Advection of the Gaussian Profile

To investigate the diffusion and the dispersion of the proposed schemes and verify the long-term stability, a Gaussian profile used in [18]
u ( x , 0 ) = e 10 x 2
is set as the initial condition for the advection Equation (31) over the [ 1 , 1 ] with the periodic boundary condition. The number of cells is N = 20 and the time step length is Δ t = 0.005 for all schemes. The numerical solution by the central and the upwind RDO-M with various M after 4 periods and 40 periods through the domain can be seen in Figure 4 and Figure 5, respectively. Overall, the maximum value and the shape can be retained better as M increases, due to the improvement of computational accuracy. The profiles advanced by the upwind the central RDO-8 coincide the exact solution very well.
The diffusion and the dispersion agree with the Fourier analysis results in Table 1 and Table 2 well. Table 1 indicates that the error of central RDO-M schemes are dominated by dispersion for all M. By contrast, Table 2 reveals that the upwind RDO-M suffers from diffusion mainly for all M with the only exception being M = 5 . Correspondingly, the numerical solution by central RDO-4 and RDO-5 are remarkably distorted by the dispersion, especially in the long-term simulation. On the other hand, the upwind RDO schemes maintains a better and meanwhile a more flat profile than the central schemes do for all choices of M due to diffusion.

5.3. Comparison of Computational Costs

From the stability analysis, we can know that the higher-order schemes enjoy better accuracy at the expense of decreasing the CFL number and hence the efficiency. Therefore, it is necessary to compare the practical computational efficiency for RDO-M to achieve the same accuracy.
Then we consider simulating the advection of the Gaussian profile in the last test. All schemes use the fourth-order Runge–Kutta method with the temporal step length Δ t = 0.95 c max Δ x for time integration, where c max denotes the largest tolerable Courant number. Furthermore, we take Δ t = 0.7 c max for the non-uniform mesh since the stability condition is more stringent in the presence of mesh non-uniformness. The target accuracy is set as E 1 = 10 5 . The number of cells N is increased gradually until the norm-1 error reaches the target accuracy. The minimal required N and the corresponding CPU time for RDO schemes are shown in Table 9. It should be remarked that the total DOFs are not M N since the model variables at each cell boundary are shared by two cells.
For both uniform and non-uniform meshes, the significant reduction can be easily observed in the number of cells, the computational storage (i.e., total DOFs) and CPU time as M and the accuracy orders increases. RDO schemes on non-uniform meshes generally require slightly more CPU time to reach the target accuracy as shown in Table 9. Hence, high-order schemes are more attractive in terms of computational efficiency.

6. Conclusions

A family of compact multi-moment schemes that use high-order derivative moments have been proposed, analyzed and tested in one dimension and two dimensions. Using the model variables at cell boundaries, the proposed schemes reconstruct the solution polynomial within a single cell and then retrieve the temporal derivatives of model variables efficiently. The central scheme also achieves excellent performance in several tests, although this scheme does not possess clear physical meaning as the upwind scheme does. This indicates that one may extend the scheme to the more complicated scenarios without considering finding the upwind direction.
The future work will focus on extending this formulation to the more generalized unstructured triangular mesh using the Argyris polynomials [24] which also have continuous gradient globally. Furthermore, designing the conservative counterpart of the current version is also on our schedule.

Author Contributions

S.L., conceptualization, methodology, software, software, validation, writing—original draft preparation; Q.L., formal analysis, writing—review and editing; H.L., supervision; J.S., project administration, resources. All authors have read and agreed to the published version of the manuscript.

Funding

This work is supported by National Natural Science Foundation of China No. 41605070.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Wang, Z.J.; Fidkowski, K.; Abgrall, R.; Bassi, F.; Caraeni, D.; Cary, A.; Deconinck, H.; Hartmann, R.; Hillewaert, K.; Huynh, H.T.; et al. High-order CFD methods: Current status and perspective. Int. J. Numer. Methods Fluids 2013, 72, 811–845. [Google Scholar] [CrossRef]
  2. Yabe, T.; Ishikawa, T.; Wang, P.; Aoki, T.; Kadota, Y.k.; Ikeda, F. A universal solver for hyperbolic equations by cubic-polynomial interpolation II. Two-and three-dimensional solvers. Comput. Phys. Commun. 1991, 66, 233–242. [Google Scholar] [CrossRef]
  3. Yabe, T.; Xiao, F.; Utsumi, T. The constrained interpolation profile method for multiphase analysis. J. Comput. Phys. 2001, 169, 556–593. [Google Scholar] [CrossRef]
  4. Tanaka, R.; Nakamura, T.; Yabe, T. Constructing exactly conservative scheme in a non-conservative form. Comput. Phys. Commun. 2000, 126, 232–243. [Google Scholar] [CrossRef]
  5. Yabe, T.; Tanaka, R.; Nakamura, T.; Xiao, F. An exactly conservative semi-Lagrangian scheme (CIP–CSL) in one dimension. Mon. Weather Rev. 2001, 129, 332–344. [Google Scholar] [CrossRef]
  6. Xiao, F.; Akoh, R.; Ii, S. Unified formulation for compressible and incompressible flows by using multi-integrated moments II: Multi-dimensional version for compressible and incompressible flows. J. Comput. Phys. 2006, 213, 31–56. [Google Scholar] [CrossRef]
  7. Ii, S.; Shimuta, M.; Xiao, F. A 4th-order and single-cell-based advection scheme on unstructured grids using multi-moments. Comput. Phys. Commun. 2005, 173, 17–33. [Google Scholar] [CrossRef]
  8. Lin, S.; Luo, Q.; Leng, H.; Song, J. Alternating Polynomial Reconstruction Method for Hyperbolic Conservation Laws. Mathematics 2021, 9, 1885. [Google Scholar] [CrossRef]
  9. Aoki, T. Interpolated differential operator (IDO) scheme for solving partial differential equations. Comput. Phys. Commun. 1997, 102, 132–146. [Google Scholar] [CrossRef]
  10. Imai, Y.; Aoki, T.; Takizawa, K. Conservative form of interpolated differential operator scheme for compressible and incompressible fluid dynamics. J. Comput. Phys. 2008, 227, 2263–2285. [Google Scholar] [CrossRef]
  11. Xiao, F.; Ii, S. CIP/Multi-moment finite volume method with arbitrary order of accuracy. arXiv 2012, arXiv:1207.6822. [Google Scholar]
  12. Toro, E.F.; Titarev, V.A. Derivative Riemann solvers for systems of conservation laws and ADER methods. J. Comput. Phys. 2006, 212, 150–165. [Google Scholar] [CrossRef]
  13. Ii, S.; Xiao, F. High order multi-moment constrained finite volume method. Part I: Basic formulation. J. Comput. Phys. 2009, 228, 3669–3707. [Google Scholar] [CrossRef] [Green Version]
  14. Deng, X.; Xie, B.; Xiao, F. Multimoment finite volume solver for euler equations on unstructured grids. AIAA J. 2017, 55, 2617–2629. [Google Scholar] [CrossRef]
  15. Cheng, L.; Deng, X.; Xie, B.; Jiang, Y.; Xiao, F. Low-dissipation BVD schemes for single and multi-phase compressible flows on unstructured grids. J. Comput. Phys. 2021, 428, 110088. [Google Scholar] [CrossRef]
  16. Cockburn, B.; Shu, C.W. TVB Runge–Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework. Math. Comput. 1989, 52, 411–435. [Google Scholar]
  17. Wang, Z.J. Spectral (finite) volume method for conservation laws on unstructured grids. basic formulation: Basic formulation. J. Comput. Phys. 2002, 178, 210–251. [Google Scholar] [CrossRef] [Green Version]
  18. Huynh, H.T. A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods. In Proceedings of the 18th AIAA Computational Fluid Dynamics Conference, Miami, FL, USA, 25–28 June 2007; p. 4079. [Google Scholar]
  19. Taylor, M.; Tribbia, J.; Iskandarani, M. The spectral element method for the shallow water equations on the sphere. J. Comput. Phys. 1997, 130, 92–108. [Google Scholar] [CrossRef]
  20. Süli, E.; Mayers, D.F. An Introduction to Numerical Analysis; Cambridge University Press: Cambridge, UK, 2003. [Google Scholar]
  21. Patera, A.T. A spectral element method for fluid dynamics: Laminar flow in a channel expansion. J. Comput. Phys. 1984, 54, 468–488. [Google Scholar] [CrossRef]
  22. Sakurai, K.; Aoki, T.; Lee, W.H.; Kato, K. Poisson equation solver with fourth-order accuracy by using interpolated differential operator scheme. Comput. Math. Appl. 2002, 43, 621–630. [Google Scholar] [CrossRef] [Green Version]
  23. Zhang, S. On the full C1-Qk finite element spaces on rectangles and cuboids. Adv. Appl. Math. Mech 2010, 2, 701–721. [Google Scholar] [CrossRef] [Green Version]
  24. Ciarlet, P.G. The Finite Element Method for Elliptic Problems; SIAM: Philadelphia, PA, USA, 2002. [Google Scholar]
Figure 1. The illustration of model variables and solution points used in RDO-4 and RDO-5 in the two-dimensional space.
Figure 1. The illustration of model variables and solution points used in RDO-4 and RDO-5 in the two-dimensional space.
Axioms 10 00295 g001
Figure 2. Contour plots of spectral radii of amplification matrices of the upwind schemes.
Figure 2. Contour plots of spectral radii of amplification matrices of the upwind schemes.
Axioms 10 00295 g002
Figure 3. Contour plots of spectral radii of amplification matrices of the central schemes.
Figure 3. Contour plots of spectral radii of amplification matrices of the central schemes.
Axioms 10 00295 g003
Figure 4. Numerical results at t = 8 (4 periods) for the advection equation with the bell function as the initial profile.
Figure 4. Numerical results at t = 8 (4 periods) for the advection equation with the bell function as the initial profile.
Axioms 10 00295 g004aAxioms 10 00295 g004b
Figure 5. Numerical results at t = 80 (40 periods) for the advection equation with the bell function as the initial profile.
Figure 5. Numerical results at t = 80 (40 periods) for the advection equation with the bell function as the initial profile.
Axioms 10 00295 g005
Table 1. Errors and accuracy orders for central RDO schemes via Fourier analysis.
Table 1. Errors and accuracy orders for central RDO schemes via Fourier analysis.
RDO-MCoarse Mesh Error ( w = π / 4 )Fine Mesh Error ( w = π / 8 )Order
M = 5 −7.88 × 10 7 + 3.16 × 10 6 i−1.24 × 10 8 + 1.00 × 10 7 i4.02
M = 6 −3.14 × 10 5 − 4.27 × 10 6 i−5.04 × 10 7 − 3.40 × 10 8 i4.97
M = 7 −1.36 × 10 7 + 2.85 × 10 7 i−5.42 × 10 10 + 2.34 × 10 9 i6.04
M = 8 −1.00 × 10 7 − 1.00 × 10 8 i−3.98 × 10 10 − 1.98 × 10 11 i6.98
M = 9 −2.42 × 10 10 + 5.31 × 10 10 i−2.34 × 10 13 + 1.07 × 10 12 i8.05
Table 2. Errors and accuracy orders for upwind RDO schemess via Fourier analysis.
Table 2. Errors and accuracy orders for upwind RDO schemess via Fourier analysis.
RDO-MCoarse Mesh Error ( w = π / 4 )Fine Mesh Error ( w = π / 8 )Order
M = 5 3.36 × 10 6 i1.02 × 10 7 i4.16
M = 6 1.17 × 10 5 i8.71 × 10 8 i6.08
M = 7 8.97 × 10 9 i1.74 × 10 11 i8.00
M = 8 −2.17 × 10 7 i−1.76 × 10 9 i5.95
M = 9 6.46 × 10 10 i1.13 × 10 12 i8.15
Table 3. cond ( H ) and the E recon as the function of M.
Table 3. cond ( H ) and the E recon as the function of M.
M4567891011
cond ( H ) 9.921.40 × 10 2 4.13 × 10 2 1.49 × 10 3 5.25 × 10 4 9.24 × 10 5 3.71 × 10 7 5.12 × 10 8
E inv 1.60 × 10 16 4.55 × 10 16 1.21 × 10 15 1.39 × 10 15 3.24 × 10 14 3.17 × 10 13 4.56 × 10 12 1.24 × 10 11
Table 4. The Courant number limit for RDO schemes.
Table 4. The Courant number limit for RDO schemes.
M 4567
c max Upwind0.450.250.240.12
Central0.690.390.330.24
Table 5. Accuracy results from grid convergence on uniform meshes for the linear advection equation.
Table 5. Accuracy results from grid convergence on uniform meshes for the linear advection equation.
RDO-MN E 1 Order E Order E 1 Order E Order
Central SchemesUpwind Schemes
M = 4 103.40 × 10 2 -5.25 × 10 2 -1.31 × 10 2 -2.02 × 10 2 -
208.21 × 10 3 2.051.30 × 10 2 2.021.70 × 10 3 2.952.65 × 10 3 2.93
303.67 × 10 3 1.985.75 × 10 3 2.015.06 × 10 4 2.987.93 × 10 4 2.98
M = 5 102.40 × 10 4 -3.71 × 10 4 -2.28 × 10 4 -3.49 × 10 4 -
201.37 × 10 5 4.132.17 × 10 5 4.101.42 × 10 5 3.952.12 × 10 5 3.97
302.70 × 10 6 4.004.24 × 10 6 4.022.95 × 10 6 4.004.12 × 10 6 3.97
M = 6 101.20 × 10 5 -1.20 × 10 5 -5.26 × 10 5 -8.13 × 10 5 -
201.79 × 10 7 4.071.79 × 10 7 6.071.68 × 10 6 4.972.64 × 10 6 4.94
301.58 × 10 8 5.991.58 × 10 8 5.992.23 × 10 7 4.993.49 × 10 7 4.99
M = 7 101.57 × 10 7 -2.42 × 10 7 -3.88 × 10 7 -5.99 × 10 7 -
206.21 ×   10 10 7.989.72 ×   10 10 7.966.24 × 10 9 5.969.68 × 10 9 5.95
302.80 ×   10 11 7.646.07 ×   10 11 6.845.48 ×   10 10 6.008.59 ×   10 10 5.97
M = 8 102.96 × 10 7 -4.58 × 10 7 -1.06 × 10 7 -1.64 × 10 7 -
204.80 × 10 9 5.957.59 × 10 9 5.918.48 ×   10 10 6.971.33 × 10 9 6.94
304.22 ×   10 10 6.006.62 ×   10 10 6.024.98 ×   10 11 6.997.82 ×   10 11 7.00
M = 9 51.89 × 10 7 -2.93 × 10 7 -1.81 × 10 7 -2.71 × 10 7 -
105.38 ×   10 10 8.468.31 ×   10 10 8.465.85 ×   10 10 8.289.24 ×   10 10 8.20
151.96 ×   10 11 8.173.07 ×   10 11 8.132.24 ×   10 11 8.053.49 ×   10 11 8.08
Table 6. Accuracy results from grid convergence on non-uniform meshes generated by (62) for the linear advection equation.
Table 6. Accuracy results from grid convergence on non-uniform meshes generated by (62) for the linear advection equation.
RDO-MN E 1 Order E Order E 1 Order E Order
Central SchemesUpwind Schemes
M = 4 103.99 × 10 2 -5.80 × 10 2 -1.39 × 10 2 -2.19 × 10 2 -
208.97 × 10 3 2.151.20 × 10 2 2.271.84 × 10 3 2.922.87 × 10 3 2.93
304.05 × 10 3 1.965.50 × 10 3 1.925.49 × 10 4 2.988.57 × 10 4 2.98
M = 5 102.45 × 10 4 -3.81 × 10 4 -2.64 × 10 4 -4.09 × 10 4 -
201.43 × 10 5 4.102.29 × 10 5 4.061.43 × 10 5 4.212.34 × 10 5 4.13
302.79 × 10 6 4.024.51 × 10 6 4.013.06 × 10 6 3.804.67 × 10 6 3.98
M = 6 101.36 × 10 5 -2.45 × 10 5 -6.04 × 10 5 -9.24 × 10 5 -
204.11 × 10 7 5.059.03 × 10 7 4.761.81 × 10 6 5.062.90 × 10 6 4.99
303.57 × 10 8 6.037.54 × 10 8 6.122.37 × 10 7 5.013.68 × 10 7 5.09
M = 7 102.05 × 10 7 -4.64 × 10 7 -4.17 × 10 7 -6.14 × 10 7 -
202.30 × 10 9 6.484.96 × 10 9 6.557.04 × 10 9 5.891.02 × 10 8 5.91
301.71 ×   10 10 6.403.58 ×   10 10 6.485.88 ×   10 10 6.129.65 ×   10 10 5.82
M = 8 103.23 × 10 7 -5.94 × 10 7 -1.41 × 10 7 -2.19 × 10 7 -
205.34 × 10 9 5.929.25 × 10 9 6.001.02 × 10 9 7.111.62 × 10 9 7.08
304.55 ×   10 10 6.077.58 ×   10 10 6.176.13 ×   10 11 6.949.62 ×   10 11 6.96
M = 9 52.07 × 10 7 -3.33 × 10 7 -2.03 × 10 7 -3.27 × 10 7 -
105.83 ×   10 10 8.479.50 ×   10 10 8.457.11 ×   10 10 8.151.12 × 10 9 8.19
152.16 ×   10 11 8.123.40 ×   10 11 8.212.50 ×   10 11 8.263.71 ×   10 11 8.41
Table 7. Accuracy results on uniform meshes for the 2D linear advection equation u t + u x + u y = 0 at t = 2 π .
Table 7. Accuracy results on uniform meshes for the 2D linear advection equation u t + u x + u y = 0 at t = 2 π .
RDO-M N x × N y E 1 Order E Order E 1 Order E Order
Central SchemesUpwind Schemes
M = 4 10 × 10 6.76 × 10 2 -1.04 × 10 1 -2.64 × 10 2 -4.08 × 10 2 -
20 × 20 1.64 × 10 2 2.042.59 × 10 2 2.013.43 × 10 3 2.955.36 × 10 3 2.93
30 × 30 7.34 × 10 3 1.991.15 × 10 2 2.011.01 × 10 3 3.001.59 × 10 3 2.99
M = 5 10 × 10 7.48 × 10 4 -1.16 × 10 3 -4.37 × 10 4 -6.75 × 10 4 -
20 × 20 4.42 × 10 5 4.086.95 × 10 5 4.062.79 × 10 5 3.974.33 × 10 5 3.96
30 × 30 8.69 × 10 6 4.011.36 × 10 5 4.025.46 × 10 6 4.028.56 × 10 6 4.00
M = 6 10 × 10 3.78 × 10 5 -5.84 × 10 5 -1.04 × 10 4 -1.60 × 10 4 -
20 × 20 5.60 × 10 7 6.088.86 × 10 7 6.043.35 × 10 6 4.955.29 × 10 6 4.92
30 × 30 4.96 × 10 8 5.987.78 × 10 8 6.004.44 × 10 7 4.986.97 × 10 7 5.00
M = 7 10 × 10 2.56 × 10 6 -3.95 × 10 6 -4.91 × 10 6 -7.59 × 10 6 -
20 × 20 1.11 × 10 8 7.841.75 × 10 8 7.824.45 × 10 8 6.796.97 × 10 8 6.77
30 × 30 4.60 ×   10 10 7.867.41 ×   10 10 7.802.88 × 10 9 6.754.67 × 10 9 6.67
M = 8 5 × 53.74 × 10 4 -5.78 × 10 4 -4.13 × 10 5 -6.38 × 10 5 -
10 × 10 6.05 × 10 6 5.959.34 × 10 6 5.953.31 × 10 7 6.965.12 × 10 7 6.96
15 × 15 5.21 × 10 7 6.048.16 × 10 7 6.011.93 × 10 8 7.013.03 × 10 8 6.97
Table 8. Accuracy results for the Burger’s equation with u ( x , 0 ) = 0.5 + sin x at t = 0.2 π .
Table 8. Accuracy results for the Burger’s equation with u ( x , 0 ) = 0.5 + sin x at t = 0.2 π .
RDO-MN E 1 Order E Order E 1 Order E Order
Central SchemesUpwind Schemes
M = 4 102.33 × 10 3 -1.04 × 10 1 -8.38 × 10 4 -4.64 × 10 3 -
207.17 × 10 4 1.702.59 × 10 2 2.011.30 × 10 4 2.691.07 × 10 3 2.12
303.23 × 10 4 1.971.15 × 10 2 2.014.22 × 10 5 2.773.98 × 10 4 2.44
M = 5 103.69 × 10 4 -1.71 × 10 3 -2.63 × 10 5 -9.20 × 10 5 -
201.95 × 10 5 4.241.53 × 10 4 3.483.41 × 10 6 2.951.35 × 10 5 2.77
302.17 × 10 6 5.422.47 × 10 5 4.497.21 × 10 7 3.834.70 × 10 6 2.61
M = 6 101.68 × 10 4 -5.82 × 10 4 -1.78 × 10 5 -1.05 × 10 4 -
205.09 × 10 6 5.043.08 × 10 5 4.241.41 × 10 6 3.661.67 × 10 5 2.65
305.46 × 10 7 5.515.27 × 10 6 4.362.48 × 10 7 4.283.64 × 10 6 3.76
M = 7 102.03 × 10 5 -1.65 × 10 4 -9.18 × 10 6 -3.19 × 10 5 -
203.63 × 10 7 5.803.18 × 10 6 5.691.42 × 10 7 6.011.08 × 10 6 4.89
303.14 × 10 8 6.042.03 × 10 7 6.791.47 × 10 8 5.609.37 × 10 8 6.02
M = 8 102.72 × 10 5 -1.33 × 10 4 -1.62 × 10 7 -9.62 × 10 7 -
207.96 × 10 8 8.421.03 × 10 6 7.021.42 × 10 8 3.511.60 × 10 7 2.59
307.20 × 10 9 5.939.95 × 10 8 5.762.17 × 10 9 4.632.60 × 10 8 4.49
M = 9 101.72 × 10 6 -1.20 × 10 5 -5.30 × 10 8 -4.72 × 10 7 -
203.20 × 10 8 5.741.76 × 10 7 6.089.41 × 10 10 5.811.26 × 10 8 5.23
305.73 ×   10 10 9.922.68 × 10 9 10.326.01 × 10 11 6.787.08 ×   10 10 7.10
Table 9. The minimal required number of cells to reach the accuracy E 1 10 5 and the corresponding resumed CPU time for RDO-M.
Table 9. The minimal required number of cells to reach the accuracy E 1 10 5 and the corresponding resumed CPU time for RDO-M.
RDO-M M = 4 M = 5 M = 6 M = 7 M = 4 M = 5 M = 6 M = 7
Mesh Type Central SchemesUpwind Schemes
UniformMinimal N1129674827271624422
Total DOFs225820114410854218613288
CPU time (s)55.80.460.250.146.310.480.370.14
Non-uniformMinimal N1187.273.451.630.0292.471.249.827.2
Total DOFs2374.4220.2154.8120.0584.8213.6149.4108.8
CPU time (s)70.460.650.370.197.930.620.470.18
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Lin, S.; Luo, Q.; Leng, H.; Song, J. Reconstructed Interpolating Differential Operator Method with Arbitrary Order of Accuracy for the Hyperbolic Equation. Axioms 2021, 10, 295. https://doi.org/10.3390/axioms10040295

AMA Style

Lin S, Luo Q, Leng H, Song J. Reconstructed Interpolating Differential Operator Method with Arbitrary Order of Accuracy for the Hyperbolic Equation. Axioms. 2021; 10(4):295. https://doi.org/10.3390/axioms10040295

Chicago/Turabian Style

Lin, Shijian, Qi Luo, Hongze Leng, and Junqiang Song. 2021. "Reconstructed Interpolating Differential Operator Method with Arbitrary Order of Accuracy for the Hyperbolic Equation" Axioms 10, no. 4: 295. https://doi.org/10.3390/axioms10040295

APA Style

Lin, S., Luo, Q., Leng, H., & Song, J. (2021). Reconstructed Interpolating Differential Operator Method with Arbitrary Order of Accuracy for the Hyperbolic Equation. Axioms, 10(4), 295. https://doi.org/10.3390/axioms10040295

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop