Application of Modal Filtering to a Spectral Difference Method
AApplication of modal filtering to a Spectral
Difference Method
Jan Glaubitz, Philipp ¨Offner, and Thomas Sonar2018
We adapt the spectral viscosity (SV) formulation implemented as a modal filter to a SpectralDifference Method (SD) solving hyperbolic conservation laws. In the SD Method we use selec-tions of different orthogonal polynomials (APK polynomials). Furthermore we obtain new errorbounds for filtered APK extensions of smooth functions. We demonstrate that the modal filteralso depends on the chosen polynomial basis in the SD Method. Spectral filtering stabilizes thescheme and leaves weaker oscillations. Hence, the selection of the family of orthogonal polyno-mials on triangles and their specific modal filter possesses a positive influence on the stabilityand accuracy of the SD Method. In the second part, we initiate a stability analysis for a linearscalar test case with periodic initial condition to find the best selection of APK polynomialsand their specific modal filter. To the best of our knowledge, this work is the first that gives astability analysis for a scheme with spectral filtering. Finally, we demonstrate the influence ofthe underlying basis of APK polynomials in a well-known test case.
MSC2010:
Keywords: hyperbolic conservation laws, high order methods, Spectral Difference Method, orthogonalpolynomials, modal filtering
In the field of Computational Fluid Dynamics (CFD), low-order methods are generally robust and reliableand therefore employed in practical calculations. The main advantage of high-order methods towards low-order ones is the possibility of considerably more accurate solutions with the same computing cost, butunfortunately they are less robust and more complicated. In recent years many researchers focus on thistopic. There has been a surge of research activities to improve and refine high-order methods as well asto develop new ones with more favorable properties. The Spectral Difference (SD) Method for simplexcells was first presented by Liu et al. [24], and later extended by Wang et al. [37, 46]. For one-dimensional,two-dimensional quadrilateral and three-dimensional hexahedral grids, the classical SD Method is identicalto the multi-domain staggered grid spectral element method proposed by Kopriva et al. [22, 23]. Furthercontributions can be found inter alia in [14,27,39]. There are many publications that show that the classicalSD Method is closely related to the Discontinuous Galerkin (DG) Method and Spectral Volume Method (orthe same), see [1, 15, 16, 18, 40, 45, 49] for example. Depending on how the degrees of freedom (DOFs) arechosen, various method implementations have different numerical properties and efficiencies. But all of themhave in common that they use piecewise continuous functions as approximation space for solutions. Here, weapply the classical orthogonal polynomials on triangles in a Spectral Difference Method and focus on filteringtechniques. We don’t use the nowadays more common Flux Reconstruction (FR) or Correction Procedurevia Reconstruction (CPR) methods, because we only want to investigate the pure influence of the orthogonalpolynomials and their natural filters. In the FR/CPR approach, a correction term is applied which works atthe interface and rectify the calculation in in every step, for details see [9, 15, 17, 34, 35, 42, 43] and referencestherein. Furthermore, the SD Method on triangular grids is not stable, i.e. numerical solutions migth exceedall boundaries where they shouldn’t. This allows us to oberve the influence of filter techniques and how theyenable us to milder such instabilities of the scheme. The SD Method combines the basic ideas of spectralmethods and finite differences. It directly reconstructs a flux polynomial based on fluxes on a given nodal With appropriate choices of correction terms, the FR framework recovers specific DG, SD, as well as spectral volume schemes. a r X i v : . [ m a t h . NA ] F e b et called flux points. Then the derivatives of the flux polynomial are used to update the solutions at thesolution points. This approach prevents the use of quadrature rules like in a normal DG or Spectral Volumeansatz . However, the SD Method has stability issues especially for high orders or if discontinuities arise inthe solution. The root of the instability is that the nonlinear flux function is represented by an insufficientamount of points. This introduces aliasing errors [20]. By increasing the dissipation we are able to masksuch an aliasing problem. We apply the same approach like in spectral methods and use spectral viscosity tostabilize the calculation, see [26,38]. As suggested in [25], the spectral viscosity can be carried out within thespectral filtering framework, resulting in an efficient computational implementation. The Spectral ViscosityMethod (SV) can be seen as a spectral method, but at each time step the numerical solution is filtered byan exponential filter which depends on the chosen set of orthogonal polynomials.Meister et al. [28] apply the Proriol-Koornwinder-Dubiner (PKD) polynomials in a DG Method and derivea relation between a modal filter for DG Methods on unstructured triangulation grids and the introductionof spectral viscosity to the scheme. The basic idea is to add a high order viscosity term, which is based onthe Sturm-Liouville Operator of the polynomials to the equation.In this article, we consider the Spectral Difference Method as described in [48], but we extend it by thegeneral classical orthogonal polynomials on triangles (APK polynomials) and their specific modal filters.With the differential operator of the APK polynomials, we show analogous to [28] the close relation betweenSV and spectral filtering in the SD Method.In a theoretical framework we prove some new upper bounds for the filtered APK expansion of smoothfunctions. Here, we generalize the theoretical results from [28], where only the properties of the filteredPKD expansion were analyzed. Our viscosity term depends on the differential operator of the chosenpolynomial set and so the exponential filter as well. Therefore the selection of the orthogonal polynomialsand their natural filter have a positive impact on the stability and accuracy of the method. By starting astability analysis as presented in [15, 39] we get a better understanding of the influence of the orthogonalpolynomials and their specific filters.This paper is organized as follows:In Section 2 we will define these considered polynomials and review some properties. We will repeat themain ideas on the Spectral Difference Method and explain our implementation in the next Section. Highorder filters as the SV Method as stabilizing technique are introduced in Section 4. Also, we prove an errorbound for the filtered APK expansion of smooth functions and transfer the SV modification to our SDMethod. In Section 5 we start with the stability analysis for our SD Method. A numerical experiment ispresented before a conclusion and an outlook for future work finishes this paper. In this Section we introduce the orthogonal polynomials under observation. We call the family of classicalorthogonal polynomials on triangles Appell-Proriol-Koornwinder polynomials (APK polynomials). In [30]the authors prove spectral convergence for the APK series. Their result gives us the theoretical foundationto use these polynomials in a spectral method.Let P α,βn ( x ) be the n -th Jacobi polynomial, T := { ( x, y ) ∈ R | x ≥ , y ≥ , x + y ≤ } be the unit triangleand h ( x, y ) := x α − y β − (1 − x − y ) γ − α − β ( α, β, γ ∈ N , γ > α + β − N = { , , · · · } ) be the weightfunction, given in this domain. For the sake of brevity, we introduce p = γ − α − β and a l = p + β + 2 l .Note that we only use α, β, γ ∈ N for simplicity. In principle, α, β, γ ∈ R +0 is possible.The polynomials A m,l ( x, y ), m, l ∈ N , defined as A m,l ( x, y ) := P α − ,a l m (1 − x ) P p,β − l (cid:18) y − x − (cid:19) (1 − x ) l (1)on T are called Appell-Proriol-Koornwinder polynomials (APK polynomials).If the triangle (cid:101) T := { ( x, y ) ∈ R | x ≥ − , y ≥ − , x + y ≤ } is used instead, then (1) transforms to (cid:101) A m,k ( x, y ) = P α − , l + γ − αm ( − x ) P γ − α − β,β − l (cid:18) y + 1)1 − x − (cid:19) (cid:18) − x (cid:19) l . The special case α = β = 1 and γ = 2 is called Proriol-Koornwinder-Dubiner polynomial (PKD polynomi-als)(see [6, 19]) . These are the classical orthogonal polynomials on triangles. Details of their properties There exists also quadrature-free implementations of DG, see for example [2]. The PKD or Dubiner polynomials are only one specific family of these polynomials. h ( x, y ) be the weight function on the triangle T . We denote by L ( T , h ) the Hilbert space with theinner product ( u, v ) := (cid:90) T h ( x, y ) u ( x, y ) v ( x, y ) dxdy which induces the weighted norm || u || L ( T ,h ) := (cid:18)(cid:90) T h ( x, y ) | u ( x, y ) | dxdy (cid:19) . We introduce by H m ( T , h ) a weighted Sobolev space . Precisely we set H m ( T , h ) := { v ∈ L ( T , h ):for each non-negative multi-index σ with | σ | ≤ m , the distributional derivative D σ v belongs to L ( T , h ) } . The space is endowed with the norm || v || H m ( T ,h ) := (cid:88) | σ |≤ m || D σ v || L ( T ,h ) . The APK polynomials are characterized by the integers m and l as can be seen in the definition. Thedegree of an APK polynomial then is m + l . Furthermore the APK polynomials are orthogonal in L ( T , h ),i.e. ( A m,l , A n,s ) = δ m,n δ l,s l + γ − α )(2( m + l ) + γ ) κ l,m with κ l,m := (cid:115) ( l + β ) p m ( m + a l ) α ( l + 1) p ( m ) α ( m + a l ) , (2)where ( ξ ) = 1 , ( ξ ) j = j (cid:89) i =0 ( ξ + i −
1) for j ∈ N is the usual Pochhammer symbol. See [36, Chapter X. p.288ff].Another important property of the APK polynomials can be found in the fact that they are solutions tothe eigenvalue problem DA m,l = λ m,l A m,l (3)where the differential operator D is given by( x − x ) ∂ ∂x + 2 xy ∂ ∂x∂y + ( y − y ) ∂ ∂y + [( γ + 1) x − α ] ∂∂x + [( γ + 1) y − β ] ∂∂y , and the eigenvalue λ m,l = ( m + l )( m + l + γ ), see [7, p.46]. Note that this eigenvalue equation is not aSturm-Liouville problem and also that the differential operator is not self-adjoint for all choices of α, β, γ .It is only self-adjoint in case of α = β = 1 and γ = 2. But in fact, the operator is, what is called, potentiallyself-adjoint in ˚ T , see [36, p.136]. Hence there exists a positive C - function g , ( x, y ) (cid:55)→ g ( x, y ), so that gD isself-adjoint in ˚ T . This will be important in Section 4, when we use the differential operator in the viscosityterm. In the following Lemma we summarize two estimates which will be needed in the proof of Theorem4.1. The proofs can be found in [30]. Lemma 2.1.
The following norm estimate holds for APK polynomials, || A m,l || L ( T ,h ) ≤ m + l + γ ) κ l,m (4)3 ith κ l,m given by (2).Furthermore let ( x, y ) ∈ ˚ T . For all l, m ∈ N the following estimate holds, | A m,l ( x, y ) | ≤ ˜ D ( x, y )(2 l + β + p ) (2( m + l ) + γ ) κ l,m , (5) where ˜ D ( x, y ) = D − x − y ) p + y + β − x + α − (1 − x ) and D < is a positive constant. The value of A m,l (1 , is | A m,l (1 , | = (cid:40)(cid:0) m + γ − αm (cid:1) if l = 00 else . (6) Remark 2.2.
Similar estimates also hold for the other edges of T and can be seen in [30]. In this paper we consider two-dimensional hyperbolic conservation laws of the form ∂∂t u ( x , t ) = −∇ x F ( u ( x , t )) , ( x , t ) ∈ Ω × R +0 , (7)where Ω ⊂ R is an open polygonal domain and u ( x , t ) ∈ R n . Furthermore, initial conditions u ( x ,
0) = u ( x ) and appropriate boundary conditions are assumed to be given.Hence, let T h be a conforming triangulation of the closure Ω of the computational domain and let P h bethe piecewise polynomial space defined by P h = { v h | τ i ∈ P N ( τ i ) , ∀ τ i ∈ T h } , where P N ( τ i ) denotes thespace of all polynomials on τ i of degree less than or equal to N .The classical Spectral Difference Method, which has been proposed by Liu [24] in 2006, can be seen asa FR or collocation method. The basic idea of this method is to discretize the right hand side of theunderlaying conservation law (7) at certain solution points x j in each cell. Then, the resulting ODEs (in t )at each x j can be solved by an arbitrary explicit time-stepping scheme. Here we used the 4-th order lowstorage Runge-Kutta scheme defined by Carpenter and Kennedy, see [5]. We say that a scheme is of order N + 1, when it is exact for u ∈ (cid:2) P N (cid:3) n . Since the derivate of the flux F is applied to update u , one needs tobe exact for F ∈ (cid:2) P N +1 (cid:3) n × . We approximate the flux F in each element τ i of our triangulation T h usingthe basis polynomials ϕ k in the following way (cid:18) F ( x , t ) F ( x , t ) (cid:19) = K F (cid:88) k =1 (cid:32) ˆ F k, ( t )ˆ F k, ( t ) (cid:33) ϕ k ( x ) . (8)The classical approach uses Lagrange polynomials L k with corresponding coefficients ˆ F k,ν ( t ) = F ν ( x k , t )for ν = 1 , x k are chosen flux points. If we want to achieve a method of order N + 1, we need thereconstruction of the solution u to lie in (cid:2) P N (cid:3) n and the reconstruction of the flux F to lie in (cid:2) P N +1 (cid:3) n × .Hence we need K s := ( N +1)( N +2)2 solution points x j and K F := ( N +2)( N +3)2 flux points x k . That meanspoints where the flux value is computed. Here we apply the classic orthogonal polynomials on triangles ϕ k = A m,l (in lexicographic order k ). In order to adapt the polynomial basis ϕ k to every triangular cell τ i , we introduce an orientation-preserving affine transformation T i which maps an arbitrary triangle τ i tothe standard triangle T of the orthogonal basis ϕ k , see Figure 1. By linearity of the differential operator,inserting (8) into (7) leads to u t ( x , t ) = − K F (cid:88) k =1 (cid:32) ˆ F k, ( t )ˆ F k, ( t ) (cid:33) · ∇ x ϕ k ( T i ( x )) . Applying the transformation and the chain rule to ∇ x ϕ k ( T i ( x )) = J T i ∇ ξ ϕ k ( T i ( x )) , igure 1: transformation of an arbitrary triangle to the standard element finally leads to the universal update scheme u t ( x , t ) = − K F (cid:88) k =1 (cid:32) ˆ F k, ( t )ˆ F k, ( t ) (cid:33) · J T i ∇ ξ ϕ k ( T i ( x )) . (9)where for each cell only the Jacobian J T i has to be stored. For the stability analysis in Section 5 we willderive benefit from the following matrix representation of (9). If we denote the vector (cid:0) u t ( x j , t ) (cid:1) j at the K s solution points x j by d u s d t ( t ) and J T i by (cid:18) ξ x ξ y η x η y (cid:19) , the universal update scheme readsd u s d t ( t ) = − (cid:16) K F (cid:88) k =1 ˆ F k, ( t ) (cid:2) ξ x ∂ ξ ϕ k ( T i ( x j )) + ξ y ∂ η ϕ k ( T i ( x j )) (cid:3)(cid:17) j − (cid:16) K F (cid:88) k =1 ˆ F k, ( t ) (cid:2) η x ∂ ξ ϕ k ( T i ( x j )) + η y ∂ η ϕ k ( T i ( x j )) (cid:3)(cid:17) j = − ξ x D ξ ˆ F ( t ) − ξ y D η ˆ F ( t ) − η x D ξ ˆ F ( t ) − η y D η ˆ F ( t ) (10)with D ξ = (cid:0) ∂ ξ ϕ k ( T i ( x j )) (cid:1) j,k , ˆ F ( t ) = (cid:0) ˆ F k, ( t ) (cid:1) k and D η = (cid:0) ∂ η ϕ k ( T i ( x j )) (cid:1) j,k , ˆ F ( t ) = (cid:0) ˆ F k, ( t ) (cid:1) k . (11)In spite of the Lagrange reconstruction approach where the fluxcoefficients are directly known we have to compute ˆ F k,ν ( t ), inevery time step. We choose an interpolation approach, i.e. weinterpolate the values of the flux function at certain flux points { x k } from the given flux F , which will be specified later on.The nodal set { x k } is chosen as the set of two dimensional Lo-batto points on a triangle, see Figure 3, proposed by Blyth andPozrikidis [4] since they are easy to implement and have goodinterpolation properties, as for instance a low condition number.The condition number of the Vandermonde matrix V in the interpolation approach depends on the param-eters α, β and γ of the APK polynomials, see Table 1. The numbers in the braces describe the parameters( α, β, γ ) of the APK-polynomials. We realize that Lagrange polynomials lead to completely bad conditionedbasis compared to any of the considered APK families. Regarding this, one should note that the coefficientsof the Lagrange reconstruction are directly given by the values of the flux at the flux points. This valueshowever are obtained from the values of the flux at the solution points by Lagrange interpolation. Additionalnumerical errors might arise.Additionally, enough flux points lie on the edges of an element to ensure global conservation [47], which isthen realized by replacing the original flux F at flux points x k at the edges by a numerical flux F num whose5 · · · · (0 . , . , .
2) 72 23 103 162 305 315 429 577 681(0 . , . ,
1) 12 17 30 42 52 71 82 100 121(1 , ,
2) 12 20 39 53 71 94 121 151 196(1 , ,
6) 29 117 391 973 2365 4484 8008 10 · (2 , ,
5) 17 44 96 161 244 450 662 853 1341(3 , ,
20) 28 131 416 929 2310 4711 10 · · (10 , ,
20) 306 2307 2 · · · · · · · (1 , ,
3) 11 33 66 99 147 209 295 409 532(2 , ,
3) 17 34 66 123 178 243 388 446 570
Table 1: condition numbers κ N normal component is computed from a numerical flux function H , i.e. F numk · n = H ( u − ( x k , t ) , u + ( x k , t ) , n ).The whole flux is used in the scheme, so that one needs a second condition to determine the numericalflux. As in [44], we enforce F num to maintain the same tangential component as the original flux, i.e. F numk · t = F ( u ( x k , t )) · t , which uniquely defines the numerical flux F numk = ( F numk, , F numk, ) T at each edgepoint x k . Now, the coefficients ˆ F k,ν ( t ) in equation (9) can be computed by solving the system of equation F ν ( t ) = V · ˆ F ν ( t ) for each component ν = 1 ,
2, where ˆ F ν = (cid:16) ˆ F k,ν ( t ) (cid:17) k , V = (cid:0) ϕ k ( x j ) (cid:1) j,k is the Vandermondematrix and ˆ F ν ( t ) = ( F k,ν ) k is given by F k,ν = (cid:40) F ν ( u ( x k , t )) , if x k ∈ ˚ τ i F numk,ν , if x k ∈ ∂τ i . (12)If we apply this to (10) the universal update scheme readsd u s d t ( t ) = − ( ξ x D ξ + ξ y D η ) V − F ( t ) − ( η x D ξ + η y D η ) V − F ( t ) (13)in a matrix representation. In Section 5 we will use this result to obtain the formd u s d t ( t ) = S u s ( t )with a matrix S for a suitable test case. In the numerical examples in Section 6 we will use the Godunovflux H ( u − , u + , n ) = min u − ≤ u ≤ u + F ( u ) · n , if u − ≤ u + , max u + ≤ u ≤ u − F ( u ) · n , else , in the scalar case for the numerical flux function. The Godunov flux is an Upwind flux. For the lineartransport equation we obtain the numerical flux H ( u − , u + , n ) = (cid:40) ( a · n ) u − , ( a · n ) ≥ , ( a · n ) u + , ( a · n ) < . We will apply this numerical flux function in the stability analysis in Section 5 as well. For a more detailedexplanation of the Spectral Difference Method and the extension we advise [47]. Yet, to at least highlight asignificant advantage of the Spectral Difference Method, we observe the experimental order of convergence(EOC) in the L ∞ -norm for the linear advection equation ∂∂t u ( x , t ) = − ∂∂x u ( x , t ) − ∂∂y u ( x , t ) with ( x , t ) ∈ [ − , × [0 , . u ( x ,
0) = sin π ( x + y ). In Table 2, this is done with respect to both thenumber of triangles (EOC(k)) and the degree of the polynomial approximation in each triangle (EOC(N)).For these approximations the parameters ( α, β, γ ) = (2 , ,
5) were chosen. Note that for unstructured6
68 272 1088 4352 h h for the size of the triangles should be used for the rate of convergence.Here for example, h was the biggest volume of the triangles. At the same time however, in this test, h wasdirectly inversely proportional to k , i.e. h = const · k .Since EOC(h) is equals to EOC(k). For the sake of simplicity, we therefore observed convergence in thenumber of triangles k and not their size h . N k L ∞ -error EOC(k) EOC(N) L -error L -error time1 68 1.226072e+00 5.234411e-01 3.483826e-01 82 68 2.999127e-01 2.03 1.106267e-01 8.451124e-02 153 68 6.865602e-02 3.63 1.461619e-02 1.290576e-02 284 68 1.358577e-02 5.63 1.654693e-03 1.770760e-03 485 68 1.583057e-03 9.63 1.353175e-04 1.657258e-04 781 272 4.154960e-01 0.78 1.638561e-01 1.079845e-01 332 272 4.481682e-02 1.37 3.21 1.738056e-02 1.300429e-02 623 272 5.600234e-03 1.80 5.12 1.203168e-03 1.050415e-03 1104 272 5.135337e-04 2.36 8.30 6.683444e-05 6.523343e-05 1905 272 2.973202e-05 2.86 12.76 3.036868e-06 3.352791e-06 3121 1088 1.141089e-01 0.93 4.514174e-02 2.951796e-02 1292 1088 6.801320e-03 1.36 4.06 2.623593e-03 1.942249e-03 2433 1088 4.438368e-04 1.82 6.73 1.007609e-04 8.647829e-05 4374 1088 1.793560e-05 2.41 11.15 3.016446e-06 2.706330e-06 7615 1088 5.468877e-07 2.88 15.64 8.022185e-08 7.697035e-08 12471 4352 3.105485e-02 0.93 1.174054e-02 7.637636e-03 5262 4352 1.458704e-03 1.11 4.41 4.029470e-04 3.134510e-04 8603 4352 4.706511e-05 1.61 8.46 8.869529e-06 7.971082e-06 17664 4352 1.320516e-05 0.22 4.41 2.170014e-07 3.488747e-07 30795 4352 2.180251e-06 -0.99 8.07 1.510489e-08 4.853751e-08 5144 Table 2: ( α, β, γ ) = (2 , , , t = 0 . s Table 2 clearly indicates the rate of convergence to be considerably higher when increasing the polynomialdegree instead of refining the triangulation. For sufficiently smooth solutions, schemes using a polynomialapproximation, in fact, often provide significant higher rates of convergence than classical ones, where lowdegrees (in particular 0) are used. Note that in Table 2, the same behavior can be observed for for the L -and L -norm. In the conservation law, discontinuities may arise in the solution. Using a series expansion P N u ( x, y ) = (cid:88) l + m ≤ Nl,m ∈ N ˜ u m,l A m,l ( x, y ); ˜ u m,l = ( u ; A m,l ) L ( T ,h ) ( A m,l ; A m,l ) L ( T ,h ) (14)to approximate the solution leads to spurious oscillations in the vicinity of discontinuities (called Gibbsphenomenon). The oscillations occur in the approximated solution because the high coefficients of the seriesexpansion turns slowly to zero, see [11, 12].The aliasing error effects stability problems in the SD Method. The root of the instabilities is that thenonlinear flux function is represented by an insufficient amount of points. This introduces aliasing errors.These stability problems display especially near the oscillations. The amplitudes increase in time exceedingthe error expected from the pure Gibbs phenomenon, see Section 6. A higher inherent dissipation is able tomask such an aliasing problem. By adding spectral viscosity to the equation we increase the dissipation.To remedy the Gibbs’ effect, different approaches can be found in the literature [8, 10, 12, 41]. A commonapproach is to use a modal filter, which appeals directly on the high-order coefficients of the series expansion.While multiplying a filter function to the high Fourier coefficients, the series loses their approximation7roperties. We prove new error bounds for the filtered APK series expansion of smooth function. Wegeneralize the results of [28] to all series expansion with classic orthogonal polynomials on triangles. Globalfiltering in each cell degrades the order of accuracy, see [28]. Therefore we apply a well-known jumpindicator [33] and transfer the SV modification to the SD Method. A modal filter appeals directly to the coefficients of the series expansion. For an integer p ≥ filter of order p as a real function σ ∈ C p − ([0 , σ (0) = 1 ,σ ( k ) (0) = 0 , ≤ k ≤ p − . (15)Additionally to the properties (15) many authors demand the following condition for a filter of order p : σ ( k ) (1) = 0 0 ≤ k ≤ p − , (16)compare [12, 41]. In [13], Hesthaven and Kirby require in addition σ ∈ C p ([0 , σ ( η ) = exp( − αη p ) , (17)where the filter strength α yields exp( − α ) in the range of the machine accuracy. We already mentionedthat there are various partial results for the approximation property of a filtered series expansion. Forinstance, in [41] Vandeven has analyzed the filtered Fourier expansions of a piecewise smooth function andlater expanded his result to filtered Chebyshev series. On the other hand, Kirby and Hesthaven investigatethe approximation property of the filtered Legendre expansion for sufficiently smooth functions in [13]. [28]expands the investigation to two-dimensional basis functions (PKD polynomials). Finally, we now completethe investigation from [28], by proving an error bound for the filtered APK partial sum. Again, we stressthat the PKD polynomials are just a special case of the APK polynomials. Here we prove the most generalcase and consider the filtered APK expansions for sufficiently smooth functions u : T → R , u σN ( x, y ) = (cid:88) l + m ≤ N σ (cid:18) l + mN (cid:19) ˜ u m,l A m,l ( x, y ) , N ≥ , with coefficients ˜ u m,l given by (14). Theorem 4.1.
Let u ∈ H k ( T , h ) ∩ C ( T ) , k ∈ N , h ( x, y ) = x α − y β − (1 − x − y ) p with α, β ∈ N and p ∈ N and let σ be a modal filter of order k − , with the additional condition σ ∈ C k − ([0 , ε )) in an interval [0 , ε ) ⊂ [0 , , ε > . Furthermore let k > max (cid:110) + α + p , + p + β , + α + β (cid:111) .Then we obtain the pointwise error bounds with constants K − K :1. If ( x, y ) ∈ ˚ T , it is | u ( x, y ) − u σN ( x, y ) | ≤ K N k − .
2. On the left edge [0 , y ] with y ∈ [0 , we obtain | u (0 , y ) − u σN (0 , y ) | ≤ K N k − α − p − , for p > β − ,K N k − α − β + 12 , for p ≤ β − .
3. On the edge [ x, with x ∈ (0 , the error is | u ( x, − u σN ( x, | ≤ K N k − − β .
4. On the hypotenuse [ x, − x ] with x ∈ (0 , pertain | u ( x, − x ) − u σN ( x, − x ) | ≤ K N k − − p . . In the point (1 , we preserve | u (1 , − u σN (1 , | ≤ K N k + α − γ − , Proof.
We start the proof similar to [30, Theorem 3.1]. Let m + l (cid:54) = 0, then( A m,l ; A m,l ) L ( T ,h ) · | ˜ u m,l | = ( u ; A m,l ) L ( T ,h ) (3) = (cid:32) u ; DA m,l λ m,l (cid:33) L ( T ,h ) holds. We use the fact that the differential operator D is potentially self-adjoint in ˚ T and that the boundaryis a set of measure zero. Due to the operator D to be potentially self-adjoint, we find a positive C -function g so that gD is self-adjoint, g is well defined and symmetric. So we get (cid:32) u ; DA m,l λ m,l (cid:33) L ( T ,h ) = (cid:32) u ; ( gD ) A m,l g · λ m,l (cid:33) L ( T ,h ) = 1 λ m,l (cid:18) g ( gD ) u ; A m,l (cid:19) L ( T ,h )recursive = (cid:32) λ m,l (cid:33) k (cid:16) D k u ; A m,l (cid:17) L ( T ,h ) and hence have shown ( A m,l ; A m,l ) L ( T ,h ) · | ˜ u m,l | = (cid:32) λ m,l (cid:33) k (cid:16) D k u ; A m,l (cid:17) L ( T ,h ) . (18)By this and u ( x, y ) = (cid:88) m,l ∈ N ˜ u m,l A m,l ( x, y )is satisfied for every point ( x, y ) ∈ T , see [30, equation (18)], we get | u ( x, y ) − u σN ( x, y ) | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) ≤ l + m ≤ Nl,m ∈ N (cid:32) − σ (cid:18) l + mN (cid:19)(cid:33) ˜ u m,l A m,l ( x, y ) + (cid:88) l + m>Nl,m ∈ N ˜ u m,l A m,l ( x, y ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:90) T h ( x , y ) (cid:2) S N ( x, y, x , y ) + R N ( x, y, x , y ) (cid:3) D k u ( x , y ) d x d y (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , for the function u , where S N ( x, y, x , y ) and R N ( x, y, x , y ) are defined by S N ( x, y, x , y ) = (cid:88) ≤ l + m ≤ Nl,m ∈ N (cid:32) − σ (cid:18) l + mN (cid:19)(cid:33) A m,l ( x , y ) A m,l ( x, y ) || A m,l || L ( T ,h ) λ km,l ,R N ( x, y, x , y ) = (cid:88) l + m>Nl,m ∈ N A m,l ( x , y ) A m,l ( x, y ) || A m,l || L ( T ,h ) λ m,l . We apply the Cauchy-Schwarz inequality to estimate the integral. Therefore we have to estimate the weighted L ( T )-norms of S N and R N for all ( x, y ). For the R N we extract the results from the proof of Theorem 3.1from [30]. We summarize the results in Table 3. Now we are able to estimate the S N . We follow the proofof [30, Theorem 3.1] and start with some inner point ( x, y ) ∈ ˚ T . It is || S N ( x, y, · , · ) || L ( T ,h ) = (cid:88) ≤ l + m ≤ Nl,m ∈ N (cid:32) − σ (cid:18) l + mN (cid:19)(cid:33) A m,l ( x, y ) || A m,l || L ( T ,h ) λ km,l (4)&(5) ≤ (cid:88) ≤ l + m ≤ Nl,m ∈ N (cid:32) − σ (cid:18) l + mN (cid:19)(cid:33) κ m,l m + l + γ ) ˜ E ( x, y ) κ l,m (2 l + β + p ) ( m + l ) k ( m + l + γ ) k { (cid:0) l + pl (cid:1) , (cid:0) l + β − l (cid:1) } || R N ( x, y, · , · ) || L ( T ,h ) < ˚ T C R N − k + [0 , y ] (cid:0) l + pl (cid:1) C R N − k + α + p + [0 , y ] (cid:0) l + β − l (cid:1) C R N − k + α + β − [ x, C R N − k + + β [ x, − x ] C R N − k + + p (1 , C R N − k − α + γ + Table 3: C R , . . . , C R ∈ R + are constants. ≤ E ( x, y ) (cid:88) ≤ l + m ≤ Nl,m ∈ N (cid:32) − σ (cid:18) l + mN (cid:19)(cid:33) m + l ) k + ( m + l + γ ) k − Using the identity (cid:88) 0) we calculate it directly. || S N (1 , , · , · ) || L ( T ,h ) = N (cid:88) m =1 (cid:32) − σ (cid:18) mN (cid:19)(cid:33) A m, (1 , || A m, || L ( T ,h ) λ km, < (cid:32) ( β ) p γ α − (1) p (cid:33) N (cid:88) m =1 (cid:32) − σ (cid:18) mN (cid:19)(cid:33) (cid:0) m + γ − αm (cid:1) m k ( m + γ ) k − < C N (cid:88) m =1 (cid:32) − σ (cid:18) mN (cid:19)(cid:33) m k ( m + γ ) k − − γ +2 α N k − − γ +2 α N k − − γ +2 α < C S N − k +3+2 γ − α N N (cid:88) m =1 (cid:32) − σ (cid:18) mN (cid:19)(cid:33) (cid:18) mN (cid:19) − k +2+2 γ − α , and after all it ensues | u (1 , − u σN (1 , | < K N − k − α + γ + . To verify Theorem 4.1 we test the approximation speed. On the triangle T , we approximate the function f ( x, y ) = sin (cid:0) π ( x + y ) (cid:1) by the filtered and unfiltered truncated series expansion with respect to the APKpolynomials for parameters ( α, β, γ ) = (1 , , 2) and (2 , , . 001 and of order 4 and strength 0 . σ ( η ) = 0 . (cid:0) πη ) (cid:1) of second order. In Figure 2 the maximal errors of these truncated series expansions are illustrated withrespect to the polynomial degree . A and E are the unfiltered APK-expansions, B, F are the filtered APK-expansions, where the exponential filter of order 2 was used, for C, G the cosine-filter was used and for D, H the exponential filter of order 4 was used. Therefore, A is comparable to E , B to F , C to G and D to H . All calculations were done in Mathematica and by its integration subroutine in order to calculate thecoefficients ˆ u m,l .As Figure 2 illustrates, the highest rate of convergence is obtained by the unfiltered APK expansions.On the other hand, the APK expansions filtered by the cosine-filter (C,G) and the exponential filter ofsecond order (B,F) show a notably slower rater of convergence. Yet, by increasing the order of the filterand applying the exponential filter of fourth order to the APK expansion (D,H), the rate of convergencesignificantly increases. Thus, for higher filter orders we get a stronger decrease in the max-error. In Table 4the position and the value of the maximum-errors are plotted for C and G . For the others it looks analogous.The maximum error lies in ˚ T and the error decreases with N . . For C K . G K . K , K , K , K K D .All of this verifies our result. Since spectral methods are known to lack sufficient dissipation, Tadmor [38] proposed the Spectral Viscosityor Super-Spectral Viscosity Method (SV Method). The main idea of the SV Method is to add a smallviscosity term to the conservation law (7). We show analogously to [25,32] that by choosing a viscosity term11 − − − − − m a x e rr o r ABCDEFGH Figure 2: A - D parameters ( α, β, γ ) = (1 , , 2) and E - H parameters (2 , , N ( x, y ) max- error ( C -case ) ( x, y ) max- error ( G )-case1 (0.331, 0.169) 0.36338 (0.331, 0179) 0.2668692 (0.269, 0.273) 0.303294 (0.198, 0.345) 0.1920263 (0.118, 0.435) 0.219701 (0.129, 0.442) 0.1329154 (0.162, 0.391) 0.149895 (0.219, 0.353 ) 0.0894925 (0.261, 0.293) 0.106279 ( 0.284, 0.289 ) 0.06292296 (0.263, 0.291) 0.0781872 (0.258, 0.316 ) 0.04607197 (0.197, 0.357) 0.0595273 (0.196, 0.378) 0.02735618 (0.220, 0.334) 0.0466554 (0.137, 0.437 ) 0.0219375 Table 4: C and G K K K K K . 496 0 . 741 0 . 975 0 . 496 0 . Table 5: K − K which depends on the differential operator of the orthogonal basis, the SV Method can be seen as a spectralmethod with modal filtering.For the SD scheme based on the APK expansions, we propose to consider high order operators of theform ( − D ) p on the reference element T . Let I N be an index set and ϕ k = A m,l the APK polynomials inlexicographic order k , see Section 3. Theorem 4.2. Let N ∈ N and { ϕ k | k ∈ I N } the APK polynomials on the triangle T . The polynomials solvethe eigenvalue problem − Dϕ k = − λ k ϕ k , (19) where D is the differential operator of the polynomials. P ϕ denotes the projection on the space of { ϕ k } . Tosolve the viscosity equation ∂∂t u N ( x , t ) + ∇ x · P ϕ f ( u N ( x , t )) = ε N ( − p +1 ( − D ) p u N ( x , t ) (20) by a splitting method is equivalent to multiply the coefficients ˆ u k with the function σ ( k ) = e − ε N ∆ tλ pk in every update step of the equation ∂∂t u N ( x , t ) + ∇ x · P ϕ f ( u N ( x , t )) = 0 . roof. We solve the equation (20) by a splitting method in two steps ∂∂t u N ( x , t ) = ε N ( − p +1 ( − D ) p u N ( x , t ) (21)and ∂∂t u N ( x , t ) + ∇ x · P ϕ f ( u N ( x , t ) = 0 (22)With u N ( x , t ) = (cid:80) k ∈ I N ˆ u k ϕ k ( x ) equation (21) implies (cid:88) k ∈ I N ∂ ˆ u k ∂t ϕ k ( x ) = (cid:88) k ∈ I N ε N ( − p +1 ˆ u k ( − D ) p ϕ k ( x ) = (cid:88) k ∈ I N − ε N ˆ u k ( t ) λ pk ϕ k ( x ) , where we applied the eigenvalue equation (19) in the last step. By comparing the coefficients we have tosolve the ordinary differential equations ∂ ˆ u k ( t ) ∂t = − ε N ˆ u k ( t ) λ pk , ∀ k ∈ I N . The solution is ˆ u k ( t ) = ce − ε N λ pk t , c ∈ R . With ∆ t := t n +1 − t n and the requirement ˆ u k ( t n +1 ) = ˆ u k ( t n ) for∆ t = 0 follows ˆ u k ( t n +1 ) = e − ε N λ pk (∆ t + t n ) = e ( − ε N λ pk ∆ t ) (cid:124) (cid:123)(cid:122) (cid:125) =: σ ( k ) ˆ u k ( t n ) . In order of us to speak of a modal filter for σ ( k ) = σ (( l, m )) = e − ε N ( l + m ) p ( l + m + γ ) p ∆ t , we have to multiplythe exponenet with N p , meaning σ (cid:18) l + mN (cid:19) = e − ε N N p ( l + mN ) p ( l + m + γN ) p ∆ t ≈ e − ε N N p ∆ t ( l + mN ) p . (23)Hence σ : [0 , → [0 , 1] can be seen as an exponential filter of order 2 p with filter strength α i := − ε N N p ∆ t .We consider (23) in detail and realize that our modal filter depends on the parameter γ of the APKpolynomials, especially if the polynomials have minor degree. So for different families of APK polynomialswe get various specific modal filters.We can prove similar results to Theorem 4.2 for every orthogonal basis { ϕ k } , if the ϕ k fulfil a comparableeigenvalue equation.For the transfer of the SV Method to our Spectral Difference Method we follow the steps according to [47].The transformation T i from τ i to T has no effect on the flux function, so that the SD update scheme of ascalar hyperbolic equation in the cell τ i with viscosity term and the solution points ξ j ∈ T reads like ∂∂t u N ( T − i ( ξ j , t )) + ∇ ξ · ˜ P N ( J T i ) T ˜ F ( u N ( T − i ( ξ j ) , t ))= ε N ( − p +1 ( − D ) p ( u N ( T − i ( ξ j ) , t )) , (24)where ˜ F is the flux function, ˜ P is the projection and ∇ ξ is the nabla operator on T . The same approach asin Theorem 4.2 yields the equivalent of (24) and the SD Method with the modal filter (23).Hence, the matrix representation (13) for the SD update scheme in cell τ i becomesd u s d t ( t ) = − ( ξ x D ξ + ξ y D η ) (cid:0) σ ( k ) ˆ F k, ( t ) (cid:1) k − ( η x D ξ + η y D η ) (cid:0) σ ( k ) ˆ F k, ( t ) (cid:1) k = − ( ξ x D ξ + ξ y D η ) M σ ˆ F ( t ) − ( η x D ξ + η y D η ) M σ ˆ F ( t )= − ( ξ x D ξ + ξ y D η ) M σ V − F ( t ) (25) − ( η x D ξ + η y D η ) M σ V − F ( t )13ith viscosity, where M σ := diag (cid:0) σ (1) , . . . , σ ( N F ) (cid:1) is a diagonal matrix. So we have to choose the size of ε N , thereby we stick on the analysis of [10]. The authors suggested in the Fourier case a viscosity strength ε N , which depends on the approximation order as ε N ∼ C P N p − , where the constant C P may be chosen by C p ≤ p (cid:88) k =1 || ∂ ku ( J T i ) T ˜ F ( u ) || L ∞ = || ( J T i ) T || ∞ p (cid:88) k =1 || ∂ ku ˜ F ( u ) || L ∞ . Here, we ignore the dependence on u and approximate || ( J T i ) T || ∞ ∼ h i by a measure of length h i of thetriangle τ i . Finally we arrive at ε iN := ch i N p − < p (cid:80) k =1 || ∂ ku ˜ F ( u ) || L ∞ h i N p − , (26)where c ∈ R is a constant. We have to select the constant c and the filter order 2 p in our numerical tests. Remark 4.3. The numerical tests from [29, 47] imply that we have to increase c if we enlarge the filterorder 2 p .Furthermore, it is a known fact that the application of a modal filter in the global domain will destroy theorder of convergence. Shock indicators are used to detect corrupted cells. There are different approachesto detect jumps in a cell, see [31, 48] for example. Here we use the shock indicator of [3]. The indicatorcompares the higher coefficients with the lower ones to detect oscillations.In general, the time step ∆ t in the explicit time-stepping scheme for the resulting ODEs (in t ) ateach solution point x j , which also appears in the filter strength α i = − ε N N p ∆ t in (23), can’t be chosenarbitrarily. For a scalar conservation law u t ( x, t ) + λu x ( x, t ) = 0 in one space dimension and a numericalscheme on a uniform grid with length h one should choose ∆ t small enough, so that a wave with propagationspeed λ can’t travel more than h in one time step. I.e. the Courant-Friedrichs-Levy condition λ ∆ t ≤ h hasto hold. The maximum max λ λ ∆ th =: CFLis called Courant number and the following example gives a first intuition how it can determined on atriangulation. Example 4.4. We observe the two dimensional advection equation u t ( x, y, t ) + u x ( x, y, t ) + u y ( x, y, t ) = 0 . Here the propagation speed in (1 , 1) direction is √ 2, which leads to the Courant number CFL = √ th , inwhich h is now a measure of length of the observed triangle.As we can see in Section 5.2 the SD Method has stability problems for all parameters α, β and γ , whichare increasing for rising order N . Hence, we can’t even compute a sufficiently numerical Courant numberand will choose the time step ∆ t := C fix ( N + 1) · hλ max (27)with a fixed value C fix and maximal propagation speed λ max . Therefore, the time step adapts an increasinginstability for rising N , the ”geometry” of the triangulation by h and by C fix = , which will be our choice,(27) coincides with the Courant number in example 4.4 for N = 0. We already mentioned, that we use the 4-th order low storage Runge-Kutta scheme defined by Carpenter and Kennedy(see [5]) in our implementation for Section 6. Stability of the Spectral Difference Method By von Neumann analysis in matrix form for a scalar advection equation in two spatial dimensions withperiodic initial condition, one can observe the linear stability of the semi-discretization. This was alreadydone by Abeele, Lacor and Wang for the classical SD Method in [39], by Huynh for the FR Method in [15]and many other authors. Here we observe the SD Method for APK polynomials on a triangular grid, whichis generated by a pattern like in Figure 3 for a scalar advection equation x axisy axis(0,0)(0, 1)(0,1) (1,0)( 1,0) Figure 3: Generating pattern u t + cos( ψ ) u x + sin( ψ ) u y = 0with ψ ∈ [0 , π ] and the periodic initial condition u init ( x, y ) = e I ( w x x + w y y ) , where w x , w y ∈ [ − π, π ] are the so called wave numbers in x - and y -direction. The exact solution is thengiven by u exact ( x, y, t ) = e I ( w x x + w y y − [ w x cos ψ + w y sin ψ ] t ) (28)for t ≥ 0. We will concentrate on the semi discretization at t = 0 and wish to obtain a relation in form of alinear ODE system d u s d t = S u s (29)of first order and with constant coefficients in the solution points. Similar to the scalar case d u d t = λu , wherethe solution is given by u ( t ) = ce λt for c ∈ R , the eigenvalues of the matrix S take the place of λ . Hence theasymptotic behavior of the solution of (29) depends on the real parts of the eigenvalues: If all eigenvaluesof S lie in the left half of the complex plane, the solution is bounded, which coincides with the behavior ofthe exact solution u exact . In this case we call the semi discretization to be stable . On the other hand, if15ny eigenvalue lies in the right half of the complex plane, the solution of (29) will blow up for t → ∞ andso might the scheme. This will happen unless the instability is mild, i.e. max λ Re( λ ) is small, and the timesteps adapt to this and are sufficiently small. In 5 . t = 0 and for the test case above. In 5 . α, β and γ . By observation of the numerical test cases, one cansee that some stabilization is necessary and we did so in Section 4 by the SV Method in terms of the SDMethod with modal filtering. Following this, in 5 . . p and α, β, γ from the APK polynomials. We already stated in Section 3 that the universal update scheme (9) at the cell τ i,j readsd u si,j ( t )d t = − ( ξ x D ξ + ξ y D η ) V − F ( t ) − ( η x D ξ + η y D η ) V − F ( t ) , (30)when evaluated at the solution points. If we look at the cell τ , and denote d u si,j ( t =0)d t for sake of simplicityby d u s d t equation (30) becomes d u s d t = − ( D ξ V − F + D η V − F ) , (31)since τ , and its neighbors are parametrized by the orientation-preserving affine transformations T , ( x, y ) = (cid:18) xy (cid:19) , T − , ( x, y ) = (cid:18) y − x (cid:19) , T , − ( x, y ) = (cid:18) − yx (cid:19) , T , ( x, y ) = (cid:18) − x − y (cid:19) . Note that F , F are given by (12), where the upwind flux is used and the flow direction was restricted to ψ ∈ [0 , π ]. Consequently F , F are independent of the values from τ , and so (31) now readsd u s d t = − cos( ψ ) D ξ V − (cid:2) M − , u f − , + M , u f , + M , − u f , − (cid:3) − sin( ψ ) D η V − (cid:2) M − , u f − , + M , u f , + M , − u f , − (cid:3) . Here the M m,ni for the second order SD Method, and so K F = (2+1)(2+2)2 = 6, are M − , = e T e + e T e + e T e , M , = e T e + e T e , M , − = e T e ,M − , = e T e , M , = e T e + e T e , M , − = e T e + e T e + e T e . For the exact computation see [47, chapter 3 . . 1] and the M m,ni for N = 3 , , u f − , and u f , − by substitutingthem by terms T − , u f , and T , − u f , . Therefore, note that u f − , , u f , and u f , − are given by the initialcondition. Hence u f , =( u init ( T − , ( ξ k , η k ))) k = ( u init ( ξ k , η k )) k =( e I ( w x ξ k + w y η k ) ) k , u f − , =( u init ( T − − , ( ξ k , η k ))) k = ( u init ( − η k , ξ k )) k =( e I ( − w x η k + w y ξ k ) ) k = ( e I ( − w x ( ξ k + η k )+ w y ( ξ k − η k )) · e I ( w x ξ k + w y η k ) ) k = T − , u f , , u f , − =( u init ( T − , − ( ξ k , η k ))) k = ( u init ( η k , − ξ k )) k =( e I ( w x η k − w y ξ k ) ) k = ( e I ( w x ( η k − ξ k ) − w y ( ξ k + η k )) · e I ( w x ξ k + w y η k ) ) k = T , − u f , with T − , := (cid:88) k e I ( − w x ( ξ k + η k )+ w y ( ξ k − η k )) e Tk e k , , − := (cid:88) k e I ( w x ( η k − ξ k ) − w y ( ξ k + η k )) e Tk e k , and we obtain d u s d t = − (cid:16) cos( ψ ) D ξ V − (cid:2) M − , T − , + M , + M , − T , − (cid:3) + sin( ψ ) D η V − (cid:2) M − , T − , + M , + M , − T , − (cid:3)(cid:17) u f , . (32)Last, we want to write u f , in terms of u s , or just u s to obtain a relation as in (29). While we reconstructedthe flux F by APK polynomials, in order to later use their natural filter given by (23), we reconstruct thesolution u by Lagrange polynomials. Hence u f , = (cid:0) K s (cid:88) k =1 u , ( x sk ) L k ( T , ( x fj )) (cid:1) K F j = E Lag u s with E Lag := (cid:0) L k ( x fj ) (cid:1) j,k ∈ R K F × K s . If we use that in (32) we finally obtaind u s d t = − (cid:16) cos( ψ ) D ξ V − (cid:2) M − , T − , + M , + M , − T , − (cid:3) + sin( ψ ) D η V − (cid:2) M − , T − , + M , + M , − T , − (cid:3)(cid:17) E Lag u s , which we want to retain in the following Lemma. Lemma 5.1. For a linear advection equation with flow direction ψ ∈ [0 , π ] and periodic initial condition u init on a triangular grid like in Figure 3, where the fully upwind is used, the universal update scheme (9)on τ , reads d u s d t = S u s with semi-discretization S given by S := − (cid:16) cos( ψ ) D ξ V − (cid:2) M − , T − , + M , + M , − T , − (cid:3) + sin( ψ ) D η V − (cid:2) M − , T − , + M , + M , − T , − (cid:3)(cid:17) E Lag . (33)We already stated in Section 4 that the universal update scheme (9) with modal filtering at the cell τ i,j reads d u si,j ( t )d t = − ( ξ x D ξ + ξ y D η ) M σ V − F ( t ) − ( η x D ξ + η y D η ) M σ V − F ( t ) , (34)when evaluated at the solution points. If we combine this with Lemma 5.1, we obtain a similar corollary forthe the filtered SD Method. Corollary 5.2. For a linear advection equation with flow direction ψ ∈ [0 , π ] and periodic initial condition u init on a triangular grid like in Figure 3, where the full upwind is used, the universal update scheme withmodal filtering on τ , reads d u s d t = S σ u s with semi discretization S σ given by S σ := − (cid:16) cos( ψ ) D ξ M σ V − (cid:2) M − , T − , + M , + M , − T , − (cid:3) + sin( ψ ) D η M σ V − (cid:2) M − , T − , + M , + M , − T , − (cid:3)(cid:17) E Lag . (35)17 .2 Stability Analysis for the SD Method We already stated that the asymptotic behavior of the solution of (29) depends on the real parts of theeigenvalues of S ≡ S ( α, β, γ ) given by (33). For linear stability, we wish them to lie in the left half of thecomplex plane. In case there are no such parameters α, β ∈ R + and γ > α + β − 1, we want at least to finda set of parameters ( α, β, γ ) with preferably small real parts. If we denote the maximal real part over alleigenvalues corresponding to a certain parameter set ( α, β, γ ) and a fixed set of test cases with respect to w x , w y , ψ by Λ ≡ Λ( α, β, γ ), this leads to the optimization problemarg min α,β ∈ R + ,γ>α + β − Λ( α, β, γ ) . (36)Note that this is a continuous and non-convex problem. Future studies could give a strict treatment of theoptimization problem (36), but since we just want to demonstrate the influence of the chosen polynomialbasis on the stability of the SD Method with and without modal filtering, we will just concentrate on certainsubranges for α, β and γ . In the following, we will concentrate on the subrange of P := { ( α, β, γ ) | α, β = 0 . , . , . . . , γ = α + β, α + β + 0 . , . . . , } due to the observation that higher parameters lead to much worse condition numbers, see 1. For a fixedparameter set ( α, β, γ ) in such a subrange we will compute all eigenvalues of S and look for the one with thegreatest real part. This will be done for several test cases with respect to the flow direction ψ ∈ [0 , π ] andthe wave numbers w x , w y ∈ [ − π, π ]. Then we will again look after the greatest real part among all thesetest cases. The resulting eigenvalue and its real part will be taken as an evidence for the stability of theunderlying parameter set ( α, β, γ ) and its corresponding polynomial basis. All of this is done in MatLab,where a descriptive pseudo code is given by Algorithm 1. The numerical results show that the maximal Algorithm 1 without filtering for α, β = 0 . . , γ = α + β : 0 . do L << − 42 ) for ψ = 0 : π : π , w x , w y = − π : π : π do (cid:46) test cases compute S v = eig(S) v re = real(v) Λ = max(v re ) if Λ > L then L = Λ return ( α, β, γ, L ) (cid:46) L is the greatest real partreal part of all eigenvalues L is independent of α, β and γ . Hence the linear stability is independent of thecorresponding basis of APK polynomials, which coincides with the following Theorem. Theorem 5.3. In Lemma 5.1 the eigenvalues of the semi-discretization S , given by (33), are independentof the basis of APK polynomials corresponding to α, β and γ .Proof. Looking at S , only the differential matrices D ξ , D η , given by (11), and the Vandermonde matrix V ,coming from the interpolation approach, depend on the underlying basis { ϕ k | ≤ k ≤ K F } of P N ( T ) .Hence the proof is done, when we can eliminate this matrices from S . Therefore note that D ξ = (cid:16) ∂ ξ ϕ k (cid:0) T , ( x sj ) (cid:1)(cid:17) K s ,K F j,k = (cid:0) ∂ ξ ϕ k ( x ) | x = x sj (cid:1) j,k = (cid:16) ∂ ξ (cid:2) K F (cid:88) i =1 ϕ k ( x Fi ) L i ( x ) (cid:3) | x = x sj (cid:17) j,k = (cid:0) K F (cid:88) i =1 ϕ k ( x Fi ) ∂ ξ L i ( x ) | x = x sj (cid:1) j,k = (cid:0) ∂ ξ L i ( x sj ) (cid:1) K s ,K F j,i · (cid:0) ϕ k ( x Fi ) (cid:1) K F ,K F i,k D Lagξ V and analogous D η = D Lagη V with D Lagξ := (cid:0) ∂ ξ L i ( x sj ) (cid:1) K s ,K F j,i , D Lagη := (cid:0) ∂ η L i ( x sj ) (cid:1) K s ,K F j,i independent of the basis { ϕ k } . Hence D ξ V − = D Lagξ , D η V − = D Lagη (37)and S = − (cid:0) cos( ψ ) D Lagξ [ M − , T − , + M , + M , − T , − ]+ sin( ψ ) D Lagη [ M − , T − , + M , + M , − T , − ] (cid:1) E Lag (38)is independent of the underlying basis { ϕ k } .Note that we have done the proof for an arbitrary basis { ϕ k } of P N ( T ) and that by (37) moreover themore general universal update scheme (13) for a general (non-linear) conservation law becomesd u s d t ( t ) = − ( ξ x D Lagξ + ξ y D Lagη ) F ( t ) − ( η x D Lagξ + η y D Lagη ) F ( t )with F ( t ) , F ( t ) independent of the basis { ϕ k } , when we use an interpolation approach to reconstruct theflux. This leads to the following corollary. Corollary 5.4. The SD Method (without filtering) with a interpolation approach is independent of theunderlying polynomial basis. Remark 5.5. While the SD Method with an interpolation approach theoretically is independent of theunderlying polynomial basis, in computation their condition properties may have an influence on the re-sulting scheme, see Table 1. By substituting D ξ by D Lagξ and D η by D Lagη and working with (38) in ourimplementation, we blind out condition properties of the underlying basis and just focus on the pure stabilityproperties of the SD Method. However, one should note that this won’t be the case in the next Subsection,where the Vandermonde matrix and its inverse will not cancel out each other. Then, to cover conditionissues as well, we focus on a proper subrange of parameters which will ensure good condition numbers. Insuch a subrange we then will optimize the linear stability of the SD Method.Table 6 shows the maximal real parts among all test cases for certain orders N . These values show N maximal real part L Table 6: maximal real parts that the SD Method on triangles is not stable, independent of the polynomial basis in the interpolationapproach. This was already stated by Van den Abeele, Lacor and Wang in [39]. As one can see in Table 6the instabilities, i.e. Re( λ ), of the method can be ’quite high’. To deal with this problem we will apply themodal filter in the SD Method in the next Section. As we proposed in Section 4, we use the SV Method seen as a spectral method with modal filtering toobtain milder instabilities in the von Neumann analysis and hence preserving the scheme to blow up atdiscontinuities. While we used D V − = D Lag in Section 5.2 to show that the SD Method is independent ofthe underlying basis this term becomes DM σ V − 19n the semi discretization S σ , given by (35), for the SD Method with modal filtering. This again leads to S σ ≡ S σ ( α, β, γ ) and the method to depend on the corresponding basis of APK polynomials. Followingsection 5.2 we use Algorithm 1 with S = S σ for a (linear) stability analysis for the SD Method with modalfiltering. But first, note that the values of the modal filter σ in the filter matrix M σ := diag (cid:0) σ , . . . , σ K F (cid:1) depend on certain constants, which were explained in Section 4.2. For the following numerical results we set h = √ , see (26) ,C fix = 12 , λ max = 1 , see (27) . For order N = 3, different orders p ∈ N and constants c ∈ N in the filter strength (26), Table 7 - 9 showsthe smallest real parts L (maximum over all test cases) and their corresponding parameters ( α, β, γ ) ∈ P . p c ( α, β, γ ) L p c ( α, β, γ ) L Table 7: maximal real parts for p = 1 and p = 2 p c ( α, β, γ ) L p c ( α, β, γ ) L Table 8: maximal real parts for p = 3 and p = 4 We can see in Tables 7 - 9 that one can obtain much milder instabilities of the method by modal filtering, p c ( α, β, γ ) maximal real part L Table 9: maximal real parts for p = 5 where also the choice of family of APK polynomials for the underlying basis of P N ( T ) has an influence.Note that in fact both p , c and every parameter α, β and γ effect the linear stability. Approximately, theparameter tuples (1 , , 6) and (2 , , 5) appear promising. However, by considering their conditional numbersin Table 1, the parameter tuple (1 , , 6) becomes unreasonable, due to its high condition numbers. Thus,we will focus on parameter tuple (2 , , 5) in the following example which will help us to find suitable choicesfor the parameter c and the filter order p . Example 5.6. First, we observe that the constant c has an influence on the linear stability by comparingthe parameter tuple (2 , , 5) for N = 3, p = 2 and different c . In this case we have L =2 . e + 00 for c = 2 , =2 . e + 00 for c = 3 ,L =1 . e + 00 for c = 4 ,L =1 . e + 00 for c = 5 ,L =1 . e + 00 for c = 6 ,L =1 . e + 00 for c = 7 ,L =1 . e + 00 for c = 8 , which indicates that by increasing parameter c the instability gets milder. This coincides with Theorem 4.2which states the equivalence of modal filtering by our natural filters and the SV Method. In this formulation c is proportional to ε N , see (26), which again is proportional to the dissipation we add to the underlyingconservation law by the viscosity term.In the same way, by comparing the parameter tuple (2 , , 5) for N = 3 and c = 8, we can observe theinfluence of the parameter p . Here L =1 . e + 00 for p = 1 ,L =1 . e + 00 for p = 2 ,L =1 . e + 00 for p = 3 ,L =1 . e + 00 for p = 4 ,L =1 . e + 00 for p = 5 , suggest p = 2 to be the best choice.Finally, we observe the influence of the parameter γ by comparing parameter tuples (2 , , γ ) for N = 3, c = 8 and p = 2. Here L =1 . e + 00 for γ = 4 ,L =1 . e + 00 for γ = 4 . ,L =1 . e + 00 for γ = 5 ,L =1 . e + 00 for γ = 5 . ,L =1 . e + 00 for γ = 6show the influence.The natural filter of (1 , , , , 5) and the shaped raised cosine filter of 8th order can beseen and compared in Figure 4.We end this Section by recommending the parameter tuple ( α, β, γ ) = (2 , , 5) with c = 8 and filter order p = 2 for the polynomial degree N = 3. The following Section will provide numerical test which buildaround this choices. After presenting our theoretical results we want to give a short numerical investigation to show that ourconclusions are justified. We consider the two dimensional Burgers’ equation u t ( x, y, t ) + u ( x, y, t ) (cid:0) u x ( x, y, t ) + u y ( x, y, t ) (cid:1) = 0in the domain [ − , with the initial condition u ( x, y ) = 14 + 12 sin (cid:0) π ( x + y ) (cid:1) and periodic boundary conditions u ( − , y, t ) = u (1 , y, t ) and u ( x, − , t ) = u ( x, , t ). At t = 0 . s twodiscontinuities arise at y = − x and y = − x . We use 1088 triangles for the spatial discretizationand present the numerical solutions at t = 0 . s for several parameter selections. Without filtering highoscillations already develop in much earlier calculations and the SD Method collapses.In Tables 10a-10b the minimum and maximum of the numerical solutions for filter parameters p and c are21 Values of the different filters for N=3 and p=2, c=8 index l f il t e r v a l ue f o r ( m + l ) / N index m polynomial degree m+l f il t e r v a l ue f o r ( m + l ) / N Values of Different Filters for N=3 Exp. filter (1,1,2) for p=2, p=8Exp. filter (2,2,5) for p=2, p=8Shaped raised cosine filter Figure 4: ( α, β, γ ) = (1 , , , , 5) with their natural filters and the shaped raised cosine filter shown in the polynomial cases ( α, β, γ ) = (1 , , 2) and (2 , , N in these test-cases is always three. For most parameters p , c and ( α, β, γ ), wewere not able to calculate in time further than 0 . s . However, for the choice p = 2 and c = 8 motivated byour stability analysis we were indeed able to calculate till 0 . 45 for both parameters (1 , , 2) and (2 , , p = 2 and c = 8 using the APK polynomials (2 , , 5) and their natural filters show significantstronger stability properties than for instance (1 , , , , 5) has a small condition number, see Table 1. But second and quiet more important are the goodproperties - using the natural filter with γ = 5 - from point of linear stability, suggested by the stabilityanalysis in Section 5. Following this observation, we wish to combine a small condition number like for(1 , , 2) and the good stability of the choice ( α, β, γ ) = (2 , , 5) with p = 2 and c = 8. Therefor we nextobserved the polynomial basis of (1 , , , , 5) instead of it’sown one. I.e we choose γ = 5 in (23) instead of the natural choice γ = 2. The resulting numerical solutions,again presented in a 2d and 3d plot, can be seen in Figure 6 and once more show an improvement.22 c min max2 8 -3.24 6.003 8 -6.55 3.674 6 -12.32 6.034 8 -11.08 4.96 (a) ( α, β, γ ) = (1 , , , t = 0 . s p c min max2 8 -1.55 2.463 8 -6.54 4.234 6 -3.10 5.984 8 -2.44 5.64 (b) ( α, β, γ ) = (2 , , , t = 0 . s ( α, β, γ ) = (1 , , , t = 0 . s ( α, β, γ ) = (2 , , , t = 0 . s Figure 5: N = 3, p = 2, C = 8 t = 0 . s t = 0 . s Figure 6: ( α, β, γ ) = (1 , , 2) with the natural filter of p = 2, c = 8 and (2 , , We started this paper by extending the SD method by more general bases of APK polynomials insteadof Lagrange or PKD polynomials and instantly observed their numerical qualification by their conditionnumbers. The classical SD method was known to be unstable on triangular elements and unfortunately wewere able to obtain the same result for every extended version of the SD method. Hence addressing theproblem of Gibb’s phenomenon, which often leads unstable schemes to blow up, we applied the well knownSV method to milder instabilities. Due to the APK polynomials to fulfill certain eigenvalue problems, wederived an equivalent but from point of computational costs much more efficient formulation of the SVmethod via modal filtering by natural exponential filters with respect to the polynomial basis. Whole newerror estimates for filtered APK extensions for smooth functions then were given. We also gave a vonNeumann analysis for linear stability for the extended SD method with modal filtering by natural filtersof APK polynomials. The idea was for the right choice of APK polynomials to may have some positiveinfluence on the stability and accuracy of the method. To the best of our knowledge, this work was thefirst to give such an analysis for a scheme with spectral filtering. This stability analysis indeed suggesteda certain choice of APK polynomials and parameters in their natural filters, which then were observed in23 short numerical investigation. The numerical solutions we obtained show significant stronger stabilityproperties than the ones resulting from classical SD methods. By combining our results from the analysisof condition numbers and linear stability for different APK polynomials, we equipped basis polynomialswith favorable condition numbers with modal filters coming from a polynomial basis with strong stabilityproperties. Once more, this led to a significant improvement of the numerical solutions.Despite this, it should also be said that Flux Reconstruction (FR) or Correction Procedure via Recon-struction (CPR) methods nowadays are getting more and more interesting, due to their strong stabilityproperties, while SD methods are receiving decreasing attention. We still analyzed modal filtering basedon SD methods to investigate the pure influence of the orthogonal polynomials and their natural filters. Inthe FR/CPR approach a correction term is applied which might cancel out effects of the filter. However,the observations made for spectral filtering in this paper target a first step to both a better understandingof the influence of spectral filtering and in future to gainfully adapt these ideas to FR/CPR methods. Im-provements might be higher CFL numbers, thus larger time steps and more efficient methods, as well asmore advanced methods in post-processing.In future researches further developments for modal filtering itself would be favorable. Due to highcondition numbers of the Vandermonde matrix for several bases of APK polynomials, we had to restrictourself to certain subranges of parameter families. However, it could be possible to find filters with evenbetter properties outside of these subranges. To bypass the problem of quite bad condition numbers in thisnew ranges, it would be promising to work with projection approaches in the series expansion, instead ofinterpolation ones. M m,ni for the SD Method of 3rd order ( K s = 6 , K F = 10): M − , = diag ( E × , × ) ,M , = diag (0 × , E × , × , E × ) ,M , − = diag (0 × , E × , × , E × , × ) ,M − , = diag (0 × , E × , × ) ,M , = diag (0 × , E × , × , E × , × , E × , × ) ,M , − = diag ( E × , × , E × , × , E × , × , E × ) .M m,ni for the SD Method of 4th order ( K s = 10 , K F = 15): M − , = diag ( E × , × ) ,M , = diag (0 × , E × , × , E × , × , E × ) ,M , − = diag (0 × , E × , × , E × , × , E × , × ) ,M − , = diag (0 × , E × , × ) ,M , = diag (0 × , E × , × , E × , × , E × , × , E × , × ) ,M , − = diag ( E × , × , E × , × , E × , × , E × , × , E × ) .M m,ni for the SD Method of 5th order ( K s = 15 , K F = 21): M − , = diag ( E × , × ) ,M , = diag (0 × , E × , × , E × , × , E × , × , E × ) ,M , − = diag (0 × , E × , × , E × , × , E × , × , E × , × ) ,M − , = diag (0 × , E × , × ) ,M , = diag (0 × , E × , × , E × , × , E × , × , E × , × , E × , × ) ,M , − = diag ( E × , × , E × , × , E × , × , E × , × , E × , × , E × ) . eferences [1] Y. Allaneau and A. Jameson. Connections between the filtered discontinuous Galerkin method and theflux reconstruction approach to high order discretizations. Computer Methods in Applied Mechanicsand Engineering , 200(49):3628–3636, 2011.[2] H. L. Atkins and C.-W. Shu. Quadrature-free implementation of discontinuous Galerkin method forhyperbolic equations. AIAA journal , 36(5):775–782, 1998.[3] G. E. Barter and D. L. Darmofal. Shock capturing with higher-order, PDE-based artificial viscosity. AIAA paper , 3823:2007, 2007.[4] M. Blyth and C. Pozrikidis. A Lobatto interpolation grid over the triangle. IMA Journal of appliedmathematics , 71(1):153–169, 2006.[5] M. H. Carpenter and C. A. Kennedy. Fourth-order 2n-storage Runge-Kutta schemes. Nasa tm ,109112:871–885, 1994.[6] M. Dubiner. Spectral methods on triangles and other domains. Journal of Scientific Computing ,6(4):345–390, 1991.[7] C. F. Dunkl and Y. Xu. Orthogonal polynomials of several variables . Number 155. Cambridge UniversityPress, 2014.[8] A. Gelb and J. Tanner. Robust reprojection methods for the resolution of the Gibbs phenomenon. Applied and Computational Harmonic Analysis , 20(1):3–25, 2006.[9] J. Glaubitz, H. Ranocha, P. ¨Offner, and T. Sonar. Enhancing stability of correction procedure viareconstruction using summation-by-parts operators II: Modal filtering, 2016. Submitted.[10] D. Gottlieb and J. S. Hesthaven. Spectral methods for hyperbolic problems. Journal of Computationaland Applied Mathematics , 128(1):83–131, 2001.[11] D. Gottlieb and S. A. Orszag. Numerical analysis of spectral methods: Theory and applications , vol-ume 26. Siam, 1977.[12] D. Gottlieb and C.-W. Shu. On the Gibbs phenomenon and its resolution. SIAM review , 39(4):644–668,1997.[13] J. Hesthaven and R. Kirby. Filtering in Legendre spectral methods. Mathematics of Computation ,77(263):1425–1452, 2008.[14] P. Huang, Z. Wang, and Y. Liu. An implicit space-time spectral difference method for discontinuitycapturing using adaptive polynomials. AIAA paper , 5255:2005, 2005.[15] H. Huynh. A flux reconstruction approach to high-order schemes including discontinuous Galerkinmethods. AIAA paper , 4079:2007, 2007.[16] H. Huynh. High-order methods including discontinuous Galerkin by reconstructions on triangularmeshes. AIAA Paper , 44:2011, 2011.[17] H. Huynh, Z. J. Wang, and P. E. Vincent. High-order methods for computational fluid dynamics: A briefreview of compact differential formulations on unstructured grids. Computers & Fluids , 98:209–220,2014.[18] A. Jameson. A proof of the stability of the spectral difference method for all orders of accuracy. Journalof Scientific Computing , 45(1-3):348–358, 2010.[19] G. Karniadakis and S. Sherwin. Spectral/hp element methods for computational fluid dynamics . OxfordUniversity Press, 2013.[20] R. M. Kirby and S. J. Sherwin. Aliasing errors due to quadratic nonlinearities on triangular spectral/hpelement discretisations. Journal of engineering mathematics , 56(3):273–288, 2006.[21] T. Koornwinder. Two-variable analogues of the classical orthogonal polynomials. Theory and applica-tions of special functions , pages 435–495, 1975. 2522] D. A. Kopriva. A conservative staggered-grid Chebyshev multidomain method for compressible flows.II. a semi-structured method. Journal of computational physics , 128(2):475–488, 1996.[23] D. A. Kopriva and J. H. Kolias. A conservative staggered-grid Chebyshev multidomain method forcompressible flow. Technical report, DTIC Document, 1995.[24] Y. Liu, M. Vinokur, and Z. Wang. Spectral difference method for unstructured grids I: Basic formula-tion. Journal of Computational Physics , 216(2):780–801, 2006.[25] H. Ma. Chebyshev–Legendre spectral viscosity method for nonlinear conservation laws. SIAM Journalon Numerical Analysis , 35(3):869–892, 1998.[26] Y. Maday, S. M. O. Kaber, and E. Tadmor. Legendre pseudospectral viscosity method for nonlinearconservation laws. SIAM Journal on Numerical Analysis , 30(2):321–342, 1993.[27] G. May and A. Jameson. A spectral difference method for the Euler and Navier-Stokes equations onunstructured meshes. AIAA paper , 304:2006, 2006.[28] A. Meister, S. Ortleb, and T. Sonar. Application of spectral filtering to discontinuous Galerkin methodson triangulations. Numerical Methods for Partial Differential Equations , 28(6):1840–1868, 2012.[29] A. Meister, S. Ortleb, T. Sonar, and M. Wirz. An extended discontinuous Galerkin and spectral differ-ence method with modal filtering. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschriftf¨ur Angewandte Mathematik und Mechanik , 93(6-7):459–464, 2013.[30] P. ¨Offner and T. Sonar. Spectral convergence for orthogonal polynomials on triangles. NumerischeMathematik , 124(4):701–721, 2013.[31] P. ¨Offner, T. Sonar, and M. Wirz. Detecting strength and location of jump discontinuities in numericaldata. Applied Mathematics , 4(12):1, 2013.[32] S. Ortleb. Ein diskontinuierliches Galerkin-Verfahren hoher Ordnung auf Dreiecksgittern mit modalerFilterung zur L¨osung hyperbolischer Erhaltungsgleichungen . kassel university press GmbH, 2011.[33] P.-O. Persson and J. Peraire. Sub-cell shock capturing for discontinuous Galerkin methods. AIAApaper , 112:2006, 2006.[34] H. Ranocha, J. Glaubitz, P. ¨Offner, and T. Sonar. Enhancing stability of correction procedure viareconstruction using summation-by-parts operators I: Artificial dissipation, 2016. Submitted.[35] H. Ranocha, P. ¨Offner, and T. Sonar. Summation-by-parts operators for correction procedure viareconstruction. Journal of Computational Physics , 311:299–328, 2016.[36] P. K. Suetin. Orthogonal polynomials in two variables , volume 3. CRC Press, 1999.[37] Y. Sun, Z. J. Wang, and Y. Liu. High-order multidomain spectral difference method for the Navier-Stokes equations on unstructured hexahedral grids. Communications in Computational Physics ,2(2):310–333, 2007.[38] E. Tadmor. Convergence of spectral methods for nonlinear conservation laws. SIAM Journal onNumerical Analysis , 26(1):30–44, 1989.[39] K. Van den Abeele, C. Lacor, and Z. Wang. On the stability and accuracy of the spectral differencemethod. Journal of Scientific Computing , 37(2):162–188, 2008.[40] K. Van den Abeele, C. Lacor, and Z. J. Wang. On the connection between the spectral volume and thespectral difference method. Journal of Computational Physics , 227(2):877–885, 2007.[41] H. Vandeven. Family of spectral filters for discontinuous problems. Journal of Scientific Computing ,6(2):159–192, 1991.[42] P. E. Vincent, P. Castonguay, and A. Jameson. Insights from von Neumann analysis of high-order fluxreconstruction schemes. Journal of Computational Physics , 230(22):8134–8154, 2011.[43] P. E. Vincent, P. Castonguay, and A. Jameson. A new class of high-order energy stable flux reconstruc-tion schemes. Journal of Scientific Computing , 47(1):50–72, 2011.2644] Z. Wang. High-order methods for the Euler and Navier–stokes equations on unstructured grids. Progressin Aerospace Sciences , 43(1):1–41, 2007.[45] Z. Wang and H. Gao. A unifying lifting collocation penalty formulation including the discontinuousGalerkin, spectral volume/difference methods for conservation laws on mixed grids. Journal of Com-putational Physics , 228(21):8161–8186, 2009.[46] Z. Wang, Y. Liu, G. May, and A. Jameson. Spectral difference method for unstructured grids II:Extension to the euler equations. Journal of Scientific Computing , 32(1):45–71, 2007.[47] M. Wirz. Ein Spektrale-Differenzen-Verfahren mit modaler Filterung und zweidimensionaler Kanten-detektierung mithilfe konjugierter Fourierreihen . Cuvillier, 2012.[48] M. Wirz. Detecting edges in high order methods for hyperbolic conservation laws. In High OrderNonlinear Numerical Schemes for Evolutionary PDEs , pages 151–167. Springer, 2014.[49] M. Yu, Z. Wang, and Y. Liu. On the accuracy and efficiency of discontinuous Galerkin, spectraldifference and correction procedure via reconstruction methods.