ENTROPY STABLE ESSENTIALLY NONOSCILLATORY METHODS BASED ON RBF RECONSTRUCTION

To solve hyperbolic conservation laws we propose to use high-order essentially nonoscillatory methods based on radial basis functions. We introduce an entropy stable arbitrary high-order finite difference method (RBF-TeCNOp) and an entropy stable second order finite volume method (RBFEFV2) for one-dimensional problems. Thus, we show that methods based on radial basis functions are as powerful as methods based on polynomial reconstruction. The main contribution is the construction of an algorithm and a smoothness indicator that ensures an interpolation function which fulfills the sign-property on general one dimensional grids. Mathematics Subject Classification. 35L65, 65M12, 65M06, 65M08, 65D05. Received August 13, 2018. Accepted February 11, 2019.

(1.1) in the sense that ∇ v upvq is symmetric positive definite and ∇ v f pupvqq is symmetric. This can be seen by introducing the entropy potential ψpvq v T ¤ f pupvqq ¡ qpupvqq to recover ∇ v f pupvqq ∇ vv ψ [28]. Note that entropy solutions satisfy ηpuq t qpuq x 0, (1.3) in smooth regions, but satisfy (1.2) at discontinuities as entropy must dissipate.

Finite difference and finite volume methods
One goal of numerical methods is to express the approximative behaviour of the physical correct solution.
Let us assume a one-dimensional grid tx i u iZ R, partitioned into cells C i rx i¡1{2 , x i 1{2 s of size ∆x. The finite difference approach is based on approximating spatial derivatives at x i as ff fx pu i q F i 1{2 ¡ F i¡1{2 ∆x Op∆x p q, (1.4) with p ¡ 0. This results in the semi-discrete scheme du i dt 1 ∆x pF i 1{2 ¡ F i¡1{2 q 0, (1.5) where the numerical flux terms F i 1{2 depend on point values tu i¡k , . . . , u i p¡k u with 0 ¤ k ¤ p ¡ 1.
On the other hand, finite volume methods work with mean valuesū i for the cells C i . By integrating (1.1) over the cells and dividing it by its size |C i | we recover (1.5) with the difference that f pupx i 1{2 qq F i 1{2 Op∆x p q. (1.6) There exist multiple high-order accurate methods to solve conservation laws, for example the MUSCL scheme introduced in [39], the ENO scheme [18] or the WENO scheme [37]. In all cases, we can apply an arbitrary time discretization technique to recover a fully discrete scheme, e.g. an SSPRK method.
In [11], Fjordholm et al. proposed an entropy stable TeCNO scheme based on polynomial reconstruction. We follow the spirit of this work and introduce a scheme based on radial basis functions (RBF). Note that for taking advantage of the RBFs on general grids in higher dimensions, we would need a MUSCL type scheme.
In Section 2 we introduce the framework of entropy stable schemes based on an entropy conservative flux and a diffusion operator [11,38]. We describe the basics of RBFs in Section 3. Furthermore, we give an explicit representation of the interpolation function for infinitely smooth RBFs in one dimension. Sections 4 and 5 contain the main contribution. In Section 4 we introduce a smoothness indicator for RBFs which is based on the generalized divided difference method and in Section 5 we prove that it fulfills a certain stability property (the sign-property) for general grids in one dimension. Section 6 combines these results to construct an arbitrarily high-order RBF based entropy stable finite difference method and a second order entropy stable finite volume one that is generalizable to general grids in higher dimensions. In the last section we demonstrate the robustness of the numerical scheme with a variety of one dimensional examples.

Entropy stable methods
The goal is to construct methods that fulfill a discrete version of (1.2), referred to as entropy stable [19,23]. As a first step, we introduce entropy conservative methods that fulfill (1.3) at the discrete level. Next, we add a dissipation term to control oscillations at discontinuities to recover an entropy stable method.

Entropy conservative methods
A finite difference method is entropy conservative if it satisfies d dt for a consistent numerical entropy flux Q i 1{2 . To construct entropy conservative methods we use Tadmor's entropy conservation condition [38] vvw T i 1{2 F i 1{2 vψw i 1{2 .

(2.2)
This condition describes a system of equations, but its solvability is not clear. For scalar conservation laws there exists a unique solution as can be summarized in the following theorem.
Theorem 2.1 (Entropy conservative schemes for scalar equations [38]). For a given entropy pair pη, qq the numerical flux

3)
defines an entropy conservative method for scalar equations with the entropy variable v and the conserved variable u. Furthermore, it is second-order accurate in smooth regions of u.
Given a numerical second order two-point flux the work of Lefloch et al. [24] combines these linearly to construct a 2pth order accurate flux on a uniform grid.
Theorem 2.2 (High-order entropy conservative fluxes [24]). Let p N and assume that α 1,p , . . . , α p,p solve the p linear equations  2 pu i¡j l , u i l q, (2.5) is consistent, 2pth order accurate and entropy conservative provided the second order two-point conservative flux F 2 fulfills (2.2).

Entropy conservative methods for shallow water equations
The shallow water equations describe a flow under the assumptions that the horizontal length scales are much larger than the vertical ones. In one space dimension the system of equations depends on the mass flow m and the fluid height h with the gravitational constant g [25]. To apply Theorem 2.2 we need to construct a second order entropy conservative scheme by solving (2.2). One choice of an entropy pair for the one-dimensional shallow water equation is An alternative second order entropy conservative flux is

Entropy conservative methods for Euler equations
The Euler equations can be recovered from the Navier-Stokes equations by neglecting the viscosity. They consist of the continuity equation, momentum equation and the conservation law for the total energy. In one dimension they are with p RρT pγ ¡1qpE ¡ 1 2 m 2 ρ q for an ideal gas with the ratio of specific heat γ [20]. For the Euler equations the thermodynamical entropy s logppq ¡ γ logpρq is different from the entropy function and the entropy flux.
One possible pair is proposed in [11] η ¡ρs γ ¡ 1 , q ¡ms γ ¡ 1 ¤ (2.11) Chandrashekar [5] proposed the kinetic energy preserving and entropy conservative (KEPEC) flux, based on the entropy variables and the potential  The KEPEC flux makes use of the logarithmic averagesρ andβ with β ρ 2p and can be written as

Entropy stable methods
Entropy conservative methods yield good results in smooth regions, but it is well-known that spurious oscillations appear close to discontinuities. Introducing artificial dissipation, depending on the size of the jump in the interface, controls these oscillations. Based on an entropy conservative schemeF j 1{2 of second order and a symmetric positive definite matrix D i 1{2 , Tadmor constructed the entropy stable numerical flux function [38] (2.14) Combining high-order conservative fluxes with dissipation terms introduces the constraint that D i 1{2 vvw i 1{2 Op∆x p q to maintain accuracy for smooth solutions.
For each cell C i we define a stencil of cells S i on which we construct an interpolation function s i pxq of order p and replace the jump vvw i 1{2 by the jump in the reconstruction xxvyy i 1{2 s i 1 px i 1{2 q ¡ s i px i 1{2 q. Thus, the method has the form with the additional condition (2.17) Then the scheme with the flux (2.15) is entropy stable.
By introducing the scaled entropy variables with the reconstructed entropy variables vï s i px i¨1{2 q, (2.17) becomes Since B i 1{2 is diagonal and semi-positive definite, this can be reformulated componentwise as for each component l. This structural property of the reconstruction is called the sign-property.

Entropy stable finite volume methods
The setting of the one dimensional finite volume method differs only slightly from the finite difference scheme, i.e. we consider cell-average valuesū i rather than point values and the definition of higher order methods changes to Nevertheless, given a 2-point second order finite difference flux F , it is also a second order accurate finite volume flux, which follows from the midpoint rule.
Since the definition of entropy conservative schemes does not change for finite volume methods, we can conclude that a second order finite difference flux that fulfills (2.2) is also a second order entropy conservative finite volume method. This can be summarized as follows.
Theorem 2.4. Every second order finite difference scheme that fulfills Tadmor's entropy conservation condition (2.2) in one space dimension is also a second order entropy conservative finite volume method.
The construction of entropy stable schemes from entropy conservative schemes proceeds as for the finite difference case, the only difference being that the interpolation is based on cell-averages instead of point values.
Thus, Lemma 2.3 holds also for finite volume methods and we recover an entropy stable finite volume method of the form (2.22) Note that the extension to higher order as in Theorem 2.2 does not work in the finite volume case and we are not aware of any method doing it.

Radial basis functions
Radial basis functions (RBF) are successfully used for scattered data interpolation. Due to their mesh-free property, they are more flexible in terms of the geometric structure of the data points. Furthermore, their application to high-dimensional problems is immediate. Following the seminal work by Duchon [8] and Micchelli [27], RBFs are successfully used in different domains.

Basic interpolation
The goal is the interpolation of a data vector f | X pfpx 1 q, . . . , f px n qq T R n , defined at a scattered set of data points X px 1 , . . . , x n q T with x i R d for some function f : R d Ñ R. The basic idea is to use one univariate continuous function φ, the radial basis function, composed with the Euclidean norm centered at the data points as the interpolation basis B tφpε x ¡ x 1 q, . . . , φpε x ¡ x n qu, (3.1) parametrized by the shape parameter ε. To simplify the notation we use The standard radial basis function approximation is written as with a polynomial p Π m¡1 pR d q, m N, the interpolation condition spx i q f px i q, (3.4) and the additional constraints ņ i1 a i qpx i q 0, for all q Π m¡1 pR d q, (3.5) with the coefficients a i R for all i 1, . . . , n. Conditions (3.4) and (3.5) can be expressed in the system of The choice of the radial basis function φ is restricted to insure the solvability of (3.6).  for all p Π m¡1 pR d q, the quadratic form ņ j,k1 is positive (non-negative).
is positive for any n pairwise different points x 1 , . . . , x n R d and c pc 1 , . . . , c n q T R n zt0u.
Note that for a positive definite function φ, the matrix A is positive definite and there exists an unique solution a to (3.6).
Some examples of positive definite functions are the inverse multiquadratics and the Gaussians (Tab. 1). Other RBFs fulfill a slightly weaker condition and are conditionally positive definite of order k, e.g. multiquadratics and polyharmonic splines.

Interpolation of cell-averages
For the finite volume method we do not consider the pointwise interpolation, but cell-averages. Let us assume a given grid of cells C 1 , . . . , C n with its average valuesū 1 , . . . ,ū n for n N. Following [1,2] we consider spxq ņ i1 a i λ ξ Ci φpx ¡ ξq ppxq, p Π m¡1 pR d q, (3.10) with the average operator of f over the cell C, λ ξ C f , such that λ Cj s ū j , for all j 1, . . . , n, To show solvability of system (3.11) it suffices to assume that φ is conditionally positive definite in a pointwise sense.
Aboiyar et al. prove in [1] that (3.11) has a unique solution if the set tλ Ci u n i1 is Π m¡1 pR d q-unisolvent. The proof closely follows the one for the pointwise evaluation in [40] plus an estimate for the positive definiteness, based on a pointwise result in [31].
Despite the simple extension of RBF-interpolation in multiple dimensions, there is a major drawback, often referred to as the Uncertainty Principle [34]. It refers to the trade-off between the well-known properties that flat infinitely smooth RBFs (ε Ñ 0) have an increasing approximation power but a decreasing numerical stability due to ill-conditioning of the interpolation matrix [7,22,36]. To overcome the ill-conditioning there are multiple propositions for choosing an "optimal" shape parameter [9,33]. Note that a continuous scaling ε αn ¡1{d causes error stagnation [4]. However, there are multiple approaches which overcome this problem: the RBF-CP [16], the RBF-QR [15], and the RBF-GA [14]. We decided to use the more flexible vector-valued rational approximation method (RBF-RA), based on the RBF-CP algorithm [42].

Explicit formula of the RBF interpolation
Let us consider the pointwise RBF interpolation problem (3.6) with the interpolation function (3.3). Furthermore, let x 1 , . . . , x n be the grid points such that x i x i 1 and n N and y 1 , . . . , y n R its values. We are looking for an RBF interpolation function (3.12) where L j for j 1, . . . , m are the Lagrange polynomials such that L j px i q δ ij and φ a conditional positive definite RBF of order m. By assuming further that m n ¡ 1, it holds y i L i pxq, Proof. From the representation of the polynomial part using Lagrange polynomials we recover a j ¡a n L j px n q, for j 1, . . . , n ¡ 1. (3.14) This yields the interpolation function with α a n , which solves the reduced interpolation problem (3.16) By the properties of the Lagrange polynomials we recover the explicit form of α and b j α y n ¡°n ¡1 (3.17b) Remark 3.5. We can express dϕ in terms of projections dϕpxq : Ψpx, x i q pId ¡ P x qpId ¡ P y qrφpx ¡ yqqs| yxi , where the operators P z is the projection of the variable z on the polynomial space of dimension n ¡1. Schaback [35] shows that Ψ is positive definite on R d ztx 1 , . . . , x n¡1 u. Thus, it is closely related to reproducing kernels and its native spaces, introduced in [35].
Remark 3.6. Note that this representation is independent of permutations of the indices. In general we can chooseỹ n y j andỹ i ty l | l $ ju.

Smoothness indicator for RBF interpolation functions
In essentially nonoscillatory (ENO)-and weighted ENO (WENO)-type methods a key component is to measure the smoothness of the interpolation function. In the polynomial ENO scheme, the highest degree divided difference plays an important role for identifying the least oscillating interpolation of a certain degree. To extend this to RBF-based interpolation we need something similar. However, the divided differences, used in the standard Newton's interpolation formula, are valid only for polynomials.

Generalized divided differences
For non-polynomial basis functions Mühlbach [29] introduces generalized divided differences, which coincide in the monomial case with the standard one. The result is based on functions f 1 , . . . , f n that form a Chebyshev system, i.e., they satisfy . . .  Theorem 4.1 (Generalized Newton's interpolation formula [30]). Let f 1 , . . . , f n form a complete Chebyshev system. Then for any f : R Ñ R and any subset G n tx 1 , . . . x n u R of cardinality n it holds where with the interpolation error in the kth step 5) and the recursively defined coefficients Based on this we can express the generalized divided differences for the basis t1, x, . . . , x N ¡2 , ϕu for N N and ϕ from Lemma 3.4 to quantify the oscillations of the interpolation function. To distinguish between the Lagrange polynomials of different degree we write L i,d j for the Lagrange polynomial of degree d such that Let the basis be given by t1, x, . . . , x N ¡2 , ϕu for N N and ϕ as defined in Lemma 3.4. We recover the generalized divided differences of the form By comparing this results with the RBF interpolation in Lemma 3.4, we observe that the last divided difference can be written as This suggests that α may be a good choice as the smoothness indicator based on the success of the classic ENO scheme.

Relation to reproducing Kernel Hilbert spaces and its norm
As mentioned above there is a close relation to native spaces of conditionally positive definite functions (see Schaback [35]). Indeed, the RBF-based basis function dϕ can be expressed in terms of the modified kernel function Ψpx, yq pId ¡ P x qpId ¡ P y qrφpx ¡ yqqs.
By analysing the norm of the interpolation function, based on the inner product of the native space, we have Lemma 4.3. Let s be a RBF-interpolation function given by (3.13). Then, it has the norm In particular, we have This lemma proposes a scaling of dϕpx N q 1{2 of our smoothness indicator.
Proof. The inner product of the native space is for f °M j1 λ j φpx, x j q and g °N k1 µ k φpx, y k q [35].
We have ps ¡ Psqpxq β dϕpxq dϕpx N q and Finally, we insert the definition of dϕ to recover is a basis of the interpolation space. In particular, we have equivalence of the norms || ¤ || φ and || ¤ || B , where Proof. From the interpolation (3.13) we recover that B is a basis of the interpolation space.

Smoothness indicator and stencil choice
Harten et al. proposed the Essentially Nonoscillatory method to control spurious oscillations at discontinuities [18]. Its principle is based on the evaluation of multiple stencils for each cell C i in which we need to reconstruct the solution. Finally, one chooses the least oscillatory reconstruction to define s i . Fjordholm et al. [12] showed the sign-property for the polynomial reconstruction method with the recursive algorithm introduced in [18] which utilizes the last divided difference as a local smoothness indicator. A sign preserving WENO reconstruction method was proposed by Fjordholm et al. [13]. In the RBF reconstruction, the highest derivative is similar to the RBF-part of the reconstruction in Lemma 4.3 and Theorem 4.2. As we shall show, the recursive algorithm from the polynomial case, combined with the smoothness indicator is sign-stable for small enough grid sizes. Numerical experiments confirm this to be true for general grids. In the next section we prove this for the second and third degree reconstructions on general grids. Note that Corollary 4.4 ensures the definition of ISpsq.
Remark 4.6. The restriction that the sign-property holds only on grids with small grid size is not a limitation. For infinitely smooth RBFs we can choose a small shape parameter to decrease the computational grid size.
Remark 4.7. The smoothness indicator (4.19) requires an impractical and computationally expensive evaluation. However, from Lemma 4.5 we recover (4.20) Thus, we have equivalence of the smoothness indicator IS with r ISpsq : β To choose the least oscillatory stencil S i for the ith cell for the RBF-reconstruction we follow Algorithm 1 which is based on the one from Harten et al. [18]. We use the notation spj, . . . , j kq that corresponds to the reconstruction on the cells C j , . . . , C j k with the interpolation points x j , . . . , x j k and its values y j , . . . , y j k .

Sign-property for 2nd and 3rd degree reconstruction
Based on the results from the previous sections we show the sign-property of the RBF interpolation for the second and third degree reconstruction, i.e. N 2, 3. This means that we deal with stencils S i of size N which represent the interpolation points for the reconstruction on cell C i . Let us name them (5.1b) where r N ¡1 ¤ 1 s N ¡1 and C j is the jth cell with its mid-point x j on which we apply the interpolation. Further, we define d N ¡1 : 1 s N ¡1 ¡ r N ¡1 ¥ 0 as the shift between the stencils. The stencils are chosen by Algorithm 1 and there are no constraints on the stencils.

Notation
For simplicity, we introduce some general notation. We assume the stencil length to be N and we name terms by the highest appearing index j that exists in the underlying stencil C j¡N 1 , . . . , C j . We also define L i j pxq ¡ Lagrange polynomial of degree N ¡ 1 such that (5.5)

Representation of the reconstructed jumps
The idea of the proof is to give a simple representation of the reconstructed jumps and show that each term has the same sign as the jump in its neighboring cells. Let us assume that we have given the stencils S i and S i 1 for the cells i and i 1 from Algorithm 1.
Theorem 5.1 (Generalized representation). The second and third degree reconstructed jump can be written in the following form with the constants and an error term The proof relies on multiple Lemmas which we now develop.
Proof. The case N 2 is immediate since the Lagrange polynomials are constant. Thus, (5.10) is ¡y j¡1 β j ¡ y j , (5.12) For N 3, we write the left hand side of (5.10) and subtract (5.14) Lemma 5.3. The reconstructed jump jR i 1{2 for the second and third degree reconstruction method can be expressed as Proof. From Lemma 3.4 and the stencils selected from (5.1) we rewrite the N th degree reconstructed jump jR i 1{2 between cell i and i 1 as The polynomial part of the reconstructed jump is By inserting γ i α i dϕ i px i q 1{2 we recover the result.
with P k j 1 as the kth degree polynomial approximation with respect to the interpolation points x j , . . . , x j 1¡k .
Proof. In the case N 2, we have A j 1 and L j 1 and recover In the case N 3, we have which can be simplified as Next, we express the last term and insert this into (5.21) where we use that Now, we are ready to prove Theorem 5.1.
Proof. (Thm. 5.1) The goal is to show the equivalence with the representation in Lemma 5.3. Therefore, we insert (5.9) into (5.7) to recover Finally, we insert the definition of C j to obtain Remark 5.5. Let us define ε j p∆xq : 1 Thus, the error εp∆xq can be written as From Lemma 5.4 we can express ε j p∆xq by (5.27)

Sign-property for small grid size
In this section we analyse the reconstructed jumps for infinitely smooth RBFs for small grid size ∆x Ñ 0.
From Theorem 5.1 we have a simple expression for the reconstructed jump. We first show that the error εp∆xq goes to zero, as the grid size goes to zero. Then, we show that each term of the remaining equation has the sign of the jump y i 1 ¡ y i . Remark 5.6. The notation Op∆x p q for ∆x Ñ 0 should be interpreted in the way that given a gridx 0 x 1 ¤ ¤ ¤ x m we analyse the terms for the grid x 0 x 1 ¤ ¤ ¤ x m with x j x j ∆x and we use Op∆x p q ¤ C∆x p , for ∆x Ñ 0. (5.28) Remark 5.7. When calculating the errors ε j we recall that dϕ j pxq ϕ j pxq ¡ P N ¡1 j ϕ j pxq. Theorem 5.8. Let φ be an infinitely smooth RBF of first or second order. Then, we have that εp∆xq Op∆x 2 q for ∆x Ñ 0 for N 2, 3 and Proof. We start by analysing the different parts in the error term ε k p∆xq. Note that as φ is a conditionally positive definite RBF φpxq hpx 2 q. Let us start with the case N 2 and a first order RBF: We further have that which allows us to conclude ε k p∆xq Op∆x 2 q.
To find a bound for εp∆xq we further need δ k {δ k 1 Op1q and β k d Op1q. The latter is clear since the reconstructed function is bounded and L i j px i q is Op1q. Further, δ k {δ k 1 Op1q results from (5.34). Using (5.36) in (5.26) we conclude that εp∆xq Op∆x 2 q.
Next, we consider the more complicated case with N 3 and a second order RBF φ. Hence, we need to analyse the following two terms: As before, we apply the Taylor expansion We write Let us calculate the coefficients a 1 and a 2 a 1 a k¡1 From standard algebra we recover that a k¡1 1 a k 1 0. The fourth order term is We repeat this for the coefficients b 1 and b 2 The results can be summarized as From this we recover  Since the error term εp∆xq vanishes, the remaining step is to prove that each term of (5.30) has the same sign as the jump. Theorem 5.9 (Sign-property of second and third degree RBF-reconstruction). Let us assume that the stencil S i and S i 1 are chosen with the Algorithm 1. Then, for infinitely smooth RBFs of first or second order it holds that sgnpC j pγ i r N ¡1 N j ¡ γ i r N ¡1 N ¡1 j qq sgnpy i 1 ¡ y i q, (5.51) for all j 0, . . . , d N ¡1 ¡ 1 and for ∆x Ñ 0.
Proof. The proof is based on a study of all possible choices of stencils, that may result from Algorithm 1: For each case we look at any inequality due to Algorithm 1 to recover the particular stencil configuration, and show for each case that (5.51) is fulfilled.
Note that jR i 1{2 0, if S i S i 1 . Hence, we do not include such cases in the analysis. Further, we use the notation A B if A B Op∆xq, for ∆x Ñ 0.
and from (5.32) it follows that  The jump can be represented as As before it holds that sgn C 0 pγ i 1 ¡ γ i q¨ sgnpy i 1 ¡ y i q and Thus, we get for the second term where we used (5.55) and (5.56).
Case 3. In the last case of the second degree reconstruction we have the stencils S i tC i , C i 1 u, S i 1 tC i 1 , C i 2 u, equivalent to the conditions The representation of the jump is As in the first case we recover from (5.32), that This completes the proof of the sign-property for the second degree reconstruction with infinitely smooth RBFs of first order for small enough grids.
The proof for N 3 can be found in Appendix A.

Entropy stable RBF-based methods
In one space dimension there is no need to deviate from the polynomial reconstruction. For unstructured grids in multiple dimensions the problem is the construction of an interpolation function. There exist a lot cell or point configurations such that the reconstruction problem is not well-defined. This issue can be relaxed by solving an overdetermined system of equations, but then we lose the exact interpolation property. The RBFinterpolation can circumvent this problem since we do not need a unisolvent set of cells or points, but only a unisolvent subset of lower order. Thus, by adding some extra cells we significantly reduce the possibility that an unsolvable configuration occurs.

RBF-TeCNOp method
Based on the theory of entropy stable schemes and the work of Fjordholm et al. [11] we introduce the RBF-TeCNOp scheme. By using Algorithm 1 with (4.21) for calculating the least oscillatory stencil, Theorem 5.9 shows that the sign-property holds for 2nd and 3rd degree reconstruction in the limit of ∆x Ñ 0. We conjecture that this result generalizes to higher order reconstructions. Thus, by combining the framework proposed in [11] with the RBF reconstruction using multiquadratics we recover an entropy stable essentially nonoscillatory RBF-based finite difference method of arbitrary high order. Furthermore, we use the RBF-RA algorithm to circumvent ill-conditioning in the reconstruction step [42].
In more detail, for constructing a pth order RBF-TeCNOp method of the form (2.15) we use an entropy conservative flux of order 2k with k rp{2s (see Thm. 2.2) and an ENO based RBF reconstruction (Algorithm 1) on the scaled entropy variables of order p with multiquadratics of order p ¡ 1.
Based on the Roe diffusion operator R|Λ|R ¡1 vuw, (6.1) with the eigenvector matrix R and the diagonal matrix of the eigenvalues Λ, evaluated at the Roe average, we are choosing R and Λ in the same way. By Merriam [26] there is a scaling of the eigenvectors such that Thus, we get the relation that has a similar structure to that of a diffusion operator (2.16). The numerical diffusion term can be written as with the scaled entropy variables (2.18). Furthermore, we choose Λ i 1{2 diagpλ 1 pu i 1{2 q, . . . , λ N pu i 1{2 qq and u i 1{2 ui ui 1 2 with the eigenvalues Λpuq of the Jacobian ∇ u f . It is important to note that the ill-conditioning of the interpolation matrix does not only affect the evaluation of the reconstruction; it also affects the calculation of the smoothness indicator which is based on the sum of the squares of the coefficients of the RBF-part of the interpolation.
From the theory we expect that the error of the interpolation with infinitely smooth RBFs decreases for smaller shape parameters. However, computations suggest that the choice of the stencil does not depend on the shape parameter. Thus, we calculate the stencil with respect to a stable shape parameter.

RBF-finite volume method
The combination of the RBF interpolation with finite volume methods works analogeously to the RBF-TeCNOp Method. Aboiyar et al. [1] combine in their work a high-order WENO approach with a polyharmonic spline reconstruction and Bigoni et al. [3] apply a high-order WENO approach to multiquadratics.
We construct an entropy stable finite volume method of second order that is essentially nonoscillatory by combining (2.22) with a second order accurate RBF interpolation that acts on the scaled entropy variables. Therefore, we are using multiquadratics with the smoothness indicator (4.21), combined with Algorithm 1 and the vector valued rational approximation to ensure a stable evaluation of the interpolation function.
We conjecture the sign-property for the RBF reconstruction on mean values that is based on Algorithm 1 which is fulfilled in the pointwise case for second and third degree reconstruction in the limit ∆x Ñ 0 (Thm. 5.9).
Under this assumption we recover a second order entropy stable finite volume (RBF-EFV2) method. This method can be generalized to general grids in higher dimensions [32].

Numerical results
In this chapter, we are evaluating the second order entropy stable finite volume (EFV2) and the TeCNOp methods with RBF reconstruction for one-dimensional problems and compare it with its original version. Note that in one dimension we do not expect to do better than the classical methods, but also not worse.
For the polynomial reconstruction we use the original algorithm from [18] to select the stencil and in the RBF case we use Algorithm 1. The EFV2 and TeCNOp methods are based on the diffusion term (6.3).
The parameters for the vector-valued rational approximation are chosen as described in [42]. Further, we choose the shape parameter ε 0.1 for all examples.

Linear advection equation
We consider the linear advection equation and obtain the second order entropy conservative flux 3) Notes. We use periodic boundary conditions and u0pxq sinp2πxq, shape parameter ε 0.1, CFL 0.5. to construct a high-order accurate scheme. We use a 5th order SSPRK method for the time discretization in the TeCNOp method [17]. For the EFV2 method we use the second order entropy conservative flux plus a third order SSPRK method in time.
The convergence results for the smooth initial conditions are shown in Table 2. The L 1 -errors are the same for the different reconstruction methods for grids of size smaller than 1{32 and their convergence rates are as expected.
In Table 3 we present a comparison of the runtime for the 5th order TeCNO scheme with RBF reconstruction using the RBF-RA algorithm, RBF reconstruction evaluated at a stable shape parameter. The RBF method based on the RBF-RA algorithm solves multiple times the same system of equations with different shape parameters to approximate the final one, explaining the difference between the first two columns. The difference between the RBF and the polynomial reconstruction comes from the fact that the RBF algorithm solves in each recursive step of Algorithm 1 a system of equations and the polynomial case calculates just the next divided difference. which is used to construct an high-order scheme. For the time discretization we use a 5th order SSPRK method [17]. Furthermore, we choose the domain r0, 1s and the initial conditions u 0 pxq sinp2πxq.

Burgers' equation
A detailed analysis of the convergence is shown in Table 2. The convergence rate is as expected and the errors of the two different methods (polynomial reconstruction and RBF reconstruction) coincide. At time t 0.3 a discontinuity emerges at x 0.5. This can be resolved accurately with vanishing oscillations (Fig. 1).
Furthermore, we observe that the difference between the reconstruction methods vanishes in the smooth part while at the shock it stays small.

Shallow water equations
For the shallow water equations (2.6) we consider the dambreak problem with the initial conditions ph 0 , m 0 q 5 p1.5, 0q if |x| ¤ 0.2 p1, 0q if |x| ¡ 0.2 , on the domain r¡1, 1s and periodic boundary conditions. We use a second order entropy stable flux (2.9) to construct a high-order flux and a third order SSPRK method for the time integration.
Fjordholm et al. [11] showed that the standard TeCNO scheme behaves similar to the ENO-MUSCL scheme. The same holds for the RBF-TeCNOp scheme and the RBF-EFV2 scheme as seen in Figure 2. The difference between the RBF methods and the polynomial scheme is around Op10 ¡6 q in the region where the discontinuity passed and much smaller in smooth regions.

Euler equations
The one-dimensional Euler equations (2.10) are a system of size three. As for the shallow water equations we use a third order SSPRK method and as a second order entropy conservative flux we use the KEPEC-flux (2.13). Further, we choose γ 1.4 which simulates a diatomic gas such as air.

Sod's shock tube problem
Sod's shock tube problem is a Riemann problem where two gases with different densities collide. A rarefaction wave emerges, followed by a contact and a shock discontinuity. The initial conditions are pρ 0 , m 0 , p 0 q 5 p1, 0, 1q where m uρ. The results at time t 2 of the RBF-TeCNOp and RBF-EFV2 methods are shown in Figure 3, clearly representing the rarefaction wave, the contact, and the shock discontinuity. Comparing the solutions obtained with polynomial reconstruction or with RBF reconstruction, we see in Figure 4 that their difference is decreasing with the refinement of the grid.

Lax shock tube problem
The Lax shock tube problem is another Riemann problem with the initial conditions pρ 0 , m 0 , p 0 q 5 p0.445, 0.698, 3.528q if x 0 p0.5, 0, 0.571q if x ¥ 0 , where m uρ. The RBF-TeCNOp methods of order three to five represent the big shock in the density sharply with just N 100 points, see Figure 5. The second order RBF-EFV2 method does not perform well for this case.   where m uρ. The RBF-TeCNOp methods of order larger than three recover the high frequency oscillations well. The RBF-EFV2 method fits the low order oscillations and the shock, but not the high frequency one due to excessive dissipation (Fig. 6).

Low density problem
The low density problem on r0, 1s is a Riemann problem that tests the ability of preserving positive density and pressure. The initial condition are with Neumann boundary conditions. In Figure 7 we observe the increasing accuracy of the RBF-TeCNOp method with increasing order p. Note that by choosing the wrong smoothness indicator, negative pressures will occur.

Two interacting blast waves
A more complex one dimensional example is the two interacting blast waves, introduced by Woodward et al. [41]. It is based on two blast waves that interact and introduce low pressures and densities. Its initial condition is pρ 0 , m 0 , p 0 q 6 9 8 9 7 p1, 0, 1000q if x 0.1 p1, 0, 0.01q if 0.1 ¤ x 0.9 p1, 0, 100q if x ¥ 0.9 . (7.12) Compared to the low density problem we add here an additional challenge caused by the collision of the two shocks. Note that the standard version of the TeCNO method always produces negative pressures or densities for the RBF and polynomial reconstruction.
Thus, we introduce a more complicated symmetric positive definite dissipation operator (2.16), similar to the one introduced by Derigs et al. [6]. The goal is to mimic the more dissipative Rusanov-type diffusion operator such that αvuw i 1{2 D i 1{2 vvw i 1{2 . The results with the new dissipation matrix approximate the correct solution and we can see an improvement with increasing order (Fig. 8) and for increasing number of points (Fig. 9).

Conclusions
We introduce a new smoothness indicator and an algorithm to choose the least oscillatory stencil based on RBF interpolation. This smoothness indicator is directly related to the RBF interpolation and it is based on the generalized divided difference method. For this ENO reconstruction we prove the sign-property in the finite difference case for the second and third order reconstruction in the limit ∆x Ñ 0 for infinitely smooth RBFs on general grids. Further, we conjecture this property for higher order schemes and for the case of the average-based interpolation based on numerical experiments. Note that the condition ∆x Ñ 0 can be replaced by the condition for the shape parameter that ε Ñ 0.
Based on this, we construct a RBF-TeCNOp method as an arbitrary high-order entropy stable finite difference method and the RBF-EFV2 method as a second order entropy stable finite volume method. Both are based on high-order entropy conservative schemes and a diffusion term which depends on the RBF-reconstruction in the scaled entropy variables. To circumvent the ill-conditioning of the local interpolation problems we apply the vector-valued rational approximation method [42].
Thus, we propose a method that has all the properties of the original TeCNO scheme [11]. It is entropy stable, high-order accurate for smooth solutions and essentially nonoscillatory near discontinuities. To show their robustness we present a range of numerical simulations in one dimension. The solutions coincide up to a small error with those obtained from the original TeCNO method.
The main drawback of the method is the expensive evaluation of the vector-valued approximation to circumvent ill-conditioning, but this problem is being considered. The advantages of RBF-based reconstructions will become much clearer for high-dimensional problems and we hope to report on this in the near future. with k i r N ¡1 N ¡ 1. By induction one proves that for l N and we recover sgnpC l q p¡1q r N ¡1 N ¡1 l . Note that this case can be characterized by d 2 1 and s 2 r 2 ¡2. From Theorem 5.8 we know and we recover where we used (A.1), (A.2) and (A.4).
Case 2. Assume the stencil S i tC i¡2 , C i¡1 , C i u, S i 1 tC i , C i 1 , C i 2 u, equivalent to the conditions |γ 2 i | |γ 2 i 1 |, |γ 3 i | |γ 3 i 1 |, 5 |γ 2 i 1 | |γ 2 i 2 |, |γ 3 i 1 | ¡ |γ 3 i 2 |, paq |γ 2 i 1 | ¡ |γ 2 i 2 |, |γ 3 i 2 | |γ 3 i 3 |. pbq (A.6) The jump can be written by For each term we calculate its sign. The first term can be understood in the same way as above and it holds for both paq and pbq in (A.6) sgn C 0 pγ 3 i 1 ¡ γ 3 i q¨ sgnpy i 1 ¡ y i q, For the second term we first assume that paq holds and compute its sign as sgn C 1 pγ 3 i 2 ¡ γ 3 i 1 q¨ sgnpγ 3 i 1 ¡ γ 3 i 2 q sgnpγ 3 i 1 q sgnpγ 2 i 1 ¡ γ 2 i q sgnpγ 2 i 1 q sgnpy i 1 ¡ y i q. For pbq we split it in two terms using (A.1) To use this in the entropy stable framework we need to symmetrize it. By assumingD 3,1 E 1 andD 3,2 uE 1 u ρ 2β we get at least the exact jump for the density and mass flow. We recover the matrix