- Email: [email protected]
Spectral penalty methods
Applied Numerical Mathematics 33 (2000) 23–41
Spectral penalty methods ✩ J.S. Hesthaven 1,2 Division of Applied Mathematics, Brown University, Providence, RI 02912, USA
Abstract We present an overview of the theory and applications of spectral penalty methods for the stable solution of initial boundary value problems. Through a number of case studies we shall illustrate that imposing boundary conditions weakly through a penalty term rather than strongly as is classically done, removes many of the problems traditionally associated with the construction and analysis of stable and accurate pseudospectral approximations. The discussion includes the treatment of complex boundary operators, the construction of integration preconditioners and the formulation of spectral methods on general unstructured grids, highlighting the versatility and flexibility introduced by the elegant splitting of the operator and the boundary conditions exploited in the penalty methods. 2000 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Spectral methods; Pseudospectral methods; Stability; Boundary conditions
1. Introduction The analysis of numerical schemes for the solution of partial differential equations is often complicated significantly by the presence of boundary conditions, which often makes it difficult to analyze such basic properties as stability and convergence for even simple cases, and renders the analysis of more realistic problems and complex boundary operators essentially impossible. This is particularly true when considering high-order or spectral methods. However, the recent increasing interest in the modeling of transient problems in science and engineering puts an emphasis on the use of such high-order methods and, in particular, schemes for which numerical stability can be established in a rigorous manner. In a traditional approach, the value of a numerical scheme for the solution of an initial boundary value problems is often judged through an analysis of the complete scheme, including the boundary conditions, after it has been formulated. This approach often implies that the analysis must be done on ✩
Expanded version of invited talk presented at the International Conference on Spectral and High Order Methods 1998 (Herzliya, Israel). 1 E-mail: [email protected] 2 Supported by DARPA/AFOSR grant F49620-96-1-0426, by DOE grant DE-FG02-95ER25239 and by NSF grant ASC9504002. 0168-9274/00/$20.00 2000 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 9 9 ) 0 0 0 6 8 - 9
24
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
a problem by problem basis and it is often complicated significantly by the modifications of the discrete operators introduced by imposing the boundary conditions exactly. These considerations are particularly true when one considers high-order and spectral methods where the existing stability theory is far from being complete. A viable approach, known as penalty methods in the context of spectral methods, towards a resolution of these problems has attracted increasing interest in recent years. The central idea behind the penalty method lies in the observation that it is only necessary to satisfy the boundary conditions asymptotically in a way consistent with the order of the scheme, i.e., there is no reason to insist on the boundary conditions being exactly obeyed if one is willing to accept a certain overall error. This approach makes the implementation of complex boundary conditions straightforward and it provides an elegant splitting between the operator and the boundary conditions such that the former maintains its delicate structure while the latter is dealt with in a simple manner. Moreover, this splitting enables the analysis of very complex schemes and provides a constructive approach to the construction of stable spectral methods for the solution of initial boundary value problems. What remains in this paper is organized as follows. In Section 2 we provide a very selective overview of the central properties of Legendre polynomials and polynomial methods. An understanding of these issues is essential to Section 3 where we present an overview of the basic elements of spectral penalty methods in different variations. This sets the stage for Section 4 where we show how the use of penalty methods can be helpful, or even essential, in resolving a number of known problems related to the use and implementation of spectral methods. We consider, as examples only, the treatment of complex boundary conditions and the associated problem of how to construct stable multi-domain schemes, the formulation of efficient preconditioners for pseudospectral operators and the construction of stable spectral methods on unstructured grids, most notably general grids within triangles and tetrahedra. Section 5 contains a few general remarks and draws the attention to recent parallel developments in alternative high-order methods, suggesting that penalty methods provide a general and convenient framework in which to formulate and analyze stable high-order methods for the solution of initial boundary value problems. 2. Polynomials and polynomial methods While the formulation of the penalty methods is independent of the actual choice of the polynomial basis, providing the basis of the spatial approximation, the discussions and subsequent analysis are most easily carried out in the framework of Legendre methods. For the sake of simplicity we shall, therefore, focus our attention on such approximations, keeping in mind that this does not pose a general restriction. We primarily concern ourselves with Legendre methods utilizing the discrete approximation IN u(x) =
N X n=0
uen Pn (x),
N 1 X uen = u(xi )Pn (xi )wi , γen i=0
(
γen = 2
(2n + 1)−1 , n < N, N −1 , n = N,
where Pn (x) represents the Legendre polynomial of order n, and (xi , wi ) signifies the set of Legendre– Gauss–Lobatto nodes and weights given as dPN 2 = 0, grid points: 1 − xi2 weights: wi = . dx xi N(N + 1)PN2 (xi ) It should be noted that the choice of the Gauss–Lobatto quadrature does not represent a constraint and alternative Gauss type quadratures equally well could be used. However, for solving boundary value
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
25
problems, the Gauss–Lobatto quadratures is the most natural choice as it includes the boundaries in the nodal set. As we shall see shortly, it is central to recall that the Gauss–Lobatto quadrature is exact as
∀u(x) ∈ P2N−1 = span Pn (x)
2N−1 n=0
Z1
:
u(x) dx =
N X
u(xi )wi .
(1)
i=0
−1
An immediate consequence of the use of the Gauss quadratures is IN u(xi ) = u(xi ), i.e., the approximation is an interpolation IN u(x) =
N X
uen Pn (x) =
n=0
N X
u(xi )li (x),
i=0
where li (x) represents the N th order Lagrange polynomial based on xi . Hence, there are two mathematically equivalent but computationally different ways of expressing the interpolation. This again implies that the central operation of computing derivatives can be done in two ways: N N X X dIN u ? 0 e = u Pn (xi ) = u(xj )Dij , dx xi n=0 n j =0
(2)
utilizing either the expansion coefficients, uen , directly or exploiting the interpolation property to derive the differentiation matrices, D, in physical space as in finite difference methods. We recall that ue0n = (2n + 1)
N X
uep =
p=n+1 p+n odd
N X
e np u ep , D
p=0
dlj (x) Dij = . dx xi
?
The = in Eq. (2) is meant as a reminder that we have yet to understand the exact relation between the two formulations. If we introduce
T
u = u(x0 ), . . . , u(xN ) ,
T
e0 , . . . , u e= u eN , u
T
e00 , . . . , u e0 = u e0N , u
e and u e 0 , respectively, we may express Eq. (2) as as the vectors of the nodal values, u, and modal values, u
du = Du. dx As a direct result of the Gauss quadrature, Eq. (1), we recover that eu e0 = D e, u
e , D = T −1 DT
i.e., the two formulations are related through a simple similarity transform as Tij =
1 Pi (xj )wj , γei
Tij−1 = Pj (xi ),
are orthogonal transformation matrices.
26
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
e is obtained directly from the polynomials [7] The entries of D
d −1 e −1 d −1 D , T ⇒ T −1 = T dx dx which represents a well known recurrence relation e= T −1 D
0 0 (2n + 1)Pn (x) = Pn+1 (x) − Pn−1 (x).
e −1 , is tri-diagonal and well-conditioned with a purely imaginary spectrum Hence, the pseudoinverse, D and the magnitude of the maximum eigenvalue is uniformly bounded by 1/π .
3. A short course on penalty methods Let us now return to the central issue of solving initial boundary value problems. We restrict, solely for simplicity, the attention to the construction of a spectral method for the solution of the scalar waveequation ∂u ∂u = , u(1, t) = f (t), (3) ∂t ∂x with x ∈ [−1, 1]. To construct the scheme we have to make a number of choices, the first one being in which space to seek solutions. In the following, we shall seek solutions of the form uN (x, t) =
N X n=0
en Pn (x) = u
N X
u(xi )li (x),
i=0
where uen are the discrete Legendre–Gauss–Lobatto expansion coefficients based on the quadrature points, xi , on which also the Lagrange polynomial, li (x), is based. The second choice to make is in which way we require the equation to be satisfied using the methods of weighted residuals. With proper choices of the test functions, one can then recover the classic Galerkin, tau and collocation methods. In the latter case one should keep in mind that there is no requirement on which points one requires the equation to be satisfied. In particular, these points need not be the same as those on which the polynomial approximation is based. As we shall see, this latter observation opens up for the construction of a number of interesting schemes for the solution of initial boundary value problems. Regardless of the specifics of the method, however, the overarching problem of all the schemes is to impose the boundary conditions in a consistent and stable manner while maintaining the accuracy of the scheme. The Galerkin and tau methods suffer from well-known difficulties for nonlinear problems while the issue of stability and discreteness of the approximation at the boundary is a concern for collocation approximations. To address these issues, we consider the following approximation of Eq. (3): ∂uN ∂uN (4) = − τ Q(x) uN (1, t) − f (t) ∂t ∂x as was first proposed in [12,13] within the context of spectral methods for evolutionary problems. A very similar approach for elliptic problems was proposed prior to that in [10,11]. A related technique has been proposed previously in the context of finite element methods as a way of enforcing global constraints
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
27
[23,28], although this was later rejected due to excessive stiffness being introduced through this global penalization. In what remains we shall restrict our attention to situations with penalization of boundary constraints only, as this is where the advantages of using the penalty methods are highlighted, and we shall, moreover, consider mainly problems with a temporal component. Eq. (4) is solved everywhere in the domain, τ plays the role of a free parameter, or a Lagrange multiplier, and Q(x) is a polynomial, that, as we shall see, characterizes the scheme. Imposing the boundary conditions weakly through a penalty term, as in Eq. (4), introduces additional degrees of freedom in the construction of the scheme and facilitates the formulation of entirely new families of schemes. We emphasize that for τ approaching infinity, the boundary condition must become strongly imposed, hence establishing consistency by equivalence with the traditional formulation. For τ being finite, the equation as well as the boundary condition is enforced at the boundary. In completing the specification of the scheme, we need to choose Q(x) and we shall moreover need to specify a constraint that allows for the specification of τ . As we shall see, τ is most often specified under the constraint of asymptotic stability in some norm but one should keep in mind that this choice is not unique, e.g., one could also require conservation as a constraint. In considering a few examples, let us begin with the classical choice of Q(x) =
(1 + x)PN0 (x) = 2PN0 (1)
1, x = 1, 0, x 6= xN ,
and require that Eq. (4) be satisfied in a collocation sense at the Legendre–Gauss–Lobatto points, i.e., the same points being used in the construction of the approximation. We note, in particular, that this choice results in solving Eq. (3) in the interior in a way similar to a traditional scheme while only the treatment of the boundary point is somewhat different from that of a classical approach. To complete the specification of the scheme, we require that the approximation be asymptotically stable in the discrete L2 -norm as
N X 1 d d uN (xi ) kuN k2N = uN (xi ) wi − τ u2N (xN )wN 6 0. 2 dt dx i=0
Utilizing quadrature and integration by parts we recover 1 d 1 1 kuN k2N = u2N (xN ) − τ u2N (xN )wN − u2N (x0 ) 6 0, 2 dt 2 2 and stability provided τ>
1 N(N + 1) > . 2wN 4
Hence, a lower finite bound on τ is established which guarantees that the semi-discrete approximation is asymptotically stable by construction. Moreover, if u(x) can be assumed smooth, the magnitude of the penalty term will decay exponentially and the scheme maintains the overall spectral accuracy. In fact, early results in [12] suggest that the penalty methods is superior to the standard method of imposing boundary conditions in terms of accuracy. While the enhanced accuracy and inherent stability appear as the main motivation of the original work [12,13] it is by considering the versatility of the weakly imposed boundary operators that one begins to appreciate the general approach.
28
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
Rather than choosing the Q(x) associated with the collocation scheme, one could consider using 0 PN+1 (x) + PN0 (x) 2 in Eq. (4) and ask that the equation be satisfied in a Galerkin sense. We shall refrain from writing the scheme directly, as it involves mass- and stiffness matrices, but rather determine τ by requiring that the scheme be stable in the L2 inner product with stability resulting for τ > 1/2. Thus, we seek polynomial solutions that satisfy the equation
Q(x) =
∂uN ∂uN P 0 (x) + PN0 (x) uN (1, t) − f (t) . = − N+1 ∂t ∂x 4 Note, however, that Q(x) is orthogonal to all polynomials of order N or less that satisfy the boundary conditions, i.e., the method is nothing else than the Legendre–Galerkin method for the wave-equation as could also have been realized by considering the Legendre–Galerkin error-equation which is identical to the penalty scheme given above, i.e., they yield the exact same polynomial solution. In the same spirit one can consider using Q(x) = PN (x), which, under the restriction of L2 stability, yields a lower bound τ > (2N + 1)/4 resulting in a stable scheme as ∂uN ∂uN 2N + 1 = − PN (x) uN (1, t) − f (t) . ∂t ∂x 4 By comparison with the error equation it becomes clear that this scheme is nothing else that a different formulation of the well-known Legendre tau method. As a final example, let us again consider the formulation of the collocation scheme, although we shall require that the approximation is based on the Chebyshev–Gauss–Lobatto nodes, yj , while the equation is sought satisfied at the usual Legendre–Gauss–Lobatto nodes, xj . As for the Legendre collocation method discussed above, we shall use
Q(x) =
(1 + x)PN0 (x) , 2PN0 (1)
and seek a solution to the equation ∂uN ∂uN (1 + x)PN0 (x) uN (1, t) − f (t) . = −τ 0 ∂t ∂x 2PN (1)
While the approximation itself is constructed using the Chebyshev nodes we shall require that the equation be satisfied at the Legendre nodes. Note that as a consequence of this, the penalty term has a global effect as it is non-zero at all the Chebyshev nodes. Nevertheless, since both sides of the equation represents a polynomial of order N specified at N + 1 nodes, it is clearly unique and we can hence obtain stability in the discrete Legendre inner product, similar to the classical case, resulting in the bound N(N + 1) . 4 Hence, stability of the Chebyshev method in the usual L2 -norm has been established using the Legendre nodes. From the perspective of computation, however, the Legendre nodes are never used. On the other hand, by comparison with the error-equation, the scheme is identical to a Legendre collocation method while the approximation is based on the nodes of Chebyshev polynomials. τ>
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
29
This latter method, known as the Chebyshev–Legendre method [9], serves as a good example of the wealth of new schemes that one arrives at by imposing the boundary condition weakly rather than strongly in combination with an increased understanding of exactly how one needs to require that the equations be satisfied. In recent years these central ideas have been further developed and have resulted in a wealth of new schemes by introducing and combining different choices of Q(x) and the way the equation is satisfied [4,8,19]. In all previous work, τ has been chosen to ensure stability in some energy norm. As a slight generalization, let us briefly consider the formulation of a penalty method for a hyperbolic system. For this purpose, we consider the problem
∂ u 0 1 ∂ u + = 0, ∂t v 1 0 ∂x v
u(−1) = f (t),
v(1) = g(t).
As the problem is hyperbolic we have
A = S −1 ΛS,
where
1
1
1 , S=√ 2 1 −1
Λ=
1
0
0 −1
,
as the characteristic decomposition, clarifying that one boundary condition is required at each boundary. The penalty scheme for this problem becomes (q N = [uN , vN ]T ) ∂q N ∂q + A N = −τ − Q− (x)SR − S T q N − f − − τ + Q+ (x)SR + S T q N − f + ∂t ∂x
where Q± (x) and τ ± take the usual meaning at the two boundaries and
R− =
1 0
,
0 0
R+ =
0 0 0 1
,
f− =
f (t) vN (−1)
,
f+ =
uN (1)
.
g(t)
Thus, the boundary conditions are imposed through the characteristic variables such that stability can be established for each decoupled equation of the individual characteristic variable. In other words, stability for the system follows directly from that of the scalar problem with straightforward extensions to non-linear problems through linearization and localization, requiring only that the solution is locally smooth, and to multi-dimensional problems for symmetrizable systems and conservation laws possessing an entropy-flux pair [13,15,18,19].
4. Beyond the basics While the basic ideas underlying the penalty method can be illustrated well through an analysis of the linear wave-equation, one has to look towards more general applications to fully appreciate the versatility and flexibility offered by this approach. In the following we shall present a number of such examples, illustrating how the penalty methods facilitate the construction of entirely new families of methods while allowing for elegant and simple analysis and solutions to old and well-known problems associated with the use of spectral methods.
30
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
4.1. Imposing complex boundary operators One of the central issues, and by many perceived as the central obstruction, in the use of spectral methods lies in the formulation of simple, yet stable, ways of imposing boundary conditions of a complicated form. To illustrate the advantages of using a penalty method for imposing complex boundary operators, let us consider the linear advection–diffusion equation [19], ∂u ∂u ∂2 u +λ = ε 2, ∂t ∂x ∂x αu(−1, t) − εβux (−1, t) = g1 (t), γ u(1, t) + εδux (1, t) = g2 (t), where (α, β, γ , δ) > 0 must be chosen to ensure well-posedness [19]. The question to address is how to impose boundary conditions of such general nature in a consistent and stable manner. While imposing the boundary conditions strongly requires modification of the operators, a process which becomes increasingly complex for more general problems and boundary operators, enforcing the conditions in a weak way only through a penalty term immediately yields ∂u ∂u ∂2 u +λ = ε 2 − τ − Q− (x) αu(−1, t) − εβux (−1, t) − g1 (t) ∂t ∂x ∂x − τ + Q+ (x) γ u(1, t) + εδux (1, t) − g2 (t) .
Here, and in what remains, we shall not distinguish between u and uN for simplicity as the correct meaning should be clear from the context. If we require that the equation be satisfied at the Legendre collocation points and take Q± (x) to be an N th order polynomial vanishing at all the interior points except at x = ±1, respectively, we may guarantee discrete L2 -stability as N−1 X 1 d u2x (xi )wi , kuk2N 6 −ε 2 dt i=1
provided only that the penalty parameters are bounded like − + τα,β 6 τ − 6 τα,β ,
where ± τa,b =
τγ−,δ 6 τ + 6 τγ+,δ ,
q i 1 h ε + 2κ ± 2 κ 2 − εκ − 12 εwN |λ| , wN εb
κ=
wN a . b
We note that the result on stability is independent on (α, β, γ , δ) and ε, i.e., it remains uniformly valid for vanishing viscosity. In [19] similar formulations for hyperbolic, parabolic and mixed problems are discussed, including stability bounds for all cases. To appreciate the accuracy and validity of the above approach, let us apply it to solve Burgers equation ∂u ∂u ∂2 u +u = ε 2, ∂t ∂x ∂x
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
31
Fig. 1. Traveling wave solution of Burgers equation computed using a Legendre collocation method employing a penalty formulation at the boundaries. Table 1 Error of the traveling wave solution using a Legendre collocation method with a penalty method at the boundary points N
L2
L∞
16
3.41E−03
3.26E−02
32
2.43E−05
3.50E−04
64
1.09E−09
2.21E−08
128
4.98E−12
7.62E−11
subject to the initial (exact) conditions
x u(x, 0) = − tanh + 1. 2ε In this case, the boundary operator of the non-linear form u2 (−1, t) − εux (−1, t) = g1 (t),
εux (1, t) = g2 (t),
ensures well-posedness and the solution represents a kink-solution traveling at a constant speed as illustrated in Fig. 1. Computing the accuracy of the solution yields the results shown in Table 1, confirming the global spectral accuracy of the scheme even though the boundary conditions are imposed only weakly. The emphasis so far has been on the construction of stable ways in which to impose open boundary conditions. However, as has been utilized in [2,14,15,18,21,26], the generalization to multi-domain solution of the problem is straightforward. Indeed, if we consider each domain as an isolated entity
32
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
Fig. 2. Viscous shock solution of Burgers equation obtained using a multi-domain Legendre collocation method with a penalty method for patching of the domains. The dots on the initial condition represents the collocation points in the 4 domain solution. In (b) we show the derivative of the solution shown in (a).
requiring boundary conditions at every time-step, the scheme developed above generalizes directly by simply passing information between the domains as prescribed through the dynamics of the equation, e.g., advective terms are always treated by upwinding. Hence, while the boundary conditions may not be available a priori in the individual domains they are available whenever needed by passing the required information from the neighboring domain. It should be noted that this line of argument, leading to a scheme local in time and space, assumes that information propagates with a finite speed as in wave-dominated problems, although accuracy and stability has been established and confirmed also for parabolic problems [2,15,18]. To illustrate the efficacy of this approach let us again consider the solution of Burgers equation subject to the boundary and initial conditions u(±1, t) = 0,
u(x, 0) = − sin(π x),
and with a small viscosity ε = 0.01/π . Rather than attempting to directly resolve the solution, which develops a very steep gradient around x = 0, we shall employ 4 subdomains of different size to accurately resolve the viscous shock. In Fig. 2 we show the computed solution, illustrating that even though the patching conditions, which are nothing else than the open boundary operator of Burgers equation ensuring well-posedness in L2 , are imposed only weakly the solution and its derivative remains smooth even for this very steep gradient (see [15] for further details). The results given in the above suggests a systematic way of developing stable schemes for the solution of initial boundary value problems. Given a particular problem, one first derives a set of open boundary conditions which ensures well-posedness of the total problem. In general, this set of boundary conditions is very complex. However, with the flexibility provided by the penalty method, these boundary conditions can be imposed in a stable and accurate way. Moreover, the generalization from a single domain solution to a multi-domain solution is straightforward and results in a global spectrally accurate, asymptotically stable and geometrically flexible formulation for the solution of the initial boundary value problems.
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
33
Fig. 3. Two-dimensional computation of recessed cavity flameholder flow at M = 0.9 and Re = 1500 based on the half channel width and the maximum mean flow. In (a) we show a segment of the spectral multi-domain grid, using order 12 in each domain and 156 domains for this problem, and in (b) we show the density and instantaneous streamlines, illustrating the recirculating flow within the cavity.
As an illustration of the versatility and robustness of the multi-domain framework outlined above, we show in Fig. 3 a simulation of the flow of compressible air inside a narrow grooved channel. The computation is performed using a multi-domain pseudospectral scheme for the solution of the unsteady compressible Navier–Stokes equations where the patching and treatment of open boundaries is done in the penalty method setting. For the compressible Navier–Stokes equations the boundary operator is extremely complex, containing components of the inviscid part from the Euler equations and elements of the normal and tangential stress and heat flux while being singular due to mass conservation, and it is difficult to construct a stable scheme in which this operator is imposed strongly. However, imposing the boundary operator weakly is straightforward and stability of the complete three-dimensional scheme has been established within a general curvilinear hexahedral, including the edges and vertices, with the scheme maintaining stability and validity for vanishing viscosity [14,15,18]. These results provide the fundamental elements for the construction of very general three-dimensional pseudospectral solvers for the solution of problems of gas-dynamics and supplies a convincing example of the efficacy of using the penalty method for imposing complex boundary operators. In recent years, penalty methods have been used to address a number of problems associated with complex boundary conditions and the construction of multi-domain schemes. In [2] a penalty formulation was developed for the heat-equation and implemented using a mixed explicit/implicit time integration scheme. A related approach was taken in [26] where a multi-domain scheme for the compressible Navier– Stokes equations is presented, utilizing a splitting of the inviscid and viscous part of equation and treating the former part by upwinding through characteristics while the dissipative part is included using a penalty formulation borrowed from earlier work on elliptic problems [10,11]. A stable scheme for the multidomain solution of the three-dimensional advection–diffusion equations is derived in [18] and in [21] a multi-domain scheme for the solution of certain problems in non-linear optics is proposed. An elegant application of the penalty method was recently presented in [29] where the basis idea was exploited to impose non-slip boundary conditions in a velocity–vorticity (v, ω) formulation of the
34
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
incompressible Navier–Stokes equations. The well-known lack of direct boundary conditions on the vorticity was addressed as ∂ω + ∇ × (ω × v) = ν∇ 2 ω − τ Q(x) ω − ∇ × v . ∂t Since the boundary conditions on the velocity field, v, are well known, the scheme imposes the unknown vorticity condition by consistency at the boundary of the computational domain. Its accuracy and stability has been confirmed through numerous experiments [29]. 4.2. Designing efficient preconditioners While the application of penalty methods for the stable and consistent treatment of complex boundary operators has received significant attention in recent years, there are a number of alternative uses of the penalty methods that seem very promising in resolving some classical problems associated with the use of spectral methods. As an example of this, let us consider the problem of preconditioning of dq u = f ⇒ Dq u = f , dx q subject to appropriate boundary conditions. For simplicity we assume that a Legendre collocation method is used and Dq signifies the qth order differentiation matrix while u and f represents the vectors of the nodal values of the solution and the forcing function, respectively. It is well known that the conditioning of Dq scales at best as O(N 2q ) where N is the order of the approximation. This result is not special to the Legendre method, or the type of collocation used, but is shared among all the classical polynomials and quadratures. Hence, if we hope to solve the above equation for a general force, f , we must consider the issue of preconditioning. Previous attention to this problem has resulted in the development of preconditioners based on a low-order approximation of the operator with examples being finite-difference and finite-element based approaches for advection and diffusion operators (see [17] for references). The efficiency of these methods varies significantly with the type of operator being considered and, in particular, the dimension of the problem. To understand the prospects of using the penalty method within this particular setting, let us consider preconditioning of the advective operator as Du = f ,
u(1) = 0.
(5)
We recall from Section 2 that e D = T −1 DT
such that
e −1 T , D−1 = T −1 D
e −1 is tri-diagonal, well conditioned and given in closed form. Hence, in the absence of the where D boundary conditions, the exact inverse, in some sense being the optimal preconditioner, is known. However, once we impose the boundary conditions strongly by altering the last row and column of D, the delicate structure of D is ruined and the construction of the inverse is no longer clear. Let us, therefore, consider the construction of a preconditioner for a collocation penalty formulation as [17]
D − τ INN u = f ,
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
35
where Iab is the diagonal matrix with unity in the rows from a to b, e.g., I0N represents the identity matrix. This suggests that since D is now unchanged the ideal preconditioner is given as e −1 T Pe−1 = T −1 D
since
Pe −1 D = T −1 I1N T .
e −1 having However, Pe−1 is singular since integration leaves a constant undetermined as manifested in D the first row being zero. To realize a solution to this problem we must recall that we, for all practical purposes, is working in a finite dimensional vector space, PN , spanned by the orthogonal polynomials. As a consequence of this we have that D : PN → PN−1 or, in other words, the last row of T DT −1 must be a zero-row—a well known fact for the spectral differentiation matrices. Since also f ∈ PN−1 e −1 , hence for consistency we are free to introduce a non-zero constant in the upper right corner of D −1 e eliminating the singularity of P while maintaining the exact integrating nature of the modified integration preconditioner, P −1 . The preconditioner is in general full, as can only be expected for global methods. However, the entries are given in explicit form and the application of the preconditioner is fast in cases where a fast transform, e.g., a fast cosine transform or a fast multipole method, is applicable. The crucial question about the efficiency of the preconditioner in combination with the penalty terms, however, still needs to be clarified. Let us recall that after preconditioning of the original problem, Eq. (5), we have
T −1 I1N T − τ P (−1) INN u = P (−1) f .
Hence, we need to consider the properties of the operator e (−1) T I N T −1 , I1N − τ D N
with a particular emphasis on its eigenvalue spectrum. The central observation to make is that utilizing the properties of the orthogonal polynomials we recover α(N + 1)−1 , 0,
i = 0, i ∈ [1, N − 2], e (−1) T I N T D N ij = −1 −1 N (2N + 1) , i = N − 1, −1 N (N + 1)−1 , i = N, which simplifies the analysis tremendously. An immediate consequence of this is that at least N − 2 of the eigenvalues are exactly one. The remaining 3 eigenvalues can be recovered from a 3 × 3 system, yielding one more unit eigenvalue while the remaining two eigenvalues depend on the actual choice of the free parameter,α, that was introduced to remove the singularity of the original preconditioner. Regardless of this parameter, however, the spectrum can be shown to be bounded independent of N provided only that τ ∼ N 2 —a situation that we have determined previously as a requirement for stability. The actual bounds on the spectrum naturally depends on the free parameters as well as the choice of τ . Making the usual choice of τ = N(N + 1)/4 for the Legendre penalty approximation of Eq. (4) we recover the spectrum on closed form as −1
s
λ0,1 = 1 +
3N + 2 , 4(2N + 1)
such that λ P
(−1)
N
D(1) − τ IN
λ2,N = 1,
"
∈ 1, 1 +
s
#
5 . 12
36
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
As has been shown in [17] the situation for the general operator is very similar. Although the spectrum in general is unbounded, albeit with a slow polynomial growth rate, the qth operator have at most 2q eigenvalues not being unity, i.e., the preconditioning supplies a very efficient matrix factorization, which can be exploited for efficient iterative solution. Among other things, the fixed number of non-unit eigenvalues guarantees convergence in a fixed number of conjugate gradient iterations independent of N . It should be noted that since the development of the preconditioner depends solely on the properties of the orthogonal polynomials equivalent developments, and results, can be obtained for all the classical orthogonal polynomials in combination with any of the Gauss quadratures. Moreover, the treatment of complex boundary operators does not pose any additional problems with the penalty formulation, as we have seen in Section 4.1, and the development of preconditioners for such problems follow the same route as outlined above. The preconditioning of mixed problems can done efficiently by using the highest order preconditioner as the inverse matrices have a bounded spectrum and hence will not affect the overall performance of approach. In a similar spirit one could extend the procedure to the multi-dimensional case by means of tensor-products. Despite all these favorable properties, there is still much room for improvement. In particular, a better understanding of which penalty formulation to choose such as to minimize the spectral radius of the preconditioned operator is needed. As we saw in Section 3 there is nothing unique about the choice of Q(x) and alternative choices to those discussed here may well result in a better overall performance. Regardless of this, the combination of the penalty method and integration preconditioners appears to offer an interesting alternative to the widespread use of low-order based approximation although more work is required to validate the performance for multi-dimensional problems. 4.3. Constructing spectral methods on unstructured grids As we discussed in Section 3 the central point as related to stability lies not in the details of the approximation of the solution but rather in the way the equation is satisfied. With this in mind one could argue that any distribution of grid points would suffice for the purpose of stability and that there certainly is no reason that one is restricted to consider the classical grid points related to quadrature point. For the purpose of accuracy, however, the position of the grid points is essential as is exemplified by the appearance of the Runge-phenomenon when using an equidistant grid for interpolation of an analytic function on the interval. Hence, even if a linear problem is considered in which case stability can be established for any set of grid points, the projection of the initial conditions onto the nodes poses a problem which requires that the Gauss quadrature points be introduced as discussed at length in [4]. The whole procedure fails for nonlinear problems where the distribution of the grid points again plays a central role. There is, consequently, little reason to consider the construction and use of stable methods on unstructured grids for one-dimensional problems as the need for grid points leading to a well behaved approximation enters at one point or another. For problems beyond one dimension the situation is entirely different. For multi-dimensional problems, polynomial approximations are usually constructed through a tensor product of one dimensional approximations. However, for domains which are not amenable to such a construction, prominent examples being the triangular and tetrahedral domain, the formulation of practical high-order schemes for the solution of initial boundary value problems on such domains becomes less than trivial and viable modal based schemes have only recently appeared (see references in [20,22,24]).
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
37
Fig. 4. Example of nodal sets for polynomial interpolation of order n within the triangle. The symmetry of the points is emphasized by n1 + 3n3 + 6n6 = Nn2 where n1 , n3 and n6 signifies the number of one, three, and six symmetries, respectively.
The realization that for linear problems the distribution of the grid points affects the accuracy but not the stability suggests that one can always construct stable schemes on general domains provided only that the polynomials are unique. To have an accurate scheme however one must have families of grid points within the general domain that lead to well behaved polynomial interpolation of analytic functions. In recent years such families have been computed for the triangular and tetrahedral domain by various techniques ([16] and references therein) with an example of such grid points within a triangular domain given in Fig. 4. In light of these recent results, and the experience gained by discussing the one-dimensional schemes, it seems that a combination of these unstructured, yet well behaved, grid points with the penalty method would allow for the construction of stable and accurate schemes for the solution of initial boundary value problems on triangular and tetrahedral domains, hence providing the essential building block for an unstructured high-order multi-domain scheme with a nodal basis. To investigate these prospects further, let us consider the linear wave equation ∂u + V · ∇u = 0, ∂t where D is taken as a general planar straightsided triangle, i.e., the metric is constant, and the boundary conditions are given as (x, y) ∈ D:
∀x ∈ δ D:
α(x)u(x, t) = g(x, t),
α(x) =
|V · n|, if V · n 6 0, 0, if V · n > 0.
We shall seek polynomial solutions of the form u(x, y, t) =
N X
u(xj , yj , t)Lj (x, y),
j =1
where we have N grid points, (xj , yj ), and Lj (x, y) represents the genuinely two-dimensional Lagrange polynomial based on these grid points. It is assumed that the grid points allow for unique interpolation. In the spirit of the penalty method we ask that the equation be satisfied in a Galerkin sense inside the domain but with a modification at the boundaries as ∂u ∀j ∈ [1, . . . , N]: + V · ∇u, Lj (x) = −Tjj Ajj u(x j , t) − g(x j , t) , ∂t
38
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
where (·, ·) represents the inner product over D. Introducing
Mij = Li (x, y), Lj (x, y) ,
Sij = V · ∇Li (x, y), Lj (x, y) ,
representing the mass- and stiffness matrix, respectively, yields the short form
M u˙ + Su = −T A u − g(x, t) , where T is diagonal matrix with all entries corresponding to internal grid points being zero and A is likewise a diagonal matrix with entries corresponding to the boundary conditions, α(n). Hence, T is the multi-dimensional representation of the penalty parameter, τ . The assumption of uniqueness of the polynomial basis implies that S is almost skew-symmetry which again has the consequence that uT Su =
1 2
T
T
T
V · na ua M a ua + V · nb ub M b ub + V · nc uc M c uc ,
where M a,b,c are mass matrices depending only on the grid points along the three edges with the outward pointing normal vectors, na,b,c , u is a vector of u(xj , yj ) and ua,b,c refers to vectors of those function values located along one of the three edges of the triangle. With this result stability 1 d T u Mu 6 0, 2 dt follows directly by requiring that along each edge the elements of T is chosen such that ua
T
1 Ma 2
− T a ua 6 0.
However, since M a,b,c are mass-matrices, and hence symmetric positive definite, this is clearly always possible and stability is guaranteed. We have consequently succeeded in constructing a stable spectral method on the form du + M −1 Su = −M −1 T A u − g(x, t) , dt where the restrictions on the grid appears solely by requiring accuracy rather than stability or, in other words, stable methods can be constructed on any distribution of grid points allowing unique polynomial interpolation. Moreover, one can show that [20,22] that M −1 S = Vx Dx + Vy Dy , where Dx,y is nothing else that the differentiation matrices with the entries ∂Lj (xi , yi ) ∂Lj (xi , yi ) y Dxij = , Dij = , ∂x ∂y and V = (Vx , Vy ), i.e., the implementation of the method is equivalent to one-dimensional collocation methods with the difference being a global correction due to the boundary conditions, similar to examples we discussed in Section 3. The construction of stable spectral methods on general unstructured grids has been discussed extensively in [20,22] where also schemes for advection–diffusion equations and hyperbolic systems on triangles and tetrahedra are given. Apart from ensuring the accuracy of the schemes, it turns out that the choice of grids plays a significant practical role. Indeed, the schemes discussed in [20,22] have a time-step restriction 1t ∼ n2 , similar to
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
39
the one-dimensional case. More importantly, however, it has been shown that provided only that the grid distributions possess a certain amount of symmetry, the differentiation matrices can be factorized to allow for fast computation of derivatives. This renders the unstructured schemes very competitive with tensor-product based alternatives while providing a spectrally accurate stable algorithm with no geometric restrictions.
5. Concluding remarks The conceptual difficulty of the penalty method lies in the observation that it is only necessary to satisfy the boundary conditions asymptotically in a way consistent with the order of the scheme, i.e., there is no reason to insist on the boundary conditions being exactly obeyed if one is willing to accept a certain overall error. Accepting this has a number of interesting and important consequences for the analysis as well as the applications of spectral methods for the solution of initial boundary value problems. As we have discussed through a number of examples, imposing the boundary conditions weakly through a penalty term provides a constructive approach to achieving stability of the scheme. Hence, rather than formulating the scheme and then addressing stability, the issue of stability is dealt with as an integral part of the formulation of the scheme. Moreover, the complexity of the boundary operators has little impact on the formulation and implementation of the schemes. The use of the penalty method, however, have even more profound consequences. By providing an elegant splitting between the operator and the boundary conditions, the analysis of the spectral methods is often greatly simplified as exemplified in the development of the integration preconditioners for the pseudospectral operators. Moreover, we discussed how the use of the penalty methods emphasize the separate importance of accuracy and stability by allowing for the construction of stable schemes on any grid supporting unique polynomial interpolation while the requirement of spectral accuracy is the main cause for restrictions on the distribution of grid points. As a last observation, it is interesting to note that developments, closely related to the spectral penalty methods, have appeared in other areas of high-order methods for initial boundary value problems in recent years. Examples of such schemes are SAT-compact finite-difference schemes [3], the summation-by-parts schemes in high-order finite-difference schemes [27] and multi-domain versions there of [5], and errorbounded finite-difference schemes [1]. Also the discontinuous Galerkin method [6] can be expressed as a penalty formulation with additional constraints imposed by the requirement of conservations [25]. The success of these parallel efforts suggests that imposing the boundary conditions weakly, through a penalty term, may well provide a constructive unifying framework for the development of stable highorder methods, hence addressing and perhaps resolving one of the central difficulties related to such schemes. References [1] S. Abarbanel, A. Ditkowski, Asymptotically stable fourth-order accurate schemes for the diffusion equation on complex shapes, J. Comput. Phys. 133 (1997) 279–288. [2] K. Black, Polynomial collocation using a domain decomposition solution to parabolic PDE’s via the penalty method and explicit/implicit time marching, J. Sci. Comput. 7 (1992) 313–338.
40
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
[3] M.H. Carpenter, S. Abarbanel, D. Gottlieb, Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: Methodology and application to high-order compact schemes, J. Comput. Phys. 111 (1994) 220–236. [4] M.H. Carpenter, D. Gottlieb, Spectral methods on arbitrary grids, J. Comput. Phys. 129 (1996) 74–86. [5] M.H. Carpenter, J. Nordström, D. Gottlieb, A stable and conservative interface treatment of arbitrary spatial accuracy, J. Comput. Phys. 148 (1999) 341–365. [6] B. Cockburn, S. Hou, C.W. Shu, TVB Runge–Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: The multidimensional case, Math. Comp. 54 (1990) 545–581. [7] E.A. Coutsias, T. Hagstrom, J.S. Hesthaven, D. Torres, Integration preconditioners for differential operators in spectral τ -methods, in: Proc. of International Conference on Spectral and High Order Methods, ICOSAHOM’95, Houston, USA, 1995, pp. 21–38. [8] L. Dettori, B. Yang, On the Chebyshev penalty method for parabolic and hyperbolic equations, M2 AN 30 (1996) 907–920. [9] W.S. Don, D. Gottlieb, The Chebyshev–Legendre method: Implementing Legendre methods on Chebyshev points, SIAM J. Numer. Anal. 31 (1994) 1519–1534. [10] D. Funaro, A multidomain spectral approximation of elliptic equations, Numer. Methods Partial Differential Equations 2 (1986) 187–205. [11] D. Funaro, Domain decomposition methods for pseudospectral approximations. Part I. Second order equations in one dimension, Numer. Math. 52 (1988) 329–344. [12] D. Funaro, D. Gottlieb, A new method of imposing boundary conditions in pseudospectral approximations of hyperbolic equations, Math. Comp. 51 (1988) 599–613. [13] D. Funaro, D. Gottlieb, Convergence results for pseudospectral approximations of hyperbolic systems by a penalty-type boundary treatment, Math. Comp. 57 (1991) 585–596. [14] J.S. Hesthaven, A stable spectral multi-domain method for the unsteady, compressible Navier–Stokes equations, in: Proc. of the 9th International Conference on Domain Decomposition, Bergen, Norway, 1996, pp. 121–129. [15] J.S. Hesthaven, A stable penalty method for the compressible Navier–Stokes equations. II. One-dimensional domain decomposition schemes, SIAM J. Sci. Comput. 18 (1997) 658–685. [16] J.S. Hesthaven, From electrostatics to almost optimal nodal sets for polynomial interpolation in a simplex, SIAM J. Numer. Anal. 35 (1998) 655–676. [17] J.S. Hesthaven, Integration preconditioning of pseudospectral operators. I. Basic linear operators, SIAM J. Numer. Anal. 35 (1998) 1571–1593. [18] J.S. Hesthaven, A stable penalty method for the compressible Navier–Stokes equations. III. Multi dimensional domain decomposition schemes, SIAM J. Sci. Comput. 20 (1999) 62–93. [19] J.S. Hesthaven, D. Gottlieb, A stable penalty method for the compressible Navier–Stokes equations. I. Open boundary conditions, SIAM J. Sci. Comput. 17 (1996) 579–612. [20] J.S. Hesthaven, D. Gottlieb, Stable spectral methods for conservation laws on triangles with unstructured grids, Comput. Methods Appl. Mech. Engrg. 175 (1999) 361–381. [21] J.S. Hesthaven, J. Juul Rasmussen, L. Bergé, J. Wyller, Numerical studies of localized wave fields governed by the Raman-extended derivative nonlinear Schrödinger equation, J. Phys. A: Math. Gen. 30 (1997) 8207–8224. [22] J.S. Hesthaven, C.H. Teng, Stable spectral methods on tetrahedral elements, SIAM J. Sci. Comput. (2000) to appear. [23] T.J.R. Hughes, W.T. Liu, A. Brooks, Finite element analysis of incompressible viscous flows by the penalty function formulation, J. Comput. Phys. 30 (1979) pp. 1–60. [24] G.E. Karniadakis, S.J. Sherwin, Spectral/hp Element Methods for CFD, Oxford University Press, Oxford, 1999. [25] M. Kirby, T. Warburton, Private communication, 1998.
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
41
[26] D. Kopriva, Multidomain spectral solution of compressible viscous flows, J. Comput. Phys. 115 (1994) 184– 199. [27] P. Olsson, Summation by parts, projections and stability. I, Math. Comp. 64 (1995) 1035–1065. [28] R. Temam, Navier–Stokes Equations, North-Holland, Amsterdam, 1979. [29] J. Tjujillo, G.E. Karniadakis, A penalty method for the vorticity–velocity formulation, J. Comput. Phys. 149 (1999) 32–58.
Spectral penalty methods ✩ J.S. Hesthaven 1,2 Division of Applied Mathematics, Brown University, Providence, RI 02912, USA
Abstract We present an overview of the theory and applications of spectral penalty methods for the stable solution of initial boundary value problems. Through a number of case studies we shall illustrate that imposing boundary conditions weakly through a penalty term rather than strongly as is classically done, removes many of the problems traditionally associated with the construction and analysis of stable and accurate pseudospectral approximations. The discussion includes the treatment of complex boundary operators, the construction of integration preconditioners and the formulation of spectral methods on general unstructured grids, highlighting the versatility and flexibility introduced by the elegant splitting of the operator and the boundary conditions exploited in the penalty methods. 2000 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Spectral methods; Pseudospectral methods; Stability; Boundary conditions
1. Introduction The analysis of numerical schemes for the solution of partial differential equations is often complicated significantly by the presence of boundary conditions, which often makes it difficult to analyze such basic properties as stability and convergence for even simple cases, and renders the analysis of more realistic problems and complex boundary operators essentially impossible. This is particularly true when considering high-order or spectral methods. However, the recent increasing interest in the modeling of transient problems in science and engineering puts an emphasis on the use of such high-order methods and, in particular, schemes for which numerical stability can be established in a rigorous manner. In a traditional approach, the value of a numerical scheme for the solution of an initial boundary value problems is often judged through an analysis of the complete scheme, including the boundary conditions, after it has been formulated. This approach often implies that the analysis must be done on ✩
Expanded version of invited talk presented at the International Conference on Spectral and High Order Methods 1998 (Herzliya, Israel). 1 E-mail: [email protected] 2 Supported by DARPA/AFOSR grant F49620-96-1-0426, by DOE grant DE-FG02-95ER25239 and by NSF grant ASC9504002. 0168-9274/00/$20.00 2000 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 9 9 ) 0 0 0 6 8 - 9
24
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
a problem by problem basis and it is often complicated significantly by the modifications of the discrete operators introduced by imposing the boundary conditions exactly. These considerations are particularly true when one considers high-order and spectral methods where the existing stability theory is far from being complete. A viable approach, known as penalty methods in the context of spectral methods, towards a resolution of these problems has attracted increasing interest in recent years. The central idea behind the penalty method lies in the observation that it is only necessary to satisfy the boundary conditions asymptotically in a way consistent with the order of the scheme, i.e., there is no reason to insist on the boundary conditions being exactly obeyed if one is willing to accept a certain overall error. This approach makes the implementation of complex boundary conditions straightforward and it provides an elegant splitting between the operator and the boundary conditions such that the former maintains its delicate structure while the latter is dealt with in a simple manner. Moreover, this splitting enables the analysis of very complex schemes and provides a constructive approach to the construction of stable spectral methods for the solution of initial boundary value problems. What remains in this paper is organized as follows. In Section 2 we provide a very selective overview of the central properties of Legendre polynomials and polynomial methods. An understanding of these issues is essential to Section 3 where we present an overview of the basic elements of spectral penalty methods in different variations. This sets the stage for Section 4 where we show how the use of penalty methods can be helpful, or even essential, in resolving a number of known problems related to the use and implementation of spectral methods. We consider, as examples only, the treatment of complex boundary conditions and the associated problem of how to construct stable multi-domain schemes, the formulation of efficient preconditioners for pseudospectral operators and the construction of stable spectral methods on unstructured grids, most notably general grids within triangles and tetrahedra. Section 5 contains a few general remarks and draws the attention to recent parallel developments in alternative high-order methods, suggesting that penalty methods provide a general and convenient framework in which to formulate and analyze stable high-order methods for the solution of initial boundary value problems. 2. Polynomials and polynomial methods While the formulation of the penalty methods is independent of the actual choice of the polynomial basis, providing the basis of the spatial approximation, the discussions and subsequent analysis are most easily carried out in the framework of Legendre methods. For the sake of simplicity we shall, therefore, focus our attention on such approximations, keeping in mind that this does not pose a general restriction. We primarily concern ourselves with Legendre methods utilizing the discrete approximation IN u(x) =
N X n=0
uen Pn (x),
N 1 X uen = u(xi )Pn (xi )wi , γen i=0
(
γen = 2
(2n + 1)−1 , n < N, N −1 , n = N,
where Pn (x) represents the Legendre polynomial of order n, and (xi , wi ) signifies the set of Legendre– Gauss–Lobatto nodes and weights given as dPN 2 = 0, grid points: 1 − xi2 weights: wi = . dx xi N(N + 1)PN2 (xi ) It should be noted that the choice of the Gauss–Lobatto quadrature does not represent a constraint and alternative Gauss type quadratures equally well could be used. However, for solving boundary value
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
25
problems, the Gauss–Lobatto quadratures is the most natural choice as it includes the boundaries in the nodal set. As we shall see shortly, it is central to recall that the Gauss–Lobatto quadrature is exact as
∀u(x) ∈ P2N−1 = span Pn (x)
2N−1 n=0
Z1
:
u(x) dx =
N X
u(xi )wi .
(1)
i=0
−1
An immediate consequence of the use of the Gauss quadratures is IN u(xi ) = u(xi ), i.e., the approximation is an interpolation IN u(x) =
N X
uen Pn (x) =
n=0
N X
u(xi )li (x),
i=0
where li (x) represents the N th order Lagrange polynomial based on xi . Hence, there are two mathematically equivalent but computationally different ways of expressing the interpolation. This again implies that the central operation of computing derivatives can be done in two ways: N N X X dIN u ? 0 e = u Pn (xi ) = u(xj )Dij , dx xi n=0 n j =0
(2)
utilizing either the expansion coefficients, uen , directly or exploiting the interpolation property to derive the differentiation matrices, D, in physical space as in finite difference methods. We recall that ue0n = (2n + 1)
N X
uep =
p=n+1 p+n odd
N X
e np u ep , D
p=0
dlj (x) Dij = . dx xi
?
The = in Eq. (2) is meant as a reminder that we have yet to understand the exact relation between the two formulations. If we introduce
T
u = u(x0 ), . . . , u(xN ) ,
T
e0 , . . . , u e= u eN , u
T
e00 , . . . , u e0 = u e0N , u
e and u e 0 , respectively, we may express Eq. (2) as as the vectors of the nodal values, u, and modal values, u
du = Du. dx As a direct result of the Gauss quadrature, Eq. (1), we recover that eu e0 = D e, u
e , D = T −1 DT
i.e., the two formulations are related through a simple similarity transform as Tij =
1 Pi (xj )wj , γei
Tij−1 = Pj (xi ),
are orthogonal transformation matrices.
26
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
e is obtained directly from the polynomials [7] The entries of D
d −1 e −1 d −1 D , T ⇒ T −1 = T dx dx which represents a well known recurrence relation e= T −1 D
0 0 (2n + 1)Pn (x) = Pn+1 (x) − Pn−1 (x).
e −1 , is tri-diagonal and well-conditioned with a purely imaginary spectrum Hence, the pseudoinverse, D and the magnitude of the maximum eigenvalue is uniformly bounded by 1/π .
3. A short course on penalty methods Let us now return to the central issue of solving initial boundary value problems. We restrict, solely for simplicity, the attention to the construction of a spectral method for the solution of the scalar waveequation ∂u ∂u = , u(1, t) = f (t), (3) ∂t ∂x with x ∈ [−1, 1]. To construct the scheme we have to make a number of choices, the first one being in which space to seek solutions. In the following, we shall seek solutions of the form uN (x, t) =
N X n=0
en Pn (x) = u
N X
u(xi )li (x),
i=0
where uen are the discrete Legendre–Gauss–Lobatto expansion coefficients based on the quadrature points, xi , on which also the Lagrange polynomial, li (x), is based. The second choice to make is in which way we require the equation to be satisfied using the methods of weighted residuals. With proper choices of the test functions, one can then recover the classic Galerkin, tau and collocation methods. In the latter case one should keep in mind that there is no requirement on which points one requires the equation to be satisfied. In particular, these points need not be the same as those on which the polynomial approximation is based. As we shall see, this latter observation opens up for the construction of a number of interesting schemes for the solution of initial boundary value problems. Regardless of the specifics of the method, however, the overarching problem of all the schemes is to impose the boundary conditions in a consistent and stable manner while maintaining the accuracy of the scheme. The Galerkin and tau methods suffer from well-known difficulties for nonlinear problems while the issue of stability and discreteness of the approximation at the boundary is a concern for collocation approximations. To address these issues, we consider the following approximation of Eq. (3): ∂uN ∂uN (4) = − τ Q(x) uN (1, t) − f (t) ∂t ∂x as was first proposed in [12,13] within the context of spectral methods for evolutionary problems. A very similar approach for elliptic problems was proposed prior to that in [10,11]. A related technique has been proposed previously in the context of finite element methods as a way of enforcing global constraints
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
27
[23,28], although this was later rejected due to excessive stiffness being introduced through this global penalization. In what remains we shall restrict our attention to situations with penalization of boundary constraints only, as this is where the advantages of using the penalty methods are highlighted, and we shall, moreover, consider mainly problems with a temporal component. Eq. (4) is solved everywhere in the domain, τ plays the role of a free parameter, or a Lagrange multiplier, and Q(x) is a polynomial, that, as we shall see, characterizes the scheme. Imposing the boundary conditions weakly through a penalty term, as in Eq. (4), introduces additional degrees of freedom in the construction of the scheme and facilitates the formulation of entirely new families of schemes. We emphasize that for τ approaching infinity, the boundary condition must become strongly imposed, hence establishing consistency by equivalence with the traditional formulation. For τ being finite, the equation as well as the boundary condition is enforced at the boundary. In completing the specification of the scheme, we need to choose Q(x) and we shall moreover need to specify a constraint that allows for the specification of τ . As we shall see, τ is most often specified under the constraint of asymptotic stability in some norm but one should keep in mind that this choice is not unique, e.g., one could also require conservation as a constraint. In considering a few examples, let us begin with the classical choice of Q(x) =
(1 + x)PN0 (x) = 2PN0 (1)
1, x = 1, 0, x 6= xN ,
and require that Eq. (4) be satisfied in a collocation sense at the Legendre–Gauss–Lobatto points, i.e., the same points being used in the construction of the approximation. We note, in particular, that this choice results in solving Eq. (3) in the interior in a way similar to a traditional scheme while only the treatment of the boundary point is somewhat different from that of a classical approach. To complete the specification of the scheme, we require that the approximation be asymptotically stable in the discrete L2 -norm as
N X 1 d d uN (xi ) kuN k2N = uN (xi ) wi − τ u2N (xN )wN 6 0. 2 dt dx i=0
Utilizing quadrature and integration by parts we recover 1 d 1 1 kuN k2N = u2N (xN ) − τ u2N (xN )wN − u2N (x0 ) 6 0, 2 dt 2 2 and stability provided τ>
1 N(N + 1) > . 2wN 4
Hence, a lower finite bound on τ is established which guarantees that the semi-discrete approximation is asymptotically stable by construction. Moreover, if u(x) can be assumed smooth, the magnitude of the penalty term will decay exponentially and the scheme maintains the overall spectral accuracy. In fact, early results in [12] suggest that the penalty methods is superior to the standard method of imposing boundary conditions in terms of accuracy. While the enhanced accuracy and inherent stability appear as the main motivation of the original work [12,13] it is by considering the versatility of the weakly imposed boundary operators that one begins to appreciate the general approach.
28
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
Rather than choosing the Q(x) associated with the collocation scheme, one could consider using 0 PN+1 (x) + PN0 (x) 2 in Eq. (4) and ask that the equation be satisfied in a Galerkin sense. We shall refrain from writing the scheme directly, as it involves mass- and stiffness matrices, but rather determine τ by requiring that the scheme be stable in the L2 inner product with stability resulting for τ > 1/2. Thus, we seek polynomial solutions that satisfy the equation
Q(x) =
∂uN ∂uN P 0 (x) + PN0 (x) uN (1, t) − f (t) . = − N+1 ∂t ∂x 4 Note, however, that Q(x) is orthogonal to all polynomials of order N or less that satisfy the boundary conditions, i.e., the method is nothing else than the Legendre–Galerkin method for the wave-equation as could also have been realized by considering the Legendre–Galerkin error-equation which is identical to the penalty scheme given above, i.e., they yield the exact same polynomial solution. In the same spirit one can consider using Q(x) = PN (x), which, under the restriction of L2 stability, yields a lower bound τ > (2N + 1)/4 resulting in a stable scheme as ∂uN ∂uN 2N + 1 = − PN (x) uN (1, t) − f (t) . ∂t ∂x 4 By comparison with the error equation it becomes clear that this scheme is nothing else that a different formulation of the well-known Legendre tau method. As a final example, let us again consider the formulation of the collocation scheme, although we shall require that the approximation is based on the Chebyshev–Gauss–Lobatto nodes, yj , while the equation is sought satisfied at the usual Legendre–Gauss–Lobatto nodes, xj . As for the Legendre collocation method discussed above, we shall use
Q(x) =
(1 + x)PN0 (x) , 2PN0 (1)
and seek a solution to the equation ∂uN ∂uN (1 + x)PN0 (x) uN (1, t) − f (t) . = −τ 0 ∂t ∂x 2PN (1)
While the approximation itself is constructed using the Chebyshev nodes we shall require that the equation be satisfied at the Legendre nodes. Note that as a consequence of this, the penalty term has a global effect as it is non-zero at all the Chebyshev nodes. Nevertheless, since both sides of the equation represents a polynomial of order N specified at N + 1 nodes, it is clearly unique and we can hence obtain stability in the discrete Legendre inner product, similar to the classical case, resulting in the bound N(N + 1) . 4 Hence, stability of the Chebyshev method in the usual L2 -norm has been established using the Legendre nodes. From the perspective of computation, however, the Legendre nodes are never used. On the other hand, by comparison with the error-equation, the scheme is identical to a Legendre collocation method while the approximation is based on the nodes of Chebyshev polynomials. τ>
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
29
This latter method, known as the Chebyshev–Legendre method [9], serves as a good example of the wealth of new schemes that one arrives at by imposing the boundary condition weakly rather than strongly in combination with an increased understanding of exactly how one needs to require that the equations be satisfied. In recent years these central ideas have been further developed and have resulted in a wealth of new schemes by introducing and combining different choices of Q(x) and the way the equation is satisfied [4,8,19]. In all previous work, τ has been chosen to ensure stability in some energy norm. As a slight generalization, let us briefly consider the formulation of a penalty method for a hyperbolic system. For this purpose, we consider the problem
∂ u 0 1 ∂ u + = 0, ∂t v 1 0 ∂x v
u(−1) = f (t),
v(1) = g(t).
As the problem is hyperbolic we have
A = S −1 ΛS,
where
1
1
1 , S=√ 2 1 −1
Λ=
1
0
0 −1
,
as the characteristic decomposition, clarifying that one boundary condition is required at each boundary. The penalty scheme for this problem becomes (q N = [uN , vN ]T ) ∂q N ∂q + A N = −τ − Q− (x)SR − S T q N − f − − τ + Q+ (x)SR + S T q N − f + ∂t ∂x
where Q± (x) and τ ± take the usual meaning at the two boundaries and
R− =
1 0
,
0 0
R+ =
0 0 0 1
,
f− =
f (t) vN (−1)
,
f+ =
uN (1)
.
g(t)
Thus, the boundary conditions are imposed through the characteristic variables such that stability can be established for each decoupled equation of the individual characteristic variable. In other words, stability for the system follows directly from that of the scalar problem with straightforward extensions to non-linear problems through linearization and localization, requiring only that the solution is locally smooth, and to multi-dimensional problems for symmetrizable systems and conservation laws possessing an entropy-flux pair [13,15,18,19].
4. Beyond the basics While the basic ideas underlying the penalty method can be illustrated well through an analysis of the linear wave-equation, one has to look towards more general applications to fully appreciate the versatility and flexibility offered by this approach. In the following we shall present a number of such examples, illustrating how the penalty methods facilitate the construction of entirely new families of methods while allowing for elegant and simple analysis and solutions to old and well-known problems associated with the use of spectral methods.
30
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
4.1. Imposing complex boundary operators One of the central issues, and by many perceived as the central obstruction, in the use of spectral methods lies in the formulation of simple, yet stable, ways of imposing boundary conditions of a complicated form. To illustrate the advantages of using a penalty method for imposing complex boundary operators, let us consider the linear advection–diffusion equation [19], ∂u ∂u ∂2 u +λ = ε 2, ∂t ∂x ∂x αu(−1, t) − εβux (−1, t) = g1 (t), γ u(1, t) + εδux (1, t) = g2 (t), where (α, β, γ , δ) > 0 must be chosen to ensure well-posedness [19]. The question to address is how to impose boundary conditions of such general nature in a consistent and stable manner. While imposing the boundary conditions strongly requires modification of the operators, a process which becomes increasingly complex for more general problems and boundary operators, enforcing the conditions in a weak way only through a penalty term immediately yields ∂u ∂u ∂2 u +λ = ε 2 − τ − Q− (x) αu(−1, t) − εβux (−1, t) − g1 (t) ∂t ∂x ∂x − τ + Q+ (x) γ u(1, t) + εδux (1, t) − g2 (t) .
Here, and in what remains, we shall not distinguish between u and uN for simplicity as the correct meaning should be clear from the context. If we require that the equation be satisfied at the Legendre collocation points and take Q± (x) to be an N th order polynomial vanishing at all the interior points except at x = ±1, respectively, we may guarantee discrete L2 -stability as N−1 X 1 d u2x (xi )wi , kuk2N 6 −ε 2 dt i=1
provided only that the penalty parameters are bounded like − + τα,β 6 τ − 6 τα,β ,
where ± τa,b =
τγ−,δ 6 τ + 6 τγ+,δ ,
q i 1 h ε + 2κ ± 2 κ 2 − εκ − 12 εwN |λ| , wN εb
κ=
wN a . b
We note that the result on stability is independent on (α, β, γ , δ) and ε, i.e., it remains uniformly valid for vanishing viscosity. In [19] similar formulations for hyperbolic, parabolic and mixed problems are discussed, including stability bounds for all cases. To appreciate the accuracy and validity of the above approach, let us apply it to solve Burgers equation ∂u ∂u ∂2 u +u = ε 2, ∂t ∂x ∂x
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
31
Fig. 1. Traveling wave solution of Burgers equation computed using a Legendre collocation method employing a penalty formulation at the boundaries. Table 1 Error of the traveling wave solution using a Legendre collocation method with a penalty method at the boundary points N
L2
L∞
16
3.41E−03
3.26E−02
32
2.43E−05
3.50E−04
64
1.09E−09
2.21E−08
128
4.98E−12
7.62E−11
subject to the initial (exact) conditions
x u(x, 0) = − tanh + 1. 2ε In this case, the boundary operator of the non-linear form u2 (−1, t) − εux (−1, t) = g1 (t),
εux (1, t) = g2 (t),
ensures well-posedness and the solution represents a kink-solution traveling at a constant speed as illustrated in Fig. 1. Computing the accuracy of the solution yields the results shown in Table 1, confirming the global spectral accuracy of the scheme even though the boundary conditions are imposed only weakly. The emphasis so far has been on the construction of stable ways in which to impose open boundary conditions. However, as has been utilized in [2,14,15,18,21,26], the generalization to multi-domain solution of the problem is straightforward. Indeed, if we consider each domain as an isolated entity
32
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
Fig. 2. Viscous shock solution of Burgers equation obtained using a multi-domain Legendre collocation method with a penalty method for patching of the domains. The dots on the initial condition represents the collocation points in the 4 domain solution. In (b) we show the derivative of the solution shown in (a).
requiring boundary conditions at every time-step, the scheme developed above generalizes directly by simply passing information between the domains as prescribed through the dynamics of the equation, e.g., advective terms are always treated by upwinding. Hence, while the boundary conditions may not be available a priori in the individual domains they are available whenever needed by passing the required information from the neighboring domain. It should be noted that this line of argument, leading to a scheme local in time and space, assumes that information propagates with a finite speed as in wave-dominated problems, although accuracy and stability has been established and confirmed also for parabolic problems [2,15,18]. To illustrate the efficacy of this approach let us again consider the solution of Burgers equation subject to the boundary and initial conditions u(±1, t) = 0,
u(x, 0) = − sin(π x),
and with a small viscosity ε = 0.01/π . Rather than attempting to directly resolve the solution, which develops a very steep gradient around x = 0, we shall employ 4 subdomains of different size to accurately resolve the viscous shock. In Fig. 2 we show the computed solution, illustrating that even though the patching conditions, which are nothing else than the open boundary operator of Burgers equation ensuring well-posedness in L2 , are imposed only weakly the solution and its derivative remains smooth even for this very steep gradient (see [15] for further details). The results given in the above suggests a systematic way of developing stable schemes for the solution of initial boundary value problems. Given a particular problem, one first derives a set of open boundary conditions which ensures well-posedness of the total problem. In general, this set of boundary conditions is very complex. However, with the flexibility provided by the penalty method, these boundary conditions can be imposed in a stable and accurate way. Moreover, the generalization from a single domain solution to a multi-domain solution is straightforward and results in a global spectrally accurate, asymptotically stable and geometrically flexible formulation for the solution of the initial boundary value problems.
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
33
Fig. 3. Two-dimensional computation of recessed cavity flameholder flow at M = 0.9 and Re = 1500 based on the half channel width and the maximum mean flow. In (a) we show a segment of the spectral multi-domain grid, using order 12 in each domain and 156 domains for this problem, and in (b) we show the density and instantaneous streamlines, illustrating the recirculating flow within the cavity.
As an illustration of the versatility and robustness of the multi-domain framework outlined above, we show in Fig. 3 a simulation of the flow of compressible air inside a narrow grooved channel. The computation is performed using a multi-domain pseudospectral scheme for the solution of the unsteady compressible Navier–Stokes equations where the patching and treatment of open boundaries is done in the penalty method setting. For the compressible Navier–Stokes equations the boundary operator is extremely complex, containing components of the inviscid part from the Euler equations and elements of the normal and tangential stress and heat flux while being singular due to mass conservation, and it is difficult to construct a stable scheme in which this operator is imposed strongly. However, imposing the boundary operator weakly is straightforward and stability of the complete three-dimensional scheme has been established within a general curvilinear hexahedral, including the edges and vertices, with the scheme maintaining stability and validity for vanishing viscosity [14,15,18]. These results provide the fundamental elements for the construction of very general three-dimensional pseudospectral solvers for the solution of problems of gas-dynamics and supplies a convincing example of the efficacy of using the penalty method for imposing complex boundary operators. In recent years, penalty methods have been used to address a number of problems associated with complex boundary conditions and the construction of multi-domain schemes. In [2] a penalty formulation was developed for the heat-equation and implemented using a mixed explicit/implicit time integration scheme. A related approach was taken in [26] where a multi-domain scheme for the compressible Navier– Stokes equations is presented, utilizing a splitting of the inviscid and viscous part of equation and treating the former part by upwinding through characteristics while the dissipative part is included using a penalty formulation borrowed from earlier work on elliptic problems [10,11]. A stable scheme for the multidomain solution of the three-dimensional advection–diffusion equations is derived in [18] and in [21] a multi-domain scheme for the solution of certain problems in non-linear optics is proposed. An elegant application of the penalty method was recently presented in [29] where the basis idea was exploited to impose non-slip boundary conditions in a velocity–vorticity (v, ω) formulation of the
34
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
incompressible Navier–Stokes equations. The well-known lack of direct boundary conditions on the vorticity was addressed as ∂ω + ∇ × (ω × v) = ν∇ 2 ω − τ Q(x) ω − ∇ × v . ∂t Since the boundary conditions on the velocity field, v, are well known, the scheme imposes the unknown vorticity condition by consistency at the boundary of the computational domain. Its accuracy and stability has been confirmed through numerous experiments [29]. 4.2. Designing efficient preconditioners While the application of penalty methods for the stable and consistent treatment of complex boundary operators has received significant attention in recent years, there are a number of alternative uses of the penalty methods that seem very promising in resolving some classical problems associated with the use of spectral methods. As an example of this, let us consider the problem of preconditioning of dq u = f ⇒ Dq u = f , dx q subject to appropriate boundary conditions. For simplicity we assume that a Legendre collocation method is used and Dq signifies the qth order differentiation matrix while u and f represents the vectors of the nodal values of the solution and the forcing function, respectively. It is well known that the conditioning of Dq scales at best as O(N 2q ) where N is the order of the approximation. This result is not special to the Legendre method, or the type of collocation used, but is shared among all the classical polynomials and quadratures. Hence, if we hope to solve the above equation for a general force, f , we must consider the issue of preconditioning. Previous attention to this problem has resulted in the development of preconditioners based on a low-order approximation of the operator with examples being finite-difference and finite-element based approaches for advection and diffusion operators (see [17] for references). The efficiency of these methods varies significantly with the type of operator being considered and, in particular, the dimension of the problem. To understand the prospects of using the penalty method within this particular setting, let us consider preconditioning of the advective operator as Du = f ,
u(1) = 0.
(5)
We recall from Section 2 that e D = T −1 DT
such that
e −1 T , D−1 = T −1 D
e −1 is tri-diagonal, well conditioned and given in closed form. Hence, in the absence of the where D boundary conditions, the exact inverse, in some sense being the optimal preconditioner, is known. However, once we impose the boundary conditions strongly by altering the last row and column of D, the delicate structure of D is ruined and the construction of the inverse is no longer clear. Let us, therefore, consider the construction of a preconditioner for a collocation penalty formulation as [17]
D − τ INN u = f ,
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
35
where Iab is the diagonal matrix with unity in the rows from a to b, e.g., I0N represents the identity matrix. This suggests that since D is now unchanged the ideal preconditioner is given as e −1 T Pe−1 = T −1 D
since
Pe −1 D = T −1 I1N T .
e −1 having However, Pe−1 is singular since integration leaves a constant undetermined as manifested in D the first row being zero. To realize a solution to this problem we must recall that we, for all practical purposes, is working in a finite dimensional vector space, PN , spanned by the orthogonal polynomials. As a consequence of this we have that D : PN → PN−1 or, in other words, the last row of T DT −1 must be a zero-row—a well known fact for the spectral differentiation matrices. Since also f ∈ PN−1 e −1 , hence for consistency we are free to introduce a non-zero constant in the upper right corner of D −1 e eliminating the singularity of P while maintaining the exact integrating nature of the modified integration preconditioner, P −1 . The preconditioner is in general full, as can only be expected for global methods. However, the entries are given in explicit form and the application of the preconditioner is fast in cases where a fast transform, e.g., a fast cosine transform or a fast multipole method, is applicable. The crucial question about the efficiency of the preconditioner in combination with the penalty terms, however, still needs to be clarified. Let us recall that after preconditioning of the original problem, Eq. (5), we have
T −1 I1N T − τ P (−1) INN u = P (−1) f .
Hence, we need to consider the properties of the operator e (−1) T I N T −1 , I1N − τ D N
with a particular emphasis on its eigenvalue spectrum. The central observation to make is that utilizing the properties of the orthogonal polynomials we recover α(N + 1)−1 , 0,
i = 0, i ∈ [1, N − 2], e (−1) T I N T D N ij = −1 −1 N (2N + 1) , i = N − 1, −1 N (N + 1)−1 , i = N, which simplifies the analysis tremendously. An immediate consequence of this is that at least N − 2 of the eigenvalues are exactly one. The remaining 3 eigenvalues can be recovered from a 3 × 3 system, yielding one more unit eigenvalue while the remaining two eigenvalues depend on the actual choice of the free parameter,α, that was introduced to remove the singularity of the original preconditioner. Regardless of this parameter, however, the spectrum can be shown to be bounded independent of N provided only that τ ∼ N 2 —a situation that we have determined previously as a requirement for stability. The actual bounds on the spectrum naturally depends on the free parameters as well as the choice of τ . Making the usual choice of τ = N(N + 1)/4 for the Legendre penalty approximation of Eq. (4) we recover the spectrum on closed form as −1
s
λ0,1 = 1 +
3N + 2 , 4(2N + 1)
such that λ P
(−1)
N
D(1) − τ IN
λ2,N = 1,
"
∈ 1, 1 +
s
#
5 . 12
36
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
As has been shown in [17] the situation for the general operator is very similar. Although the spectrum in general is unbounded, albeit with a slow polynomial growth rate, the qth operator have at most 2q eigenvalues not being unity, i.e., the preconditioning supplies a very efficient matrix factorization, which can be exploited for efficient iterative solution. Among other things, the fixed number of non-unit eigenvalues guarantees convergence in a fixed number of conjugate gradient iterations independent of N . It should be noted that since the development of the preconditioner depends solely on the properties of the orthogonal polynomials equivalent developments, and results, can be obtained for all the classical orthogonal polynomials in combination with any of the Gauss quadratures. Moreover, the treatment of complex boundary operators does not pose any additional problems with the penalty formulation, as we have seen in Section 4.1, and the development of preconditioners for such problems follow the same route as outlined above. The preconditioning of mixed problems can done efficiently by using the highest order preconditioner as the inverse matrices have a bounded spectrum and hence will not affect the overall performance of approach. In a similar spirit one could extend the procedure to the multi-dimensional case by means of tensor-products. Despite all these favorable properties, there is still much room for improvement. In particular, a better understanding of which penalty formulation to choose such as to minimize the spectral radius of the preconditioned operator is needed. As we saw in Section 3 there is nothing unique about the choice of Q(x) and alternative choices to those discussed here may well result in a better overall performance. Regardless of this, the combination of the penalty method and integration preconditioners appears to offer an interesting alternative to the widespread use of low-order based approximation although more work is required to validate the performance for multi-dimensional problems. 4.3. Constructing spectral methods on unstructured grids As we discussed in Section 3 the central point as related to stability lies not in the details of the approximation of the solution but rather in the way the equation is satisfied. With this in mind one could argue that any distribution of grid points would suffice for the purpose of stability and that there certainly is no reason that one is restricted to consider the classical grid points related to quadrature point. For the purpose of accuracy, however, the position of the grid points is essential as is exemplified by the appearance of the Runge-phenomenon when using an equidistant grid for interpolation of an analytic function on the interval. Hence, even if a linear problem is considered in which case stability can be established for any set of grid points, the projection of the initial conditions onto the nodes poses a problem which requires that the Gauss quadrature points be introduced as discussed at length in [4]. The whole procedure fails for nonlinear problems where the distribution of the grid points again plays a central role. There is, consequently, little reason to consider the construction and use of stable methods on unstructured grids for one-dimensional problems as the need for grid points leading to a well behaved approximation enters at one point or another. For problems beyond one dimension the situation is entirely different. For multi-dimensional problems, polynomial approximations are usually constructed through a tensor product of one dimensional approximations. However, for domains which are not amenable to such a construction, prominent examples being the triangular and tetrahedral domain, the formulation of practical high-order schemes for the solution of initial boundary value problems on such domains becomes less than trivial and viable modal based schemes have only recently appeared (see references in [20,22,24]).
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
37
Fig. 4. Example of nodal sets for polynomial interpolation of order n within the triangle. The symmetry of the points is emphasized by n1 + 3n3 + 6n6 = Nn2 where n1 , n3 and n6 signifies the number of one, three, and six symmetries, respectively.
The realization that for linear problems the distribution of the grid points affects the accuracy but not the stability suggests that one can always construct stable schemes on general domains provided only that the polynomials are unique. To have an accurate scheme however one must have families of grid points within the general domain that lead to well behaved polynomial interpolation of analytic functions. In recent years such families have been computed for the triangular and tetrahedral domain by various techniques ([16] and references therein) with an example of such grid points within a triangular domain given in Fig. 4. In light of these recent results, and the experience gained by discussing the one-dimensional schemes, it seems that a combination of these unstructured, yet well behaved, grid points with the penalty method would allow for the construction of stable and accurate schemes for the solution of initial boundary value problems on triangular and tetrahedral domains, hence providing the essential building block for an unstructured high-order multi-domain scheme with a nodal basis. To investigate these prospects further, let us consider the linear wave equation ∂u + V · ∇u = 0, ∂t where D is taken as a general planar straightsided triangle, i.e., the metric is constant, and the boundary conditions are given as (x, y) ∈ D:
∀x ∈ δ D:
α(x)u(x, t) = g(x, t),
α(x) =
|V · n|, if V · n 6 0, 0, if V · n > 0.
We shall seek polynomial solutions of the form u(x, y, t) =
N X
u(xj , yj , t)Lj (x, y),
j =1
where we have N grid points, (xj , yj ), and Lj (x, y) represents the genuinely two-dimensional Lagrange polynomial based on these grid points. It is assumed that the grid points allow for unique interpolation. In the spirit of the penalty method we ask that the equation be satisfied in a Galerkin sense inside the domain but with a modification at the boundaries as ∂u ∀j ∈ [1, . . . , N]: + V · ∇u, Lj (x) = −Tjj Ajj u(x j , t) − g(x j , t) , ∂t
38
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
where (·, ·) represents the inner product over D. Introducing
Mij = Li (x, y), Lj (x, y) ,
Sij = V · ∇Li (x, y), Lj (x, y) ,
representing the mass- and stiffness matrix, respectively, yields the short form
M u˙ + Su = −T A u − g(x, t) , where T is diagonal matrix with all entries corresponding to internal grid points being zero and A is likewise a diagonal matrix with entries corresponding to the boundary conditions, α(n). Hence, T is the multi-dimensional representation of the penalty parameter, τ . The assumption of uniqueness of the polynomial basis implies that S is almost skew-symmetry which again has the consequence that uT Su =
1 2
T
T
T
V · na ua M a ua + V · nb ub M b ub + V · nc uc M c uc ,
where M a,b,c are mass matrices depending only on the grid points along the three edges with the outward pointing normal vectors, na,b,c , u is a vector of u(xj , yj ) and ua,b,c refers to vectors of those function values located along one of the three edges of the triangle. With this result stability 1 d T u Mu 6 0, 2 dt follows directly by requiring that along each edge the elements of T is chosen such that ua
T
1 Ma 2
− T a ua 6 0.
However, since M a,b,c are mass-matrices, and hence symmetric positive definite, this is clearly always possible and stability is guaranteed. We have consequently succeeded in constructing a stable spectral method on the form du + M −1 Su = −M −1 T A u − g(x, t) , dt where the restrictions on the grid appears solely by requiring accuracy rather than stability or, in other words, stable methods can be constructed on any distribution of grid points allowing unique polynomial interpolation. Moreover, one can show that [20,22] that M −1 S = Vx Dx + Vy Dy , where Dx,y is nothing else that the differentiation matrices with the entries ∂Lj (xi , yi ) ∂Lj (xi , yi ) y Dxij = , Dij = , ∂x ∂y and V = (Vx , Vy ), i.e., the implementation of the method is equivalent to one-dimensional collocation methods with the difference being a global correction due to the boundary conditions, similar to examples we discussed in Section 3. The construction of stable spectral methods on general unstructured grids has been discussed extensively in [20,22] where also schemes for advection–diffusion equations and hyperbolic systems on triangles and tetrahedra are given. Apart from ensuring the accuracy of the schemes, it turns out that the choice of grids plays a significant practical role. Indeed, the schemes discussed in [20,22] have a time-step restriction 1t ∼ n2 , similar to
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
39
the one-dimensional case. More importantly, however, it has been shown that provided only that the grid distributions possess a certain amount of symmetry, the differentiation matrices can be factorized to allow for fast computation of derivatives. This renders the unstructured schemes very competitive with tensor-product based alternatives while providing a spectrally accurate stable algorithm with no geometric restrictions.
5. Concluding remarks The conceptual difficulty of the penalty method lies in the observation that it is only necessary to satisfy the boundary conditions asymptotically in a way consistent with the order of the scheme, i.e., there is no reason to insist on the boundary conditions being exactly obeyed if one is willing to accept a certain overall error. Accepting this has a number of interesting and important consequences for the analysis as well as the applications of spectral methods for the solution of initial boundary value problems. As we have discussed through a number of examples, imposing the boundary conditions weakly through a penalty term provides a constructive approach to achieving stability of the scheme. Hence, rather than formulating the scheme and then addressing stability, the issue of stability is dealt with as an integral part of the formulation of the scheme. Moreover, the complexity of the boundary operators has little impact on the formulation and implementation of the schemes. The use of the penalty method, however, have even more profound consequences. By providing an elegant splitting between the operator and the boundary conditions, the analysis of the spectral methods is often greatly simplified as exemplified in the development of the integration preconditioners for the pseudospectral operators. Moreover, we discussed how the use of the penalty methods emphasize the separate importance of accuracy and stability by allowing for the construction of stable schemes on any grid supporting unique polynomial interpolation while the requirement of spectral accuracy is the main cause for restrictions on the distribution of grid points. As a last observation, it is interesting to note that developments, closely related to the spectral penalty methods, have appeared in other areas of high-order methods for initial boundary value problems in recent years. Examples of such schemes are SAT-compact finite-difference schemes [3], the summation-by-parts schemes in high-order finite-difference schemes [27] and multi-domain versions there of [5], and errorbounded finite-difference schemes [1]. Also the discontinuous Galerkin method [6] can be expressed as a penalty formulation with additional constraints imposed by the requirement of conservations [25]. The success of these parallel efforts suggests that imposing the boundary conditions weakly, through a penalty term, may well provide a constructive unifying framework for the development of stable highorder methods, hence addressing and perhaps resolving one of the central difficulties related to such schemes. References [1] S. Abarbanel, A. Ditkowski, Asymptotically stable fourth-order accurate schemes for the diffusion equation on complex shapes, J. Comput. Phys. 133 (1997) 279–288. [2] K. Black, Polynomial collocation using a domain decomposition solution to parabolic PDE’s via the penalty method and explicit/implicit time marching, J. Sci. Comput. 7 (1992) 313–338.
40
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
[3] M.H. Carpenter, S. Abarbanel, D. Gottlieb, Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: Methodology and application to high-order compact schemes, J. Comput. Phys. 111 (1994) 220–236. [4] M.H. Carpenter, D. Gottlieb, Spectral methods on arbitrary grids, J. Comput. Phys. 129 (1996) 74–86. [5] M.H. Carpenter, J. Nordström, D. Gottlieb, A stable and conservative interface treatment of arbitrary spatial accuracy, J. Comput. Phys. 148 (1999) 341–365. [6] B. Cockburn, S. Hou, C.W. Shu, TVB Runge–Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: The multidimensional case, Math. Comp. 54 (1990) 545–581. [7] E.A. Coutsias, T. Hagstrom, J.S. Hesthaven, D. Torres, Integration preconditioners for differential operators in spectral τ -methods, in: Proc. of International Conference on Spectral and High Order Methods, ICOSAHOM’95, Houston, USA, 1995, pp. 21–38. [8] L. Dettori, B. Yang, On the Chebyshev penalty method for parabolic and hyperbolic equations, M2 AN 30 (1996) 907–920. [9] W.S. Don, D. Gottlieb, The Chebyshev–Legendre method: Implementing Legendre methods on Chebyshev points, SIAM J. Numer. Anal. 31 (1994) 1519–1534. [10] D. Funaro, A multidomain spectral approximation of elliptic equations, Numer. Methods Partial Differential Equations 2 (1986) 187–205. [11] D. Funaro, Domain decomposition methods for pseudospectral approximations. Part I. Second order equations in one dimension, Numer. Math. 52 (1988) 329–344. [12] D. Funaro, D. Gottlieb, A new method of imposing boundary conditions in pseudospectral approximations of hyperbolic equations, Math. Comp. 51 (1988) 599–613. [13] D. Funaro, D. Gottlieb, Convergence results for pseudospectral approximations of hyperbolic systems by a penalty-type boundary treatment, Math. Comp. 57 (1991) 585–596. [14] J.S. Hesthaven, A stable spectral multi-domain method for the unsteady, compressible Navier–Stokes equations, in: Proc. of the 9th International Conference on Domain Decomposition, Bergen, Norway, 1996, pp. 121–129. [15] J.S. Hesthaven, A stable penalty method for the compressible Navier–Stokes equations. II. One-dimensional domain decomposition schemes, SIAM J. Sci. Comput. 18 (1997) 658–685. [16] J.S. Hesthaven, From electrostatics to almost optimal nodal sets for polynomial interpolation in a simplex, SIAM J. Numer. Anal. 35 (1998) 655–676. [17] J.S. Hesthaven, Integration preconditioning of pseudospectral operators. I. Basic linear operators, SIAM J. Numer. Anal. 35 (1998) 1571–1593. [18] J.S. Hesthaven, A stable penalty method for the compressible Navier–Stokes equations. III. Multi dimensional domain decomposition schemes, SIAM J. Sci. Comput. 20 (1999) 62–93. [19] J.S. Hesthaven, D. Gottlieb, A stable penalty method for the compressible Navier–Stokes equations. I. Open boundary conditions, SIAM J. Sci. Comput. 17 (1996) 579–612. [20] J.S. Hesthaven, D. Gottlieb, Stable spectral methods for conservation laws on triangles with unstructured grids, Comput. Methods Appl. Mech. Engrg. 175 (1999) 361–381. [21] J.S. Hesthaven, J. Juul Rasmussen, L. Bergé, J. Wyller, Numerical studies of localized wave fields governed by the Raman-extended derivative nonlinear Schrödinger equation, J. Phys. A: Math. Gen. 30 (1997) 8207–8224. [22] J.S. Hesthaven, C.H. Teng, Stable spectral methods on tetrahedral elements, SIAM J. Sci. Comput. (2000) to appear. [23] T.J.R. Hughes, W.T. Liu, A. Brooks, Finite element analysis of incompressible viscous flows by the penalty function formulation, J. Comput. Phys. 30 (1979) pp. 1–60. [24] G.E. Karniadakis, S.J. Sherwin, Spectral/hp Element Methods for CFD, Oxford University Press, Oxford, 1999. [25] M. Kirby, T. Warburton, Private communication, 1998.
J.S. Hesthaven / Applied Numerical Mathematics 33 (2000) 23–41
41
[26] D. Kopriva, Multidomain spectral solution of compressible viscous flows, J. Comput. Phys. 115 (1994) 184– 199. [27] P. Olsson, Summation by parts, projections and stability. I, Math. Comp. 64 (1995) 1035–1065. [28] R. Temam, Navier–Stokes Equations, North-Holland, Amsterdam, 1979. [29] J. Tjujillo, G.E. Karniadakis, A penalty method for the vorticity–velocity formulation, J. Comput. Phys. 149 (1999) 32–58.