Implicit summation by parts operators for ﬁnite difference approximations of ﬁrst and second derivatives

Implicit ﬁnite difference approximations are derived for both the ﬁrst and second derivatives. The boundary closures are based on the banded-norm summation-by-parts framework and the boundary conditions are imposed using a weak (penalty) enforcement. Up to 8th order global convergence is achieved. The ﬁnite difference approximations lead to implicit ODE systems. Spectral resolution characteristics are achieved by proper tuning of the internal difference stencils. The accuracy and stability properties are demonstrated for linear hyperbolic problems in 1D and the 2D compressible Euler equations. © 2022 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
It is well-known that higher-order methods (as compared to first and second order accurate methods) capture transient phenomena more efficiently since they allow a considerable reduction in the degrees of freedom for a given error tolerance.In particular, high-order finite difference methods (HOFDMs) are ideally suited for problems of this type, cf. the pioneering paper by Kreiss and Oliger [1].
For long-time simulations, it is imperative to use finite difference approximations that do not allow growth in time if the PDE does not allow growth-a property termed time stability [2].Achieving stable HOFDM has received considerable past attention [3][4][5][6][7][8].
SBP operators found in literature (see for example [9][10][11][27][28][29][30][31]) are essentially central FD stencils closed at the boundaries with a careful choice of one-sided difference stencils, to mimic the underlying integration-by-parts formula in a discrete norm.For first derivative finite difference approximations we distinguish between collocated grid and staggered grid SBP operators.Staggered grids SBP operators were recently derived in [32,22].It is well-known [33] that staggered grid finite difference approximations (SGFDA) have a much smaller interior dispersion error compared to collocated grid finite difference approximations (CGFDA).SGFDA are popular for hyperbolic wave equations (acoustic wave equation, Maxwell's equations or the elastic wave equation) on Cartesian geometries.For more general IBVP [34], and/or non-Cartesian geometries [35], staggering makes grid-interpolation necessary, potentially introducing significant numerical errors.It is an open question whether SGFDA are truly efficient when subjected to grid-interpolation.In the present study, we restrict the analysis to CGFDA.
SBP operators may be further categorized by the structure of their norm: a) diagonal, b) diagonal interior with block boundary closures, c) fully banded.The early papers [36,10] report (collocated) first-derivative SBP operators up to 6thorder, based on diagonal-norms and block-norms, while in [37] a 4th-order accurate (collocated) first derivative Padé type banded SBP operator is derived.
For a linear IBVP with smooth data (here referring either to physical data or the underlying curvilinear grid) an SBP-SAT (or an SBP-Projection) approximation constructed from diagonal-norm SBP operators yields a time stable approximation.This is also one of the main reasons why diagonal-norm SBP operators are the most widely used in practice.It is wellknown [38,39] that the usage of block-norm (or banded-norm) SBP operators might lead to instability when applied to IBVP mapped onto curvilinear grids.A remedy to this problem is to add a sufficient amount of artificial dissipation (AD), as shown in for example [39].This is something we will address in the present study.
Unfortunately, however, traditional (here meaning that they are constructed for equidistant grids) diagonal-norm SBP operators (see [9][10][11]) suffer from reduced order of accuracy near boundaries.It has been shown that diagonal-norm SBP operators where all FD stencils have accuracy at least p are associated with a quadrature rule of degree at least 2p − 1.
In an effort to improve the accuracy of first-derivative diagonal-norm SBP operators, [19] introduced a type of boundary optimized operators, of orders up to eight, that was later extended to orders up to twelve in [30].In the construction of the operators, the authors allowed for a few non-equispaced gridpoints near the boundaries.Although the optimized 2pth-order operators were closed with pth order accurate boundary stencils just like their traditional counterparts, the leading order error constants were reduced by several orders of magnitude.SBP-SAT or SBP-projection approximations based on diagonalnorm (or block-norm) SBP operators lead to explicit ODE systems.The traditional and optimized SBP operators all have minimal stencil width in the interior.This is ideal considering the truncation error of the spatial discretization.However, a small truncation error alone does not guarantee an accurate solution of wave-dominated problems where high-frequency content is present (such as in computational aeroacoustics (CAA) [40]).A small truncation error implies that well-resolved waves converge fast, but is not a reliable indicator of the resolution characteristics of shorter length scales, i.e., on coarser grids.A more accurate measure of the propagation characteristics (of FD schemes) is given by the dispersion relation (in Fourier space).Compared to the traditional (explicit) finite difference approximations, the implicit schemes provide a better representation of the shorter length scales (indicated by the dispersion relation).By proper tuning of the interior stencils it is possible to achieve spectral-like (SL) properties for implicit schemes.
Tam and Webb [41] introduced the concept of Dispersive Relation Preserving (DRP) FD schemes, where the main idea was to sacrifice optimal truncation error in favor of a more accurate approximation of the exact dispersion relation, better suited for CAA.Starting from an explicit 4th-order accurate (internal) 7-point stencil, the DRP stencil was derived by tuning the free parameter such that the integrated dispersion relation error is minimized.An extension of the DRP concept for diagonal-norm SBP operators is presented in [42].The dispersion relation is strictly speaking only relevant for periodic solutions and is performed in Fourier space.However, this strategy is shown to yield accurate interior FD stencils applied to bounded domains.How to tune the parameters in the non-central boundary stencils in the presence of non-periodic boundary closures is a much more difficult problem, in particular concerning stability.An attempt to expand the DRP concept for non-periodic boundaries is presented in [29].Explicit FD schemes however have a rather limited ability to capture the exact dispersion relation, unless a close to a full stencil width is used, i.e., the spectral limit.
It is well-known [43,3,44] that implicit FD schemes have much better resolution characteristics than corresponding explicit FD schemes.Implicit FD stencils with minimal stencil width are referred to as Padé (or compact) schemes, see for example [3].It is found that it is not possible to achieve more than first order accurate SBP boundary closures for Padé approximations on equidistant grids, with the exception of a 4th-order SBP operator presented in [37].In this paper, it is shown that a slightly wider stencil approximation is required to achieve an accurate implicit SBP boundary closure, for orders higher than 4. In the present study we derive implicit SBP operators for orders up to twelve, for both first and second derivatives.We also derive novel 4th-and 6th-order accurate implicit SBP operators with SL properties, following the ideas presented in [3].Recently, high-order accurate implicit FD methods have been constructed using a Galerkin framework [45][46][47][48], referred to as Galerkin Difference Methods.In Figs.1-3 we compare the dispersion properties for a few implicit FD schemes, presented in this study, in [3] and in [45].
For nonlinear IBVP the exclusive usage of (collocated) first derivative SBP operators in combination with SAT or Projection does not guarantee time stability.For nonlinear hyperbolic IBVP (such as the Euler equations) the addition of AD is most often necessary in order to damp spurious oscillations.How to add AD in the framework of SBP-SAT approximations was first analyzed in [49] and later applied in [50,51].The addition of AD following the recipe in [49] makes the boundary treatment and AD independent of each other.In the present study we have constructed AD compatible with the implicit first derivative SBP operators.
In Section 2 the SBP definitions are presented, including details of the necessary steps in the derivation of the novel implicit SBP operators.The dispersion properties for the implicit SBP operators are presented in Section 3. In Section 4 the convergence and dispersion properties of the novel implicit SBP operators are verified and compared against the corresponding explicit SBP operators.Verification of accuracy and stability by numerical long-time simulations of an analytic 2D Euler vortex is presented in Section 5. Section 6 summarizes the work.The coefficients of the implicit first-and second-derivative SBP operators are included online as supplementary data for this paper, as described in Section 2.3.

Definitions
Here we will introduce some necessary notation and definitions used throughout the present study.We start with the continuous definitions and then go through the semi-discrete summation-by-parts definitions.Let u, v ∈ L 2 [x l , x r ] where u T = (u (1) . . .u (d) ) and v T = (v (1) . . .v (d) ) are real-valued vector functions with d scalar components.The inner product is defined as with the corresponding norm u 2 ≡ (u, u) .

Summation by parts
SBP operators (see for example [9][10][11]27,19,39,28,31]) are essentially central finite difference stencils closed at the boundaries with a careful choice of one-sided difference stencils to mimic the underlying integration-by-parts formula in a discrete norm.
Remark.SBP operators with non-central finite-difference stencils in the interior were introduced in [20], referred to as upwind SBP operators.Upwind SBP operators are useful for nonlinear problems due to their inherent artificial dissipation.In the present study, we focus exclusively on implicit SBP operators with central finite difference stencils in the interior.Implicit upwind SBP operators are not treated in the present study, although we derive SBP preserving AD to accompany the implicit first derivative SBP operators.
In the present paper, the SBP operators are addressed by the accuracy of the interior finite difference stencil.Finite difference SBP operators may be further categorized by the structure of their norm: a) diagonal, b) diagonal interior with block boundary closures, c) fully banded.For stability reasons (see for example [52,39,16,[53][54][55]) the diagonal-norm SBP operators are the preferred choice in practice.In the present study the focus is on the derivation of implicit (i.e., fully banded) SBP operators on equidistant grids.The domain x l ≤ x ≤ x r is discretized using m + 1 grid points: For an equidistant grid we have ( The following vectors and matrices will be used frequently: ( Note that e T 0 x = x 0 and e T m x = x m .The vectors d 0 and d m are first derivative finite difference stencils at x 0 = x l and x m = x r , and are building blocks of the second derivative SBP operator.
The discrete norm is subject to the following constraint: Definition 2.2.A symmetric positive definite matrix H defines a discrete norm if 1 T H1 = x r − x l independently of the number of gridpoints m + 1.
Similarly, SBP approximations of second order derivatives can be defined as: approximating ∂ 2 /∂ x 2 using a pth order accurate stencil in the interior, is said to be a pth order second-derivative SBP operator if the symmetric positive definite matrix H defines a discrete norm and M = M T ≥ 0.
The matrices B and B D are defined in (3).Traditional second derivative SBP operators were derived in [11,27].
The semi-discrete analysis (when d > 1, i.e., for non-scalar problems) in the present study will utilize the Kronecker product, where M is a p + 1 × q + 1 matrix and N is an m + 1 × n + 1 matrix.Two useful rules for the Kronecker product are The approximate solution vector is given by v (2)   . . .
Here v ( j) denotes the discrete solution vector of the j-th component.Let I d denote the unit matrix of size d × d.To simplify notation we introduce For scalar equations (i.e., when d = 1) we will simply write

Construction of SBP operators
For the reader interested in the construction of the implicit SBP operators, the procedure will be outlined in some detail.First some definitions are needed, before motivating the necessary order of steps.Let x q be the projection of the polynomial x q q !onto the discrete gridpoints, i.e., a vector denoted x.The following two definitions are central to the present study.Definition 2.4.Let e (q) = Hx q−1 − Q + B 2 x q be the qth-order error vector.We say that D 1 is pth order if e (q) vanishes for q = 1 . . .p, in the interior and for q = 1 . . .s, at the boundaries where s ≥ p/2.For diagonal-norm SBP operators s = p/2.Definition 2.5.Let e (q) = Hx q−2 − (−M + B D) x q be the qth-order error vector.We say that D 2 is pth order if e (q) vanishes for q = 1 . . .p + 1, in the interior and for q = 1 . . .s, at the boundaries where s ≥ p/2.For diagonal-norm SBP operators It is important to understand that the formal boundary accuracy alone does not dictate the expected convergence rate of the numerical approximation.The expected convergence rate can be shown by a careful error analysis (not presented in the present study).This is summarized in the following remark: Remark.In [56] it is shown that a stable approximation of an IBVP involving derivatives up to order q yields a convergence rate of order q + b, where b is the order of accuracy at the boundaries.This implies that the convergence rate drops to (b + 1)th order for first order hyperbolic problems, and (b + 2)th order for parabolic problems and second order hyperbolic problems, where b denotes the boundary accuracy.For diagonal-norm SBP operators b = p/2, and the expected convergence rate p/2 + 1.The numerical convergence studies presented in Sections 4 and 5 indicate that the expected convergence rates are obtained.
For implicit finite difference approximations of first-and second-derivatives, the minimum stencil width for a given order of accuracy (often referred to as a compact scheme) leads to something referred to as a Padé approximation, first introduced in [43].On an equidistant grid, it is not possible to derive SBP operators based on the original Padé approximations, strictly following the Definition 2.1 and Definition 2.3.In [37] a 4th-order accurate (collocated) first derivative Padé type banded first derivative SBP operator is derived, but the SBP property is proven using a different structure (as compared to Definition 2.1).
In the present study we derive implicit SBP operators for orders up to twelve, for both first and second derivatives.We also derive implicit SBP operators with SL properties, following the ideas presented in [3].These new implicit SBP operators are not Padé approximations since the stencil widths are not minimal.In [57,58] implicit 4th and 6th order accurate first derivative finite-difference operators with a modified SBP property are presented.These modified SBP operators are of Padé type, but stability proofs for general 2D hyperbolic problems using the energy method is problematic.

Procedure
The first step when building SBP operators is to make sure that the symmetry requirements for the various matrices involved are met, and at the same time include sufficiently many unknowns in the boundary closures to fulfill the accuracy requirements.When building the SL SBP operators, additional unknowns are added, necessary to minimize the dispersion error of the interior stencils.The SBP operators are derived using Maple, a symbolic math software (although any software with symbolic math tools can be used).
The structure of Q is given by the ansatz: The structure of 1 h H is given by the ansatz: The structure of h M is given by the ansatz: Let p denote the interior order of accuracy.In the present study, implicit SBP operators (D 1 and D 2 ) of orders p = 6, 8, 10, 12 are derived.We also derive a set of SL operators of orders 4 and 6.The unknowns in Q , M and H are tuned to obey the accuracy requirements for both D 1 and D 2 .Here it is imperative that H is positive definite, and that both D 1 and D 2 are built on the same norm H for a given order of accuracy.There are "free" parameters left after imposing the accuracy and symmetry conditions.The number of free parameters grow with the order of accuracy.When deriving the SL SBP operators we introduce additional free parameters for the interior stencil, that are tuned to minimize the dispersion error.Building the SL SBP operators is a significantly more difficult task, as compared to the more standard implicit SBP operators.We will explain the additional requirements in Section 3.
A provably efficient strategy [39,20] when tuning these "free" parameters is to minimize the leading order truncation errors e (q) (for q larger than a number given by the boundary accuracy of the SBP operator) defined in Definition 2.4 and Definition 2.5.When tuning these "free" parameters it is imperative that we do not violate the requirements that H is positive definite and that M is positive semi-definite. The

Artificial dissipation
An efficient AD (here disregarding strong nonlinear features such as shocks) should have (at least) the following four characteristics, 1.The linear stability properties of the numerical approximation should not be destroyed when the AD is added.2. It should efficiently damp spurious oscillations.3. It should not severely reduce the convergence properties of the non-dissipative numerical scheme.4. It should not introduce excessive stiffness (as compared to the non-dissipative numerical scheme).
The AD derived in this paper is not intended for strong nonlinear features such as shocks, that requires efficient limiters (see for example [59]).For nonlinear problems with relatively "smooth" solutions, the addition of AD is often necessary to stabilize the numerical solution.In Section 5 we study the effect of AD for a smooth solution to the compressible 2D Euler equations.
For each of the derived SBP operators we have derived an SBP preserving AD (see [49] for more details).The AD is sometimes also beneficial for linear problems with non-smooth parameters, and to stabilize non-diagonal SBP operators when applied to IBVP mapped onto curvilinear grids, as discussed in Section 4.2.2.The AD is added to the Right Hand Side (RHS), and is denoted H −1 S, where the matrix S = S T ≤ 0. The symmetric matrix S for an implicit 2pth order accurate SBP operator is basically S = βh 2p+2 D T p+1 D p+1 , where D p+1 is a consistent approximation of the (p + 1)th order derivative, except at a few (roughly p/2) boundary stencils, that are set to zero (in order to reduce stiffness).β = −1/max(S), where max(S) is the largest parameter in S.This particular choice of S is to guarantee efficiency according to the list above.For the pth order accurate S L operators S is given by S = γ h 2p+4 D T p+2 D p+2 and γ < 0 is chosen to stabilize the Euler vortex problem, presented in Section 5.
The needed amount of AD is somewhat problem dependent.For certain problems you will need less or more AD, and it is up to the user to tune the amount by a positive parameter α multiplying S. The AD is then denoted αH −1 S. In Fig. 7 we present some numerical results where we have tuned α to stabilize the 12th order implicit SBP operator applied to a problem mapped onto a curvilinear grid.For more details about SBP preserving AD we refer to [49,51,20].The exact form of the AD is shown in the MATLAB files (see next section).

SBP operators
The implicit SBP operators are available online as MATLAB function files at https://github .com/ylvarydin /imlicit _sbp.The first derivative SBP operator is given by The SBP preserving AD is given by H −1 S where S = S T ≤ 0. The MATLAB functions return the differentiation matrices M, Q , S, the boundary derivative B D, the norm matrix H , the gridpoint vector x, and the grid spacing h.Inputs to the functions are the grid size m and the width L of the domain.

Dispersion relation
For wave propagation on periodic domains, the propagation characteristics are governed by the (exact) dispersion relation in frequency and wave number space.To analyze the resolution characteristics of the FD schemes, we therefore derive the corresponding discrete dispersion relation.The difference between the exact and discrete dispersion relation quantifies the resolution characteristics of the FD schemes.This analysis (see for example [3,41]) is strictly speaking only valid for periodic domains, but it still says a lot about the interior stencil approximation on bounded domains.It also provides a way to optimize the FD schemes, to better mimic the exact dispersion relation.Such implicit FD schemes are referred to as SL and was first derived in [3].The construction of SL schemes with SBP boundary closures is one of the main results in the present study.
The action of a derivative of order n on a pure Fourier mode e i κ x , results in (i κ) n e i κ x .The second derivative for example gives −κ 2 e i κ x .Consider the same Fourier mode on a grid over [−1, 1] with grid spacing h.The Fourier mode defined on the grid is given by ûT = [e i κ x 1 , e i κ x 2 , • • • , e i κ x N ], assuming periodicity.It is convenient to introduce a scaled wavenumber ω = κ h, where ω ∈ [0, π].The Fourier mode for the wavenumber highest frequency that can exist on the grid).The action of the first derivative on a pure Fourier mode e i κ x , results in ω i h e i κ x , whereas the explicit second order FD approximation would give sin (ω) i h e i κ x .It makes sense to compare the exact value (ω) against the corresponding discrete approximation (sin (ω)).The dispersion relation for the second order FD approximation is given by ω (ω) = sin (ω), where the exact relation is ω (ω) = ω.The dispersion property (see [3] for details) for first derivative (central) FD schemes are given by Remark.For explicit FD schemes h 0 = 1 and h j = 0, for j > 0.
The action of the second derivative on a pure Fourier mode e i κ x , results in ω 2 −1 h 2 e i κ x , whereas the explicit second order FD approximation would give 4 sin 2 ( ω 2 ) −1 h 2 e i κ x .The dispersion relation for the second order FD approximation is given by ω (ω) = 4 sin 2 ( ω 2 ), where the exact relation is ω (ω) = ω 2 .The dispersion property (see [3] for details) for second derivative (central) FD stencils are given by In Figs. 1 and 2 the dispersion properties for explicit and implicit first-and second-derivative difference approximations (i.e., the interior SBP stencils) are presented.

Spectral-like SBP operators
When deriving the SL SBP operators we introduce a wider (than necessary) stencil for a given order of accuracy, i.e., we increase n i and r i .This introduces free parameters in the interior stencil (defined by the parameters h j , q j and m j ) in Q , M and H .In the present study we construct two different SL operators, 4th-and 6th-order accurate (denoted SL 4 and SL 6).
The procedure to construct SL SBP operators is more involved.The first step is to tune the free interior parameters to significantly reduce the dispersion error.For the SL 4 operators we set n i = 4 and r i = 5.This leads to 5 free parameters in The dispersion relation is given by (4).For the corresponding second derivative D 2 there are 2 free parameters tuned such that ω (ω k ) = ω 2 k , where ω 1 = 2.0, ω 2 = 2.6.The dispersion relation is given by (5).
For the SL 6 operators we set n i = 5 and r i = 6.This leads to 5 free parameters in D 1 that are tuned such that ω (ω k ) = ω k , where ω 1 = 1.4,ω 2 = 1.8, ω 3 = 2.1, ω 4 = 2.4, ω 5 = 2.6.For the corresponding second derivative D 2 there are 2 free parameters tuned such that ω (ω k ) = ω 2 k , where ω 1 = 2.0, ω 2 = 2.6.The most difficult part of the procedure is to find an accurate SBP boundary closure, compatible with the optimized interior stencil, in particular for the second derivative.It is necessary to increase the number of boundary stencils (r b ) in D 2 to accomplish this.In Figs. 1 and 2 the dispersion relations for explicit and implicit interior stencils (for the first-and second-derivative SBP operators) are presented.Here we have also included the dispersion relations for the implicit p = 5 approximations presented in [45] (here denoted Banks 6), the 4th order accurate explicit DRP presented in [41] (here denoted Tam), and the explicit SBP DRP scheme (Remez(2, 9) presented in [42] (here denoted Linders).In Fig. 3 we compare the newly constructed SBP operators SL 4 and SL 6 against previously derived SL FD stencils derived in [3] (here denoted Lele 4).A close up of the dispersion errors are presented in Fig. 4.

Numerical computations in 1D
In this section, we will verify the convergence and dispersion properties of the newly derived implicit SBP operators by numerical simulations in 1D.We will verify the resolving capability of the SL SBP operators, by performing periodic long-time simulations on coarse grids.To verify the accuracy of the SBP-SAT boundary closures, we include convergence studies for the SBP-SAT approximations of the wave-equation on first and second order form, comparing the implicit SBP-SAT approximations against the corresponding explicit SBP-SAT approximations.

The advection equation in 1D
We start with the advection equation u t = −u x on a periodic domain x ∈ [0, 1], initialized with a narrow Gaussian profile given by: .
We use 61 gridpoints (i.e., a coarse grid) and integrate the solution, using periodic operators (i.e., using the interior SBP stencils), until t = 500 (500 periods) using a 6th-order accurate Runge-Kutta method [60] with time step dt = 0.2 h.Here  we compare the 12th order implicit and explicit schemes against SL 6 and the explicit 6th-order scheme.The results are presented in Fig. 5.

Wave equation on first order form
To verify the convergence properties of the implicit first derivative SBP operators we consider the 1D wave equation on first order form, where Multiplying the first equation in ( 6) by u T , integrating by parts and adding the transpose leads to The eigenvalues of A are ±1 and one BC at each boundary should be specified.There are many possibilities but here the focus is on homogeneous Dirichlet BC.One set of well-posed BC is for example The semi-discrete approximation of ( 9) is given by where e (1) = (1 0) and e 0, m are given by (3).The semi-discrete problem using the SBP-SAT method reads (1)   m , (11) where τ T l, r = (τ (1) l, r τ (2)   l, r ) are the penalty parameter vectors at the left and right boundary, respectively.The following lemma (first presented in [25]) is relevant to the present study: Lemma 4.1.The scheme (11) = −1, and τ (1)   l, r ≤ 0 hold.
Table 2 log(l 2 −errors) and convergence rates (q) using the explicit SBP operators applied to the wave equation.
m denote the number of gridpoints.Proof.Multiplying ( 11) by v T H and adding the transpose yields, m , where we have used Definition 2.1.Time stability follows if τ (2)   l = 1, τ (2)   r = −1 and τ (1)   l, r ≤ 0, such that Remark.The free SAT parameters τ (1)  l, r = τ (1) in Lemma 4.1 affect the numerical efficiency.A careful numerical study shows (although not presented here) that the optimal value (referring to spectral radius and accuracy) is given by τ (1) = −1.This is the value chosen in the numerical computations of (11).
The analytic solution to ( 6) and ( 9) with the above provided initial data is given by, after the Gaussian pulses have been reflected at the boundaries, i.e., when t ∈ [L − L/4, L + L/4].In the numerical simulations r * = 0.1.The numerical approximations are integrated to t = 1.8 using a 6th-order Runge-Kutta method, with a time step dt = 0.1 h.The waves reach the boundaries around t = 0.9.In Tables 2 and 3, the convergence results for the explicit SBP operators are compared to the corresponding implicit SBP operators.The results for the spectral-like operators are found in Table 4.In Fig. 6 a long-time simulation is presented, comparing the explicit and implicit 12th order SBP-SAT approximations of the 1D wave equation.The numerical approximations are integrated to t = 101.8(i.e., after 51 reflections against the boundaries) using a 6th-order Runge-Kutta method, with a time step dt = 0.2 h, using 71 gridpoints.This test is primarily done to verify the accuracy of the boundary closures of the SBP operators.
Remark.The discretization error for the internal (SBP) difference stencil is negligible compared to the discretization errors of the boundary stencils, for the standard explicit SBP operators.

Spectral radius
The time step restriction (or the CFL number) when employing an explicit time-integrator, depends on the spectral radius of the solution matrix, here given by Table 3 log(l 2 −errors) and convergence rates (q) using the implicit SBP operators applied to the wave equation.m denote the number of gridpoints.Remark: Machine precision reached at m = 801 for 10th and 12th order methods.log(l 2 − errors) and convergence rates (q) using the two SL SBP operators applied to the wave equation.m denote the number of gridpoints.      m .
A time stable choice of the penalty parameters is given by τ (2)   l = 1, τ (2)   r = −1 and τ (1)   l, r = 0, resulting in energy conservation (see proof of Lemma 4.1), which means that the eigenvalues to R reside on the imaginary axis.In Table 5 we present the spectral radius of h R when employing the explicit and implicit SBP operators.As expected the spectral radius grows with the order of accuracy, in particular for the implicit SBP operators.

Stability on curvilinear grids
Here we study the robustness (or stability) of the implicit SBP operators when applied to the hyperbolic system (6) mapped onto a curvilinear grid x = x(ξ ), where J = ξ −1 x is the Jacobian matrix.The semi-discrete SBP-SAT approximation (11) for the mapped problem with addition of AD reads, (1)  0 + τ r ⊗ H −1 e m v (1)  m + α H S , (14) where ¯J = I 2 ⊗ J and J is the projection of the Jacobian matrix onto the diagonal, and S = I 2 ⊗ S is the AD derived for each of the implicit SBP operators.
Multiplying ( 14) by v T H and adding the transpose yields, for the energy conservative choice of the penalty parameters, i.e., by setting τ (2)   l = 1, τ (2)   r = −1 and τ (1)   l, r = 0. Obtaining an energy estimate requires that the following holds such that the LHS in ( 15) can be written as Hence, time stability on curvilinear grids can only be guaranteed if the norm H is diagonal.For implicit SBP operators, H are banded matrices, and time stability cannot be guaranteed.However, in [39] it is shown that ( 16) is not a necessary condition for stability, when applying SBP operators with non-diagonal norms to IBVP on curvilinear grids.But sometimes AD is needed to achieve a time stable approximation.
As we will show, the possible instability (i.e., positive eigenvalues) vanish with increased grid resolution.To see this we split H J into a symmetric and a skew-symmetric part leading to Since (H J) S ymm is symmetric and positive definite we identify the term destroying the energy estimate in (15), i.e., v T (H J) Skew v t .This term is of order O(h).Hence we expect that an eventual instability (i.e., positive eigenvalues) will disappear with grid refinement, i.e., letting h decrease.This is in fact what we see in Fig. 7. Another observation [39] is that the instability on a given grid is more likely to appear, the higher we go in order of accuracy.Hence, for a given, relatively coarse grid, and high order of accuracy, we can not guarantee a time stable solution when employing the implicit SBP operators.However, the AD term 2α v T S v in the RHS of (15), is constructed to add efficient artificial damping without destroying the accuracy of the approximation.
In the numerical experiments we use the following mapping The solution matrix, including AD is given by, (1)  0 + τ r ⊗ H −1 e m v (1)  m + α H S .
Let ρ(h R) denote the spectral radius of h R and η the largest real part of the spectrum (eigenvalues to h R).For a time stable approximation it is necessary that η ≤ 0. (If η ≤ h, for some constant the approximation is still stable, but not time stable.) Remark.The conclusion from various numerical approximations of IBVP on curvilinear grids (not presented here) is that all implicit SBP operators are time stable without the addition of AD, i.e., with α = 0, with the exception of the 12th order accurate SBP operator.This is of course not a strict proof, but rather a strong indication that all (except the 12th order) implicit SBP operators yield time stable approximations, when applied to well-posed (energy conservative) IBVP on curvilinear grids.
In Fig. 7 we present the spectral radius ρ(h R) and η (the largest eigenvalue to h R) employing the 10th and 12th order accurate implicit SBP operators, for the case with no AD (α = 0) and with addition of AD.We test this for different levels of grid refinements.For the 12th order case α = 20 is needed to stabilize the approximation on all grids.We also see that the 12th order approximation becomes time stable with α = 0, for sufficiently fine grids (here around m = 250).Setting α = 20 introduces some stiffness.The value α = 1 is the upper limit, such that the addition of AD does not interfere with the CFL number (as compared to the approximation without AD).The 10th order implicit SBP operator yields time stable approximations on all grids, without AD.This is also true for all of the other implicit SBP operators presented in this paper (although we have omitted the numerical results here).

Wave equation on second order form
Consider the wave equation, Here we consider both periodic and Neumann boundary conditions.The Gaussian profile θ 1 (x, t) is given by (12).The width of the Gaussian profile is set to r * = 0.1.In the first numerical test we verify the dispersion properties of the interior second derivative finite difference stencils.We use 41 gridpoints (i.e., a coarse grid) and integrate the solution, using periodic operators (i.e., using the interior second derivative SBP stencils), until t = 500 (500 periods) using a 6th-order accurate Runge-Kutta method [60] with time step dt = 0.2 h.Here we compare the 12th order implicit and explicit schemes against SL 6 and the explicit 6th-order scheme.The results are presented in Fig. 8.To verify the convergence properties of the second derivative SBP operators we now consider a bounded domain, and impose homogenous Neumann boundary conditions.The semi-discrete problem using the SBP-SAT method (see [61,27] for more details) reads The analytic solution to (18) when imposing Neumann boundary conditions is given by, u(x, t) = 1 2 θ (1) (x, L −t) − 1 2 θ (2) (x, L − t), after the Gaussian pulses have been reflected at the boundaries, i.e., when t ∈ [L − L/4, L + L/4].In the numerical simulations r * = 0.1.The numerical approximations are integrated to t = 2.2 using a 6th-order Runge-Kutta method, with a time step dt = 0.1 h.The waves reach the boundaries around t = 0.9.In Tables 6 and 7, the convergence results for Fig. 8. Comparing the dispersion properties for the periodic second order wave equation.Integrated to t = 500 using 41 gridpoints.Analytic solution given by red solid line.The numerical solutions are given by the blue dotted lines.Table 6 log(l 2 − errors) and convergence rates (q) using the explicit D 2 SBP operators applied to the wave equation.m denote the number of gridpoints.Remark: Machine precision reached at m = 801 for 10th order method.the explicit second derivative SBP operators are compared to the corresponding implicit SBP operators.The results for the spectral-like operators are found in Table 8.

Euler vortex in 2D
To test the dispersion and convergence properties of the novel SL SBP operators for a more challenging nonlinear problem, we consider an Euler-vortex that satisfies the 2D Euler equations, under the assumption of isentropy.Here, we solve the problem in a periodic setup with an SBP-SAT coupling of the left and right boundaries.The domain size is 12 × 12.
Details concerning how to discretize this problem with the SBP-SAT method can be found in [50,39].This cyclic SBP-SAT coupling makes it possible to simulate a long-time simulation of the Euler-vortex problem, where the vortex crosses multi-Table 7 log(l 2 − errors) and convergence rates (q) using the implicit D 2 SBP operators applied to the wave equation.m denote the number of gridpoints.Remark: Machine precision reached at m = 401 for 12th order method.exp ( f (x, y, t)) where f (x, y, t) = 1− ((x−x 0 )−t) 2 +y 2 r 2 * and x 0 is the initial position of the vortex (in the x-direction).The vortex is introduced into the computational domain by using the analytic solution as boundary data and initial data.In the present study the Mach number is set to M = 0.5, and the nondimensional effective radius r * = 4.The vortex is initialized in the middle of the domain, i.e., x 0 = 0 and y 0 = 0.This means that the vortex is centered at the SBP-SAT interface (coupling right and left boundary) at t = 6.The numerical approximations are integrated to t = 6 using the classical 4th-order Runge-Kutta method, with a time step dt = 0.1 h.In Fig. 9 the initial data (pressure contours), and the numerical solution at t = 6 are presented, using the SL 6 operator and 60 2 gridpoints.To stabilize the vortex we add SBP preserving AD, as defined in Section 2.2.2.Here we set α = 1.We add identical amount of AD for the explicit and implicit SBP-SAT approximations, for each order of accuracy (4 and 6).For more details how to add AD to the 2D Euler equations we refer to [50,39].The exact form of the AD is shown in the MATLAB files, that can be downloaded (see Section 2.3).In Table 9, the convergence results for the SL 4 and SL 6 SBP-SAT approximations are compared to the explicit 4th-and 6th-order accurate SBP-SAT approximations.
In Fig. 10 we present the l 2 -error as a function of time in a long-time simulation (up to t = 1000), comparing SL 4, SL 6 against the corresponding explicit 4th-and 6th-order SBP operators using 60 2 gridpoints.The vortex is initially centered at (x, y) = (0, 0) and integrated using RK4 with a time step dt = 0.1 h.

Conclusions
The main focus has been to construct novel banded-norm (implicit) first-and second-derivative SBP operators up to order 12.We have also constructed spectral-like implicit SBP operators by optimizing the dispersion properties of the         internal discretization.The boundary conditions are treated with the SAT technique.Numerical computations show that the novel implicit SBP operators have superior dispersion properties, compared to the corresponding explicit diagonal-norm SBP operators.In particular we show that the spectral-like SBP operators are superior for long-time simulations on coarse grids.This was also the main motivation of constructing them.More work is still needed on implementing these implicit SBP operators for large 2D and 3D implementations, and further analyze the efficiency compared to utilizing explicit SBP operators.In a coming study we will explore the possibility of combining SL properties and optimal boundary closures, by allowing non-equidistant grids close to the boundaries.

Fig. 3 .
Fig.3.Comparing the dispersion properties for "spectral-like" first derivatives (left) and second derivatives (right).Here we also include the 12th order explicit and implicit (interior SBP stencils) as relevant references.Banks 6 is the p = 5 approximation presented in[45], Lele 4 is a 4th order SL approximation presented in[3].

Fig. 4 .
Fig. 4. Comparing the dispersion error for SL first derivatives (left) and second derivatives (right).

Fig. 5 .
Fig. 5. Comparing the dispersion properties for the periodic advection equation.Integrated to t = 500 using 61 gridpoints.Analytic solution given by red solid line.The numerical solutions are given by the blue dotted lines.(For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.Comparing the long-time accuracy of the explicit (left) and implicit (right) SBP-SAT approximations of the 1D wave equation.Integrated to t = 101.8using 71 gridpoints.The implicit SBP-SAT solution is very close to the exact solution.

Fig. 7 .
Fig.7.The spectral radius ρ(h R) and the largest eigenvalue to h R, when employing the 10th and 12th order accurate implicit SBP operators, for the case with no AD (α = 0) and with addition of AD.For the 12th order case α = 20.

2 , 2 ρ = 1 − 2 (γ − 1 )M 2 8π 2 p
where is the nondimensional circulation, r * the nondimensional effective radius and (r, ) the polar coordinates.The strength of the vortex is governed by the parameter and the size of the vortex is governed by the parameter r * .The analytic solution, in the nondimensionalized primitive variables: density (ρ), velocity in x-direction (u), velocity in y-direction (v) and pressure (p), in a fixed frame of reference (x, y, t) becomes,u = 1 − y 2π √ p ∞ r 2 * exp f (x,y,t) v = ((x−x 0 )−t) 2π √ p ∞ r 2 * exp f (x,y,t) ∞ r 2 *

Fig. 9 .
Fig. 9.The setup for the Euler vortex problem.The vortex is initially (left subfigure) centered at (x, y) = (0, 0) and integrated using RK4 to t = 6 (right subfigure).Here presenting the pressure contours.Employing the SL 6 operator and 60 2 gridpoints.At t = 6 the vortex has propagated with the free-stream to the interface between the left and right boundary.

Table 9 Convergence
2D Euler vortex t end = 6, with the addition of AD.

Definition 2.1. A
difference operator D 1 = H −1 Q + B 2 approximating ∂/∂ x, using a pth order accurate interior stencil, is said to be a pth order first derivative SBP operator if the symmetric positive definite matrix H defines a discrete norm, and

Table 1
The characteristics of the new implicit SBP operators.
matrices ( Q , M and H ) building the SBP operators have different number of boundary stencils.For H and Q the number of boundary stencils are denoted n b , whereas the width of the interior stencils are 2 n i + 1.For M the number of boundary stencils are denoted r b , whereas the width of the interior stencils are 2 r i + 1.We denote the boundary accuracy of D 1 and D 2 by b 1 and b 2 respectively.In Table1the values for b 1 , b 2 , n i , n b , r i and r b are listed for the derived SBP operators.

Table 5
The spectral radius of the solution matrix R multiplied by h for the explicit and implicit SBP operators.

Table 8
log(l 2 − errors) and convergence rates (q) using the two SL D 2 SBP operators applied to the wave equation.m denote the number of gridpoints.
This will allow us to verify both the long-time stability behavior and the boundary accuracy, employing the novel SL SBP operators.The analytic vortex solution is steady in the frame of reference moving with the free-stream.Let p ∞ = (γ M 2 ) −1 denote the nondimensional background pressure, where M is the Mach number and γ = c p /c v = 1.4 (in air).The scaled vortex has the velocity field