A SPACE-TIME DISCONTINUOUS GALERKIN METHOD FOR LINEAR HYPERBOLIC PDE’S WITH HIGH FREQUENCIES

The main purpose of this paper is to describe a space-time discontinuous Galerkin (DG) method based on an extended space-time approximation space for the linear first order hyperbolic equation that contains a high frequency component. We extend the space-time DG spaces of tensorproduct of polynomials by adding trigonometric functions in space and time that capture the oscillatory behavior of the solution. We construct the method by combining the basic framework of the space-time DG method with the extended finite element method. The basic principle of the method is integrating the features of the partial differential equation with the standard space-time spaces in the approximation. We present error analysis of the proposed spacetime DG method for the linear first order hyperbolic problems. We show that the new space-time DG approximation has an improvement in the convergence compared to the space-time DG schemes with tensor-product polynomials. Numerical examples verify the theoretical findings and demonstrate the effects of the proposed method.


Introduction
In computational acoustics, the medium frequency regime and multiscale wave propagation governed by the wave equation have been gained a constant interest in last decades.When multiscale wave propagation presents a high frequency component, developing an e¢ cient numerical methods for these classes of problem is a challenging task.Some example of high frequency problems include the high-intensity focused ultrasound (HIFU) treatment of cancer [1], coupled atomistic continuum modeling in nanomaterials [2] and tunneling in quantum mechanics [3].The reason for ine¢ ciency of the existing methods is that the standard numerical methods such as the …nite element (FEM) or discontinuous Galerkin (DG) methods based on semi-discrete approach require a very …ne mesh in the discretization 214 ŞUAYIP TOPRAKSEVEN in both space and time, and this leads to huge computational cost and makes the numerical methods ine¢ cient.Moreover, these methods based on semi-discrete approach may not suitable for multiscale approximations in the temporal domain.These issues on the standard numerical techniques have lead to high order methods that solve wave propagation phenomena in the time domain.One promising approach that has gained considerable popularity is space-time approximation spaces in which the time domain is also discretized.In these methods, two approaches have been proposed during the last decades.The …rst approach is called the time continuous space-time Galerkin methods (TCG) that do not require continuity in time.This approach generalizes the semi-discrete discretization to time domain with continuous time functions.The detailed explanations of such methods are given in [4].The drawback of these methods is high computational cost because of discretization of whole domain.The second approach is based on space-time discontinuous Galerkin methods that use standard polynomials spaces to discretize the problem in space and time, while temporal domains are divided into time slab and discontinuities and jumps are allowed in time.In each slab, TCG method is applied and the next slab uses the information from the previous slab.This second approach is more robust and e¢ cient than the …rst one.The wave equation can be discretized by a space-time setting in two ways.One way is to discretize the wave equation directly in a one …eld formulation with only one unknown as [5] and [6] .The second way is to convert the second order equation to a system of …rst order equations as done in [7] and [8].Using this second formulation, a priori and a posteriori error estimates have been proved in [9] using linear interpolation.This approach clearly increases the unknowns in the resulting systems.Error estimates to prove convergence of the methods have been derived by French [6] and Hughes and Hulbert [5].In the latter work, Galerkin least-squares stabilization terms are added for convergence analysis.In [6], the weighted inner product is included for the stability.A space-time DG method in which discontinuities and jumps are allowed both in space and time have been developed in [10] and recently proposed in [11] and [12] with discontinuous Petrov-Galerkin method in temporal domain for linear hyperbolic systems.Furthermore, many applications require boundary movement such as Stefan problems and water waves.In such problems, the mesh points also move in order to capture boundary movement.These movements in mesh points make the numerical scheme in e¢ cient or need more complicated numerical discretization.In this case, it is natural to consider the space-time discontinuous Galerkin approach.Analysis and survey of space-time DG method for hyperbolic and parabolic conservation laws on time dependent domains are explained in details in [13].Recently, space-time methods have become popular for the time dependent problems discussed in [14] and [15].An application of this method to the compressible Vaiver-Stokes equations is discussed in [17].Space-time DG method for the advection-di¤usion equation has been given in details in [18] and [19].This method also has been successively applied to nonconservative hyperbolic PDEs as models for dispersed multiphase ‡ows in [20].Furthermore, space-time DG methods have been proposed for the nonlinear water waves in [21] as well.
The medium or high frequency in wave propagation has been dealt with high order numerical methods including the ultra-weak variational method [22] that is a special case of the Tre¤tz-DG formulation for the wave equation [23] and the discontinuous enrichment method [24].In these methods, the approximation space is enriched by the solution of the equation under consideration.In [24], the DG space is extended by solutions of the homogeneous di¤erential equation that capture the high frequency in the solutions.In the same direction, an enriched space-time FEM for the …rst-order hyperbolic systems with discontinuities in both space and time has been studied by Chessa and Belytschko [25].This enriched space-time approach is based on the extended FEM studied in [26].These methods are based on the partition of unity approach developed in [27].Motivated by these approaches, in this paper we propose a high-order accurate space-time DG method that is wellsuited for …rst order linear hyperbolic problem with high frequency components.We construct the extended space-time DG space by enriching space-time DG space with the trigonometric functions in space and time.These trigonometric-function spaces intuitively capture the high frequency solutions and should be used to the highly oscillatory problems.This extended space time DG method is an extension to an extended DG method presented in [28].We will show global convergence in error estimates.Our error analysis based on the DG method proposed by [29].
The outline of this paper is as follows: Section 2 describes the mathematical analysis, formulations, and an introduction to a space-time DG method for scalar hyperbolic linear equation with high-frequency components.The basic properties of the proposed space-time DG method, the geometry of the space-time domain and elements and the space-time formulation of the problem have been explained and discussed and general solution form for linear hyperbolic equation with high frequency components is also given in Section 2. In Section 3, we introduce preliminaries and notations and recall some basic facts on DG methods for linear hyperbolic equations.We present our extended DG method for linear hyperbolic equation with high frequency components and our special interpolation operators have been given in Section 4. Stability and error analysis are given in Section 5.In Section 6 numerical example is given to show that our theoretical results agree with numerical results.Finally we explain some conclusion and future direction in Section 7.

Space-Time Formulation with Trigonometric Functions
The basic principle of an extended DG method is to enrich the DG space by special functions that are, generally, the solution of homogeneous di¤erential equation.The linear hyperbolic equation has the solution of the form h(x t) in one dimension.In Section 2.2, we show that if the initial condition has a high frequency component, then the homogeneous di¤erential equation will also have a high frequency functions.This observation suggests that the enrichment shape functions consist of the polynomials and the trigonometric functions in the space E := spanfsin(x t); cos(x t)g.A similar idea has been proposed for the wave equation in [30] and Tre¤tz DG method in [23].
Here, LU := @U @t + @U @x , Q T = (0; T ] and c 0 and c 1 are constants with 2 (0; 1] and U denotes a scalar quantity, and t represents time with T the …nal time.This problem has been chosen purely for its simplicity.This analysis can be easily extended to more general hyperbolic and scalar conservation law problems. We propose a space-time DG method based on extended DG approximation space for the equation (1).In this method, we directly consider the domain Q T R 2 in which spatial and temporal variables are not distinguished and a point x 2 Q T has coordinates (x 0 ; x 1 ) with x 0 = t representing a time variable and x = x 1 space variable.Thus, we de…ne the space-time domain as the open domain Q T R 2 .For space-time discretization, we need space-time slabs and elements.To do this, we partition the time interval I = (0; T ] into an ordered time levels 0 = t 0 < t 1 < < t N = T .Let I n = (t n ; t n+1 ) so that I = [ n I n with the time length t = t n+1 t n .Let (t n ) denote the space-time domain at the time level t = t n .Then, we de…ne space-time slabs as Q n T = Q T \ I n .We divide further (t n ) into non-overlapping spatial elements K n and similarly we divide the spatial domain (t n+1 ) into elements K n+1 .We then connect the elements K n and K n+1 to obtain space-time element K n by using linear interpolation in time.We also describe the tessellation of the space-time slab T n h = [ n K n and all space-time elements By @K we denote the boundary of the space-time element K.These space-time elements can be mapped to reference element (square or rectangle) by a suitable map, e.g., see [13] for construction of such a map. Figure 1 show a sketch of the space-time slab in Q T .
In this paper, we require c and g are slowly varying smooth functions with bounded derivatives of many orders while f has the high frequency components.For instance, if f (z) = cos(!z), then the solution has the form of: where the frequency ! is a large number in absolute value.We further assume the functions S(x; t) and R(x; t) are slowly-varying functions of x and t in the sense that they have many derivatives all of which have norms that are moderately sized in space.We can assume that the forcing term g(x; t) has also frequency components and we can show a similar solution form to (2).However, in this case extending the DG space with trigonometric functions is not easy task since we should extend the approximation space in all characteristic lines.For example, if we let g(x; t) = sin( x) + cos( t) with ; >> 1, then the solution form looks like U (x; t) = S 1 (x; t) sin( (x t)) + S 2 (x; t) cos( t) + R(x; t) cos(!(x t)) so that S 1 ; S 2 and R do not have high frequency component.Therefore, enriching the space-time DG space is not easy job in this simple example.As an application of this phenomena, we consider high-intensity focused ultrasound (HIFU) treatment of cancer that uses sound wave.Tumors in body tissues are destroyed when HIFU is focused onto them.The initial condition in partial di¤erential equation will generally help to determine high frequency shape to destroy tumor.In Figure 2, high frequency components in the initial condition determine acoustic pressure (high frequency shape) that heats and destroys the tumor.
We de…ne an interpolation based on the assumption (2) in the error analysis of the proposed method.Hence, we prove this assumption in the next subsection.
2.2.General Form of the Solutions.In this section, we give the explicit solution form of the following problem: 8 > < > : The variable change to characteristic lines helps transform the PDE to an in…nite set of ODE's.Let x = t + x 0 where x 0 2 R and de…ne Then, we …nd that or, with C(t) = c(t + x 0 ; t) and G(t) = g(t + x 0 ; t) we have We multiply this equation by an integrating factor ~ and …nd This …rst order linear BVP (in Ũ ) can now be solved and we …nd that So, now we unwind the variable change to produce a solution for U .Note that x 0 = x t and, thus, letting we have If we now assume that g and c are slowly varying smooth functions with bounded derivatives of many orders while f has the high frequency components; that is, say, then, we have U (x; t) = S(x; t) + R(x; t) cos(!(x t)): This proves the assumption (2).

Preliminaries and Notations
The Sobolev space, W m;p (K), for a domain K, consists of functions with m derivatives in the L p (K) norm.We will use the following notation for Sobolev space semi-norms and norms for 1 p < 1 where D v = @ j j v @x 1 1 : : : @x n n with j j = P n k=1 k for k 0; k = 1; : : : ; n, and x = x 1 1 : : : x n n and the standard Lebesgue space L p (K) norms and kuk L 1 (K) = ess sup x2K jv(x)j: For simplicity, we occasionally denote k:k 0;2;K by k:k.We will primarily be working with the Hilbert space H m (K) = W m;2 (K).Let P r (K) be the space of polynomials with degree q in K.We also assume there is an interpolation operator [ [31], Theorem 4.4.4]h : W q+1;p (K) !P r (K) for which k(I h )vk `;p;K Ch r `kvk r;p;K ; (0 ` r q + 1): and h = for 2 P r (K).Moreover, if v 2 C( ) then h v is continuous on as well.
The inverse inequality [ [31], Theorem 4.5.11]for functions 2 P r (K) states that there exists C > 0, which is independent of h, so that k k `;p;K Ch m `+1=p 1=r k k m;r;K ; m `; 1 p 1; and 1 r 1: (7) The arithmetic-geometric mean inequality states that for scalars a and b, where C = 1=(4 ) and > 0.
To be able to easily present our results and compare with previous works, we follow the paper by C. Johnson and J. Pitkaranta [29].Given a piecewise smooth function v write v n (:) = v (:; nh) and the approximate solution u is computed successively on the strips S n = fx 2 : (n l)h < t < nhg; n = l; :::; N so that ku n U n k is the error on each time level t = nh.
Let n K = (n K x ; n K t ) represent the outward pointing unit normal vector on @K with space coordinate n K x and the time coordinate n K t .Let := (1; ) and @Q T := .The in ‡ow boundary is de…ned In an element K, its in ‡ow boundary @K_ and its out ‡ow @K + = @Kn@K_ is similarly de…ned by Space-time DG space is then de…ned as where P r (K) denotes the space of polynomials of maximum degree at most r in (x; t).Functions in V h are allowed to be discontinuous at discrete time level.For K 2 T h , and a piecewise smooth function v, we de…ne the jump operator by when x 2 @K_ E, interior faces, and [u](x) := u(x) when x 2 _ \ @K_.The jump of v across @K_n de…ned similarly by where v + K the trace of v on @K taken from within the element K and v K is the exterior trace of u.Note that the sign of the jump depends on the direction of the ‡ow.The average of a function u is de…ned by We de…ne the equivalent space-time DG method for (1) by summing over where and where = (1; ) and n K is the outward unit normal to @K.Note that, for di¤erentiable functions u and v, we have the following integration by parts formula Equivalently, using this formula we may write this as where h := S K @K n .Now, we have the Galerkin orthogonality relation by replacing u by the exact solution U in ( 11) Let us recall the DG method for (1).Given a …nite element partitioning J h := fKg of Q T , we look for a solution u de…ned on Q T such that for all K 2 J h and uj K 2 P r (K) so that Z where M K := jn K : j.As shown in [29], u is uniquely determined by (15) and it is possible to compute u successively on each K starting at the in ‡ow boundary _ where u is given.Then it is possible to …nd the numerical solution u successively on one time level after another computing space-time element by element by starting for each strip on the left.The order of elements on which u will be computed is shown in Figure 3. Thus, given g and u on in ‡ow boundary, we can solve u locally in each K as shown below.For detailed proof, we refer the reader to [29].
Below, we denote by C , a positive constant which may take di¤erent values on di¤erent occurrences.Lemma 1. [29] Assume that g 2 L 2 (Q T ) and f 2 L 2 ( _) are given in (1).Then u is determined by (15) and the following local stability holds for each K Let us introduce a norm jjj:jjj h for the error analysis : Then, we have the following a generalization of the Poincare-Friedrichs inequality [31]:

Space-time Discontinuous Galerkin Discretization
In this section, we discuss the space-time DG method based on the extended space-time approximation spaces and that combines the framework of space-time DG with XFEM for the linear hyperbolic equation.We enrich the space-time DG space by adding a range of Fourier-series components to handle the high-frequency terms in the exact solution.As shown in Section 2.1, if the initial condition has only one high frequency component, then the solution form given by (2).In the same direction, if we assume that the initial condition has L high frequency components, that is, f (z) = P L `=1 cos(! `z) with !`>> 1, then we wil have the solution form of Thus, our "enriched space" X r h (Q T ) consists of the functions of the form (a `cos(n `(x t)) + b `sin(n `(x t)))g; where n 1 , : : : , n L are integers and s as well as the a i 's and the b i 's are all elements of the space-time DG space V h (9).More precisely, these s, a `and b `functions are tensor-product of piecewise discontinuous polynomials of degree at most r in x and t variables.Note that these functions are allowed to be discontinuous at the nodal points both in space and time and continuous in each element.This enriched space provides good approximations to the solutions of (1) if the range of high frequencies are known a priori.
For a high-frequency component of ( 17), we have, by using a simple trigonometric identities, R(x; t)cos(!`y) = R(x; t)cos(n `y + (! ` n `)y) = R(x; t) cos((!` n `)y) cos(n `y) R(x; t) sin((!` n `)y) sin(n `y); where y := x t and n `is an integer and can be chosen between n 0 and n L with 0 !` n ` 1.The key idea is that the functions `(x; t) = R(x; t) cos((!` n `)(x t)); and `(x; t) = R(x; t) sin((!` n `)(x t)) oscillate slowly since their frequencies are small and can be well approximated by functions in V h .
We now directly approximate the form (2) using interpolation.Let [( h `)(x; t)cos(n `(x t)) + ( h `)(x; t)sin(n `(x t))] : i : So, using (6), there is constant C, independent of h and !`, we have Also, since U is continuous, it follows that U A is continuous as well.
We remark that if U 2 P r (Q T ), that is U is a polynomials, then our special interpolation agrees with the interpolation operator h so that in this case we have Furthermore, we have, by the trace inequality, To construct our space-time approximation solution on the extended space, we perform the space-time discretization of the linear hyperbolic equation ( 1).Thus, we de…ne the space-time DG scheme : First, we show that the discrete problem ( 23) is stable from which the existence and uniqueness of the problem follows.Then we prove the bilinear form a(:; :) is coercive and continuous.For short reference, we take M K = M .
Lemma 2. The solution u to the problem (23) satis…es the following stability estimates jjjujjj Proof.Taking v = u + Lu when = Ch for some constant C in (23) we have, for each K 2 T h a(u; u + Lu) = (Lu + cu; u + Lu) Now, using Green's formula we have and since c c 0 > 0, we get Since every side of interior element boundary @K + agrees with a side of @K 0 _ for an neighbour element K 0 , we have (26) and consequently if we take = Ch for some constant C with h 1 and using the fact that Thus we obtain Now we estimate the right-hand side of (25).Applying the Cauchy-Schwarz and the arithmetic-geometric inequalities we get Using the fact that u_ = f on _ and combining ( 27) and ( 28), the desired result follows.
In particular this estimate shows the uniqueness and existence of a solution to (23) Now we prove the improved stability estimate.
Lemma 3. The following improved stability holds for = Ch with suitable constant Proof.From ( 27) and the de…nition of the norm (3) , the result easily follows.

Error Analysis of the Space-time DG Method
We can now state and prove the basic global error estimate for our space-time DG method (23).
where C does not depend on ! and h.
Proof.Let U A 2 X r h be the special interpolation of U de…ned by (19).Let us write := U U A ; := u U A ; e = : Using Lemma 3 with u = e and = Ch and the orthogonality property (14) with v = and the fact that e_ = 0 on _, we have In order to bound T 1 , we …rst prove that a(e; e) = jej 2 h : which proves the result (32).We next bound the term T 2 .To end this, we bound the bilinear form a(e; L ).So using the Cauchy-Schwarz and the arithmetic-geometric inequalities we have [e] 2 jM j ds + h Z @K_ 2 jM j ds ; thus we …nd that Cha(e; L ) h 2 [e] 2 jM j ds + h 2 Z @K_ 2 jM j ds ; and 2 jM j ds: Using the bounds (20) and (22) and the fact that the number of elements is O(h 2 ) we can bound the right hand-side by Finally inserting (32) and ( 33) into (31), it follows that or jjjejjj h C!h r+ 1 2 : This …nishes the proof of the error estimate (30).

Numerical Results
In this section, we will demonstrate some numerical experiments to verify our theoretical …ndings.Let us consider the scalar hyperbolic equation.We take Q T = (0; T ] = [0; 2 ] (0; T ] and the initial condition u(x; 0) = sin(!x);!= 100: We choose the boundary conditions so that the exact solution is given by u(x; t) = sin(!(xt)): In this example, the initial condition has only one high frequency component so we take L = 1 and n 1 !.For simplicity, we consider here only structured discretizations in space and time and choose h = 2 `; `= 3; 4; 5; 6; 7 and we use linear tensor-product polynomials, that is, r = 1.Thus we have 2 shape functions in the (unenriched) space-time DG space, and we have 6 shape functions (2 unenriched and 4 enriched) in the extended space-time DG space for the reference element I 2 = (0; 1) (0; 1).Matrix integrals are all done on a reference element by using 10 Gauss-Lobatto points numerical integration.The most simple shape functions of maximum degree r in the reference element can be given by ( 0 ; 1 ) = r0 0 r1 1 ; r = r 0 + r 1 : These shape functions give better conditioned mass and sti¤ness matrices, and make the computations relatively easier.De…ne the transformation G n K ( 0 ; 1 ) = (x; t); where with 0 and 1 are linear …nite element shape functions that are the images of 1 to the elements K n and K n+1 , respectively, by a suitable mapping.An example of such a mapping, F n K can be given as where x k (K n ) is the vertices of the element K n and ( 1 ) is the standard linear …nite element shape functions de…ned on I 2 .Thus, the space-time tessellation consists of the union of all the partitioning of the space-time slabs.For more detailed discussions of such mappings and other basis functions, see [13] and [14].
The numerical results are shown in Table 1 and Table 2.The observed convergence rates (OCR) of the proposed method in L 2 and the energy norms are given at T = 1 and T = 2: The observed convergence rate R 1 in L 2 norm is computed by the formula R 1 = log(jje 2h jj=jje h jj)= log(2) and the observed convergence rate R 2 in the energy norm is computed by the formula R 2 = log(jjje 2h jjj h =jjje h jjj h )= log( 2) where e h = u U is the error on the mesh.It is known that optimal convergence is observed only by using suitable chosen meshes.The loss of order h 1=2 in the order of convergence of L 2 norm is still under discussion, e.g., see [32].In practice, the optimal convergence h r+1 is achieved when polynomials of degree at most r used even if there is no uniform requirement on the chosen meshes.See [33] for the computational results for conforming triangulations for an example of this issue .Thus, typically there is a gap of order h 1=2 between computed convergence rate and the optimal convergence rate in DG methods.2. The errors and the order of convergence of the spacetime DG for the …rst order polynomial approximation (r = 1) at T = 1 and T = 2 in the energy norm.
The results in Table 2 clearly indicate that the numerical results are in good agreement with the theoretical …ndings and show that the proposed method convergences with the expected (r + 1=2)-th order of convergence when the polynomial space of order r is used without any mesh re…nement.

Conclusion
In this paper, we presented a space-time discontinuous Galerkin method for the scalar hyperbolic problems that contain high frequency components.We extend the space-time approximation space with trigonometric functions to capture the oscillatory behavior of the solutions.We applied discontinuous Galerkin methodology in both space and time and derived a stable space-time DG scheme.Thus, the method can be seen as a space-time framework of extended DG method.The key feature of the method is that it uses the solutions of PDE under consideration.Furthermore, the choice of DG space enriched by the solutions of the governing di¤erential equation enables an e¢ cient evaluation of integral terms.The proposed method here performs well when compared to standard space-time DG method.With conventional space-time DG method, one needs to re…ne the mesh size to get an acceptable accuracy for high frequency component.This leads to the computational costs in each space-time slab for solving the resulting system.We showed optimal a priori error estimates in a mesh dependent space-time DG norm.Additionally, we gave a numerical experiments to verify the theoretical …ndings.An extension of the analysis for an extended space-time solutions for the linear hyperbolic problems or conservation laws in two and three dimensional computational domains will be considered in the future.

2. 1 .
Problem Statement.In this paper, we consider a scalar hyperbolic equation in an open domain with boundary @

Figure 1 .
Figure 1.Space-time slab in space-time domain Q T .On the right, the rectangular mesh is an example of structured discretizations in space and time.

Figure 2 .
Figure 2. High frequency sound waves are concentrated on body tissues and tumor heats up and dies.

Figure 3 .
Figure 3.The order of the space-time elements on which u is computed.

By Green's formula for each K 2 (+ 2kek 2 :
Le; e) K = e )e + jM j ds o Now using the identity (26) with u = e along with the fact that e = 0 on _ we obtain that 2a(e; e) = X ) 2 2e + e + (e ) 2 jM j ds o + Z + (e ) 2 M ds + 2kek 2 ;

Table 1 .
The errors and the order of convergence of the spacetime DG for the …rst order polynomial approximation (r = 1) at T = 1 and T = 2 in L 2 norm.h jjje h jjj h;T =1 R 2 jjje h jjj h;T =2