Bernardo Cockburn

Discontinuous Galerkin methods

Key words discontinuous Galerkin methods, finite element methods

This paper is a short essay on discontinuous Galerkin methods intended for a very wide audience. We present the discontinuous Galerkin methods and describe and discuss their main features. Since the methods use completely discontinuous approxima-tions, they produce mass matrices that are block-diagonal. This renders the methods highly parallelizable when applied to hyperbolic problems. Another consequence of the use of discontinuous approximations is that these methods can easily handle irregular meshes with hanging nodes and approximations that have polynomials of different degrees in different elements. They are thus ideal for use with adaptive algorithms. Moreover, the methods are locally conservative (a property highly valued by the computational fluid dynamics community) and, in spite of providing discontinuous approximations, stable, and high-order accurate. Even more, when applied to non-linear hyperbolic problems, the discontinuous Galerkin methods are able to capture highly complex solutions presenting discontinuities with high resolution. In this paper, we concentrate on the exposition of the ideas behind the devising of these methods as well as on the mechanisms that allow them to perform so well in such a variety of problems.

Introduction

The discontinuous Galerkin (DG) methods are locally conservative, stable, and high-order accurate methods which can easily handle complex geometries, irregular meshes with hanging nodes, and approximations that have polynomials of different degrees in different elements. These properties, which render them ideal to be used with hp-adaptive strategies, not only have brought these methods into the main stream of computational fluid dynamics, for example, in gas dynamics [START_REF] Bassi | High-order accurate discontinuous finite element solution of the 2D Euler equations[END_REF][START_REF] Baumann | A discontinuous hp-finite element method for the solution of the Euler equation of gas dynamics[END_REF][START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin finite element method for conservation laws V: Multidimensional systems[END_REF], compressible [START_REF] Bassi | A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations[END_REF][START_REF] Lomtev | Simulations of viscous supersonic flows on unstructured hp-meshes[END_REF][START_REF] Lomtev | A discontinuous Galerkin method for the Navier-Stokes equations[END_REF][START_REF] Lomtev | Spectral/hp methods for viscous compressible flows on unstructured 2D meshes[END_REF] and incompressible [START_REF] Baumann | A discontinuous hp-finite element method for the Navier-Stokes equations[END_REF][START_REF] Cockburn | Local discontinuous Galerkin methods for the Oseen equations[END_REF][START_REF] Cockburn | Local discontinuous Galerkin methods for the Stokes system[END_REF] flows, turbomachinery [START_REF] Bassi | A high-order accurate discontinuous finite element method for inviscid and viscous turbomachinery flows[END_REF], magneto-hydrodynamics [START_REF] Warburton | A discontinuous Galerkin method for the viscous MHD equations[END_REF], granular flows [START_REF] Gremaud | On the computation of steady hopper flows I. Stress determination for coulomb materials[END_REF][START_REF] Gremaud | Simulation of gravity flow of granular materials in silos[END_REF], semiconductor device simulation [START_REF] Chen | Quantum hydrodynamic simulation of hysteresis in the resonant tunneling diode[END_REF][START_REF] Chen | Mixed-RKDG finite element methods for the 2-D hydrodynamic model for semiconductor device simulation[END_REF], viscoplastic crack growth and chemical transport [START_REF] Carranza | An adaptive discontinuous Galerkin model for coupled viscoplastic crack growth and chemical transport[END_REF], viscoelasticity [START_REF] Baaijens | Application of low-order discontinuous Galerkin methods to the analysis of viscoelastic flows[END_REF][START_REF] Bahhar | Galerkin discontinuous approximation of the transport equation and viscoelastic fluid flow on quadrilaterals[END_REF][START_REF] Fortin | New approach for the finite element method simulation of viscoelastic flows[END_REF], and transport of contaminant in porous media [START_REF] Aizinger | Local discontinuous Galerkin method for contaminant transport[END_REF][START_REF] Cockburn | Some extensions of the local discontinuous Galerkin method for convection-diffusion equations in multidimensions[END_REF][START_REF] Cockburn | Approximation of the velocity by coupling discontinuous Galerkin and mixed finite element methods for flow problems, Special issue: Locally Conservative Numerical Methods for Flow in Porous Media[END_REF][START_REF] Dawson | Local discontinuous Galerkin methods for problems in contaminant transport, in: Discontinuous Galerkin Methods[END_REF], but have also prompted their application to a wide variety of problems for which they were not originally intended like, for example, Hamilton-Jacobi equations [START_REF] Hu | The effect of the lest square procedure for discontinuous Galerkin methods for Hamilton-Jacobi equations[END_REF][START_REF] Hu | A discontinuous Galerkin finite element method for Hamilton-Jacobi equations[END_REF][START_REF] Lepsky | Analysis of the discontinuous Galerkin method for Hamilton-Jacobi equations[END_REF], second-order elliptic problems [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF][START_REF] Babuška | A discontinuous hp-finite element method for diffusion problems: 1-D analysis[END_REF][START_REF] Brezzi | Discontinuous Galerkin approximations for elliptic problems[END_REF][START_REF] Castillo | An a priori error analysis of the local discontinuous Galerkin method for elliptic problems[END_REF][START_REF] Cockburn | Superconvergence of the local discontinuous Galerkin method for elliptic problems on Cartesian grids[END_REF][START_REF] Cockburn | The local discontinuous Galerkin method for time-dependent convection-diffusion systems[END_REF][START_REF] Oden | A discontinuous hp-finite element method for diffusion problems[END_REF][START_REF] Rivière | Improved energy estimates for interior penalty, constrained and discontinuous Galerkin methods for elliptic problems. Part I[END_REF], elasticity [START_REF] Engel | Continuous/discontinuous finite element approximations of fourth-order elliptic problems in structural and continuum mechanics with applications to thin beams and plates, and strain gradient elasticity[END_REF][START_REF] Hansbo | Discontinuous finite element methods for incompressible and nearly incompressible elasticity by use of Nitsche's method[END_REF], and Korteweg-deVries equations [START_REF] Shu | A local discontinuous Galerkin method for KdV type equations[END_REF][START_REF] Shu | Local discontinuous Galerkin for partial differential equations with higher order derivatives[END_REF].

An introduction to DG methods can be found in the short monograph [START_REF] Cockburn | Discontinuous Galerkin methods for convection-dominated problems[END_REF]. A history of their development and the state of the art up to 1999 can be found in [START_REF] Cockburn | The development of discontinuous Galerkin methods, in: Discontinuous Galerkin Methods[END_REF]. Finally, a fairly complete and updated review is given in [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF]. This paper is a short essay on DG methods which differs from the above mentioned references in that it is intended for a wider audience and focuses on the exposition of the ideas behind the devising of these methods as well as in the mechanisms that allow them to perform so well in such a variety of problems.

Let us briefly carry out this program in what is perhaps the simplest situation, namely, that of approximating the solution of an ordinary differential equation (ODE). So, consider the following model initial-value problem:

d dt u(t) = f (t) u(t), t ∈ (0, T ), u(0) = u 0 , (1)
and suppose that we want to compute an approximation u h to u on the interval (0, T ) by using a DG method. To do that, we first find a partition of the interval (0, T ), { t n } N n=0 and set I n = (t n , t n+1 ) for n = 0, . . . , N -1. Then we look for a function u h which, on the interval I n , is the polynomial of degree at most k n determined by requiring that -

I n u h (s) d dt v(s) ds + u h v t n+1 t n = I n f (s) u(s) v(s) ds, (2)
for all polynomials v of degree at most k n . To complete the definition of the DG method, we still need to define the quantity u h . Note that the method establishes a link between the values of u h in different intervals only through u h . Since for the ODE, the information travels "from the past into the future", it is reasonable to take u h as follows:

u h (t n ) = ⎧ ⎨ ⎩ u 0 , if t n = 0, lim ↓0 u h (t n -) otherwise. (3)
This completes the definition of the DG method.

In this simple example, we already see the main components of the method, namely, (i) The use of discontinuous approximations u h , (ii) The enforcing of the ODE on each interval by means of a Galerkin weak formulation, and (iii) The introduction and suitable definition of the so-called numerical trace u h . The choice of the numerical trace is perhaps the most delicate and crucial aspect of the definition of the DG methods as it can affect their consistency, stability, and even accuracy. The simple choice we have made is quite natural for this case and gives rise to a very good method; however, other choices can also produce excellent results. Next, we address the question of how to choose the numerical trace u h .

Let us begin with the problem of the consistency of the DG method. As it is typical for most finite element methods, the method is said to be consistent if we can replace the approximate solution u h by the exact solution u in the weak formulation [START_REF] Alotto | Discontinuous finite element methods for the simulation of rotating electrical machines[END_REF]. We can immediately see that this is true if and only if u = u.

Next, let us consider the more subtle issue of the stability of the method. Our strategy is to begin by obtaining a stability property for the ODE (1) which we will then try to enforce for the DG method (2) by a suitable definition of the numerical trace u h . If we multiply the ODE by u and integrate over (0, T ), we get the equality

1 2 u 2 (T ) -1 2 u 2 0 = T 0 f (s) u 2 (s) ds,
from which the L ∞ -stability of the solution follows. The result we have to obtain for the DG method is a similar equality.

To do that, it is enough to set v = u h in the weak formulation [START_REF] Alotto | Discontinuous finite element methods for the simulation of rotating electrical machines[END_REF], integrate by parts and add over n. We get

N -1 n=0 - 1 2 u 2 h + u h u h t n+1 t n = 1 2 u 2 h (T -) + Θ h (T ) - 1 2 u 2 0 = T 0 f (s) u 2 h (s) ds,
where

Θ h (T ) = - 1 2 u 2 h (T -) + N -1 n=0 - 1 2 u 2 h + u h u h t n+1 t n + 1 2 u 2 0 .
Note that if Θ h (T ) were a non-negative quantity, the above equality would imply the stability of the DG method. In other words, the stability of the DG method is guaranteed if we can define the numerical trace u h so that Θ h (T ) ≥ 0. Setting

u h (t) = u 0 , t<0,
and using the notation

{u h } = 1 2 u - h + u + h , [[u h ]] = u - h -u + h , u ± h (t) = lim ↓0 u h (t ± ) ,
we rewrite Θ h (T ) as follows:

Θ h (T ) = - 1 2 u 2 h (T -) + - 1 2 u 2 h (T -) + u h (T ) u h (T -) + N -1 n=1 - 1 2 [[u 2 h ]] + u h [[u h ]] (t n ) -- 1 2 u 2 h (0 + ) + u h (0) u h (0 + ) + 1 2 u 2 0 = u h (T ) -u h (T -) u h (T -) + N -1 n=1 (( u h -{u h }) [[u h ]]) (t n ) -( u h (0) -u 0 ) u h (0 + ) + 1 2 [[u h ]] 2 (0),
where we used the simple identity

[[u 2 h ]] = 2 {u h } [[u h ]]. It is now clear that if we take u h (t n ) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ u 0 , if t n = 0, ({u h } + C n [[u h ]]) (t n ), if t n ∈ (0, T ), u h (T -), if t n = T,
where C n ≥ 0, we would have, setting

C 0 = 1/2, Θ h (T ) = N -1 n=0 C n [[u h ]] 2 (t n ) ≥ 0,
just as we wanted.

We thus see that it is possible to define the numerical trace u h to enforce the stability of the method. Note that the choice C n ≡ 1/2 corresponds to the numerical trace we chose at the beginning, namely, u h (t) = u h (t -). Note that the resulting DG methods are not only stable but consistent for all values of C n ≥ 0, as the condition u = u is satisfied. Moreover, in our search for stability, we found, in a very natural way, that the numerical trace u h (t) can only depend on both traces of u h at t, that is, on u h (t -) and on u h (t + ).

Finally, let us address the issue of how the choice of the coefficients C n can affect the accuracy of the method. It can be proven that if we take C n = 1/2, the order of the method at the points t n is 2k + 1, [START_REF] Lesaint | On a finite element method for solving the neutron transport equation[END_REF]. However, if we simply take C n = 0, the order is 2k + 2, [START_REF] Delfour | Discontinuous Galerkin methods for ordinary differential equations[END_REF]. Of course, in this case we must sacrifice the ability of solving for u h interval by interval since, if C n = 1/2, we must solve for u h in the whole computational domain (0, T ) at once. Thus, if we insist in having a DG method that is consistent, stable and can solve for the approximate solution interval by interval, we must take

C n ≡ 1/2!
Next, we want to emphasize three important properties of the DG methods that do carry over to the multi-dimensional case and to all types of problems. The first is that the approximate solution of the DG methods does not have to satisfy any interelement continuity constraint. As a consequence, the method can be highly parallelizable (when dealing with time-dependent hyperbolic problems) and can easily handle different types of approximations in different elements which renders it ideal for use with hp-adaptivity.

The second is that the DG methods are locally conservative. This is a reflection of the fact that the method enforces the equation element-by-element and of the use of the numerical trace. In our simple setting, this property reads

u h | t n+1 t n = I n f (s) u(s) ds,
and is obtained by simply taking v ≡ 1 in the weak formulation [START_REF] Alotto | Discontinuous finite element methods for the simulation of rotating electrical machines[END_REF]. This a much valued property in computational fluid dynamics.

The third property, which has not been properly discussed in the available literature on DG methods, is the strong relation between the residuals of u h inside the intervals and its jumps across inter-interval boundaries. To uncover it, let us integrate by parts in (2) to get

I n d dt u h (s) v(s) ds + ( u h -u h ) v| t n+1 t n = I n f (s) u(s) v(s) ds,
or, equivalently,

I n R(s) v(s) ds = (u h -u h ) v| t n+1 t n ,
where R denotes the residual d dt u hf u h . If we now take v = 1 and use the definition of the numerical trace u h , we obtain

I n R(s) ds = [[u h ]](t n ).
In other words, the jump of u h at t n , [[u h ]](t n ), is nothing but the integral of the residual over the interval I n . This means that if the ODE has been very well approximated in the interval I n , the jump [[u]](t n ) is very small. On the other hand, if the ODE was not well approximated therein, the jump [[u h ]](t n ) is big and the DG method becomes more dissipative as can be seen directly from the stability identity

1 2 u 2 h (T -) + 1 2 N -1 n=0 [[u h ]] 2 (t n ) - 1 2 u 2 0 = T 0 f (s) u 2 h (s) ds.
Thus, roughly speaking, the jumps [[u h ]] act as dampers that stabilize the DG method whenever the ride becomes bumpy and the ODE cannot be well approximated. This renders the method very stable without degrading its accuracy. In multi-dimensions, the residual is a more complicated linear functional of the jumps, but the dissipative mechanism just described remains essentially the same. This mechanism is shared with the so-called stabilized finite element methods for convection and for second-order problems.

In what follows, we consider DG methods for a variety of problems. In each of them, we study the properties of consistency, stability, and accuracy of the method, and discuss the relevance of the properties of local conservativity and stabilization by the jumps. We also comment on some computational issues and establish various links with other numerical schemes. In Sect. 2, devoted to linear problems, we consider (i) the transport equation, (ii) the wave equation, (iii) the Poisson equation, and (iv) the Oseen system. In Sect. 3, we consider the more difficult non-linear hyperbolic problems. We end the paper in Sect. 4, with some concluding remarks.

Linear problems 2.1 The transport equation

In this sub-section, we consider DG methods for the transport equation

u t + ∇ • (a u) = 0 in R d × (0, T ), u(t = 0) = u 0 on R d .
We consider only the discretization of this equation in space; the full discretization will be studied when when deal with non-linear hyperbolic problems in Sect. 3.

The objective here is to examine three properties that are especially relevant in this case. The first is that when the DG methods are strongly related to classical finite volume methods like the up-winding and the Lax-Friedrichs methods. The second is that when polynomials of high degree are used, the DG method can achieve high-order accuracy while remaining highly parallelizable. The third is that the artificial viscosity of the method is given by the size of the jumps which, in turn, are associated with the residual inside the elements. As a consequence, as the polynomial degree of its approximate solution increases, the artificial viscosity diminishes even in the presence of discontinuities.

The DG methods. To discretize the transport equation in space by using a DG method, we first triangulate the domain R d ; let us denote by T h such triangulation. We then seek a discontinuous approximate solution u h which, in each element K of the triangulation T h , belongs to the space V (K). There is no restriction in how to choose the space V (K), although a typical choice is the space of polynomials of degree at most k, P k (K). We determine the approximate solution on the element K by weakly enforcing the transport equation as follows:

K (u h ) t v - K a u h • ∇v + ∂K a u h • n v ds = 0, (4)
for all v ∈ V (K). To complete the definition of the DG method, it only remains to define the numerical trace a u h .

To do that, we proceed as in the ODE considered in the Introduction and begin by obtaining a stability result for the problem under consideration. We thus multiply the transport equation by u, and integrate over the space and time to get

1 2 R d u 2 (x, T ) dx + 1 2 T 0 R d ∇ • a(x) u 2 (x, t) dx dt = 1 2 R d u 2 0 (x) dx.
From this equation, a stability result immediately follows if we assume, for example, that -∇ • a ≤ L. Now, let us mimic the above procedure for the DG method under consideration. Taking v = u h in the weak formulation defining the approximate solution u h , and adding on the elements K, we get

1 2 R d u 2 h (x, T ) dx + 1 2 T 0 R d ∇ • a(x) u 2 h (x, t) dx dt + T 0 Θ h (t) dt = 1 2 R d u 2 h,0 (x) dx,
where

Θ h (t) = K∈T h - 1 2 K ∇ • (a u h )(x, t) dx + ∂K a u h (x, t) • n u h (x, t) ds .
Next, we investigate if it is possible to define the numerical trace a u h in such a way as to render Θ h non-negative. At this point, it is convenient to introduce the following notation. Let x be a point on the set e = ∂K + ∩ ∂K -and let n± denote the unit outward normal to ∂K ± at the point x. Let u ± h (x) denote the value lim ↓0 u h xn ± and set

{u h } = 1 2 u + h + u - h , [[u h ]] = u - h n -+ u + h n + .
Finally, let E h denote the set of sets e = ∂K + ∩ ∂K -for all K + and K -∈ T h .

We are now ready to rewrite Θ h in a suitable way. Indeed, dropping the argument t, we have

Θ h = K∈T h ∂K a u h • n u h - 1 2 a u 2 h • n ds = e∈E h e [[ a u h u h - 1 2 a u 2 h ]] ds = e∈E h e a u h [[u h ]] - 1 2 [[a u 2 h ]] ds = e∈E h e ( a u h -a {u h }) • [[u h ]] ds since 1 2 [[u 2 h ]] = {u h } [[u h ]]
. Thus, if we take

a u h = a {u h } + C[[u h ]],
we get

Θ h = e∈E h e C [[u h ]] • [[u h ]] ds ≥ 0,
if C is a non-negative definite matrix. This completes the definition of the DG space discretization and shows that it is possible to define the numerical trace as to render the method consistent and stable.

Examples of the DG methods.

There are two main examples of DG methods in this case. The first uses the following choice for its artificial viscosity coefficient:

C = 1 2 | a • n | Id. This implies that the numerical trace is a u h (x) = a lim ↓0 u h (x -a),
which is nothing but the classical up-winding numerical flux.

The second example is when we take C = 1 2 | a | Id. For this choice, we have

a u h = a {u h } + 1 2 | a |[[u h ]],
which is the so-called Lax-Friedrichs numerical flux.

Some properties. (i)

From the two examples above, we see that the DG methods are strongly related to finite volume methods. Indeed, the method of lines, that is, the discretization in space, for the up-winding scheme and the local Lax-Friedrichs scheme coincide with the corresponding DG method under consideration when the local space V (K) is taken to be the space of constant functions.

The DG methods, like finite volume methods, can easily handle complex computational domains. Also like finite volume methods, they have the property of being locally conservative, that is,

K (u h ) t dx + ∂K a u h • n ds = 0,
provided that the local space V (K) contains the constant functions. Indeed, this property is easily obtained from eq. ( 4) by simply taking the test function v to be a constant.

Unlike finite volume methods, however, DG achieve with ease high-order accuracy. Indeed, a theoretical order of convergence of k +1/2 can be proven simply by requiring that the local spaces V (K) contain all polynomials of degree at most k. Moreover, this is achieved while keeping a high degree of locality since to evolve the degrees of freedom of the approximate solution u h in an element, only the degrees of freedom of u h in the immediate neighbors are involved.

(ii) This last property and the fact that the mass matrix is block diagonal, and hence easily invertible, renders the DG methods extremely parallelizable when they are discretized in time by, for example, an explicit Runge-Kutta method.

We give an example taken from [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF]. In Table 1 below, we see the solution time and total execution time for the twodimensional linear problem with initial condition u(x, y, 0) = sin(πx) sin(πy) and periodic boundary conditions. We can see that as we progress from 1 to 256 processors while keeping the work per processor constant, a remarkable parallel efficiency is achieved. Polynomials of degree two and an explicit time-marching scheme were used.

u t + u x + u y = 0, in (-π, π) 2 × (0, T ),
(iii) Next, we show that the dissipation of the DG methods is given by the jumps of their approximate solution. This can be immediately seen when we realize that the DG method has a higher rate of dissipation of the energy, which in this case is nothing but the square of the L 2 -norm, than the exact solution of the transport equation. The extra rate of dissipation for the DG method is given by

Θ h = e∈E h e C [[u h ]] • [[u h ]] ds.
In the literature for monotone finite difference schemes for hyperbolic problems, the above term, when C = ν Id, is introduced by what could be considered to be a term modeling a viscosity effect, with ν being the viscosity coefficient, artificially inserted to render the scheme stable. That is why it is also called artificial viscosity. We thus see that the artificial viscosity of the DG method solely depends on the jumps of their approximate solution. Moreover, the jumps and the local residual (u h ) t + ∇ • (a u h ), which we denote by R, are strongly related. To see this, it is enough to carry out a simple integration by parts in the definition of the approximate solution:

K R v = ∂K (a u h • n v -a u h • n v) dx.
Note that in the case of the up-winding flux, we get

K R v = ∂K - a • [[u h ]] ds, where ∂K -= {x ∈ ∂K : a(x) • n(x) ≤ 0}.
In other words, the residual of u h in K is linearly related to the jump of u h on its inflow boundary ∂K -. A similar, but more complicated relation holds for general DG methods.

This indicates that the artificial viscosity of the method depends on the polynomial degree of the approximate solution. If the solution is very smooth, as the polynomial degree increases, we expect to have a smaller residual, hence a smaller jump and hence a smaller artificial viscosity. This is indeed what happens as the method can be proven to be of order k + 1/2 uniformly in time in the L 2 -norm when polynomials of degree k are used and the exact solution is smooth. Moreover, the artificial viscosity or, better, the dissipation of the DG methods, does decrease as the polynomial degree increases. This can be seen in Fig. 1 where we compare the exact solution of the transport equation

u t + u x = 0, in (0, 1) × (0, T = 100), with initial condition u(x, 0) = ⎧ ⎨ ⎩ 1, if x ∈ (.4, .6), 0, otherwise,
and periodic boundary conditions, with the DG approximations obtained with polynomial approximations of degree k = 0, 1, 2.

To march in time, a Runge-Kutta method of order k + 1 was used. Fig. 1 Effect of the polynomial degree on the dissipation of the DG method. The exact solution u (solid line) is contrasted against the approximate solution obtained on a mesh of 160 elements with piecewise-constant (dash-point line), piecewiselinear (dotted line) and piecewise-quadratic (dashed line) approximations.

The wave equation

In this subsection, we want to show how to extend the definition of DG methods to linear hyperbolic systems. To do that, we consider DG methods for the wave equation

u tt -c 2 Δu = 0, in R d × (0, T ),
where, for simplicity, we the take speed c to be constant, and rewrite it as a first-order hyperbolic problem, namely, as

U t + ∇ • F (U ) = 0, in R d × (0, T ), where U = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ q 1 q 2 . . . q d u ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , F(U ) = -c ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ u 0 . . . 0 0 u . . . 0 . . . . . . . . . . . . 0 0 . . . u q 1 q 2 . . . q d ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .
We want to show that the properties of local conservativity, stability, dissipation, and parallelizability of DG methods for the transport equation also hold for the above model hyperbolic problem. We also want to establish a natural link between the DG discretization of the wave equation and the DG discretization of the Laplacian operator.

The DG methods. To discretize the wave equation in space by using a DG method, we proceed as follows. After triangulating the domain R d , for each element K of the triangulation T h , we choose the local space U(K) which U h | K belongs to. Note, once again, that there is no constraint on how to take the space U(K); a typical choice, however, is

U(K) = P k (K) × . . . × P k (K).
Then, we determine the approximate solution on the element K by weakly enforcing the conservation law as follows:

K (U hi ) t V i - K F ij (U h ) V i,j + ∂K F ij n j V i dx = 0, for all V ∈ U(K).
Here, we adopt Einstein summation convention; note that V i,j denotes ∂ x j V i . To complete the definition of the DG method, it only remains to define the numerical trace F ij .

To do that, we begin by obtaining the stability result for the problem under consideration. We multiply the hyperbolic system by U , integrate over the space and time to get

1 2 R d U 2 (x, T ) dx = 1 2 R d U 2 (x, 0) dx,
where U 2 = U i U i . Now, let us mimic the above procedure for the DG method under consideration. Taking V = U h and adding on the elements K, we get

1 2 R d U 2 h (x, T ) dx + T 0 Θ h (t) dt = 1 2 R d U 2 h (x, 0) dx,
where

Θ h (t) = K∈T h - K F ij (U h )(x, t) U hi,j (x, t) dx + ∂K F ij (x, t) n j U hi (x, t) dx .
Next, we show can define the numerical trace F ij in such a way that the above quantity is non-negative. Dropping the argument t, we get

Θ h = - K F ij U hi,j dx + ∂K F ij n j U hi dx = K∈T h ∂K F ij n j U hi -c u q j n j dx, since F ij U hi,j = (c u q j ) ,j . Then, setting [[V i ]] j = V - i n j -+ V + i n j + , we get Θ h = e∈E h e F ij [[U hi ]] j -c [[u q j ]] j dx = e∈E h e F ij -F ij • [[U hi ]] j dx, since c [[u q j ]] j = c {u} [[q j ]] j + q j [[u]] j = c {u} δ ij [[q i ]] j + q j [[u]] j = F ij [[U hi ]] j .
Thus, if we take,

F ij = F ij + C ijk [[U hk ]] ,
we obtain

Θ h = e∈E h e C ijk [[U hk ]] [[U hi ]] j ds ≥ 0, provided C ijk is non-negative definite.
Examples of DG methods. We give three examples of DG methods for this case. We only have to determine their numerical trace F ij . Two of them are widely known in computational fluids dynamics as the up-winding and Lax-Friedrichs numerical fluxes. The third has been recently discovered in the context of DG methods for second-order elliptic problems. To write them down, we use the following notation:

[[q]] = [[q j ]] j .
They are as follows:

(i) F ij = F ij + |c| 2 [[q]] δ ij + |c| 2 [[u]] j δ i(d+1) (up-winding), (ii) F ij = F ij + |c| 2 [[q i ]] j + |c| 2 [[u]] j δ i(d+1) (Lax-Friedrichs), (iii) F ij = F ij + (C 22 [[q]] -C 12 • [[u]]) δ ij + (C 12 ) j [[q]] + C 11 [[u]] j δ i(d+1) .
The last numerical trace is, of course, a generalization of the up-winding numerical trace. It is not difficult to see that the dissipation term Θ h for each of these numerical traces is given by

(i) Θ h = e∈E h e |c| 2 [[q]] 2 + |c| 2 [[u]] 2 ds (up-winding), (ii) Θ h = e∈E h e |c| 2 [[q i ]] • [[q i ]] + |c| 2 [[u]] 2 ds (Lax-Friedrichs), (iii) Θ h = e∈E h e C 22 [[q]] 2 + C 11 [[u]] 2 ds.
Note that the dissipation of the first two methods is proportional to the speed of propagation |c|. Note also that the DG method with Lax-Friedrichs numerical trace is more dissipative than the DG method with the up-winding numerical trace since

[[q i ]] • [[q i ]] ≥ [[q]] 2 ,
where the equality holds only if the tangential component of q is continuous. Finally, note that if C 22 = C 11 = |c|/2, the third method produces the same amount of dissipation than the DG method using the up-winding numerical trace. This indicates that, for that method, the vector-parameter C 12 does not have any stabilizing effect; it could be used, however, to enhance the accuracy of the method. Some properties. It is not difficult to extend to the case of linear hyperbolic systems, what was discussed for the transport equation. Let us just briefly point out that the DG methods for those systems (i) are related to finite volume methods like the up-winding and the Lax-Friedrichs methods and always are locally conservative, (ii) are high-order accurate (in fact, of order (k + 1/2) when polynomials of degree k are used) while remaining highly parallelizable when discretized in time with explicit methods, and (iii) have a dissipative mechanism that solely depends on the jumps which, in turn, are strongly linked with the residuals inside the elements.

A point we are particularly interested in stressing here is that working with the wave equation allows us to see, in a very natural way, how to discretize the Laplace operator -Δ by using DG methods. Indeed, since the hyperbolic system for the wave equation can be rewritten as

q t -c∇u = 0, u t -c∇ • q = 0,
if we eliminate the term u t from the second equation and formally replace q t by q, we get q -c∇u = 0, -c∇ • q = 0, which is nothing but a rewriting of the equation c 2 Δu = 0. Thus, we can establish a one-to-one correspondence between the numerical traces used to discretize the hyperbolic system for the wave equation and those used to discretize the Laplace operator. This shows that switching from hyperbolic problems to elliptic ones is, in fact, does not entail a dramatic change as far as DG discretizations are concerned.

Second-order elliptic problems

In this sub-section, we consider DG methods for the model elliptic problem

-Δu = f in Ω, u= 0 on ∂Ω,
where Ω is a bounded domain of R d . Following what was done for the wave equation, to define them, we rewrite our elliptic model problem as

q = ∇u, -∇ • q = f in Ω, u= 0 on ∂Ω.
If when dealing with hyperbolic problems, we emphasized the DG methods are a generalization of finite volume methods, here we are going to show that DG methods are in fact mixed finite element methods. We also emphasize the fact that also in this context, DG methods are locally conservative methods ideally suited for adaptivity. Finally, we show that the dissipation mechanism of the DG methods, which is associated to the idea of artificial diffusion in the framework of hyperbolic equations, is associated to the idea of penalization of the discontinuity jumps in this context.

The DG methods. A DG numerical method is obtained as follows. After discretizing the domain Ω, the approximate solution (q h , u h ) on the element K is taken in the space Q(K) × U(K) and is determined by requiring that

K q h • v dx = - K u h ∇ • v dx + ∂K u h v • n ds, K q h • ∇w dx - ∂K w q h • n ds = K fw dx, for all (v, w) ∈ Q(K) × U(K).
Note that now we have two numerical traces, namely, u h and q h , that remain to be defined.

To do that, we begin by finding a stability result for the solution of the original equation. To do that, we multiply the first equation by q and integrate over Ω to get

Ω |q| 2 dx - Ω q • ∇u dx = 0.
Then, we multiply the second equation by u and integrate over Ω to obtain

- Ω ∇ • q u dx = Ω f u dx.
Adding these two equations, we get

Ω |q| 2 dx = Ω f u dx.
This is the result we sought. Next, we mimic this procedure for the DG method.

We begin by taking v = q h in the first equation defining the DG method and adding on the elements K to get

Ω |q h | 2 dx - K∈T h - K u h ∇ • q h dx + ∂K u h q h • n ds = 0.
Next, we take w = u h in the second equation and add on the elements to obtain

K∈T h K q h • ∇u h dx - ∂K u h q h • n ds = Ω fu h dx.
Adding the two above equations, we find that

Ω |q h | 2 dx + Θ h = Ω f u dx,
where

Θ h = - K∈T h - K ∇ • (u h q h ) dx + ∂K ( u h q h • n + u h q h • n) ds .
It only remains to show that we can define consistent numerical traces u h and q h that render Θ h non-negative. Since,

Θ h = K∈T h ∂K (u h q h • n -u h q h • n -u h q h • n) ds = e∈E h e [[ u h q h -u h q h -u h q h ]] ds = e∈E ih e ([[u h q h ]] -u h [[q h ]] -[[u h ]] • q h ) ds + ∂Ω (u h q h -u h q h • n -u h q h • n) ds = e∈E ih e (({u h } -u h ) [[q h ]] + [[u h ]] • ({q} -q h )) ds + ∂Ω (u h (q h -q h ) • n -u h q h • n) ds,
it is enough to take, inside the domain Ω,

q h = {q h } + C 11 [[u h ]] + C 12 [[q h ]], u h = {u h } -C 12 • [[u h ]] + C 22 [[q h ]],
and on its boundary,

q h = q h -C 11 u h n, u h = 0,
to finally get

Θ h = e∈E ih e C 22 [[q h ]] 2 + C 11 [[u h ]] 2 ds + ∂Ω C 11 u 2 h ds ≥ 0,
provided C 11 and C 22 are non-negative. Note how the boundary conditions are imposed weakly through the definition of the numerical traces. This completes the definition of the DG methods.

Some properties. (i)

Let us show that to guarantee the existence and uniqueness of the approximate solution of the DG methods, the parameter C 11 has to be greater than zero and the local spaces U(K) and Q(K) must satisfy the following compatibility condition:

u h ∈ U(K) : K ∇u h v dx = 0 ∀ v ∈ Q(K) then ∇u h = 0.
Indeed, the approximate solution is well defined if and only if, the only approximate solution to the problem with f = 0 is the trivial solution. In that case, our stability identity gives

Ω |q h | 2 dx + e∈E ih e C 22 [[q h ]] 2 + C 11 [[u h ]] 2 ds + ∂Ω C 11 u 2 h ds = 0,
which implies that q h = 0, [[u h ]] = 0 on E ih , and u h = 0 on ∂Ω, provided that C 11 > 0. We can now rewrite the first equation defining the method as follows:

K ∇u h v dx = 0, ∀ v ∈ Q h ,
which, by the compatibility condition, implies that ∇u h = 0. Hence u h = 0, as wanted.

(ii) When all the local spaces contain the polynomials of degree k, the orders of convergence of the L 2 -norms of the errors in q and u are k and k + 1, respectively; see [START_REF] Castillo | An a priori error analysis of the local discontinuous Galerkin method for elliptic problems[END_REF] when C 11 is of order O(h -1 ). Superconvergence in q has been proven and numerically observed in Cartesian grids, with a special choice of the numerical fluxes and equal-order elements with Q k -polynomials; see [START_REF] Cockburn | Superconvergence of the local discontinuous Galerkin method for elliptic problems on Cartesian grids[END_REF].

(iii) Next, we show that DG methods are in fact mixed finite element methods. To see this, let us begin by noting that the DG approximate solution (q h , u h ) can be also be characterized as the solution of

a(q h , v) + b(u h , v) = 0, -b(w, q w ) + c(u h , w) = F (w), for all (v, w) ∈ Q h × U h where Q h = {v : v ∈ Q(K) ∀ K ∈ T h }, U h = {w : w ∈ U(K) ∀ K ∈ T h },
and

a(q, r) := Ω q • r dx + E i C 22 [[q]][[r]] ds, b(u, r) := K∈T K u ∇ • r dx - E i ({u} + C 12 • [[u]]) [[r]] ds, c(u, v) := E ih C 11 [[u]] • [[v]] ds + ∂Ω C 11 uv ds, F (r) := Ω fv dx.
As a consequence, the corresponding matrix equation has the form

A -B t B C Q U = 0 F ,
which is typical of stabilized mixed finite element methods. As it is well known, those methods are not well defined unless the 'stabilizing' form c(•, •), usually associated with residuals, is introduced. For DG methods, the 'stabilizing' form c(•, •) solely depends on the parameter C 11 and the jumps across elements of the functions in U h . This is why we could think that this form stabilizes the method by penalizing the jumps, C 11 being the penalization parameter; thus, what was interpreted to be the artificial dissipation coefficient in hyperbolic problems can now be thought of a penalization parameter. Finally, let us emphasize that, for DG methods, penalizing the jumps is also a way of introducing stabilization by using residuals. Indeed, just as for the hyperbolic case, the residuals are related to the jumps. To see this, set R 1 = q h -∇u h and R 2 = -∇ • q hf and use the weak formulation of the DG method and the definition of its numerical trace to get

K R 1 • v dx = ∂K - 1 2 n -C 12 • [[u h ]] + C 22 [[q h ]] v • n ds, K R 2 w dx = ∂K - 1 2 n + C 12 [[q h ]] + C 11 [[u h ]] • n w ds, for all (v, w) ∈ Q(K) × U(K).
(iv) Let us now comment on an essential difference between the DG and the classical mixed methods. To solve the system associated to classical mixed methods, namely,

A -B t B C Q U = 0 F ,
we can try to eliminate Q from the equations to obtain

BA -1 B t + C U = F.
Since the matrix A is not easily invertible, due to continuity constraints on the approximation q h , we can hybridize the method; see [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF]. We then obtain a new system of the form

⎛ ⎜ ⎝ A -B t -D t B C 0 D 0 0 ⎞ ⎟ ⎠ ⎛ ⎜ ⎝ Q U Λ ⎞ ⎟ ⎠ = ⎛ ⎜ ⎝ 0 F 0 ⎞ ⎟ ⎠ ,
where Λ is the vector of degrees of freedom associated to the so-called Lagrange multiplier λ h . As is now well-known, the new vectors of degrees of freedom Q and U actually define the same approximation (q h , u h ) as the original mixed method. Moreover, since now A is block-diagonal, both Q and U can now be easily eliminated to obtain an equation for the multiplier only, namely,

E Λ = H ,
where E and H are given by

E = DA -1 A -B t BA -1 B t + C -1 B A -1 D t , H = -DA -1 B t BA -1 B t + C -1 F. (5)
In the case of DG methods, the above procedure it not necessary as the matrix A is block diagonal when C 22 is identically zero. In this case, the DG method remains well-posed and the matrix A can be made equal to the identity, if a suitable basis is used. This allows us to easily eliminate q h from the equations. Finally, note that, unlike classical and stabilized mixed methods, this can be achieved even if polynomials of different degrees are used in different elements. This renders DG methods ideal for adaptivity.

(v) The methods we have presented are locally conservative. As we saw in the hyperbolic case, this is a reflection of the form of the weak formulation and the fact that the definition of the numerical traces on the face e does not depend on what side of it we are. More general DG methods define the approximate solution by requiring that

K q h • r dx = - K u h ∇ • r dx + ∂K u h,K r • n K ds, K q h • ∇v dx = K fv dx + ∂K v q h,K • n K ds, for all (r, v) ∈ Q(K) × U(K).
In this general formulation, the numerical traces u h,K and q h,K can have definitions that might depend on what side of the element boundaries we are. Hence they are not locally conservative. This is the case for the numerical fluxes in u of the last four schemes in Table 2 taken from [START_REF] Arnold | Discontinuous Galerkin methods for elliptic problems[END_REF]. The acronym LDG stands for local DG, and IP for interior penalty.

Finally, let us point out that, in that table, the function α r ([u h ]) is a special stabilization term introduced by Bassi and Rebay [START_REF] Bassi | A high-order accurate discontinuous finite element method for inviscid and viscous turbomachinery flows[END_REF] and later studied by Brezzi et al. [START_REF] Brezzi | Discontinuous Galerkin approximations for elliptic problems[END_REF]; its stabilization properties are equivalent to the one originally presented. In [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF], a complete study of these methods is carried out in a single, unified approach.

The Oseen system

In this sub-section, we show how to discretize the Oseen system, namely,

-Δu + (a • ∇) u + ∇p = f , ∇ • u = 0 in Ω, u = 0 on ∂Ω,
where Ω is a bounded domain of R d and ∇ • a = 0.

The objective here is to put to work what we have learned about DG-space discretizations for convection and diffusion operators.

Table 2 Some DG methods and their numerical fluxes.

Method q e,K u h,K

Bassi-Rebay [START_REF] Bassi | A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations[END_REF]

{q h } { u h } LDG [38] {q h } + C 11 [[u h ]] -C 12 [[q h ]] {u h } + C 12 • [[u h ]] DG [22] {q h } + C 11 [[u h ]] -C 12 [[q h ]] {u h } + C 12 • [[u h ]] + C 22 [[q h ]] Brezzi et al. [20] {q h } -α r ([[u h ]]) {u h } IP [45] {∇u h } + C 11 [[u h ]] {u h } Bassi-Rebay [12] {∇u h } -α r ([[u h ]]) {u h } Baumann-Oden [15] {∇u h } { u h } -n K • [[u h ]] NIPG [70] {∇u h } + C 11 [[u h ]] {u h } -n K • [[u h ]] Babuška-Zlámal [7] C 11 [[u h ]] u h | K Brezzi et al. [20] -α r ([[u h ]]) u h | K
The DG methods. To define a DG method for the Oseen system, we begin by rewriting it as a first-order system,

σ i = ∇u i , -∇ • σ i + ∇ • (a u i ) + ∂ i p = f i , 1 ≤ i ≤ d, ∇ • u = 0 in Ω, u = 0 on ∂Ω.
where u i denotes the i-th component of the velocity u. To discretize the above equations, we take the approximate solution (σ h , u h , p h ) on the element K to be in the space S(K) d × U(K) d × P(K) and determine it by imposing that, for 1 ≤ i ≤ d, and for all (τ , v, w) ∈ S(K) × U(K) × P(K),

K σ i h • τ dx = - K u i h ∇ • τ dx + ∂K u σ,i h τ • n ds, K (σ i h • ∇v -u i h a • ∇v -p h ∂ i v) dx + ∂K -σ i h • nv + a u i h • nv + p h v n i ds = K f i v dx, - K u h • ∇q dx + ∂K u p,h • n q ds = 0.
It remains to find numerical traces that ensure that the method is well defined and stable. As usual, we begin by finding the stability equality for the continuous case. To do that, we first multiply the equation defining σ i by σ i sum over i and integrate over Ω to get

Ω |σ| 2 dx = Ω σ i • ∇u i dx,
where |σ| 2 = σ i • σ i ; we are using the Einstein summation convention. Now, we multiply the second equation of the Oseen system by u i , sum over i and integrate over Ω to get

Ω -∇ • σ i u i - 1 2 ∇ • a|u| 2 + u • ∇p dx = Ω f • u dx,
where |u| 2 = u i u i . Adding the last two equations, and integrating by parts, we get

Ω |σ| 2 dx - Ω p ∇ • u dx = Ω f • u dx,
and using the incompressibility of the velocity u,

Ω |σ| 2 dx = Ω f • u dx.
Next, we mimic the above procedure for the DG method. First, we take τ = σ i h , add over i and then over K, to get

Ω |σ h | 2 dx = - K∈T h K u i h ∇ • σ i h dx + K∈T h ∂K u σ,i h σ i h • n ds.
Now, take v = u i,h , add over i and then over K to obtain

K∈T h K σ i h • ∇u i h - 1 2 a • ∇|u h | 2 -p h ∇ • u dx + K∈T h ∂K -σ i h • nu i h + a u i h • nu i h + p h u h • n ds = Ω f • u h dx.
Adding the last two equations, we get

Ω |σ h | 2 dx + K∈T h K ∇ • (σ i h u i h ) - 1 2 a • ∇|u h | 2 -p h ∇ • u dx + K∈T h ∂K -u σ,i h σ i h • n -σ i h • nu i h + a u i h • nu i h + p h u h • n ds = Ω f • u h dx.
Finally, we take q = p h , add over K and add the result to the above equation to get

Ω |σ| 2 dx + Θ h = Ω f • u dx,
where

Θ h = Θ 1 h + Θ 2 h ,
and

Θ 1 h = K∈T h K ∇ • σ i h u i h - 1 2 a |u h | 2 -p h u h dx = e∈E h e σ i h u i h - 1 2 a |u h | 2 -p h u h dx,
and

Θ 2 h = K∈T h ∂K -u σ,i h σ i h • n -σ i h • nu i h + a u i h • nu i h + p h u h • n + u p,h • n p h ds = eE h e [[ -u σ,i h σ i h -σ i h u i h + a u i h u i h + p h u h + u p,h p h ]] ds.
Adding and rearranging terms, we get

Θ h = e∈E ih e {u i h } -u σ,i h [[σ i h ]] + ({σ i h } -σ i h ) [[u i h ]] -{a u i h } -a u i h [[u i h ]] -({p h } -p h ) [[u h ]] -{u h } -u p,h [[p h ]] ds + ∂Ω -u σ,i h σ i h + (σ i h -σ i h ) u i h - 1 2 a u i h -a u i h u i h -(p h -p h ) u h + u p,h p h • n ds.
We now take the numerical traces as to render Θ h non-negative. On the interior of the domain, we take

σ i h = {σ i h } -C 11 [[u i h ]] + C 12 [[σ i h ]], a u i h (x) = a lim ↓0 u i,h (x -a) , u σ,i h = {u i h } -C 12 • [[u i h ]],
and, on the boundary,

σ i h = σ i h -C 11 u i h n, u σ,h = 0.
The numerical traces associated with the incompressibility constraint, u p,h and p h , are defined by using an analogous recipe.

In the interior of Ω, we take

u p,h = {u h } + D 11 [[p h ]] -D 12 [[u h ]], p h = {p h } + D 12 • [[p h ]],
and on the boundary, we take

u p,h = 0, p h = p + h .
With the above choice, we get

Θ h = e∈E ih e (C 11 + |a • n|) |[[u h ]]| 2 + D 11 [[p h ]] 2 ds + 1 2 ∂Ω |a • n||u h | 2 ds ≥ 0,
whenever the parameters C 11 and D 11 are non-negative. This completes the definition of the DG method for the Oseen system. Some properties. (i) Let us show that the above DG method is well defined when the stabilization parameters C 11 and D 11 are positive and when the two following compatibility conditions on the local spaces are satisfied:

u h ∈ U(K) : K τ • ∇u h dx = 0 ∀ τ ∈ S d (K) then ∇u h = 0,
and

p h ∈ P(K) : K v • ∇p h dx = 0 ∀ v ∈ U d (K) then ∇p h = 0.
Thus, we must show that when f = 0, the only solution is u h = σ i h = 0 and p h is a constant. From the stability result, we get that

Ω |σ| 2 dx + e∈E ih e (C 11 + |a • n|) |[[u h ]]| 2 + D 11 [[p h ]] 2 ds + 1 2 ∂Ω |a • n||u h | 2 ds = 0, which implies that σ i h = 0, [[u i h ]] = 0, u i h = 0 on ∂Ω, and [[p h ]] = 0.
With this information, the first equation defining the DG method becomes

K τ • ∇u i h dx = 0 ∀ τ ∈ S d ,
which, by the first compatibility condition implies that ∇u i h = 0 on each element K. Since u h is continuous and zero in the border, we must have that u h = u. Finally, the second equation defining the DG method becomes

K v • ∇p h dx = 0 ∀ v ∈ U d (K),
which, by the second compatibility condition, implies that ∇p h = 0 on each element K. Since p h is continuous, this shows that p h is a constant.

(ii) In obtaining error estimates, one of the main issues is the so-called inf-sup condition. It turns out that, when D 11 > 0, it is possible to circumvent that condition and obtain optimal error estimates. Indeed, if we take S(K) = U(K) = P(K) to be the space of polynomials of degree k, the L 2 -norm of the error in p and σ i converge with order k, and the L 2 -norm of the velocity with order k + 1. Moreover, if polynomials of degree k -1 are used to approximate the pressure p and the stress tensor σ i , the above mentioned orders of convergence remain the same. This method is, however, not more efficient; see [START_REF] Cockburn | Local discontinuous Galerkin methods for the Stokes system[END_REF].

The fact that the inf-sup condition can be circumvented is typical of stabilized numerical methods. In the case of DG methods, this happens because it is possible to obtain a much weaker, generalized inf-sup condition. Moreover, in some cases, this inf-sup condition can be proven even if we take D 11 = 0. In such a case, the conditions on the velocity and pressure spaces are more restrictive, of course; see [START_REF] Schötzau | hp-DGFEM for incompressible flows[END_REF][START_REF] Toselli | hp-discontinuous Galerkin approximations for the Stokes problem[END_REF].

(iii) Let us end this section by pointing out that the method is locally conservative and that all the local residuals can be easily computed solely in terms of the jumps of the approximate solution. Moreover, the method is robust with respect to the Reynolds number. To support this claim, we quote a result from [START_REF] Cockburn | Local discontinuous Galerkin methods for the Oseen equations[END_REF] about the Kovasznay flow with different Reynolds numbers. We use quadrilateral meshes generated by consecutive refinements. In each refinement step, each grid cell is divided into four similar cells by connecting the edge midpoints. Therefore, grid level L corresponds to a mesh-size h L = 2 1-L . All the unknowns are approximated with tensor product polynomials of degree k ≥ 1. The discontinuity and pressure stabilization functions σ and δ are chosen of the order O(1/h) and O(h), respectively. Fig. 2 shows robustness of the discretization with respect to the Reynolds number.

Non-linear hyperbolic conservation laws

In all the previous section, we have considered DG-space discretizations for linear problems ranging from ODEs to the Oseen system. In computational fluid dynamics, the next natural step would be to consider the incompressible case. However, this topic is currently being vigorously addressed by several researchers; see the pioneering work [START_REF] Baker | Piecewise solenoidal vector fields and the Stokes problem[END_REF][START_REF] Karakashian | A nonconforming finite element method for the stationary Navier-Stokes equations[END_REF].

In this section, we switch to the numerical solution of non-linear hyperbolic equations. One of the main applications is to devise a locally conservative, stable and high-order accurate method to discretize the Euler equations of gas dynamics, as this is the bottleneck to deal with the compressible Navier-Stokes equations.

The RKDG methods: Introduction

The DG methods we consider are called the Runge Kutta discontinuous Galerkin (RKDG) methods, [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: The multidimensional case[END_REF][START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One dimensional systems[END_REF][START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for scalar conservation laws II: General framework[END_REF][START_REF] Cockburn | The Runge-Kutta local projection P 1 -discontinuous Galerkin method for scalar conservation laws[END_REF][START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin finite element method for conservation laws V: Multidimensional systems[END_REF]. To describe these methods, we use the simple model problem of the non-linear hyperbolic scalar conservation law

u t + ∇ f (u) = 0.
The RKDG methods are obtained in three steps:

Step 1: The DG space discretization. First, the conservation law is discretized in space by using a DG method. A discontinuous approximate solution u h is sought such that when restricted to the element K, it belongs to the finite dimensional space U(K). It is defined by imposing that, for all v h ∈ U(K),

K (u h ) t v h dx - K f (u h ) • ∇v h dx + ∂K f (u h ) • n K v h ds = 0.
Here, as we have seen in the previous section, the proper definition of the numerical trace f (u h ) (called, in this kind of problems, the approximate Riemann solver or also the numerical flux) is essential for the stability and convergence of the method.

Step 2: The RK time discretization. Then, we discretize the resulting system of ordinary differential equations, d dt u h = L(u h ), by using special explicit Runge-Kutta (RK) methods:

1. Set u (0) h = u n h . 2. For i = 1, . . . , K compute the intermediate functions:

u (i) h = i-1 l=0 α il w il h , w il h = u (l) h + β il α il Δt n L h u (l) h . 3. Set u n+1 h = u K h .
The distinctive feature of these RK methods is that their stability follows from the stability of the mapping u (l)

h → w il h defining the intermediate steps. Unfortunately, this mapping is not necessarily stable and, as a consequence, the method requires another component to enforce stability.

Step 3: The generalized slope limiter. This component is nothing but the so-called generalized slope limiter ΛΠ h . This non-linear projection operator, is devised in such a way that if u (l) h = ΛΠ h v h for some function v h , then the mapping u

(l) h → w il h is stable.
Thus, we incorporate the generalized slope limiter in the above time-marching algorithm as follows:

1. Set u (0) h = u n h . 2. For i = 1, . . . , K compute the intermediate functions:

u (i) h = ΛΠ h i-1 l=0 α il w il h , w il h = u (l) h + β il α il Δt n L h u (l) h . 3. Set u n+1 h = u K h .

This is the general form of the RKDG methods

In what follows, we elaborate on each on the above points and pay special attention to how to enforce the stability of the method. This is by far more delicate than what was done in the previous section because, as we can see, the stability of the method relies on the stability of the forward Euler step u

(l) h → w il h .
Unfortunately, this operator is always unstable in the L 2 -norm. Nevertheless, we are going to show that it is possible to find a weaker stability property for the forward Euler and then enforce it with the generalized slope limiter to obtain a stable method. The remarkable fact is that this can be achieved while maintaining the high-order accuracy of the method.

The DG-space discretization

With respect to the DG-space discretization, we do not have to say too many things that are different from the DG-space discretization of the transport equation. However, let us very briefly emphasize, that the method is locally conservative, and has order k + 1/2 when polynomials of degree k are used. It is highly parallelizable since its mass matrix is block diagonal and since, as it is typical of DG methods, the numerical trace f solely depends on the traces of u h on both sides of the inter-element boundary. The classical examples are the following:

(i) The Godunov flux:

f G (a, b) = ⎧ ⎨ ⎩ min a≤u≤b f (u), if a ≤ b, max b≤u≤a f (u), otherwise.
(ii) The Engquist-Osher flux:

f EO (a, b) = b 0 min(f (s), 0) ds + a 0 max(f (s), 0) ds + f (0).
(iii) The Lax-Friedrichs flux:

f LF (a, b) = 1 2 [f (a) + f (b) -C (b -a)], C = max inf u 0 (x)≤s≤sup u 0 (x) |f (s)|.
When piecewise constant approximate solutions are taken and when the forward Euler method is used to march in time, a monotone scheme is obtained. Monotone schemes are not only very stable but converge to the physically relevant solution, the so-called entropy solution. Unfortunately, monotone schemes are at most first order accurate. By using high-degree polynomials in a DG-space discretization, the accuracy of the scheme is raised. However, a RK time-discretization that that has the same accuracy in time and renders the scheme stable has to be devised.

The RK time discretization

The RK method we consider is required to satisfy the following conditions:

(i) If β il = 0 then α il = 0. (ii) α il ≥ 0. (iii) i-1 l=0 α il = 1.
Note that, by the first property, we can express the RK method in terms of the functions w il h . Note also that if we assume that, for some semi-norm | • |, we have that h → w il h implies the stability of the RK method! We must also keep in mind that, the RK method being explicit, we must ensure that the round-off errors are not amplified. For DG-space discretizations using polynomials of degree k and a (k + 1)-stage RK method of order k + 1 (which give rise to a (k + 1)-st order accurate method), a von Neumann stability analysis for the one-dimensional linear case f (u) = c u gives us the stability

w il h ≤ u (l) h , then u (i) h = i-1 l=0 α il w il h , ≤ i-1 l=0 α il w il h ,
|c| Δt Δx ≤ 1 2 k + 1 .
This is a condition to be respected, since otherwise the round-off errors will be amplified.

The stability of the step u h → w

h = u h + δ L h (u h )
It is not difficult to prove that the forward Euler step is not stable in the L 2 -norm, except in the case in which polynomials of degree 0 are used. If polynomials of degree k > 0 are used, the Euler step can be rendered stable if Δt/Δx is proportional to (Δx) p (k), where p(k) > 0; for example, p(1) = 1/2, see [START_REF] Chavent | The local projection P 0 P 1 -discontinuous-Galerkin finite element method for scalar conservation laws[END_REF]. This is clearly an unacceptable situation for hyperbolic problems. This negative result prompted the search of weaker norms, or semi-norms, for which the forward Euler step would be stable. To do that, recall that, when using polynomials of degree zero, the DG methods are nothing but monotone schemes, which are stable methods for the L ∞ norm in several space dimensions and in the total-variation semi-norm (in one space dimension). Thus, the idea is to see if these stability properties remain invariant when high-order degree polynomials are used.

To address this issue as clearly as possible, let us restrict ourselves to the one dimensional case. In this case, the forward Euler step for the method reads

I j (w h -u h ) δt v h dx - I j f (u h ) (v h ) x dx + f (u h ) v h x j+1/2 x j-1/2 = 0.
Taking v h ≡ 1, we obtain,

w j -u j /δ + f u - j+1/2 , u + j+1/2 -f u - j-1/2 , u + j-1/2 /Δ j = 0, (6)
where u j denotes the mean of u h on the interval I j , u - j+1/2 denotes the limit from the left and u + j+1/2 the limit from the right. When the approximate solution is piecewise-constant, we obtain a monotone scheme for small enough values of | δ | and, as a consequence, we do have that the scheme is total variation diminishing (TVD), that is, that

| w h | T V (0,1) ≤ | u h | T V (0,1) , where | u h | T V (0,1) ≡ 1≤j≤N | u j+1 -u j |
is the total variation of the local means. For approximate solutions that are not necessarily piecewise-constant, the above result still holds provided that the following sign conditions are satisfied:

sign u + j+1/2 -u + j-1/2 = sign u j+1 -u j , sign u - j+1/2 -u - j-1/2 = sign u j -u j-1 .
Since these conditions are not necessarily satisfied it is necessary to enforce them by means of what will be called a generalized slope limiter, ΛΠ h . We see that, unlike what happened in the previous section, the introduction of the numerical traces is not enough to guarantee the stability of the DG method. For non-linear hyperbolic problems, the use of a generalized slope limiter is indispensable, as was originally shown for the so-called high-resolution methods. See also [START_REF] Cockburn | Devising discontinuous Galerkin methods for non-linear hyperbolic conservation laws[END_REF] for a motivation of the introduction of this operator.

The generalized slope limiter

For piecewise-linear functions v h | I j = v j + (xx j ) v x,j , we define u h = ΛΠ 1 h (v h ), [START_REF] Osher | Convergence of generalized MUSCL schemes[END_REF], as follows:

u h | I j = v j + (x -x j ) m v x,j , v j+1 -v j Δ j /2 , v j -v j-1 Δ j /2 ,
where the minmod function m is defined by Note that this projection is non-linear and can be rewritten as follows:

u - j+1/2 = v j + m v - j+1/2 -v j , v j -v j-1 , v j+1 -v j (7) u + j-1/2 = v j -m v j -v + j-1/2 , v j -v j-1 , v j+1 -v j . ( 8
)
Next, we define a generalized slope limiter ΛΠ h for general discontinuous functions. Let us denote by v 1 h the L 2 -projection of v h into the space of piecewise-linear functions. We then define u h = ΛΠ h (v h ) on the interval I j , as follows:

(i) Compute u - j+1/2 and u + j-1/2 by using ( 7) and ( 8). (ii

) If u - j+1/2 = v - j+1/2 and u + j-1/2 = v + j-1/2 set u h | I j = v h | I j . (iii) If not, take u h | I j equal to ΛΠ 1 h (v 1 h ).
Since the above generalized slope limiter enforces the sign conditions, the forward Euler step is stable, that is, we have the following result.

Proposition 1 (The TVBM property). Suppose that for j = 1, . . . , N

| δ | | f (a, •) | Lip Δ j+1 + | f (•, b) | Lip Δ j ≤ 1/2.
Then, if u h = ΛΠ h v h , for some v h , then | w h | T V (0,1) ≤ | u h | T V (0,1) .

The stability of the RKDG method

For this method, we have the following stability result.

Theorem 2 (TVBM-stability of the RKDG method). Let each time step Δt n satisfy the following CFL condition:

max il β il α il Δt n | f (a, •) | Lip Δ j+1 + | f (•, b) | Lip Δ j ≤ 1/2. ( 9
)
Then we have

| u n h | T V (0,1) ≤ | u 0 | T V (0,1) ∀ n ∈ N.

Extensions

We must point out that with the above generalized slope limiter, there is loss of accuracy when the exact solution displays critical points. It is possible, however, to modify the above limiter to completely overcome this difficulty. It is also possible to extend the above results to the multi-dimensional case if we use the L ∞ -norm of the local means instead of their total variation. Finally, let us point that, although there are no stability proofs for DG methods applied to non-linear hyperbolic problems, the methods work extremely well.

To show this, we show some computational results from [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin finite element method for conservation laws V: Multidimensional systems[END_REF]. We display the contours of the density for the so-called double Mach reflection problem. We see that both the strong shocks as well as the contacts are well approximated and that no spurious oscillations appear when we go from linear to quadratic approximations. Moreover, the contact discontinuities seem to be better approximated by using quadratic polynomials. In fact, as argued in [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin finite element method for conservation laws V: Multidimensional systems[END_REF], even though the use of higher degree polynomials entails a more restrictive CFL condition, the enhanced quality of the approximation more than off-sets the increase of cost per mesh point.

Concluding remarks

In this paper, we have displayed the main ingredients (discontinuous approximations, element-by-element Galerkin weak formulations, and the numerical traces) and properties (local conservativity, high-order accuracy, and dissipation or stabilization through the jumps) of the DG methods as applied to a wide variety of problem ranging from ODEs to non-linear hyperbolic systems. We have also shown some of the particularities of these methods when applied to different problems: High accuracy and parallelizability for hyperbolic problems, easy elimination of the auxiliary unknowns for elliptic problems, and robustness with respect to the Reynolds numbers for the Oseen problem.

In [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF], some of the challenges that the development of these methods present are described. Let us end this paper by pointing out that a very important issue that has not been touched in this paper is the ease with which these methods can handle hp-adaptive strategies [16-18, 43, 44, 47-50, 55-58, 74, 75] and can be coupled with already existing methods [START_REF] Alotto | Discontinuous finite element methods for the simulation of rotating electrical machines[END_REF][START_REF] Cockburn | Approximation of the velocity by coupling discontinuous Galerkin and mixed finite element methods for flow problems, Special issue: Locally Conservative Numerical Methods for Flow in Porous Media[END_REF][START_REF] Perugia | On the coupling of local discontinuous Galerkin and conforming finite element methods[END_REF]. This is already having an important impact on the development of hp-adaptive algorithms and will also facilitate the handling of multi-physics models.

Fig. 2

 2 Fig. 2 Scaled L 2 -errors in u with k = 2 for different Reynolds numbers.

Fig. 3

 3 Fig. 3 Example of slope limiters: The MUSCL limiter (left) and the less restrictive ΛΠ 1 h limiter (right). Displayed are the local means of u h (thick line), the linear function u h in the element of the middle before limiting (dotted line), and the resulting function after limiting (solid line).

m (a 1 , a 2 , a 3 s

 123 min 1≤n≤3 | a n | if s = sign(a 1 ) = sign(a 2 ) = sign(a 3 ), 0 otherwise.

Fig. 4

 4 Fig. 4 Euler equations of gas dynamics: Double Mach reflection problem. Isolines of the density around the double Mach stems. Quadratic polynomials on squares Δx = Δy = 1/240 (top); linear polynomials on squares Δx = Δy = 1/480 (middle); and quadratic polynomials on squares Δx = Δy = 1/480 (bottom).

Table 1

 1 Scaled parallel efficiency. Solution times (without I/O) and total execution times measured on the nCUBE/2.

	Number of Work (W ) Solution time Solution parallel Total time Total parallel
	processors		(secs.)	efficiency	(secs.)	efficiency
	1	18,432	926.92	-	927.16	-
	2	36,864	927.06	99.98 %	927.31	99.98 %
	4	73,728	927.13	99.97 %	927.45	99.96 %
	8	147,456	927.17	99.97 %	927.58	99.95 %
	16	294,912	927.38	99.95 %	928.13	99.89 %
	32	589,824	927.89	99.89 %	929.90	99.70 %
	64 1,179,648	928.63	99.81 %	931.28	99.55 %
	128 2,359,296	930.14	99.65 %	937.67	98.88 %
	256 4,718,592	933.97	99.24 %	950.25	97.57 %

  by the positivity property (ii),

	i-1
	≤	α il u	(l) h	,	by the stability assumption,
	l=0
	≤ max 0≤l≤i-1	u	(l) h	,	by the consistency property (iii),
	which readily implies that | u n h | ≤ | u 0 h |, ∀ n ≥ 0, by a simple induction argument. In other words, the stability of the Euler forward step u (l)

Acknowledgements

The author would like to thank Prof. Dr. Ronald H. W. Hoppe and Prof. Dr. H.-J. Bungartz for the kind invitation to deliver a plenary lecture at the Annual Meeting of GAMM at Augsburg in March of 2002. The content of paper is based on that lecture. The author would also like to thank Ryuta Suzuki for discussions leading to a better presentation of the material in this paper.


Professor Firdaus E Udwadia
email: fudwadia@usc.edu

Aaron D Schutte
email: schutte@usc.edu

An Alternative Derivation of the Quaternion Equations of Motion for Rigid-Body Rotational Dynamics

This note provides a direct method for obtaining Lagrange's equations describing the rotational motion of a rigid body in terms of quaternions by using the so-called fundamental equation of constrained motion.

Introduction

In this note, a general formulation for rigid body rotational dynamics is developed using quaternions, also known as Euler parameters. The use of quaternions is especially useful in multibody dynamics when large angle rotations may be involved since their use does not cause singularities to arise, as it occurs when using Euler angles.

The equations of rotational motion in terms of quaternions appear to have been first obtained by Nikravesh et al. ͓1͔. The authors in this particular paper developed an approach to model and simulate constrained mechanical systems in which the components ͑bodies͒ in the mechanical system are connected by nonredundant holonomic constraints. Since they use unit quaternions to parametrize the angular coordinates of a rigid body, they provide many useful unit quaternion identities and show their relation to the components of the angular velocity of the rigid body. To develop a suitable set of equations of motion involving quaternions, they utilize Lagrange's equation. Of course the components of the unit quaternion must be of a unit norm and they are not all independent, and so the unit norm requirement must be imposed as a constraint on the system. They employ the Lagrange multiplier method to deal with this characteristic, and in so doing, they arrive at a mixed set of ordinary differential equations and an algebraic equation, which comprises a differential algebraic equation ͑DAE͒. The accelerations are found by inverting the DAE mass matrix, and the Lagrange multiplier is explicitly found. The generalized torques in this framework are given in a four-vector format, and its relation to a physically applied torque in the bodyfixed coordinate frame is provided by computing the virtual work of a force located at an arbitrary point of the rigid body and utilizing quaternion identities.

In a later paper, Morton ͓2͔ obtains both Hamilton's and Lagrange's equations in terms of quaternions. However, neither of these sets of equations is developed solely through the use of Hamiltonian or Lagrangian dynamics, as one might have expected. Morton's starting point for each of these sets of equations-the backbone for all of his derivations-is the Newtonian equations of rotational motion that describe the rate of change in the angular momentum in terms of quaternions ͑Eqs. ͑49͒ and ͑57͒ in Ref. ͓2͔͒. To obtain Hamilton's equations, Morton compares the rate of change in the angular momentum, expressed in terms of quaternions using Newton's equations and what is obtained using the Hamiltonian formalism ͑Eqs. ͑57͒ and ͑62b͒ are compared in Morton's paper ͓2͔͒. He thus obtains, through this somewhat long and involved procedure, the crucial connection between the physically applied torque vector and the generalized quaternion torques that are necessary to complete his Hamilton's equations. In deriving Lagrange's equations, Morton likewise starts with the Newtonian equations of motion in terms of quaternions ͑Eq. ͑49͒ in Morton's paper ͓2͔͒ and compares these equations with what is obtained by carrying out the Lagrange derivative of the kinetic energy of rotation ͑Morton's Eqs. ͑79͒ and ͑82͒ ͓2͔͒. Results obtained through comparisons between the Hamiltonian formulation and the Newtonian formulation, regarding the connection between the physically applied torque and the quaterion torque, are then utilized in the Lagrange derivative to identify the Lagrange multiplier that is needed to enforce the constraint that the quaternion must have a unit norm. Hence, the derivation of the Lagrangian equations of rotational motion in terms of quaternions appears to be carried out via a mixed line of thinking-some of it Newtonian, some Hamiltonian, some Lagrangian,-and the final Hamilton and Lagrange equations are both obtained basically through a comparison with those obtained by the Newtonian approach.

Having described in detail the two principal methods used to date in arriving at the requisite equations of motion, in this paper, we develop the Lagrange equations for rotational motion using a direct Lagrangian approach and follow a simple and straightforward line of reasoning rooted solely within the framework of Lagrangian dynamics. The development does not require the use of Lagrange multipliers, and it yields a positive definite mass matrix involving only the quaternion components. After the unconstrained equations of motion are obtained, the so-called fundamental equation is directly employed to get the final equations describing rotational motion. The ease, conceptual simplicity, and clarity, with which the equations are derived, are striking when compared with the derivation of Morton ͓2͔. No appeal to Newtonian mechanics is made, and no comparisons with results from it are made to arrive at the requisite equations. Of special importance is the simple derivation provided herein of the connection between the physically applied torque and the generalized quaternion torque. Unlike in Ref. ͓2͔, it is obtained completely within the framework of Lagrangian mechanics from simple arguments related to virtual work, and unlike in Ref. ͓1͔, it is achieved in a few simple steps that require less appeal to quaternion identities and algebraic manipulations. The approach used in the derivation is new and relies on some simple recent results dealing with the development of the equations of motion of constrained mechanical systems obtained by Udwadia and Kalaba ͓3,4͔. In addition, it provides new insights that were unavailable before.

Lagrange's Equations for Rigid Body Dynamics Using Quaternions

Consider a rigid body that has an angular velocity R 3 , with respect to an inertial coordinate frame. The components of this angular velocity with respect to a coordinate frame fixed in the body and whose origin is located at the body's center of mass are denoted by 1 , 2 , and 3 . We can express the four-vector = ͓0, 1 , 2 , 3 ͔ T in terms of quaternions by the relation

= 2Eu ˙= -2E ˙u ͑1͒
where the unit quaternion u ª ͓u 0 , u 1 , u 2 , u 3 ͔ T = ͓u 0 , u T ͔ T = ͓cos / 2,e T sin / 2͔ T , and the orthogonal matrix

E = ΄ u 0 u 1 u 2 u 3 -u 1 u 0 u 3 -u 2 -u 2 -u 3 u 0 u 1 -u 3 u 2 -u 1 u 0 ΅ ª ΄ u T --- E 1 ΅ ͑2͒
By a unit quaternion, we mean that

N͑u͒ ª u T u = 1 ͑3͒
Geometrically, e R 3 is any unit vector defined relative to an inertial coordinate frame and is the rotation about this unit vector.

Without loss of generality, we assume that the body-fixed coordinate axes are aligned along the principal axes of inertia of the rigid body whose principal moments of inertia are J i , i = 1,2,3. To formulate the equations of motion of the rotating rigid body, which may be subjected to a generalized torque four-vector ⌫ u ͑part of which may be derived from a potential͒, we begin by considering its kinetic energy

T = 1 2 T J ˆ = 1 2 T J = 2u ˙TE T JEu ˙= 2u T E ˙TJE ˙u ͑4͒
The inertia matrix of the rigid body is represented by J = diag͑J 1 ,J 2 ,J 3 ͒ and the 4 ϫ 4 diagonal matrix J = diag͑J 0 ,J 1 , J 2 ,J 3 ͒, where J 0 is any positive number.

We now proceed to obtain Lagrange's equation of motion for the rotating rigid body. We begin by assuming at first that each component of the four-vector u is independent. This assumption is tantamount to taking the generalized virtual displacements ␦u i , i = 0,1,2,3 to be all independent of one another. After we obtain these equations under this assumption, we will then impose on them the required unit quaternion constraint ͑Eq. ͑3͒͒. The constrained equations of motion are then the equations of rotational motion of the body. Thus, Lagrange's equation becomes

d dt ͩ ץT ץu ˙ͪ - ץT ץu = ⌫ u ͑5͒
The generalized quaterion torque four-vector ⌫ u is the torque that would exist if all the components of u were actually independent. Noting that ץT / ץu ˙=4E T JEu ˙, we obtain

d dt ͩ ץT ץu ˙ͪ = 4E T JEu ¨+ 4E ˙TJEu ˙+ 4E T JE ˙u ˙͑6͒ Also, since ץT ץu = ץ ץu 2u T E ˙TJE ˙u = 4E ˙TJE ˙u = -4E ˙TJEu ˙͑7͒
Lagrange's equation, assuming that the components of the unit quaternion are all independent, then becomes

4E T JEu ¨+ 8E ˙TJEu ˙+ 4E T JE ˙u ˙= ⌫ u ͑8͒
Using Eq. ͑2͒ and noting that E ˙u ˙= ͓ N͑u ˙͒ 0 3ϫ1 ͔ , the last member on the left-hand side of Eq. ͑8͒ can be simplified to

4E T JE ˙u ˙= 4͓u E 1 T ͔ ͫ J 0 0 0 J ˆͬͫ N͑u ˙͒ 0 3ϫ1 ͬ = 4J 0 N͑u ˙͒u ͑9͒
where N͑u ˙͒ = u ˙Tu ˙. Hence, we can rewrite Eq. ͑8͒ as

4E T JEu ¨+ 8E ˙TJEu ˙+ 4J 0 N͑u ˙͒u = ⌫ u ͑10͒
Note that the matrix M u =4E T JE multiplying the four-vector u ¨in Eq. ͑10͒ is symmetric and positive definite. In order to not distract us from our line of thinking, we shall later on show how we express the generalized quaterion torque four-vector ⌫ u in terms of the physically applied torque three-vector, ⌫ ˆB = ͓⌫ 1,B ⌫ 2,B ⌫ 3,B ͔ T , applied to the body, where ⌫ 1,B , ⌫ 2,B , and ⌫ 3,B denote the components of the physically applied torque about the body-fixed axes.

The equation of motion in Eq. ͑10͒ presumes that the components of the four-vector u are all independent of each other, which is indeed not the case in actuality, because these four components are constrained by Eq. ͑3͒. The explicit equation of motion of the constrained system is now written directly using the fundamental equation of constrained motion ͓3͔. To do this, we first differentiate the constraint equation ͑3͒ twice with respect to time to yield

A u u ¨ª u T u ¨= -N͑u ˙͒ ª b u ͑11͒
The explicit equation of motion is then directly and simply given by ͓4͔

4E T JEu ¨+ 8E ˙TJEu ˙+ 4J 0 N͑u ˙͒u = ⌫ u + M u 1/2 ͑A u M u -1/2 ͒ + ͑b u -A u a u ͒ ͑12͒
where a u is the acceleration of the unconstrained system obtained from Eq. ͑10͒ and is given by

a u = E T J -1 Eͩ ⌫ u 4 -2E ˙TJEu ˙-J 0 N͑u ˙͒uͪ ͑13͒
The superscript "+" in Eq. ͑12͒ denotes the Moore-Penrose inverse of the matrix A u M u -1/2 . We next simply proceed to assemble the various quantities required to determine the last member on the right-hand side of Eq. ͑12͒; this is done by computing A u a u , then b u -A u a u , and lastly M u -1/2 ͑A u M u -1/2 ͒ + . We begin by noting that

E T J -1 E = ͓u E 1 T ͔ ΄ 1 J 0 0 0 J ˆ-1 ΅ ͫ u T E 1 ͬ= uu T J 0 + E 1 T J ˆ-1 E 1 ͑14͒
From Eq. ͑2͒, we find that Eu = ͓ u T E 1 ͔u=͓ u T u E 1 u ͔=͓ 1 0 3ϫ1 ͔, so that using Eq. ͑14͒, we obtain the relations

E T J -1 Eu = u J 0 , u T E T J -1 E = u T J 0 , and u T E T J -1 Eu = 1 J 0 ͑15͒
where we have used the fact that N͑u͒ = 1. Hence, we find that

A u a u = u T E T J -1 Eͩ ⌫ u 4 -2E ˙TJEu ˙-J 0 N͑u ˙͒uͪ = u T ⌫ u 4J 0 -2 u T E ˙TJEu J0 -N͑u ˙͒ = u T ⌫ u 4J 0 + T J 2J 0 -N͑u ˙͒ = u T ⌫ u 4J 0 + T J ˆ 2J 0 -N͑u ˙͒ ͑ 16͒
where in the second and third equalities above, we have used Eqs. ͑1͒ and ͑15͒. This yields

b u -A u a u = -N͑u ˙͒ - u T ⌫ u 4J 0 - T J ˆ 2J 0 + N͑u ˙͒ = - u T ⌫ u 4J 0 - T J ˆ 2J 0 ͑17͒
Lastly, since the matrix A u = u T is a nonzero row vector, we have

M u 1/2 ͑A u M u -1/2 ͒ + = A u T ͑A u M u -1 A u T ͒ -1 = 4u͑u T E T J -1 Eu͒ -1 = 4J 0 u ͑18͒
where Eq. ͑15͒ is again used in the last equality. Using Eqs. ͑17͒ and ͑18͒, Eq. ͑12͒ reduces to

4E T JEu ¨+ 8E ˙TJEu ˙+ 4J 0 N͑u ˙͒u = ͑I -uu T ͒⌫ u -2͑ T J ˆ͒u ͑19͒ or alternatively to 4E T JEu ¨+ 8E ˙TJEu ˙+ 4J 0 N͑u ˙͒u = ͑I -uu T ͒⌫ u -2͑ T J͒u ͑20͒
In Eqs. ͑19͒ and ͑20͒, I is the 4 ϫ 4 identity matrix. The term 2͑ T J͒ on the right-hand side in Eq. ͑20͒ arises because of the constraint N͑u͒ = 1, and is twice the kinetic energy of rotation of the rigid body. Equation ͑20͒ also informs us-and this does not seem to be widely recognized-that the generalized torque must appear in the form ⌫ ˆu ª ͑Iuu T ͒⌫ u . Furthermore, because I = ͑Iuu T ͒ + uu T , postmultiplying this relation by ⌫ u , we find that

⌫ u = ͑I -uu T ͒⌫ u + u͑u T ⌫ u ͒ ͑ 21͒
Thus, u T ⌫ u is the component of the generalized torque four-vector ⌫ u in the direction of the unit vector u, and ⌫ ˆu is the orthogonal projection of ⌫ u in the plane normal to u.

The generalized acceleration u ¨is now explicitly obtained as

u ¨= E T J -1 Eͩ-2E ˙TJEu ˙-J 0 N͑u ˙͒u + ͑I -uu T ͒⌫ u 4 - ͑ T J ˆ͒u 2 ͪ ͑22͒ Since EE ˙T = 1 2 ͫ 0 -T ˜ͬ ͑23͒
where

˜= ΄ 0 -3 2 3 0 -1 -2 1 0 ΅ ͑24͒
the first term on the right-hand side of Eq. ͑22͒ can be expressed as

2E T J -1 EE ˙TJEu ˙= 1 2 E T J -1 ͫ 0 -T ˜ͬͫ 0 J ˆ ͬ = 1 2 ͩ-͑ T J ˆ͒u J 0 + E 1 T J ˆ-1 ˜J ˆͪ ͑25͒
By Eq. ͑15͒, the second and fourth terms on the right-hand side of Eq. ͑22͒ become

J 0 N͑u ˙͒E T J -1 Eu = N͑u ˙͒u ͑26͒
and

E T J -1 E ͑ T J ˆ͒u 2 = T J ˆ 2 E T J -1 Eu = ͑ T J ˆ͒u 2J 0 ͑27͒
respectively. Finally, using Eq. ͑15͒ again, the third term on the right-hand side of Eq. ͑22͒ simplifies to

E T J -1 E ͑I -uu T ͒⌫ u 4 = ͓E T J -1 E -͑E T J -1 Eu͒u T ͔⌫ u 4 = ͩE T J -1 E -uu T J 0 ͪ ⌫ u 4 = E 1 T J ˆ-1 E 1 ⌫ u 4 ͑28͒
where the last equality follows from Eq. ͑14͒. Thus, Eq. ͑22͒ simplifies to

u ¨= - 1 2 E 1 T J ˆ-1 ˜J ˆ -N͑u ˙͒u + E 1 T J ˆ-1 E 1 ⌫ u 4 ͑29͒
Note the arbitrary positive scalar J 0 has vanished, and we see that Eq. ͑29͒ is obtained directly in a simple straightforward manner. In Eq. ͑29͒, u and ⌫ u are four vectors, is a three-vector, and E 1 is the 3 ϫ 4 matrix defined in Eq. ͑2͒.

Our final task is to express the generalized torque four-vector ⌫ ˆu = ͑Iuu T ͒⌫ u on the right-hand side of Eq. ͑20͒, in terms of the physically applied torques on the body. Consider the four-vector

⌫ B = ͓0 ͉ ⌫ ˆB T ͔ T = ͓0, ⌫ 1,B , ⌫ 2,B , ⌫ 1,B ͔ T ͑30͒
The virtual work done under an infinitesimal, virtual rotation ␦ = ͓0,␦ 1 , ␦ 2 , ␦ 3 ͔ T , in which ␦ 1 , ␦ 2 , and ␦ 3 are the virtual "displacements" about the body-fixed axes, is then given by ⌫ B T ␦.

On the other hand, the virtual work done by the generalized quaternion torque vector ⌫ ˆu is given by ⌫ ˆu T ␦u, where the virtual displacement of the corresponding quaternion is ␦u. Equating these two virtual work expressions, we obtain

⌫ B T ␦ = ⌫ ˆu T ␦u = ⌫ u T ͑I -uu T ͒␦u ͑31͒
Here, we set = ͓0, ˙1 , ˙2 , ˙3͔ T , where ˙1, ˙2, and ˙3 are the components of the angular velocity about the body-fixed axes, and by Eq. ͑1͒, we have

΄ 0 ˙1 ˙2 ˙3 ΅ = 2E ΄ u ˙0 u ˙1 u ˙2 u ˙3 ΅ ͑32͒ so that ␦ = 2E␦u ͑33͒
However, since u T u ˙= 0, the virtual displacement ␦u in Eq. ͑33͒ must satisfy the relation

u T ␦u = 0

͑34͒

Solving for ␦u in Eq. ͑34͒, we get ͓5͔ ␦u = ͑Iuu T ͒␦w ͑35͒

where the four-vector ␦w is any arbitrary, infinitesimal column vector. In view of Eqs. ͑33͒ and ͑35͒, we find that Eq. ͑31͒ becomes

⌫ B T 2E͑I -uu T ͒␦w = ⌫ u T ͑I -uu T ͒␦w = ⌫ ˆu T ␦w ͑36͒ because ͑I -uu T ͒ is idempotent. Since ␦w is arbitrary, we have 2͑I -uu T ͒E T ⌫ B = ͑I -uu T ͒⌫ u = ⌫ ˆu ͑37͒
The left-hand member in Eq. ͑37͒ can be simplified and written as

2͑I -uu T ͒E T ⌫ B = 2E T ⌫ B = 2͓u E 1 T ͔ ͫ 0 ⌫ ˆB ͬ = 2E 1 T ⌫ ˆB ͑38͒
where the first equality follows because u T E T ⌫ B = ͓1 ͉ 0 1ϫ3 ͔ ϫ͓0 ͉ ⌫ ˆB T ͔ T = 0, and the last from u T E 1 T = 0. Hence, relation ͑37͒ becomes

⌫ ˆu = 2E 1 T ⌫ ˆB = 2E T ⌫ B ͑39͒
Thus, the equation of motion ͑20͒ describing the rotational motion of the system becomes

4E T JEu ¨+ 8E ˙TJEu ˙+ 4J 0 N͑u ˙͒u = 2E T ⌫ B -2͑ T J͒u ͑40͒
where ⌫ B is the four-vector containing the components of the body torque about the body-fixed axes. Furthermore, using the first equality in Eq. ͑39͒, premultiplication by E 1 yields

⌫ ˆB = 1 2 E 1 ⌫ ˆu = 1 2 E 1 ͑I -uu T ͒⌫ u = 1 2 E 1 ⌫ u ͑41͒ since E 1 E 1 T = I. Using Eq. ͑41͒ in Eq. ͑29͒, we get u ¨= - 1 2 E 1 T J ˆ-1 ˜J ˆ -N͑u ˙͒u + E 1 T J ˆ-1 ⌫ ˆB 2 ͑42͒
Equations ͑40͒ and ͑42͒ represent the requisite Lagrange equations that describe the rotational motion of a rigid body in terms of quaternions and the applied torques in the body-fixed coordinate frame.

We end with some observations and remarks, which provide additional insights.

1. We observe that Eq. ͑40͒ is the same as Eq. ͑12͒, for in getting to Eq. ͑40͒, all we have done is compute the various entities on the right-hand side of Eq. ͑12͒. From a computational stand-point, we could just as well as have used Eq. ͑12͒ along with Eq. ͑39͒ directly. 2. Were we to have replaced ⌫ u in Eqs. ͑5͒ and ͑10͒ by ͑I uu T ͒⌫ u , we would, accordingly, need to replace the term ⌫ u on the right-hand side of Eq. ͑20͒ also by ͑Iuu T ͒⌫ u . But since ͑Iuu T ͒ is idempotent, this leaves Eq. ͑20͒ unchanged. This shows that rotational dynamics only involves the component ͑Iuu T ͒⌫ u of the generalized quaternion torque ⌫ u . Alternately, we conclude that the component of the fourvector ⌫ u along u, namely, u T ⌫ u , does not play any role in the rotational dynamics of a rigid body ͑see Eq. ͑21͒͒. 3. The simplicity of the approach developed herein becomes apparent were we to be mainly interested in computing the rotational response of a rigid body subjected to an impressed torque. Then, from the abovementioned remarks, we see that instead of Eq. ͑10͒, one could have started with the unconstrained equation of motion given by

M u u ¨ª 4E T JEu ¨= -8E ˙TJEu ˙-4J 0 N͑u ˙͒u + 2E T ⌫ B ͑43͒
where ⌫ B is the physically impressed torque four-vector on the rigid body. To arrive at Eq. ͑43͒, we have replaced ⌫ u in Eq. ͑10͒ by ͑Iuu T ͒⌫ u and used Eq. ͑39͒. Let us denote the right-hand side of Eq. ͑43͒ by

Q ª -8E ˙TJEu ˙-4J 0 N͑u ˙͒u + 2E T ⌫ B ͑44͒
Direct application of the fundamental equation now gives the rotational equation of motion of a rigid body ͑subjected to the impressed torque ⌫ B ͒ in a single step as

4E T JEu ¨= Q + M u 1/2 ͑A u M u -1/2 ͒ + ͑b u -A u M u -1 Q͒ ͑45͒
where the various quantities are defined as before, and Q is given in Eq. ͑44͒. For computational purposes, the righthand side of Eq. ͑45͒ can be numerically computed directly, and there is no need to simplify it any further as we did before in going from Eqs. ͑12͒-͑40͒.

Conclusion

This note provides a simple and direct route for obtaining Lagrange's equation describing rigid body rotational motion in terms of quaternions. The derivation is carried out in a uniform manner without any appeal to the notion of Lagrange multipliers or Newtonian mechanics by using the fundamental equation of constrained motion. A new and simple way of establishing the connection between the physically applied torque and the generalized quaternion torque is presented. Besides providing new insights hereto unavailable, this derivation reveals the explicit Lagrange's equation in a very transparent manner, an aspect previous derivations have had difficulty attaining.