Alternating Polynomial Reconstruction Method for Hyperbolic Conservation Laws

: We propose a new multi-moment numerical solver for hyperbolic conservation laws by using the alternating polynomial reconstruction approach. Unlike existing multi-moment schemes, our approach updates model variables by implementing two polynomial reconstructions alternately. First, Hermite interpolation reconstructs the solution within the cell by matching the point-based variables containing both physical values and their spatial derivatives. Then the reconstructed solution is updated by the Euler method. Second, we solve a constrained least-squares problem to correct the updated solution to preserve the conservation laws. Our method enjoys the advantages of a compact numerical stencil and high-order accuracy. Fourier analysis also indicates that our method allows a larger CFL number compared with many other high-order schemes. By adding a proper amount of artiﬁcial viscosity, shock waves and other discontinuities can also be computed accurately and sharply without solving an approximated Riemann problem.


Introduction
In this paper, we aim to solve the hyperbolic conservation law ∂u ∂t where u and f (u) can be either scalars or vectors. The classical higher-order numerical schemes such as the finite difference (FD) and the finite volume (FV) frameworks usually involve only one degree of freedom (DOF) per cell. That means to achieve the kth order discretization of the first-order spatial derivative, a stencil of at least k + 1 cells is required. The compactness of a scheme can help simplify the polynomial interpolation in unstructured meshes and reduce communication costs between computational nodes on modern supercomputers. In this paper, the "compact" is used specifically for numerical schemes with a minimal stencil, which excludes so-called Padê schemes [1,2] since its approximation of the derivative involves solving a diagonal system of the full difference stencil.
A natural thought for the design of a compact high-order scheme is involving more than one DOF per cell so as to construct higher-order polynomials under the same stencil width. A series of multi-moment schemes is developed based on this idea. The constrained interpolation profile (CIP) scheme [3,4], which uses the semi-Lagrangian method to advance in time, individually updates the physical variable and its derivatives on each cell by the governing equations in different forms. The interpolated differential operator (IDO) method [5] is similar to CIP but uses an explicit time-advancing approach. For the IDO method, a stencil containing three adjacent cells suffices to approximate the first-order spatial derivative with fifth-order accuracy in one dimension, which makes it possible to obtain comparably accurate results with spectral methods in direct numerical simulation [6]. For IDO, updating the derivative-value variables induces the problem of computing the time derivative of u x , which involves f uu . However, for the system of conservation laws, f uu has a fairly complicated form and can be difficult to compute. Goodrich et al. proposed a similar Hermite method [7] on staggered grids.
Neither CIP nor IDO are conservative schemes, which means they are not ideal for shock-capturing problems and long-term simulations. To address the non-conservativeness, so-called CIP-conservative semi-Lagrangian (CIP-CSL) schemes [8,9] and the IDO scheme in conservative form (IDO-CF) [10,11] add the cell-averaged value as an extra independent model variable to evolve. The cell-averaged variable is advanced by the flux-form conservation law and is hence exactly conserved, while other kinds of variables are updated using the differential form of conservation laws. However, these conservative CIP algorithms need to find polynomials that match cell-averaged variables across cells, which can be time-consuming when encountered with non-uniform meshes. A framework termed the multi-moment finite volume method (MM-FVM) has been proposed on the basis of CIP-CSL. Compared with CIP-CSL, MM-FVM can construct the high-order interpolation on a local base and is more flexible when treating unstructured grids, as shown in [12]. Nonetheless, MM-FVM uses the semi-Lagrangian method in the time integration of point-based variables, which is difficult to implement in the multi-dimensional case.
Another class of high-order conservative schemes based on local flux reconstructions also uses a compact stencil and has become increasingly attractive in recent decades. Examples include the well-known discontinuous Galerkin (DG) [13][14][15], spectral volume [16,17], and spectral difference [18] methods and, more recently, the flux reconstruction/correction via procedure (FR/CPR) method [19,20]. These methods evolve k DOFs in each cell for a kth order spatial accuracy (in one dimension) and compute the numerical flux on cell boundaries to represent the interaction between adjacent cells. Thus, polynomial interpolation across cells is not needed, which makes these methods compact and flexible in handling unstructured meshes. The multi-moment constrained finite volume (MCV) method proposed by Ii and Xiao [21] can be viewed as a combination of CIP and the flux reconstruction approach. The MCV method introduces high-order spatial derivatives of numerical flux and reconstructs the local flux function more directly. MCV is generalized as multi-moment constrained flux reconstruction (MMC-FR) [22]. However, the nonlinear aliasing phenomenon [23,24] in treating nonlinear flux functions, which causes instability and accuracy loss, is one typical demerit for these methods based on local flux reconstruction. This paper presents a novel conservative multi-moment approach termed the AltPoly method, which uses the alternating polynomial reconstruction as the time integration method. AltPoly divides the computational domain into cells and adopts point-based and cell-averaged values as model variables. The point-based variables including the physical variable and its spatial derivatives are updated together, which is different from the conventional multi-moment methods that update different types of model variables according to governing equations in different forms. AltPoly then solves a constrained least-squares problem to reconstruct the solution within the cell and update the model variables, which preserves conservativeness exactly. Solving Riemann problems on cell boundaries is not needed. Furthermore, AltPoly can also handle non-uniform meshes in one dimension efficiently with the aid of the coordinate transform.
The remainder of this paper is organized as follows: Section 2 describes the formulation of AltPoly and constructs artificial viscosity for discontinuous problems. Section 3 studies the accuracy and stability of AltPoly by Fourier analysis. Section 4 validates the numerical convergence rate and capability of solving some widely used benchmark problems involving both the scalar equation and Euler systems in 1D. Finally, a short summary in Section 5 ends this paper.

Formulation of AltPoly Method
This section gives the formulation of the AltPoly method for a one-dimensional conservation law. For concreteness and ease of comprehension, we only present the detailed formulation for AltPoly with 2 point-value variables and one extra cell-averaged variable per cell, which is termed AltPoly-3. As for the general AltPoly-r with r model variables within a cell, the formulation is similar to AltPoly-3 and is hence briefly illustrated in Appendix A.
To begin with, we first introduce some notation. The bold font lower-case letter denotes a column vector, e.g., a, u. The transpose a is denoted by a T . e k denotes the kth standard unit vector in a space in proper dimensions. [a k ] m k=1 denotes a column vector in the form of [a 1 , . . . , a m ] T . Bold font capital letters such as A denote matrices. The submatrix of A constructed from rows {r 1 , r 2 , . . . , r k } and columns {c 1 , c 2 , . . . , c l } is denoted as A[r 1 , r 2 , . . . , r k ; c 1 , c 2 , . . . , c l ]. a denotes the 2-norm of the vector a, while A for a matrix A is defined by A = sup x =0 Ax x .

Spatial Discretization and Model Variables
Suppose the spatial domain [x min , x max ] is divided into N non-overlapping cells, among which the j-th cell is I j := [x j−1/2 , x j+1/2 ], j = 1, . . . , N. It should be noted that the cells do not necessarily have a uniform length. In the j-th element, we consider computing the time integration of three model variables: u, u x at the middle point of the element, and additionally the element-integrated valueū j = 1

Updating Model Variables
Denote the model variables at time level t = t 0 as u This part illustrates the procedure to obtain u where ∆t is the time step length. In a nutshell, this procedure consists of 3 steps:
Computing solution values on selected solution points based on reconstructed polynomials and updating solution values and cell-averaged variables via the Euler forward method; 3.
Reconstructing the solution polynomial that matches updated solution values under the constraint of cell-averaged variables to preserve conservativeness.
Details are presented in the following and algorithm illustrations are presented in Figure 1 for convenience of illustration.

Step 1: Hermite Interpolation across Cells
In the first step, we find the interpolation polynomial h j− 1 Figure 1a. This step does not involve cell-averaged variablesū j . This interpolation polynomial is of the third order since there are 4 conditions in total. Before constructing h j− 1 2 (x), we introduce a local coordinate for x ∈ [x j , x j+1 ]as is usually done in the finite element method [25]. The local coordinate can unify the interpolation templates so as to reduce the computational cost. The polynomial h j− 1 2 (x) on the local coordinate s is denoted byh j− 1 2 (s) := h j− 1 2 (x(s)). Notice that ∂h ∂s and ∂h ∂x can be connected by the following rule: Then the coefficients ofh j− 1 2 (s) = ∑ 3 i=0 a i s i can be determined by the following system: The matrix form is where In practice, the inverse of H is computed and stored in advance so as to reduce the computation cost in solving (4).

Step 2: Updating Solution Values on Solution Points
Based on Step 1, we can now approximate and update the solution in cell I j using polynomials h j− 1 2 (x) and h j+ 1 2 (x). For AltPoly-3, four symmetric solution points x j,k = x j + ξ k ∆x j /2, k = 1, . . . , 4 on I j are needed in this step with ξ k ∈ [−1, 1]. In this paper, we choose ξ 1 = −1, ξ 2 = −1 + σ, ξ 3 = 1 − σ and ξ 4 = 1 as solution points with σ = 0.03. To update the solution values on these solution points, we need to compute to approximate u and u x . For convenience, we transfer the solution points ξ k into local coordinate s(ξ k ) for k = 1 . . . , 4. Combining the conservative law (1) we can then compute u j,k by the Euler forward method: where the function f u stands for ∂u . The cell-averaged variableū (t 1 ) j is updated by the flux form of the conservation law: where u j+ 1 2 is the function value at the cell boundary x j+1/2 computed by When the computing grid is uniform, (10) is simply u j+1/2 =h j+ 1 2 (0). Figure 1b describes this step. Remark 1. The Euler time integration used for updating solution values on solution points is simple and easy to implement. However, the semi-Lagrangian time integration scheme [26] might be robust with a larger time step.

Step 3: Updating Model Variables
The cell-averaged variableū j has been updated (9). As for the point-value model variables, the four updated solution values on I j suffice to reconstruct a Lagrange interpolation polynomial l j (x), so as to obtain the updated model variables u j . However, such a direct routine may cause the mismatch of u j andū j . Similar to Step 1, we introduce a local coordinate ξ ∈ [−1, 1]: This is an over-determined linear system since it contains 5 conditions and 4 unknowns. However, (13) must hold for ensuring the conservation law. In other words, we need to find the optimal coefficients {b i } i=1,...,4 to fit (12) under the constraint (13). This is a typical constrained least-squares problem. To solve this, we rewrite the constraint (13) as and plug it into (12). Then the constrained least-squares problem is cast into the standard least-squares problem as where Then the solution of (15) is where A † is the generalized inverse or the Moore-Penrose inverse [27] of non-square matrix A.
From the polynomial expression of l j (x), we can now retrieve the point-value model variables for the next time level: So far, (9), (17) and (18) together give the updating rule for all model variables. This polynomial reconstruction of l j (x) based on solving constrained squares problem (13) can be efficiently implemented if we compute and store A † in advance.

Runge-Kutta Time Integration
Up to this point, we have only presented AltPoly using the Euler forward method, which is merely of first-order temporal accuracy. Since the presented AltPoly is not in the semi-discretization form, the Runge-Kutta algorithm cannot be directly applied for AltPoly. Therefore, a slight modification of the Runge-Kutta algorithm is adopted here.
To proceed, we introduce the following notation. Let u be a compounded list of all model variables, and its j-th entry is a vector u j = u j , (u x ) j ,ū j . Denote the result of AltPoly upon u after one step with length ∆t of the temporal forward Euler method by E (u, ∆t) and the corresponding increment part is denoted by Then, one time step of the Runge-Kutta procedure for AltPoly is where

Artificial Viscosity
To tackle problems that involve discontinuities such as shocks and contact discontinuities, we add artificial viscosities to the proposed method. The critical step is depicting the smoothness of a cell and determining whether the discontinuity exists. Our shock sensor borrows the idea from the sub-cell shock capturing for DG [28].
Within each cell, we define the following truncated cell average based on the moments for AltPoly-r: where 2k is the largest even number smaller than r. It can be seen thatū trun j utilizes the even-order partial differential point-value moments and is an approximation ofū j with the error of magnitude O ∆x 2k+2 j . Therefore, we expect that the solution is continuous and thenū trun j is close enough toū j , and hence declaim that a discontinuity exists in I j if ū j −ū trun j exceeds some threshold. Based on this, the following smoothness indicator is defined: where u range = u max − u min and the parameter = 10 −10 is added to avoid dividing zero. This smoothness indicator has a simple form and can be calculated directly from the model variables. Then the amount of viscosity at x j+ 1 2 is calculated by the following smooth function, where ϑ j = log 10 θ j , the parameters ϑ ∼ (∆x j ) 2k ,ν = ∆x j /4, and κ is large enough to maintain a sharp but smooth shock profile. So far, we have obtained the viscosity at the middle of each cell. Next we use the linear interpolation to construct the viscosity function between the middle points of two adjacent cells. Then, to solve the problem with this viscosity function, one only needs to replace (8) and (9), respectively.

Fourier (von Neumann) Analysis
This part gives the Fourier analyses of stability and accuracy for AltPoly in solving a linear advection equation. We only give a detailed analysis for AltPoly-3 for concreteness. AltPoly-r with other choices of r can be analyzed following a similar routine.

Formulation of the Amplification Matrix
Consider the linear scalar equation Suppose the computational domain is [0, L] with a uniform mesh spacing ∆x. The initial condition is periodic: u init (x) = e iwx/∆x where i denotes the imaginary unit, and The initial values of model variables on the cell I j = [x j−1/2 , x j+1/2 ] and its neighbor cells are given by The discrete Fourier transform of the series u j is u k is not 0 only when k = k under the given initial condition. Hence, we only need to consider the Fourier coefficientû k and denote it asû in the following text for simplicity.û x andû are defined similarly for series (u x ) j and ū j .
The key procedure in our analysis is rewriting the AltPoly-3 formulation in Section 2 into the matrix form: j ], and S denotes the amplification matrix connectingû (t 1 ) at the next time level t = t 1 = t 0 + ∆t.
The coefficient vector a j±1/2 of the reconstructed polynomialh(s) in Step 1 of AltPoly-3 is According to (26), we have d j+1/2 = e iw d j−1/2 and a j+1/2 = e iw a j−1/2 according to the periodicity of the spatial profile. According to (10), the function values at the cell boundaries x j±1/2 are where e k denotes the k-th standard unit vector. Now we turn to Step 2 of AltPoly-3. As the solution points x j,1 and x j,2 are in the domain of h j−1/2 (x), we transfer them into the local coordinate in [x j−1 , x j ]. Let s 1 = s(x j,1 ) = 0 and s 2 = s(x j,2 ) = α. Then we can derive where s 3 = −α and s 4 = 0.
Combining a j+1/2 = e iw a j−1/2 , the updating of solution values on solution points can be expressed in the following matrix form: where the items of B, C are, respectively: with 0 0 defined as 1, and W is a diagonal matrix with 1, 1, e iw , e iw on the diagonal line.
Since f u = 1, the updated solution values here are computed by where the CFL number c = |∂ u f (u)|∆t/∆x = ∆t/∆x for this linear problem. As for Step 3, the updated cell average over the jth cell is Solving a constrained least-squares problem we have From b we can easily obtain point-value variables of the next time level: Gathering (28), (30), (35), (37) and (38), we finally obtain the amplification matrix S in (27): which is the basis for the accuracy and the stability analysis as follows.

Accuracy Analysis
At the beginning, we assumed AltPoly-3 only had third-order spatial accuracy since all reconstruction polynomials are of the fourth-order, which has third-order spatial accuracy for first-order derivatives. Instead, however, we observed fourth-order spatial accuracy in numerical tests. Analyzing only truncation error is not enough to explain this phenomenon. We believe that the cell-averaged momentū j and the constrained polynomial reconstruction help improve the spatial accuracy.
The amplification matrix S can be divided into two parts, where S 0 and S 1 are independent from c. Letû = [1, iw/2, 1 iw (e iw/2 − e −iw/2 )]. The principal eigenvalue of S 0 is 1 (see Appendix B for details). We denote the corresponding eigenvector as ψ, namely and we make the last items of ψ andû equal, namely ψ 3 = 1 iw (e iw/2 − e −iw/2 ). With the assistance of Wolfram Mathematica 12, we decompose the following terms into Taylor series ψ =û + ε 1 w 4 e 1 + O(w 5 ) and where the scalars ε 1 , ε 2 , and the vectors v 1 , v 2 are independent with w.
Suppose we choose a proper c that ensures the stability, namely S i ≤ M for i ∈ N and some positive constant M. Divide the time domain [0, T] uniformly with step length ∆t and denote n = T/∆t. It is obvious that the exact solution at t = T for (53) with the initial condition u init (x) = e iwx/∆x is u exact = e iw(x−T)/∆x = e −incw u init (x). Hence the error of the numerical solution by AltPoly-3 is However, According to the bound result (A16) in Appendix B, S k e 2 ≤ µ k 1 + (k − 1)cµ k 2 + O(w) with some positive constants µ 1 , µ 2 ∈ (0, 1). Then one can deduce that Remember that w ∼ ∆x, c = ∆t/∆x, and n = T/∆t, then we have Combining (44)-(47), we can see that the error order of AltPoly-3 is To ensure stability, the order of magnitude of ∆t is usually not higher than that of ∆x, namely c ∼ O(1). Therefore, the truncated error (48) can be simply reduced to O(∆t + ∆x 4 ). That means AltPoly-3 using the Euler time scheme has fourth-order spatial accuracy and first-order temporal accuracy.
When using RK4 as time integration strategy (20), the amplification matrix S RR4 for this linear equation can be expressed as where Expanding S RK4 , we have One can easily show that S RK4 has fifth-order truncated temporal accuracy and then the total error becomes

Stability Analysis
To obtain stable time integration, the CFL number c should be carefully chosen such that the magnitude of the principal eigenvalue of the amplification matrix is not greater than 1 for all wavenumbers k. The principal eigenvalue of the amplification matrix S RK4 for AltPoly with various moments can be seen in Figure 2. It is remarkable that all schemes have a fairly large range of c except for AltPoly-5. In particular, AltPoly-3, which has fourth-order spatial accuracy, is stable even with c > 1, which is better than many other existing explicit high-order schemes using a compact stencil. As more moments are used, the stability condition becomes more stringent. For AltPoly-6, the tolerated CFL number is approximately c ∈ [0.2, 0.4]. However, the maximum principal eigenvalue is around 1 + 10 −2 when c ≤ 0.2 for AltPoly-6, which means short-term time integration can still maintain stability.

Numerical Results
In this part, we solve several types of one-dimensional hyperbolic law to illustrate the performance of the proposed method. We choose the four-stage Runge-Kutta method (20) for the time discretization.
For all numerical experiments, AltPoly-r with the moment number r varying from 3 to 6 are tested on uniform meshes. In addition, we check the accuracy of the proposed method on non-uniform meshes. Non-uniform grids are generated by adding a 40% random perturbation upon the element boundary of the uniform mesh. Specifically, suppose x j+1/2 and ∆x are separately the cell boundary and cell length for the uniform grid, then the element boundary generated for the non-uniform grid is where ζ j are independently identically distributed random variables, Unif(−1, 1) denotes the uniform distribution over [−1, 1], and β is the amount of perturbation, which is set as 0.2 in our experiments. Special attention might be paid to specify the initial values of model variables. This is trivial if the initial profile is analytically given. Otherwise, these unknowns can be specified by implementing high-order interpolations and numerical integration.
Since the temporal accuracy is only fourth-order, we take ∆t ∼ (∆x) (2r−2)/4 for AltPoly-r in the accuracy tests.

Linear Scalar Equation
In this subsection, we test the performance of AltPoly schemes on the following scalar equation [29]: We first check the accuracy of the proposed method through grid refinement tests. The initial condition is given as u(x, 0) = sin(πx) with a periodical boundary condition on x ∈ [0, 2π). The exact solution is u(x, t) = sin(x − t). We refine the grid and record numerical errors after one period (t = 2).
Two typical norms for errors are adopted here, i.e., E 1 = 1 N ∑ N j=1 u j − u true j , and The errors and the convergence rates are presented in Table 1. We can see that the spatial accuracy order is 2r − 2 for the AltPoly-r scheme, which matches our theoretical analysis in Section 3. Additionally, the order of accuracy is also maintained well for non-uniform meshes. Next, we use a complicated initial condition from [30] to evaluate the ability of the proposed method in handling both continuous and discontinuous profiles. This initial condition contains four piece-wise functions: where the involved functions are with the constants set as α = 10, a = 0.5, z = −0.7, δ = 0.005, and β = log 2/(36δ 2 ). As for the initial profile, we suggest using the routine used in the discontinuous Galerkin method. A smaller time step can help reduce the oscillations around discontinuities. We take ∆t = 0.3∆x, 0.1∆x, 0.075∆x, and 0.01∆x, respectively, for AltPoly-3, AltPoly-4, AltPoly-5, and AltPoly-6.
The computational domain is divided into 200 elements. We first compute this problem without adding artificial viscosity, for which higher-order methods can produce better results in both continuous and discontinuous regions. The results at t = 2 are illustrated in Figures 3 and 4 separately. Specifically, AltPoly-3 and AltPoly-4 generate significant oscillations around x = −0.4 and x = −0.2. These oscillations are reasonable as our method is based on the smoothness assumption and the polynomial does have its limitations on presenting discontinuities. However, it is interesting that spurious oscillations are significantly reduced and the transitions around strong discontinuities are quite sharp for AltPoly-5 and AltPoly-6. After adding a proper amount of artificial viscosity, the proposed method is able to smooth the artificial oscillations for AltPoly method of various orders while not losing high accuracy in the smooth regions.

Nonlinear Burgers Equation
This subsection considers the 1D inviscid Burgers equation The initial condition is u(x, 0) = 0.5 + sin(x) for x ∈ [0, 2π) with the periodical boundary condition.
Although the initial profile is smooth, shock waves exist when t is large enough since the flux function is nonlinear. To evaluate the order of accuracy, we calculate the numerical solution before the shock wave occurs and compute the errors in respect of grid refinement. Table 2 lists the errors at t = 0.5 when the solution is still smooth. The numerical convergence rate results are consistent with our theoretical analysis.
The numerical results for the Burgers equation at t = 1.5 are shown in Figure 5 when a shock has already emerged. To handle the shock discontinuity, we add a proper amount of artificial viscosity according to Section 2.4. We can see that our schemes sup-press the oscillations well around the discontinuity while still maintaining a sharp and accurate profile.

Euler Equations
We then use our scheme to solve the one-dimensional Euler system: where ρ refers to the density, u the velocity, m = ρu the momentum, e the total energy, and p the pressure: with γ = 1.4 for a perfect dry gas. Just like in the scalar case, we introduce point-valued and cell-averaged moments for each variable. The cell-averaged moments are updated by the flux-form Equation (56). In Step 2 of AltPoly, solution values on auxiliary points are updated using the following non-conservative form equations: where Neither characteristic decomposition nor Riemann solver are needed here. We first check the order of accuracy for this nonlinear system by simulating advection of the density perturbation. The initial conditions are ρ(x 0 ) = 1 + 0.2 sin x, v(x, 0) = 1, and p(x, 0) = 1 for x ∈ [0, 2π) with the periodic boundary condition. The exact solution is a traveling wave: ρ(x, t) = 1 + 0.2 sin(x − t), v(x, t) = 1 and p(x, t) = 1. The error results of the density are listed in Table 3 and a satisfactory convergence rate is obtained. We then consider the classical shock tube problem, which is described by Euler equations with the following Riemann-type initial distribution: Two initial conditions are tested here. The first one is known as Sod's problem and the initial condition is defined as The second one is Lax's problem, which starts from the following initial condition: We compute these shock tube problems with 100 cells using the AltPoly schemes equipped with the artificial viscosity. Generally, a certain degree of spurious oscillation still exists but is in a tolerant range. As shown in Figure 6 and 7, all AltPoly schemes can depict shocks sharply using only one or two cells. The higher-order schemes (AltPoly-5 and AltPoly-6) capture the contact discontinuities with only one or two cells, while more cells are needed for the lower-order AltPoly schemes.
To show the accuracy of our smoothness indicator (23), Figure 8 illustrates the cells where the artificial viscosity takes effect on the x-t plane. Generally, the number of cells recognized as discontinuous for lower-order schemes is significantly smaller than that for higher-order schemes. This might be because lower-order schemes have a larger dissipation and the smoothness indicator is less accurate than that of high-order schemes. It is remarkable that the artificial viscosity vanishes in most areas and is in effect only around the shock wave. The contact discontinuity does not trigger the smoothness indicator as it can be smoothed by the artificial viscosity in the first few time steps.

Conclusions
We have presented a high-order conservative solver termed AltPoly for 1D hyperbolic conservation laws. This method can be categorized as a multi-moment scheme since it adopts both point-values at the middle of the cell and cell-averaged values as the model variables. AltPoly updates the model variables via alternately implementing two polynomial reconstructions. The point-based variables are used to reconstruct a high-order accurate approximation of the solution via Hermite interpolation. The cell-averaged variables updated by the flux-form equations serve as the constraint to guarantee numerical conservation. The solution within a cell is corrected by polynomial reconstruction via solving a constrained least-squares problem. In other words, our method directly reconstructs the solution to maintain the conservation, which is different from existing schemes based on the local flux reconstruction. Fourier analysis shows the accuracy and the admissible range of the CFL number for the linear scalar equation. Further, our scheme can be extended to other kinds of PDEs such as parabolic equations, advection-diffusion equations, etc., which will be our future work.
where A † denotes the generalized inverse or Moore-Penrose inverse [27] of A. From b we can obtain the updated point-based model variables.
In practice, we find that A † has a large condition number when r ≥ 5 and we calculate it via symbolic computation using mathematical software to reduce the error. Throughout this paper, we used the symmetrical solution points − 1, −1 + σ 1 , −1 + σ 2 , . . . , −1 + σ r−2 , 1 − σ r−2 , 1 − σ r−1 , . . . , 1 − σ 1 , 1 , for AltPoly-r with the choices of {σ i } shown in Table A1. Through numerical experiments, we found that smaller {σ i } can generally induce a more stable numerical scheme, which allows for a larger CFL number. However, when {σ i } approaches 0, the condition number of A tends to infinity, which will cause a loss of accuracy. Based on numerical trials, the choices of {σ i } are determined as a compromise between stability and accuracy.
We can see that the magnitudes of all entries in S 0 [1 : 2; 1 : 2] are smaller than 1 when w is small enough. In addition, from (39) we know that S 0 e 3 = e 3 . (A11) Therefore the principal eigenvalue of S 0 is 1 as long as w is small enough. Now we try to bound S n e 1 and S n e 2 . We first assume that S i ≤ M for any i ≥ 1, which naturally holds as long as the scheme is stable. Then we have Since the magnitude of c = ∆t/∆x is not greater than O(1) and S i ≤ M, then A 1 can be estimated by As for A 2 , the following bound can be deduced based on (A12),