Analysis of the SBP-SAT Stabilization for Finite Element Methods Part I: Linear Problems

In the hyperbolic community, discontinuous Galerkin (DG) approaches are mainly applied when finite element methods are considered. As the name suggested, the DG framework allows a discontinuity at the element interfaces, which seems for many researchers a favorable property in case of hyperbolic balance laws. On the contrary, continuous Galerkin methods appear to be unsuitable for hyperbolic problems and there exists still the perception that continuous Galerkin methods are notoriously unstable. To remedy this issue, stabilization terms are usually added and various formulations can be found in the literature. However, this perception is not true and the stabilization terms are unnecessary, in general. In this paper, we deal with this problem, but present a different approach. We use the boundary conditions to stabilize the scheme following a procedure that are frequently used in the finite difference community. Here, the main idea is to impose the boundary conditions weakly and specific boundary operators are constructed such that they guarantee stability. This approach has already been used in the discontinuous Galerkin framework, but here we apply it with a continuous Galerkin scheme. No internal dissipation is needed even if unstructured grids are used. Further, we point out that we do not need exact integration, it suffices if the quadrature rule and the norm in the differential operator are the same, such that the summation-by-parts property is fulfilled meaning that a discrete Gauss Theorem is valid. This contradicts the perception in the hyperbolic community that stability issues for pure Galerkin scheme exist. In numerical simulations, we verify our theoretical analysis.


Introduction
In recent years, significant efforts have been made to construct and develop high-order methods for the solution of hyperbolic balance laws, and most of the common methods are either based on finite difference (FD) or finite element (FE) approaches.In the FE framework, one favorable, if not the most favorable scheme, seems to be the discontinuous Galerkin (DG) method introduced by Reed and Hill [1] because of its good stability properties [2,3,4,5].In the stability proofs, so called summation-by-parts (SBP) operators are used [6,7,8,9,5,10].SBP operators originate in the FD framework [11] and lead to a way to demonstrate stability similar to the one in the continuous analysis [12,13,14].Together with SBP operators, Simultaneous Approximation Terms (SATs) that impose the boundary conditions weak are applied.The SBP-SAT technique is powerful and universally applicable.Certainly, one of the reasons for the popularity of DG is that the numerical solution is allowed to have discontinuities at the element boundaries, and, since non-linear hyperbolic problems are supporting shocks (like Riemann problems), this property is believed to be desirable.Another reason, maybe the most important, is that DG methods leads to block diagonal mass matrices which are easy to invert.The difference between a DG approach and continuous Galerkin (CG), besides the structure of the mass matrix, is that in CG the approximated solution is forced to be continuous also over the element boundaries and this restriction seems to be quite strong also in terms of stability where a common belief in the research community is that a pure CG scheme is notoriously unstable1 .Therefore, stabilization terms have been developed and are frequently applied to remedy this issue [15,16,17].Even there exist some preliminary stability results [18,19,20] including the procedure at the boundary where the main idea is to switch the norm of the trial space.However, these results may seem forgotten in the community.In this paper, we also focus on the stability property of a pure Galerkin scheme, but follow a different approach.Our preliminary idea/thought is: If one considers the DG method with one element, the method is stable through the investigation done in the literature mentioned above.There is nothing that says that the approximation space must be a broken polynomial space, the only thing that is needed is that the trial and test function must have some kind of regularity within the elements, so that the divergence theorem (or SBP techniques) can be applied.Continuity is enough.Hence, one can see a CG method as a DG one, with only one element (the union of the simplex) with an approximation space made of polynomials with continuity requirement between the simplex.Hence, what is the difference between these two approaches?The answer to this question points to the procedure at the boundary.In the stability proofs, the usage of SATs is essential.However, up-to-our-knowledge SATs have never been used together with a pure CG scheme and this is the topic of this paper.We divide the paper as follows: In the second section, we shorty introduce the continuous Galerkin scheme which is used and investigated in the following.Next, we introduce and repeat the main idea of the SAT procedure from the FD framework and extend it to the Galerkin approach.We show that the determination of the boundary operators is essential.In section 4, we focus on the eigenvalue analysis of the spatial operators and derive conditions from the continuous setting to build adequate boundary operators in the discrete framework.We give some recipes which will be used in section 5 to support our analysis in numerical experiments.Finally, we conclude and discuss future work.

Continuous Galerkin Scheme
In this section, we shortly introduce the pure continuous Galerkin scheme (CG) as it is also known in the literature [21,16,11].We are interested in the numerical approximation of a hyperbolic problem with suitable initial and boundary conditions.Later, we will focus on the boundary condition, but for the explanation of CG this is not important.The domain Ω is split into subdomains Ω h (e.g triangles/quads in two dimensions, tetrahedrons/hex in 3D).We denote by K the generic element of the mesh and by h the characteristic mesh size.Then, the degrees of freedom (DoFs) σ are defined in each K: we have a set of linear forms acting on the set P k of polynomials of degree k such that the linear mapping ) is one-to-one.The set S denote the set of degrees of freedom in all elements.The solution U will be approximated by some element from the space V h defined by A linear combination of basis functions ϕ σ ∈ V h will be used to describe the numerical solution As basis functions we are working either with Lagrange interpolation where the degrees of freedom are associated to points in K or Bézier polynomials.
To start the discretisation, we apply a Galerkin approach and multiply with a test function V h and integrate over the domain.This gives Using the divergence theorem, we get By choosing V h = ϕ σ for any σ ∈ S and inserting (3), we obtain a system of equations: that in practice we compute using a quadrature rule: where represents the quadrature rules for the volume and surface integrals.
In this paper, we are considering only linear problems, i.e. the flux is linear in U , but may depend on the spatial coordinate.In all the numerical experiences, we will make the spatial dependency simple enough (i.e.typically polynomial in x), so that it will always be possible to find a standard quadrature formula that make the formula exact.In other words, in this paper we proceed such that ( 6) is always exactly reproduced, unless it is specified.Using a matrix formulation, we obtain the classical FE framework: where U h denotes the vector of degrees of freedom, F is the approximation of div f and M is a mass matrix2 .In case of continuous elements, this matrix is sparse but not block diagonal, contrary to the discontinuous Galerkin methods.It is well-known that the continuous Galerkin scheme suffers from stability issues.Therefore, it is common to add stabilization terms to the scheme as for example in [17].However, we follow a different approach in this paper and will renounce these classical stabilisation techniques.In order to do this, we need more known results from the literature, which we will briefly repeat here.
3 Weak Boundary Conditions

SATs in SBP-FD framework
Weak boundary conditions implemented using simultaneous approximation terms (SATs) was originally developed in the finite difference (FD) framework.Together with summation-by-parts (SBP) operators it provides a powerful tool for proofs of semidiscrete (L 2 ) stability of linear problems by the energy method, see [14,12,23] for details.
Here, we present a short introductory example of the SBP-SAT technique as it presented [22,14].Consider the linear advection equation where u in is the initial condition and b is the boundary data in L 2 that is only define on the inflow part of ∂[0, 1] = {0, 1}.In other words, if a > 0, the b is only set for x = 0, and if a < 0, this will be for x = 1 only.
A semi-discretisation of ( 8) is given in terms of SBP operators.
with u = (u 0 , u 1 , , • • • , u N (t)) T are the coefficients of u and similarly for u in .The coefficients correspond to the degrees of freedom in the finite element setting and are used to express the numerical solution (3) in the grid points.M is a symmetric positive definite mass matrix which approximates the usual L 2 scalar product.The term D is a difference matrix.This is exemplified below as Instead of having an extra equation on the boundary like in (8), the boundary condition is enforced weakly by the term S = (S 0 , 0, • • • , S N ) T which is called the SAT.We will focus on this operator more precisely and demonstrate how it should be selected to guarantee stability for (9).
Definition 3.1.The scheme (9) is called strongly energy stable if holds.The term K(t) is bounded for any finite t and independent from u in , b 0 and the mesh.
Remark 3.2.The definition 3.1 is formulated in terms of the initial value problem (8) where only one boundary term is fixed.Indeed, extensions are straightforward.If an additional forcing function is considered at the right hand side of (8), we have to include the maximum of this function in (11) in the spirit of b 0 , for details see [14].
As established in [24], we can prove the following: Let a + = max{a, 0} and a − = min(a, 0), b 0 = b(0, t) and b N = b(1, t).If S 0 = τ 0 a + (u 0 − b 0 ) with τ 0 ≤ − 1 2 , and S N = τ a − N (u N − b N ) with τ N ≤ − 1 2 then the scheme (9) is energy stable.Proof.Multiplying (9) with u T M yields Transpose (13) and adding both equations together leads to Further, we obtain from ( 12) Looking at the cases a > 0 and a < 0, if This shows that boundary operator S can be chosen in such way that it guarantees stability for the SBP-SAT approximation of (8).Next, we will apply this technique in the Galerkin framework.

SATs in the Galerkin-Framework
Instead of working with SBP-FD framework we consider now for the approximation of (8), a Galerkin approach.In [6,5] it is presented that also the discontinuous Galerkin spectral element method satisfies a discrete summation-by-parts (SBP) property and can be interpreted as an SBP-SAT scheme with diagonal mass matrix.As we described already in section 2, the differences between the continuous and discontinuous Galerkin approach are the solution space (2) and structure of the mass matrix (7) which is not block diagonal in CG.However, the approach with SAT terms can still be used to ensure stability also in case of CG but one has to be precise, as we will explain in the following.Let us step back to the proof of Proposition 3.3 and have a closer look on it.Essential in the proof is condition (12).Let us focus on this condition for a FE based discretisation of (8) as described also in [22].We approximate (8) now with u h (x, t) = N j=0 u h j (t)ϕ j (x) where ϕ j are basis functions and u h j are the coefficients.Let us assume that ϕ j are Lagrange polynomials where the degrees of freedoms are associated to points in the interval.Introducing the scalar product let us consider the variational formulation of (8) with test function ϕ i and inserting the approximation yields and finally where In matrix formulation ( 14) can be written as it is also described in [22].Let us check (12).Therefore, we consider If the boundaries are included in the set of degrees of freedom, then we obtain Up to this point exact integrals are considered but the same steps are valid if a quadrature rule is applied with sufficient accuracy, especially (16) has to hold exactly.If this is the case, we can apply the proof of proposition 3.3 and demonstrate: Proposition 3.4.If the Galerkin method described above together with a SAT approach is applied to (8).
is used then the described scheme is energy stable.The weak formulation of the problem writes: For simplicity, we consider the case a > 0 only.The SAT techniques adds a penalty term into the approximation (8) on the right side.We focus now on the energy.Therefore, we multiply also with u h instead of ϕ i and rearrange the terms.We obtain for the semi-discretization ( 14): where we used the fact that By following the steps of the proof of proposition 3.3 and using the above considerations we get the final result.
In the derivation above, we restricted ourselves to one-dimensional problems using Lagrange interpolations.Nevertheless, this shows that a continuous Galerkin method is stable if the boundary condition is enforced by a proper penalty term.For the general FE semi-discretization of ( 7) it is straightforward to the procedure.For a general linear problem (scalar or systems) the formulation (7)3 can be written with penalty terms as where Π is the boundary operator which includes the boundary conditions.A U h represents the flux and Q 1 the spatial operator.The matrix Π is usually sparse and we can formulate the following.
Theorem 3.5.Apply the general FE semidiscretization (17) together with the SAT approach to a linear equation and let the mass matrix M of (17) be symmetric.If the boundary operator αΠ together with the discretization Q 1 can be chosen such that has only non-positive eigenvalues ∆ , then the scheme is energy stable.
Proof.We use the energy approach and multiply our discretization with U h instead of ϕ i and add the transposed equation using M T = M .We obtain Remark 3.6.This theorem yields directly conditions for a FE based stable scheme.If (18) is not fulfilled, stabilization terms have to be added.However, no internal stabilization terms are necessary when (16) holds.For this result to hold, a number of requirements are needed.The distribution of the degrees of freedom should be suitable for the problem and the mesh and the quadrature rule must be chosen to guarantee exactness of all the calculations.In the numerical test, we will present an example of what happens if the quadrature rule is not sufficiently exact.
Furthermore, in case of a non-conservative formulation of the hyperbolic problem or in case of variable coefficients a skew-symmetric formulation must be applied in the way as it is described in [22,25,26] .
If the implementation of the continuous Galerkin method is done in such way that the condition (18) are fulfilled, then by applying the SAT technique the method is stable only through our boundary conditions.To our opinion, this is contrary to common belief about continuous Galerkin methods.The only stabilizing factor needed is a proper implementation of boundary conditions.For the linear scalar case, the proof is given in (3.4).In the following, we will extend this theory to more general cases.

Estimation of the SAT-Boundary Operator
As described before, a proper implementation of the boundary condition is essential for stability.Here, we give a recipe how these SAT boundary operators can be chosen to get a stable CG scheme for different types of problems.First, we are considering a scalar equation in 2D and transfer the eigenvalue analysis for the spatial operator from the continuous to the discrete setting.Then, we extend our investigation to 1D systems.Using again the continuous setting, we develop estimations for Π and transfer the results to the finite element framework.Finally, we extend the investigation to the system case in two dimensions which will also be used in the numerical section later.

Eigenvalue Analysis
In the following part, we derive the conditions on the boundary operators and perform an eigenvalue analysis to get well posedness in the continuous setting.Next, the results are transformed to the discrete framework to guarantee stability of the discrete scheme.

Continuous Setting
Consider the initial boundary value problem Without loss of generality, it is enough to consider the homogeneous boundary conditions and we consider the spatial operator considered in the subspace of functions for which Bu = 0.This operator will be dissipative if u, Du > 0.
Using the Gauss-Green theorem, we obtain The operator will be dissipative if ∂Ω a n u 2 ds > 0. Looking at this question amounts to looking at the spectral properties of the operator D. The question rises: How do we guarantee this condition?This is the role of the boundary conditions, i.e. when a n ≤ 0, we need to impose u = 0.For outflow, i.e. ∂Ω out we get a n > 0 and using this information, we directly obtain We do not go further into details about well posedness, but we recommend [23,27,28] for details.Now, we transfer our analysis to the discrete framework and imitate this behavior discretely.

Discrete Setting
We have to approximate the spatial operator D and the boundary condition (B.C), i.e.Du + B.C which is approximated an operator of the form M −1 (Q − αΠ)u.The term M −1 Q u approximates Du where αΠu is used to describe Bu weakly.M denotes the mass matrix4 and M −1 Q the derivative matrix, in the context of SBP [14].Here, the projection operator Π works at the boundary points.
Looking at the dissipative nature of M −1 Q u amounts to study its spectrum.The eigenvalue problem is We denote by ũ * ,T , the adjoint of ũ, multiply with ũ * ,T M and obtain We transpose (24) and add both equations together.This results in The boundary terms (BT) correspond to ∂Ωout a n u 2 ds with a properly chosen Π.The matrix is positive semi-definite.However, the remaining terms in the boundary term makes Re(λ) > 0, i.e. the eigenvalues for the spatial operator have a strictly positive real parts only.
Next, we estimates the boundary operators such that theorem 3.5 is fulfilled and a pure CG scheme is stable.
We start with the continuous energy analysis and derive the estimate above.Afterwards, we translate the result to the discrete FE framework as done for the scalar one-dimensional case.

Continuous Energy Analysis for 1D Systems
First, we extend our result to the 1D system case.The problem under consideration is given by: where A is a n × n symmetric matrix, and we assume L 0 and L 1 to be linear.If n 0 is the number of incoming characteristics at x = 0 and n 1 the number of incoming characteristics at x = 1, we know that the rank of L 0 (resp L 1 ) is n 0 (resp n 1 ).The system (26) admits an energy: if we multiply ∂U ∂t + A ∂U ∂x = 0 by U T , we first get The energy To understand the role of boundary conditions, we follow what is usually done for conservation laws, we consider the weak form of ( 26): let ϕ be a regular vector function in space and time.We multiply the equation by ϕ T , integrate and get: In order to enforce the boundary conditions weakly, we modify this relation by: The operators Π 0 and Π 1 are selected in such a way that 1.For any t, the image of the boundary operator L 0 is the same as the image of Π 0 L 0 , i.e. there is no loss of boundary information, and the same applies for Π 1 and L 1 .
A solution to this problem is given by the following: First, let A = XΛX T where Λ are the eigenvalues of A and X is the matrix which rows are the right eigenvectors of A. Secondly, we have X T X = I and choose: (28) where Λ − are the negative eigenvalues only.Here we have introduced the operator R 0 and R 1 which are just L 0 and L 1 written using characteristic variables, and not U .We will write in the sequel.
Step 1: Strong Implementation First, we consider again the strong implementation of the problem.
For simplicity, we look only at what occurs for x = 0. Let A = XΛX T where Λ are the eigenvalues of A and X is the matrix which rows are the right eigenvectors of A. Then, we obtain with W + = X T U + are the ingoing waves and they have the size of the positive eigenvalues Λ + .Analogously, W − = X T U − are the outgoing waves with size of Λ − .A general homogeneous boundary condition is W + = RW − , since with that form we get if the matrix in the bracket is negative semidefinite.
Step 2: Weak Implementation Assume now that we have chosen an R such that The energy is given (remember ∂Ω = {0, 1}) Second, we consider everything at the boundary Since ∂Ω = x=0 + x=1 , we focus only on the x = 0 part.We define Π such that U T Π = (W + ) T Π and get for the integrands Collecting the terms, we obtain We must select Π such that the matrix W B is negative definite.Now, let us use the fact that we have the strong condition (31).By adding and subtracting, we obtain We re-order Q: We choose Π = αΛ + with α scalar and get We choose α such that the matrix is negative semidefinite.With α = −1, it is G has the eigenvalues 0 and 2 and we obtain stability thanks to (36) and (31).

Extension to 2D Symmetric Systems
Next, we will extend our investigation to the general hyperbolic system where A, B ∈ R m×n are the Jacobian matrices of the system, the matrix L n ∈ R q×m and the vector G n ∈ R q are known, n is the local outward unit vector, q is the number of boundary conditions to satisfy.We assume that A, B are constant and the system (37) is symmetrizable.It exists a symmetric and invertible matrix P such that for any vector n = (n x , n y ) T the matrix is symmetric with A n = An x + Bn y .Using the matrix P , one can introduce new variables V = P −1/2 U .The original variable can be expressed as U = P 1/2 V and the original system (37) will become Multiplying (38) form the left by P −1/2 we obtain the system which is symmetric since P −1/2 AP 1/2 n x + P −1/2 BP 1/2 n y = P −1/2 B n P −1/2 .Focusing on the boundary treatment of the problem and using the weak formulation, we get for the system (37): where δ ∂Ω is Dirac distribution on ∂Ω 5 and Π is our boundary projection operator.

Energy balance
Again, we start by considering the energy balance for the weak formulation (40) in the continuous setting.We define the global energy of the solution of (37) by where we take into account the symmetrizability of the system (37).We multiply (40) and integrate over Ω leads to We reformulate the left-hand side of (42).
For stability, the energy does not increase, i.e. dE dt ≥ 0 which is guaranteed if the integral term of the left-hand side of ( 44) is non-negative.Therefore, 1 2 has to be fulfilled.Inequality (45) imposes the restrictions on the choice of the projection operator Π n .Switching to the discrete FE setting, the information of (45) will be used to determine Π n in the following.We will give a concrete example in 5.3.

Numerical Simulations
Here, we demonstrate that a pure Galerkin scheme is stable if the boundary procedure is done via the SAT approach.This means we impose the boundary conditions weakly and use an adequate boundary operator as developed in section 4. We focus on several different examples and analyze different properties in this context (error behavior, eigenvalues, etc.).As basis functions, we use Bernstein or Lagrange polynomials of different orders (second to fourth order) on triangular meshes.Nearly, no differences can be seen.The time integration is done via strong stability preserving Runge-Kutta methods of second to fourth order, see [29] for details.We use the same order for space and time discretization.

Two-Dimensional Scalar Equations
We consider a two-dimensional scalar hyperbolic equation of the form where a = (a, b) is the advection speed and Ω the domain.In this subsection, the initial condition is given by It is a small bump with height one located at (x 0 , y 0 ).We consider homogeneous boundary conditions G n ≡ 0, we do further let the boundary matrix M n be the identity matrix.The boundary conditions reads M n U = U = 0 for (x, y) ∈ ∂Ω, t > 0 which means that the incoming waves are set to zero.

Linear advection
In our first test, we are considering the linear advection equation in Ω = [0, 1] 2 .The advection speed is assumed to be constant.The components of the speed vector a are given by (a, b) T = (1, 0) and so the flux is given by f (U ) = aU with a = (1, 0).We get inflow / outflow conditions on the left / right boundaries and periodic boundary condition on the horizontal boundaries.In our first test, we use Bernstein polynomials and a fourth order pure Galerkin scheme.The boundary operators are computed using the technique developed in section 4 where the positive eigenvalues are set to zero and the negative ones are used in the construction of Π.For the time discretization we apply strong stability preserving Runge-Kutta (SSPRK) scheme with 5 stages and fourth order with CFL=0.3.We use 1048 triangles.In figure (1), we plot the results from a pure Galerkin scheme with SAT terms at three times.Clearly, the scheme is stable, also at the boundary.Next, we check the real parts of the eigenvalues of our problem using formula (18) for different orders, different bases (Bernstein and Lagrange) and different meshes.For the calculation of the eigenvalues of ( 18), we use a PETSc routine [30,31] which can calculate up to 500 eigenvalues 6 .Have in mind that in contrast to DG and multi-block FD setting, the mass matrix in the pure Galerkin scheme is not block diagonal.Therefore, we can not split the eigenvalue calculation to each block matrix and have to consider the whole matrix M and therefore Q .Every coefficient of the numerical approximation belongs to one degree of freedom and we will obtain the same number of eigenvalues as number of DOFs are used.However, for the coefficients which belong to DOFs inside the domain, in all calculations we obtain zero up to machine precision.To get useful results, we decrease the number of elements in the following calculations and provide only the most negative and positive ones in tables 1-3 where we give results with and without the application of the SAT boundary operators.We see from tables 1-3 that the the boundary operator decreases the negative eigenvalues and forces the positive ones to zero (up to machine precision).For third and fourth order, we print only the case using Bernstein polynomials.The applications of Lagrange polynomials lead only to slightly bigger amounts of positive and negative eigenvalues of the Q + Q T operator (i.e maximum eigenvalue is 0.11713334374388217 for P 3).However, the results are similar after applying the SAT procedure, we obtain only negative or zero eigenvalues.We also mention that for higher degrees and more DOFs, we may strengthen the SAT terms to guarantee that the eigenvalues are negative and /or forced to zero and that we do not encounter stability issues because our procedure is too weak.All of our investigations are in accordance with the analysis done in subsection   4.1 and means that that everything is in order and consistent with the theory.All of our calculations demonstrate that a pure Galerkin scheme is stable if a proper boundary procedure is used.
Remark 5.1.Finally, we have done a couple of additional simulations changing both, the domain Ω (circles, pentagons, etc.) and the speed vector including also some horizontal movement.All of our calculations have remained stable if the boundary approach from section 4 was used .This is in contradiction of a common belief in the research community about continuous Galerkin schemes.But what are the reasons for this belief?In our opinion, one of the major issues is that the exactness of the quadrature rule is chosen too low and without artificial stabilization terms the continuous Galerkin scheme collapses, and the corresponding Q matrix does not become almost skew-symmetric.

Inexactness of the Quadrature Rule
To support our statement, we provide the following example.We consider the same problem as before, but in the Galerkin scheme we lower the accuracy of our quadrature rule so that the mass matrix is not evaluated exactly anymore (up to machine precision).We decrease the CFL number to 0.01 for stability reasons.However, as it is shown in figure 2 the scheme crashes after some time even with this super low CFL number.In pictures 2c , the structure of the bump can still be seen, but, simultaneously, the minimum value is ≈ −2.996 and the maximum value is around 2.7.Additional time steps will lead to a complete crash of the test.Here, again nothing has changed from the calculations before except that the quadrature rule is not exact anymore in the calculation of M , which leads to an erroneous Q matrix, which cannot be stabilized with the SAT boundary treatment.

Linear Rotation
In the next test we are considering an advection problem with variable coefficients.The speed vector has components a = 2πy, b = −2πx.
The problem is defined on the unit disk D = {(x, y) ∈ R 2 | x 2 + y 2 < 1}.The small bump rotates in the clockwise direction in a circle around zero.In figure 3a the initial state is presented where figure 3b shows the used mesh.We apply an unstructured triangular mesh with 932 triangles.In the second calculation 5382 triangles are used.The time integration is again done via a SSPRK54 scheme with CFL=0.2.A pure continuous Galerkin scheme with Bernstein polynomials is used for the space discretisation.The boundary operator is estimated via the approach presented in 4.1.3.In figure (4), the results are presented after two rotations of the bump.This test again verifies more that our scheme remains stable only through our boundary procedure.
We compute this problem up to ten rotations for different orders.We observe that all of our calculation remain stable both using Lagrange or Bernstein polynomials as can be seen for example in figure 5. Finally, we analyze the error behavior and calculate the order.In the first step, figure 6 shows the error behaviors for second-fourth order in space and in time after t = 0.01 with CFL=0.2.This shows that figure 6 represents mainly the space discretisation, i.e. the pure continuous Galerkin scheme.Here, we see that the space discretization has indeed the desired accuracy order.To include the time integration effects, we determine again the errors, but now after one rotation, see figure 7. We recognize a slight decrease of the order.The investigation of the decreased order of accuracy is not the main focus of this paper, where we focus on the stability properties of the pure continuous Galerkin scheme.In applications, we will consider the unsteady version of (47) ∂U ∂t + AU x + BU y = 0 with boundary conditions that will be detailed in the next part of this section.The aim is to look for a steady solution of this system, and hence to develop a time marching approach.With α ∈ R, the matrix cos αA + sin αB reads 2 , 0, 0).Through P , we can calculate P −1 , P 1/2 and P −1/2 without any problems.Remark 5.2.Since the system is symmetrizable, the eigenvectors are orthogonal for the quadratic form associated to P , i.e. for eigenvectors r i = r j hold P r i , r j = 0, where •, • denotes the scalar product.

The boundary conditions
The physical boundary condition follows from Maxwell's kinetic accommodation model.We have −αθ + s x n x + s y n y − αR nn βt x s x + βt y s y + t x n x R xx + (t x n y + t y n x )R xy + t y n y R yy weakly.Using this approach, we derive a suitable boundary operator from the continuous setting to guarantee that the discrete energy inequality is always fulfilled.We present a recipe how these operators can be constructed, in detail, for scalar two-dimensional problems, one-dimensional systems, and two-dimensional systems.In numerical experiments, we verify our theoretical analysis.Furthermore, in one test, we demonstrate the importance of the used quadrature rule.It has to be exact, if not, the Galerkin scheme suffers from stability issues.If stability can be reached only by enforcing the proper dissipative boundary conditions, there is no free meal: this procedure is very sensitive to any numerical error, like roundoff error or quadrature error.
We think and hope that through our results the common opinion about continuous Galerkin and its application in CFD problems changes, modulo the restriction we have already described.At least this result is interesting from the theoretical point of view, and emphasizes the positive role that the boundary conditions may have.In the companion paper [34], we consider and analyze non-linear (e.g.entropy) stability of continuous Galerkin schemes.Here, the SAT approach will also be important and some approximation for the boundary operators will be developed.

Figure 1 : 4 -
Figure 1: 4-th order scheme in space and time

Figure 2 : 4 -
Figure 2: 4-th order scheme in space and time

Figure 3 : 4 -
Figure 3: 4-th order scheme in space and time

= 2 Figure 9 :== R 1 ,
Figure 9: Results for the problem and t = 2, irregular mesh, 3rd order scheme in space and time