K Veroy

C Prud'homme

D V Rovas

A T Patera

A Posteriori Error Bounds for Reduced-Basis Approximation of Parametrized Noncoercive and Nonlinear Elliptic Partial Differential Equations

We present a technique for the rapid and reliable prediction of linear-functional outputs of elliptic partial differential equations with affine parameter dependence. The essential components are (i ) rapidly convergent global reduced-basis approximations -(Galerkin) projection onto a space W N spanned by solutions of the governing partial differential equation at N selected points in parameter space; (ii ) a posteriori error estimation -relaxations of the error-residual equation that provide inexpensive yet sharp bounds for the error in the outputs of interest; and (iii ) off-line/on-line computational procedures -methods which decouple the generation and projection stages of the approximation process. The operation count for the on-line stage -in which, given a new parameter value, we calculate the output of interest and associated error bound -depends only on N (typically very small) and the parametric complexity of the problem.

In this paper we develop new a posteriori error estimation procedures for noncoercive linear, and certain nonlinear, problems that yield rigorous and sharp error statements for all N . We consider three particular examples: the Helmholtz (reduced-wave) equation; a cubically nonlinear Poisson equation; and Burgers equation -a model for incompressible Navier-Stokes. The Helmholtz (and Burgers) example introduce our new lower bound constructions for the requisite inf-sup (singular value) stability factor; the cubic nonlinearity exercises symmetry factorization procedures necessary for treatment of high-order Galerkin summations in the (say) residual dual-norm calculation; and the Burgers equation illustrates our accommodation of potentially multiple solution branches in our a posteriori error statement. Numerical results are presented that demonstrate the rigor, sharpness, and efficiency of our proposed error bounds, and the application of these bounds to adaptive (optimal) approximation.

Introduction

The optimization, control, and characterization of an engineering component or system requires the prediction of certain "quantities of interest," or performance metrics, which we shall denote outputsfor example deflections, heat transfer rates, or drags. These outputs are typically expressed as functionals of field variables associated with a parametrized partial differential equation which describes the physical behavior of the component or system. The parame-ters, which we shall denote inputs, serve to identify a particular "configuration" of the component. We thus arrive at an implicit input-output relationship, evaluation of which demands solution of the underlying partial differential equation.

Our goal is the development of computational methods that permit rapid and reliable evaluation of this partial-differential-equation-induced input-output relationship in the limit of many queries -that is, in the design, optimization, control, and characterization contexts. Our particular approach is based on the reduced-basis method, first introduced in the late 1970s for nonlinear structural analysis, [START_REF] Almroth | Automatic choice of global shape functions in structural analysis[END_REF][START_REF] Noor | Reduced basis technique for nonlinear analysis of structures[END_REF] and subsequently developed more broadly in the 1980s and 1990s. 2-4, 12, 13, 17 The reduced-basis method recognizes that the field variable is not, in fact, some arbitrary member of the infinite-dimensional solution space associated with the partial differential equation; rather, it resides, or "evolves," on a much lowerdimensional manifold induced by the parametric dependence.

The reduced-basis approach as earlier articulated is local in parameter space in both practice and theory. [START_REF] Fink | On the error behavior of the reduced basis technique for nonlinear finite element approximations[END_REF] As a result, the computational improvements -relative to conventional (say) finite element approximation -are often quite modest. [START_REF] Porsching | Estimation of the error in the reduced basis method solution of nonlinear equations[END_REF] Our work [START_REF] Machiels | Output bounds for reduced-basis approximations of symmetric positive definite eigenvalue problems[END_REF][START_REF] Maday | A blackbox reduced-basis output bound method for noncoercive linear problems[END_REF][START_REF] Maday | Global a priori convergence theory for reduced-basis approximation of singleparameter symmetric coercive elliptic partial differential equations[END_REF][START_REF] Prud | Reduced-basis output bounds for approximately parametrized elliptic coercive partial differential equations[END_REF][START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF][START_REF] Veroy | A posteriori error estimation for reduced-basis approximation of parametrized elliptic coercive partial differential equations: "convex inverse" bound conditioners[END_REF] differs from these earlier efforts in several important ways: first, we develop global approximation spaces; second, we introduce rigorous a posteriori error estimators; and third, we exploit off-line/on-line computational decompositions (see [START_REF] Balmes | Parametric families of reduced finite element models: Theory and applications[END_REF] for an earlier application of this strategy.) These three ingredients allow usfor the restricted but important class of "parameteraffine" problems -to reliably decouple the generation and projection stages of reduced-basis approximation, thereby effecting computational economies of several orders of magnitude.

In this paper we develop new a posteriori error estimation procedures for noncoercive linear, and certain nonlinear, problems that -unlike our earlier "asymptotic" techniques [START_REF] Maday | A blackbox reduced-basis output bound method for noncoercive linear problems[END_REF][START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF] -yield rigorous error statements for all N . We consider three particular examples: the Helmholtz (reduced-wave) equation (Section 2); a cubically nonlinear Poisson equation (Section 3); and Burgers equation (Section 4) -a model for incompressible Navier-Stokes. The Helmholtz (and Burgers) example introduce our new lower bound constructions for the requisite inf-sup (singular value) stability factor; the cubic nonlinearity exercises symmetry factorization procedures necessary for treatment of highorder Galerkin summations in the (say) residual dualnorm calculation; and the Burgers equation illustrates our accommodation of potentially multiple solution branches in our a posteriori error statement. Numerical results are presented that demonstrate the rigor, sharpness, and efficiency of our proposed error bounds, and the application of these bounds to adaptive (optimal) approximation.

2 Noncoercive Linear Problems: Helmholtz Equation

Preliminaries

We consider a suitably regular domain Ω ⊂ R d , 1 ≤ d ≤ 3, with boundary ∂Ω. We then introduce a Hilbert space Y with associated inner product, ( • , • ) Y , and induced norm, • Y . We shall assume that H 1 0 (Ω) ⊂ Y ⊂ H 1 (Ω), where H 1 (Ω) ≡ {v ∈ L 2 (Ω), ∇v ∈ (L 2 (Ω)) d }, H 1 0 ≡ {v ∈ H 1 (Ω) | v| ∂Ω = 0}, and L 2 (Ω) is the space of square-integrable functions over Ω. We shall further assume that

( • , • ) Y = ( • , • ) H 1 (Ω) , • Y = • H 1 (Ω) , (1)
where

(w, v) H 1 (Ω) ≡ Ω ∇w • ∇v + wv, ∀ w, v ∈ H 1 (Ω) , v H 1 (Ω) ≡ Ω |∇v| 2 + v 2 , ∀ v ∈ H 1 (Ω) .
(2) More general inner products and norms can (and should) be considered, as discussed in Section 2.4.2.

We shall denote by Y the dual space of Y . For a g ∈ Y , the dual norm is given by

g Y = sup v∈Y g(v) v Y . (3)
If we introduce the "representation" operator

Y : Y → Y such that, for any g ∈ Y , (Yg, v) Y = g(v) , (4)
then g Y = Yg Y ; (5)
this is simply a statement of the Riesz representation theorem.

We now introduce our parametrized bilinear form. We first define a parameter set D µ ⊂ R P , a typical point in which -our input P -tuple -shall be denoted µ; we can then define, for any µ ∈ D µ , our bilinear form a( • , • ; µ) : Y × Y → R. We shall assume that a satisfies a continuity and inf-sup condition for all µ ∈ D, as we now state more precisely.

It shall prove convenient to state our hypotheses in terms of a "supremizing" operator T µ : Y → Y . In particular, for any given µ ∈ D µ , and any w ∈ Y ,

(T µ w, v) Y = a(w, v; µ), ∀ v ∈ Y ; (6)
it is readily shown that

T µ w = arg sup v∈Y a(w, v; µ) v Y . (7)
Furthermore, if we define the inf-sup (singular value) and continuity constants as

β(µ) ≡ inf w∈Y sup v∈Y a(w, v; µ) w Y v Y (8)
and

γ(µ) ≡ sup w∈Y sup v∈Y a(w, v; µ) w Y v Y , (9)
then,

β(µ) = inf w∈Y σ(w; µ) , (10)
γ(µ) = sup w∈Y σ(w; µ) , (11)
where

σ(w; µ) ≡ T µ w Y w Y . ( 12
)
Our assumptions are then: for some positive constant

ε s , ε s ≤ β(µ) ≤ γ(µ) < ∞, ∀ µ ∈ D µ . We next define the bilinear form b( • , • ; µ) : Y × Y → R as b(w, v; µ) = (T µ w, T µ v) Y , ∀ w, v ∈ Y . (13)
We then introduce the eigenproblem:

Given µ ∈ D µ , find χ i (µ) ∈ Y, λ i (µ) ∈ R, i = 1, . . . , ∞, such that b(χ i (µ), v; µ) = λ i (µ)(χ i (µ), v) Y , ∀ v ∈ Y , ( 14
) χ i (µ) Y = 1 . (15)
We shall, for convenience, assume that the spectrum is discrete (in actual practice we require only that the first few modes belong to the discrete component). In that case, we may assume that

b(χ i (µ), χ j (µ); µ) = λ i (µ)(χ i (µ), χ j (µ)) Y = λ i (µ)δ ij , (16)
where δ ij is the Kronecker-delta symbol; that 0 < λ 1 (µ) ≤ λ 2 (µ) ≤ • • • ; and that Y = span {χ i (µ), i = 1, . . . , ∞}. Note that, from ( 10)-( 14), β(µ) = λ 1 (µ); furthermore, γ(µ) is an upper bound for the spectrum.

We shall make the further assumption that a is "affine in the parameter" in the sense that, for some finite Q,

a(w, v; µ) = Q q=1 Θ q (µ) a q (w, v) , (17)
where Θ : D µ → R Q are differentiable parameterdependent coefficient functions, and the a q : Y × Y → R, 1 ≤ q ≤ Q, are parameter-independent bilinear forms. We define, for future reference,

D qp = max µ∈D µ ∂Θ q ∂µ p (µ) , (18)
for 1 ≤ q ≤ Q, 1 ≤ p ≤ P . Furthermore, we assume that the a q are continuous in the sense that there exist positive finite constants Γ q , 1 ≤ q ≤ Q, such that

|a q (w, v)| ≤ Γ q |w| q |v| q ; (19) here | • | q : H 1 (Ω) → R are seminorms that satisfy Q q=1 |v| 2 q 1/2 ≤ C 1/2 Y v Y , ∀ v ∈ Y , (20)
where C Y is a finite constant. Finally, it directly follows from ( 6) and ( 17) that, for any w ∈ Y , T µ w ∈ Y may be expressed as

T µ w = Q q=1 Θ q (µ) T q w , (21)
where, for any w ∈ Y , T q w, 1 ≤ q ≤ Q, is given by

(T q w, v) Y = a q (w, v), ∀ v ∈ Y . (22)
Note that the operators T q : Y → Y are independent of the parameter µ.

Problem Formulation

Weak Statement

We introduce an output functional ∈ Y and "data" functional f ∈ Y . Our weak statement of the partial differential equation is then: Given

µ ∈ D µ , find s = (u(µ)) , (23)
where

u(µ) ∈ Y satisfies a(u(µ), v; µ) = f (v), ∀ v ∈ Y . (24)
In the language of the introduction, s(µ) is our output, and u(µ) is our field variable.

In actual practice, we shall replace ( 23)-( 24) with a truth approximation: Given µ ∈ D µ , find

s N (µ) = (u N (µ)) , where u N (µ) ∈ Y N ⊂ Y satisfies a(u N (µ), v; µ) = f (v), ∀ v ∈ Y N , (25)
and Y N is a finite element approximation subspace. We assume that N is chosen sufficiently large that s N (µ) and u N (µ) may be effectively equated with s(µ) and u(µ), respectively. We shall thus distinguish between Y N and Y only in our discussion of computational complexity. (Note that issues associated with a possible continuous component to the spectrum of ( 14) may be addressed by considering Y as the limit of Y N , N → ∞.)

Reduced-Basis Approximation

The focus of the current paper is a posteriori error estimation. We shall thus take our reduced-basis approximation as given. In particular, we assume that we are provided with a reduced-basis approximation to u(µ), u N (µ) ∈ W N , where

W N = span {ζ n ≡ u(µ n ), 1 ≤ n ≤ N } , (26)
S N = {µ 1 ∈ D µ , . . . , µ N ∈ D µ }, and u(µ n ) satisfies (24) (in practice, (25)) for µ = µ n . It follows that u N (µ) may be expressed as

u N (µ) = N n=1 u N n (µ) ζ n . (27)
The reduced-basis approximation to the output s(µ), s N (µ), is given by s N (µ) = (u N (µ)).

For the purposes of this paper, we shall consider only standard Galerkin projections: a(u N (µ), v; µ) = f (v), ∀ v ∈ W N . However, the discrete inf-sup parameter associated with the latter may not be "good," with corresponding detriment to the accuracy of u N (µ) and hence s N (µ). More sophisticated minimumresidual [START_REF] Maday | A blackbox reduced-basis output bound method for noncoercive linear problems[END_REF][START_REF] Rovas | Reduced-Basis Output Bound Methods for Parametrized Partial Differential Equations[END_REF] and in particular Petrov-Galerkin 7, 18 approaches restore (guaranteed) stability, albeit at some additional complexity and cost.

Error Estimation: Objective

We now wish to develop a posteriori error bounds ∆ N (µ) and ∆ s N (µ) such that

e(µ) H 1 (Ω) ≤ ∆ N (µ) , (28)
and

|s(µ) -s N (µ)| = | (e(µ))| ≤ ∆ s N (µ) , (29)
where e(µ) ≡ u(µ) -u N (µ). For the purposes of this paper, we shall focus on the H 1 (Ω) bound, ∆ N (µ), in terms of which ∆ s N (µ) can be expressed as Y Y ∆ N (µ); the latter may be significantly improved by the introduction of adjoint techniques. [START_REF] Giles | Superconvergent lift estimates through adjoint error analysis[END_REF][START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF] It shall prove convenient to introduce the notion of effectivity, defined (here) as

η N (µ) ≡ ∆ N (µ) e(µ) H 1 (Ω) . ( 30
)
Our certainty requirement (28) may be stated as η N (µ) ≥ 1, ∀ µ ∈ D µ . However, for efficiency, we must also require η N (µ) ≤ C η , where C η ≥ 1 is a constant independent of N and µ; preferably, C η is close to unity, thus ensuring that we choose the smallest Nand hence most economical -reduced-basis approximation consistent with the specified error tolerance.

A Posteriori Error Estimation

Error Bound

We assume that we are given a β(µ) such that, for the given inner product ( • , • ) Y ≡ ( • , • ) H 1 (Ω) (which in our previous papers [START_REF] Prud | Reduced-basis output bounds for approximately parametrized elliptic coercive partial differential equations[END_REF][START_REF] Veroy | A posteriori error estimation for reduced-basis approximation of parametrized elliptic coercive partial differential equations: "convex inverse" bound conditioners[END_REF] would be denoted a "bound conditioner"),

β(µ) ≥ β(µ) ≥ (1 -τ ) ε s , ∀ µ ∈ D µ , (31)
where τ ∈ ]0, 1[ . We then define our error bound as

∆ N (µ) ≡ Yr( • ; µ) Y β(µ) , (32)
where

r(v; µ) = f (v) -a(u N (µ), v; µ), ∀ v ∈ Y , (33)
is the residual associated with u N (µ). Note it follows from (24) that (33) may be restated as

a(e(µ), v; µ) = r(v; µ), ∀ v ∈ Y , (34)
where we recall that e(µ) ≡ u(µ) -u N (µ). We can then state Proposition 1 For the error bound ∆ N (µ) of (32), the effectivity satisfies

1 ≤ η N (µ) ≤ γ(µ) (1 -τ ) ε s , ∀ µ ∈ D , (35)
for all N ∈ N.

Proof It follows from ( 4), (6), and (34) that

Yr( • ; µ) Y = T µ e(µ) Y . (36)
Furthermore, from (12) we know that

e(µ) Y = T µ e(µ) Y σ(e(µ); µ) , (37)
and hence from (1), (30), (32), (36), and (37)

η N (µ) = σ(e(µ); µ) β(µ) . ( 38
)
The result then directly follows from (10), (11), (31), and (38).

We note that our proof (or bound) does not exploit any special properties of e(µ) (or u N (µ)).

It remains to develop our lower bound construction, β(µ), and to demonstrate that both β(µ) and Yr( • ; µ) Y may be computed efficiently (that is, in complexity independent of N ).

Inf-Sup Lower Bound Construction

Many of the most obvious eigenvalue approximation concepts are not relevant here, since we require a lower, not upper, bound. We thus develop a construction particularly suited to our context.

We assume that we are given a set of J parameter points, L J ≡ {µ 1 ∈ D µ , . . . , µ J ∈ D µ }, and associated set of polygonal regions R µ j ,τ , 1 ≤ j ≤ J, where

R µ,τ ≡ {µ ∈ D µ |B µ q (µ) ≤ τ C Y β(µ), 1 ≤ q ≤ Q} , (39) and
B µ q (µ) = Γ q P p=1 D qp |µ p -µ p | ; (40)
we further assume that

J j=1 R µ j ,τ = D µ . (41)
We then define J : D µ → {1, . . . , J} such that, for a given µ, R µ J (µ) ,τ is that region (or a selected region) which contains µ.

Our lower bound is then:

Given µ ∈ D µ , β(µ) = β(µ J (µ) ) -C Y B µ J (µ) max (µ) , (42)
where

B µ max (µ) = max q∈{1,...,Q} B µ q (µ) (43)
for B µ q (µ) defined in (40). We can now state Proposition 2 The construction β(µ) of (42) satisfies the inequality (31).

Proof We first note that, for given µ and µ ∈ D µ , and any w ∈ Y ,

T µ w = T µ w + Q q=1 (Θ q (µ) -Θ q (µ)) T q w ,
where we have appealed to (21). It thus follows from (10), (12), and the triangle inequality that

β 2 (µ) ≥ inf w∈Y σ(w; µ) - Q q=1 (Θ q (µ) -Θ q (µ)) T q w Y w Y 2 . ( 44
)
It is then a simple matter to show, from ( 22), ( 18), (19), the Cauchy-Schwarz inequality, and (20), that

Q q=1 (Θ q (µ) -Θ q (µ)) T q w Y ≤ C Y B µ max (µ) w Y ;
(45) it thus follows from ( 10), (11), and ( 44), (45) that

β 2 (µ) ≥ inf t∈[β(µ),γ(µ)] t -(C Y B µ max (µ)) 2 . (46)
Recall that ( 46) is valid for any µ and µ ∈ D µ . We now choose, for any given µ, µ = µ J (µ) . We thus obtain

β 2 (µ) ≥ inf t∈[β(µ J (µ) ),γ(µ J (µ) )] t -(C Y B µ J (µ) max (µ)) 2 .
(47) However, from (39) and (43) we know that

C Y B µ J (µ) max (µ) ≤ τ β(µ J (µ) ) ; ( 48
)
the infimizer of (47) is thus β(µ J (µ) ), yielding

β(µ) ≥ β(µ J (µ) ) -C Y B µ J (µ) max (µ) . (49)
The desired result, (31), immediately follows from (49), ( 42), (48), and β(µ) ≥ ε s .

It should be clear that our bound of |Θ(µ) -Θ(µ)| is rather crude; a more careful treatment of this term, leading to correspondingly larger regions R µ j ,τ , and hence smaller J, is described elsewhere. 10 It may appear paradoxical to combine a linear approximation -in particular, with O(1) error -to β(µ) with an exponentially convergent reduced-basis approximation -with very small error -to u(µ). In fact, it is not inconsistent:

|β(µ) -β(µ)|/β(µ) ∼ O(1)
is manifested as a 100% "error in the error" -and is acceptable; in contrast, |s(µ) -

s N (µ)|/s(µ) ∼ O(1)
is manifested as a 100% error in the solution -and is thus clearly unacceptable. In essense, the equation for the error e(µ), (34), permits relaxations -and hence rigorous yet inexpensive bounds -that can not be directly applied to the original equation for u(µ), (24).

Offline/Online Computational Procedure

Summary. The central computational aspect of our reduced-basis approach is an offline/online computational decomposition which separates the requisite calculations into two distinct stages. The complexity of the offline -or preprocessing -stage will depend on N (large), N , Q, and J; however, the complexity of the online stage -in which, given a new value of µ, we evaluate s N (µ) and ∆ s N (µ) -will depend only on N , Q, and J. The absence of N dependence in the online stage translates, in many cases, into real-time response.

It is simple to show 15 that, for our Galerkin approximation, the online cost to evaluate

u N n (µ), 1 ≤ n ≤ N , and s N (µ) is O(QN 2 ) + O(N 3
) and O(N ), respectively. We develop here similar estimates for ∆ N (µ) (and hence ∆ s N (µ)). In particular, we shall show that, in the online stage, Yr( • ; µ) Y and β(µ) may be calculated in only O(Q 2 N 2 ) and O(P log J) operations, respectively. We also briefly address the associated offline complexity.

Calculation of the Dual Norm of the Residual . We first invoke ( 17) and ( 27) to write the residual (33) as

r(v; µ) = f (v) - Q q=1 N n=1 Θ q (µ)u N n (µ)a q (ζ n , v) . (50)
It then follows from linear superposition that

Yr( • ; µ) = ẑ00 - Q q=1 N n=1 Θ q (µ)u N n (µ)ẑ qn (51)
where

(ẑ 00 , v) Y = f (v), ∀ v ∈ Y , and
(ẑ qn , v) Y = a q (ζ n , v), ∀ v ∈ Y , (52)
for 1 ≤ q ≤ Q, 1 ≤ n ≤ N . We now insert the expression for Yr( • ; µ) of (51) into the definition of the • Y norm to obtain

Yr( • ; µ) 2 Y = (ẑ 00 , ẑ00 ) Y -2 Q,N q,n=1 Θ q (µ)u N n (µ)(ẑ qn , ẑ00 ) Y + Q,N q,n=1 Q ,N q ,n =1 Θ q (µ)Θ q (µ)u N n (µ)u N n (µ)(ẑ qn , ẑq n ) Y .
The offline/online decompostion is now clear.

In the offline stage, we evaluate and store (ẑ 00 , ẑ00 ) Y , (ẑ qn , ẑ00 ) Y , and

(ẑ qn , ẑq n ) Y , 1 ≤ q, q ≤ Q, 1 ≤ n, n ≤ N : the cost is O(N Q 2 N 2 ) operations.
Then, in the online stage, to compute Yr( • ; µ) Y -given any new value of µ -we need only perform summations over Q and N : the cost is O(Q 2 N 2 ) operations. (In many cases, domain decomposition may be exploited to further reduce the Q dependence -from quadratic to linear. [START_REF] Solodukhov | Reduced-Basis Methods Applied to Locally Non-Affine Problems[END_REF] ) As required, the complexity of the online stage is independent of N .

Calculation of the Inf-Sup Lower Bound . The offline/online decomposition for β(µ) follows directly from the construction. In the offline stage, we must solve J "inf-sup" eigenproblems -calculate λ 1 (µ) from ( 14) for µ = µ j , 1 ≤ j ≤ J: the cost is O(N t J) for some exponent t ≥ 1. In the online stage, we need only determine J (µ) and then evaluate (42): the cost for the former is O(P log J) for a suitable indexing of the R µ j ,τ , 1 ≤ j ≤ J; the cost for the latter is O(P Q). We explore the dependence of J on D µ , ε s , τ , and P in Sections 2.4 and 2.5 in the context of particular examples.

We implicitly assume in our estimates above that we know, a priori , how to choose L J such that (41) will be satisfied. In practice, that will not be the case, and hence additional inf-sup eigenproblems ( 14) will need to be solved. However, for purposes of constructing the sample L J , we may exploit a reduced-basis approximation to (14); [START_REF] Rovas | Reduced-Basis Output Bound Methods for Parametrized Partial Differential Equations[END_REF] this inexpensive yet accurate surrogate greatly reduces the design costs.

Improvements

Model Helmholtz Problem: P = 1

To motivate our improvements, we consider a very simple Helmholtz problem:

Ω ⊂ R d , Y = H 1 0 (Ω), P = 1, and a(w, v; µ) = Ω ∇w • ∇v -µ ω 2 Ω wv . (53)
This model represents harmonic forcing of a "pinned" membrane at frequency ω = √ µ. We may also identify

Q = 2, Θ 1 (µ) = 1, Θ 2 = µ, a 1 (w, v) = Ω ∇w • ∇v, a 2 (w, v) = -Ω wv; furthermore, for |w| 2 1 = Ω |∇w| 2 , |w| 2 2 = Ω w 2 , we readily calculate Γ 1 = 1, Γ 2 = 1, and C Y = 1.
Recall that the latter are for the case in which (

• , • ) Y = ( • , • ) H 1 (Ω) , • Y = • H 1 (Ω) .
It is clear that (53) exhibits resonances at σ 1 , σ 2 , . . ., where 0 < σ 1 ≤ σ 2 < . . . are the real positive eigenvalues of a

1 ( • , • ) relative to -a 2 ( • , • ) (for example, for Ω = ]0, 1[ , σ k = k 2 π 2 ).
It can be further shown that the spectrum of ( 14) is indeed discrete, and that β(µ) is piecewise linear -a "sawtooth." More quantitatively, we obtain

β(µ) = |σ k -µ| σ k + 1 (54)
for µ in a neighborhood of σ k . For purposes of illustration, we shall consider

D µ = [σ * (1 -Λ), σ * (1 -ε s )]
for some small Λ > ε s and for some σ * ≡ σ j * 1. It is then readily shown from (39), (40), and (54) that J -the number of regions R µ j ,τ required to satisfy our "coverage" constraint, (41)is given by

J(σ * , Λ, ε s , τ ) ∼ σ * ln( Λ εs ) τ ( 55
)
for sufficiently large σ * . Recall that τ ∈ ]0, 1[.

Bound Conditioner

The first evident problem with (55) is the dependence on σ * -in particular, since σ * large will often be the case of interest. This can be remedied by better choice of our bound conditioner. To wit, we now define (

• , • ) Y = ( • , • ) μ, • Y = • μ, where ( • , • ) μ ≡ Ω ∇w • ∇v + μ Ω wv , (56)
and

• 2 μ ≡ (•, •) μ; we require that μ ≥ 1.
In what follows subscript μ refers to quantities (previously introduced but now) defined relative to the (•,

•) Y = (•, •) μ and • Y = • μ inner product and norm, respectively.
We first note that e(µ) H 1 (Ω) ≤ e(µ) μ (since μ ≥ 1). We next note that, from the definition of the dual norm, ( 4)-( 5),

Yr( • ; µ) μ2 ≤ Yr( • ; µ) μ1 (57)
for μ2 ≥ μ1 . Finally, it can be shown that (54) must now be replaced with

β μ(µ) = |σ k -µ| σ k + μ . (58)
Assuming (for our purposes here) that β μ(µ)/ βμ (µ) does not depend on μ, it follows from (32), (57), and (58) that ∆ N,μ=σ * ≤ 2∆ N,μ=1 . Therefore, if we bound e(µ) H 1 (Ω) by ∆ N,μ=σ * -which from Proposition 1 is a bound for e(µ) μ=σ * and hence e(µ) H 1 (Ω) -rather than by ∆ N,μ=1 (as before), the effectivity should increase by no more than roughly a factor of two.

However, there will be a significant decrease in the requisite J. In particular, we can now take

|w| 2 1 = Ω |∇w| 2 , |w| 2 2 = μ Ω w 2 , in terms of which Γ 1 = 1, Γ 2 = 1/μ,
and C Y = 1. It then follows from (39), ( 40) and (58) that

J μ=σ * (σ * , Λ, ε s , τ ) ∼ 2 ln( Λ εs ) τ ; ( 59
)
we have successfully eliminated the σ * dependence, at little detriment to the effectivity. The reason is simple: χ 1 (µ) is high-wavenumber, and thus we may add a significant L 2 contribution to our bound conditioner without adversely affecting the inf-sup parameter; this additional L 2 term does, however, significantly improve our continuity constants -on which our lower bound construction is critically dependent. These arguments apply to Helmholtz problems generally; however, for larger ranges of frequency, we will need different bound conditioners for different subdomains of D µ -if we wish to retain the D µ -independence of J.

Deflation

The second debilitating aspect of (55) is the -ln(ε s ) dependence. This, too, can be eliminated, albeit at slightly increased effort. The strategy is deflation: remove the most dangerous components of the errorthose near χ 1 (µ), χ 2 (µ), . . ., -thereby increasing the effective inf-sup parameter. The latter should improve our bounds and effectivity; but, more importantly, it will remove the ε s dependence from (55) -our regions will be generally larger, and will not shrink to zero as we approach resonances (or, at most, except very near resonances).

To begin, we define the "trial" and "test" spaces

U M (µ), V M (µ), with U M (µ) = span{φ i (µ), 1 ≤ i ≤ M }, V M (µ) = span{ξ i (µ), 1 ≤ i ≤ M }. We next introduce the correction δ D M (µ) ∈ U M (µ) through the Petrov-Galerkin projection, a(δ D M (µ), v; µ) = r(v; µ), ∀ v ∈ V M (µ) . ( 60
)
We can then define the "deflated" reduced-basis approximation as

u D N,M (µ) = u N (µ) + δ D M (µ) , (61)
with corresponding output s D N,M (µ) = (u D N,M (µ)). The associated residual is now given by

r D (v; µ) ≡ f (v) -a(u D N,M (µ), v; µ), ∀ v ∈ Y ; (62)
and thus

a(e D (µ), v; µ) = r D (v; µ), ∀ v ∈ Y , (63)
where the deflated error is defined as e D (µ) ≡ u(µ)u D N,M (µ).

For the purposes of this paper, we shall consider a particular set of spaces U M (µ), V M (µ): for a given µ (which shall ultimately depend on µ), we set

U M (µ) = span {χ 1 (µ), . . . , χ M (µ)} , (64)
V M (µ) = span {T µ χ 1 (µ), . . . , T µ χ M (µ)} . ( 65
)
We then note that, from (61)-( 63), (6), and ( 13),

b(e

D (µ), v; µ) = r(T µ v; µ)-a(δ D M (µ), T µ v), ∀v ∈ Y . ( 66
(e D (µ), v; µ) = 0, ∀ v ∈ U M (µ) , (67)
and therefore e D (µ) ∈ Z M (µ) where

Z M (µ) = {v ∈ Y | b(v, χ m (µ); µ) = 0, 1 ≤ m ≤ M }.
We conclude from ( 12), (13), and (67) that

e D (µ) Y ≤ T µ e D (µ) Y β M (µ) , where β M (µ) ≡ inf w∈Z M (µ) σ(w; µ) . (68)
It follows from ( 4), (63), and ( 6) that

Yr D ( • ; µ) Y = T µ e D (µ) Y
and therefore e D (µ) Y ≤ ∆ N,M (µ) where

∆ N,M (µ) ≡ Yr D ( • ; µ) Y βM (µ) ,
and βM (µ) is a lower bound for β M (µ). It remains to construct βM (µ). Revisiting the arguments in the proof of Proposition 2, we note that

β 2 M (µ) ≥ inf w∈Z M (µ) σ(w; µ) - Q q=1 (Θ q (µ) -Θ q (µ)) T q w Y w Y 2 . ( 69
)
It then follows from (68), ( 11), (69), and (45), that for any µ and µ ∈ D µ ,

β 2 M (µ) ≥ inf t∈[β M (µ),γ(µ)] t -(C Y B µ max (µ)) 2 . It thus follows that, if µ ∈ R µ,τ M where R µ,τ M ≡ {µ ∈ D µ | B µ q (µ) ≤ τ C Y λ M +1 (µ), 1 ≤ q ≤ Q} , ( 70
) then βM (µ) ≡ λ M +1 (µ) -C Y β µ max (µ) (71) satisfies β M (µ) ≥ βM (µ) ≥ (1 -τ ) λ M +1 (µ);
here we have invoked (16) to deduce that β(µ) = λ M +1 (µ). This result is much improved over (31) since λ M +1 (µ) (M ≥ 1) will, generically, remain much larger than ε s -indeed, O(1) -for all µ.

In terms of our model problem of Section 2.4.1, we completely eliminate the ε s dependence of (55) (or (59)) -J is now O(1). For this model problem a simple deflation M = 1 suffices. However, for more complicated problems -that exhibit mode-crossing -M = 2 or even M = 3 (near where, say, three modes intersect) higher-order deflations will prove beneficial; we shall observe such a case in Section 2.5.

It remains to address two issues concerning δ D M (µ). First, we must show that δ D M (µ) H 1 (Ω) ≤ Const e(µ) H 1 (Ω) . This is readily demonstratedfor a Const typically considerably better that ε -1 s thanks to the choice of perforce stable Petrov-Galerkin spaces, (64), (65). Second, we must show that δ D M (µ) can be calculated efficiently online. This is indeed the case. The matrix associated with the left-hand side of (60) can be expressed, for the spaces (64), (65), as

a(χ i (µ), T µ χ j (µ); µ) = Q q=1 Θ q (µ)a q (χ i (µ), T µ χ j (µ)) = Q q=1 Q q =1
Θ q (µ) Θ q (µ) a q (χ i (µ), T q χ j (µ)) , (72) which may thus be formed in O(Q 2 M 2 ) operations; here we have invoked ( 17) and (21). Similar arguments indicate that the right-hand side of (60) may be formed in O(Q 2 M N ) operations. We conclude that the additional online complexity is negligible (e.g., compared to calculation of Yr( • ; µ) Y ). †

Numerical Results

Model Helmholtz Problem: P = 2

We consider a model problem which illustrates a potential application of our methodology: real-time, reliable solution of partial differential equations in the service of non-destructive evaluation (and adaptive mission design) "in the field."

Our domain Ω o (L o

c ) is a rectangular membrane [0, 3 2 ] × [0, 1] with a "crack" Γ o crack of length L o c , as shown in Figure 1. We assume that the boundary of the memrane is "pinned" except on the "stress-free" crack. The membrane is forced over a patch Ω o in at frequency ω; the response is measured over the patch Ω o out . Our problem statement is then: Given L o c and † The offline expense will be increased somewhat, not so much due to the λ 2 (µ j ), χ 2 (µ j ) (say for M = 1) -in particular, since J will now be much smaller -but rather due to the JQ 2 (M 2 + M N )N operations required for the inner products associated with the deflation correction (60). Relatedly, the online storage will also increase, by JQ 2 (M 2 + M N ). This effect can be very significantly reduced if we replace χ 1 (µ), χ 2 (µ), . . . in U M (µ), V M (µ) with a reduced-basis approximation to these quantities; [START_REF] Rovas | Reduced-Basis Output Bound Methods for Parametrized Partial Differential Equations[END_REF] however, there will be a concomitant slight increase in online cost.

ω 2 , find

s(ω 2 , L o c ) = 1 |Ω o out | Ω o out u o (ω 2 , L o c ) (73)
where

u o (ω 2 , L o c ) ∈ X o (L o c ) satisfies Ω o (L o c ) ∇u o (ω 2 , L o c ) • ∇v -ω 2 u o (ω 2 , L o c ) v = 1 |Ω o in | Ω o in v, ∀ v ∈ X o (L o c ) ; (74) here X o (L o c ) = {H 1 (Ω o (L o c )) | v| ∂Ω o (L o c )\Γ o crack = 0}.
Clearly, an elastic plate (and more realistic outputs) would be a much more relevant model; our methodology directly applies to this case as well.

We

now map Ω o (L o c ) to a reference domain indepen- dent of L o c , Ω ≡ Ω o (L o c = 1
2 ), through piecewise-affine subdomain co-ordinate tranformations. The resulting equations can then be cast in the desired form ( 23), (24), and ( 17 out map to invariant images in Ω, the reference domain.) We do not give here the detailed expressions for either the Θ q (µ) or a q ( • , • ), 1 ≤ q ≤ Q.

We briefly comment on the structure of this model problem. For a given L o c , the bilinear form is, apart from several scaling factors, identical to (53) of Section 2.4.1. We conclude that, for any given (ω 2 , L o c ), the spectrum of ( 14) is discrete; and furthermore that, for given

L o c , β(ω 2 , L o c
) is a piecewise-linear function of ω 2 . However, if we now permit L o c to also vary, the behavior of β(ω 2 , L o c ) is no longer trivial. In particular, unlike in Section 2.4, we can no longer characterize β(ω 2 , L o c ) in terms of a few (more generally, denumerable) "resonance" eigenvalues -our lower bound constructions are now required. We show in Figure 2 a contour plot of β(ω 2 , L o c ) over D µ ; the dark lines indicate "excluded regions" that contain resonanceslines (more generally, P -1 dimensional manifolds) along which β(µ) vanishes. (The resonances will be eliminated if we incorporate damping or radiation into our model; however, the inf-sup parameter may still remain small.)

The Inf-Sup Lower Bound

For this first set of tests we do not consider a complete "paving" of D µ ; rather L J ≡ L TEST J comprises 25 points µ 1 , . . . , µ J equidistributed along a segment (25, 0.4)(50, 0.6) in D µ . We shall consider the following cases:

I • Y = • H 1 (Ω) , no deflation; II • Y = • μ=ω 2 min =25 , no deflation; III • Y = • μ=ω 2 min =25 , deflation, M = 1; IV • Y = • μ=ω 2 min =25 , deflation, M = 2; V • Y = • μ=ω 2 min =25
, deflation, M = 3. For the definition of • μ and of deflation see Sections 2.4.2 and 2.4.3, respectively. In all cases, we choose τ = 3/4.

We plot in Figure 3 the polygons R µ j ,τ for cases I, II, III, and V. Note we shift vertically the regions for cases I (highest), II, and III for purposes of easy comparison. (In actual fact, for given j, the center of R µ j ,τ for all cases (I, II, II, and V) is µ = µ j ; in Figure 3, only R µ j ,τ for case V is honestly (vertically) located.) First, we observe that the "correct" bound conditioner (I → II) considerably increases the size of the regions; furthermore, this effect will be even more dramatic for higher frequency ranges. Second, we observe that some deflation (II → III) further improves the situation; and sufficient deflation (III → V) greatly improves the situation, in particular as we approach resonance. Note that although IV performs better than III, only with V do we have sufficient deflation in the sense that all dangerous modes are neutralized -it is clear from Figure 3 that three modes are "active" near the end of our segment (25, 0.4)(50, 0.6). Increasing M beyond 3 has little effect as all modes now appear "far away."

In short, a combination of "tuned" bound conditioners and sufficient deflation greatly increases the size of our regions -in particular at high frequencies and near resonance, respectively -such that we can expect J to be roughly independent of D µ and ε s . We show in Figure 4, for case V, a (more or less) complete paving of D µ by regions R µ j ,τ , 1 ≤ j ≤ J = 625. Note we anticipate that J can be further reduced by roughly 5-10 based on less conservative polygons which exploit the monotonicity of the Θ q (µ). 10 The result of Figure 4 is disappointing in one aspect. In actual fact, β(µ) varies significantly (locally) only in the one direction perpendicular to the P -1 dimensional "resonance" manifolds. Our construction, even with deflation, remains isotropic, and thus J ∼ (Const > 1) P . In contrast, an anisotropic construction would be relatively insensitive to P . Unfortunately, at present, we see no way to capture this anisotropy efficiently while maintaining the requisite lower bound property.

Error Bounds and Effectivity

We consider here a point µ TEST which lies within a region R µ j ,τ for all cases I, II, III, IV, and V. This point, µ TEST = (42.95, 0.55), is quite close to a resonance.

We present in Tables 1 and2 e (D) (µ TEST ) H 1 (Ω) / u(µ TEST ) H 1 (Ω) and η N (µ TEST ), respectively, as a function of N for Cases I, II, and V. Note S N , our reduced-basis approximation sample, is chosen here log-randomly over D µ ; see Sections 3 and 4 for discussion of better optimal/adaptive alternatives for sample design. As expected from the arguments of Section 2.4.2, the bound conditioner has little (detrimental) effect on the effectivity. And, as expected from the arguments of Section 2.3.3, deflation has a modest (respectively, significant) positive effect on the error (respectively, effectivity).

e (D) (µ TEST ) H 1 (Ω) / u(µ TEST ) H 1 (Ω) N I II V 4 2.7 × 10 -1 2.7 × 10 -1 8.4 × 10 -2 8 9.6 × 10 -2 9.6 × 10 -2
3.0 × 10 -2 12 6.2 × 10 -2 6.2 × 10 -2 5.8 × 10 -2 16 9.0 × 10 -3 9.0 × 10 -3 8.5 × 10 -3 20 3.5 × 10 -4 3.5 × 10 -4

3.5 × 10 -4 24 2.7 × 10 -4 2.7 × 10 -4 2.7 × 10 -4 28 1.2 × 10 -4 1.2 × 10 -4

1.2 × 10 -4 32 2.1 × 10 -5 2.1 × 10 -5

2.1 × 10 -5 36 8.9 × 10 -6 8.9 × 10 -6

8.9 × 10 -6 40 3.0 × 10 -6 3.0 × 10 -6

3.0 × 10 -6

Table 1 The normalized error as a function of N for cases I, II, and V.

η N (µ) We close with two related remarks. First, even if the regions R µ j ,τ include a resonance, this does not imply that the error (or error bound) remains finite as we approach the resonance. Second, round-off errors will become increasingly important, and ultimately dominant, in the very immediate vicinity of resonances; in particular, as we approach extremely close to a resonance, we may observe effectivities below unity. The reason is clear: our error bound assumes (67); however, in finite precision, this condition will be violated -and the resulting "round-off" error amplified by 1/β(µ). Exact orthogonalization recovers the theoretical result -η N (µ) ≥ 1; in more realistic models, damping will provide the necessary "cut-off."

N I II V 4

Cubically Nonlinear Poisson Problem

Preliminaries

We consider a suitably regular domain Ω ⊂ R 2 with boundary ∂Ω = ∂Ω R ∪ ∂Ω N . We set Y = H 1 (Ω), and 2)). The dual space of Y will be denoted Y , with norm defined as in (3) or, equivalently, (4), (5).

( • , • ) Y = ( • , • ) H 1 (Ω) and • Y = • H 1 (Ω) (as defined in (
We next define our parameter set (µ

1 , µ 2 ) ∈ D µ ≡ [µ min 1 > 0, µ max 1 ] × [µ min 2 > 0, µ max 2 ] ⊂ R P =2 + . Then, for any µ ∈ D µ , a( • , • ; µ) : Y × Y → R is given by a(w, v; µ) = a L (w, v; µ) + a NL (w, v), ∀ w, v ∈ Y , (75)
where

a L (w, v; µ) = µ 1 Ω ∇w • ∇v + µ 2 ∂ΩR wv , (76)
and

a NL (w, v) = Ω w 3 v (77)
represent a "Poisson-Robin" operator and cubic nonlinearity, respectively. We note that a L is continuous,

sup w∈Y sup v∈Y a L (w, v; µ) w Y v Y ≤ γ(µ) < ∞, ∀ µ ∈ D µ , (78)
symmetric, and coercive,

0 < ε s ≤ α(µ) = inf w∈Y a L (w, w; µ) w 2 Y , ∀ µ ∈ D µ . (79)
Furthermore, a L ( • , • ; µ) depends affinely on the parameter µ -the expansion (17) applies with Q = 2 and Θ q = µ q . In fact, our treatment of this section applies to any continuous, coercive, affine a L ( • , • ; µ).

Problem Formulation

Weak Statement

We introduce an output functional ∈ Y and "data" functional f ∈ Y ; for our model problem we take (v) = f (v) = ∂ΩN v. Our weak statement of the partial differential equation is then: Given µ, find

s = (u(µ)) , (80)
where

u(µ) ∈ Y satisfies a(u(µ), v; µ) = f (v), ∀ v ∈ Y . (81)
In the language of the introduction, s(µ) is our output, and u(µ) is our field variable. It can be shown ? that (75), (81) admits a unique solution.

As for our Helmholtz problem, in actual practice we replace s(µ) and u(µ) with corresponding "truth" Galerkin approximations s N (µ) and u N (µ), respectively (see Section 2.2.1).

Reduced-Basis Approximation

The focus of the current paper is a posteriori error estimation. We shall thus take our reduced-basis approximation as given. In particular, we assume that we are provided with a reduced-basis approximation to u(µ), u N (µ) ∈ W N , where

W N = span {ζ n ≡ u(µ n ), 1 ≤ n ≤ N } , (82)
S N = {µ 1 ∈ D µ , . . . , µ N ∈ D µ }, and u(µ n ) satisfies (81) for µ = µ n . It follows that u N (µ) may be expressed as

u N (µ) = N n=1 u N n (µ) ζ n . (83)
The reduced-basis approximation to the output s(µ), s N (µ), is given by s N (µ) = (u N (µ)).

For our equations ( 75), (81), standard Galerkin projection is the best choice -as we shall see, there are no stability issues -and we thus select this (simplest) option:

a(u N (µ), v; µ) = f (v), ∀ v ∈ W N .

Error Estimation: Objective

As for Helmholtz, we wish to provide an a posteriori error bound ∆ N (µ) for e(µ) H 1 (Ω) such that the effectivity, (30), satisfies

1 ≤ η N (µ) ≤ C η , (84)
for C η independent of N and µ -and preferably close to unity. Error bounds for |s(µ) -s N (µ)| may also be developed.

A Posteriori Error Estimation

Error Bound

We assume that we are given an α :

D µ → R + such that α(µ) ≥ α(µ) ≥ (1 -τ ) ε s , ∀ µ ∈ D µ , (85)
for given τ ∈ ]0, 1[ . We then define our error bound as

∆ N (µ) ≡ Yr( • ; µ) Y α(µ) , (86)
where

r(v; µ) ≡ f (v) -a(u N (µ), v; µ), ∀v ∈ Y , (87)
is the residual. We can then state Proposition 3 For the error bound ∆ N (µ) of ( 86),

e(µ) H 1 (Ω) ≤ ∆ N (µ), ∀ µ ∈ D µ , (88)
for all N ∈ N.

Proof We know from (75), (76), (77) that (87) may be written as

a L (e(µ), v; µ)+ Ω u 3 (µ)-u 3 N (µ) v = r(v; µ), ∀v ∈ Y , (89)
where e(µ) ≡ u(µ) -u N (µ). We now take v = e(µ), note that

Ω u 3 (µ) -u 3 N (µ) (u(µ) -u N (µ)) ≥ 0 , (90)
and invoke ( 79) and ( 85) to obtain α(µ) e(µ) 2 Y ≤ r(e(µ); µ) .

The desired result then follows from ( 5) and (91) (recall that

• Y ≡ • H 1 (Ω) ).
We do not include here a uniform upper bound for the effectivity; however, it is clear from (89) that, as r( • ; µ) Y → 0, η N (µ) ≤ γ(µ)/(1 -τ )ε s .

Coercivity Lower Bound Construction

Our approach to the inf-sup lower bound, described in Section 2.3.2, can also be adapted to general coercive problems. 10 For our purposes here, however, we consider a simple variant that exploits the monotonicity of α(µ).

In particular, it can be shown that, for a L ( • , • ; µ) as defined in (76), α(µ 1 ) ≤ α(µ 2 ) for µ 2 ≥ µ 1 (in the sense of each component); the proof follows directly from the Rayleigh quotient definition, (79). Thus, given a sample L J ≡ {µ 1 = µ min , . . . , µ J }, we may define a lower bound as

α(µ) ≡ max {j∈{1,...,J} | µ j ≤µ} α(µ j ) . (92)
The best distribution of points -to minimize J given τ and our requirement (85) -is logarithmic. Further details on these and related bound conditioners for coercive problems may be found elsewhere. [START_REF] Veroy | A posteriori error estimation for reduced-basis approximation of parametrized elliptic coercive partial differential equations: "convex inverse" bound conditioners[END_REF] 3.3.3 Offline/Online Computational Procedure Summary. Our nonlinear problem admits an offline/online decomposition quite similar to that for linear problems. The key new issue is the higher order summations that perforce arise within the Galerkin context. Our focus here will be on efficient (or as efficient as possible) treatment of these new terms.

By way of summary, the online complexity to calculate u N n (µ), 1 ≤ n ≤ N , and subsequently s N (µ) = (u N (µ)), scales as K iter (2N 2 + N 3 + N 4 ); here K iter is the number of Newton iterations required to solve for u N n (µ), 1 ≤ n ≤ N . The N 4 dependence -which arises due to the nonlinearity -is not pleasant, but for reasonably small N , not debilitating. In constrast, the online complexity to calculate ∆ N (µ) scales, to leading order for "large" N , as C 6 N 6 ; the N 6 dependencewhich arises due to the nonlinear contribution to the dual norm of the residual -is now more daunting. Fortunately, C 6 = 1/72, and thus, again for modest N , efficiency is preserved.

We focus our attention here on the dual norm calculation, and in particular on the origin of the constant C 6 .

Calculation of the Dual Norm of the Residual . We first invoke (75), ( 76),(77), and (83) to write

r(v; µ) = f (v) - N n=1 u N n (µ) a L (ζ n , v; µ) - N n,n ,n =1 u N n (µ)u N n (µ)u N n (µ) Ω ζ n ζ n ζ n v .
(93) It thus follows from (76), (93) and linear superposition that

Yr( • ; µ) = ẑL 00 -µ 1 N n=1 u N n (µ)ẑ L 1n -µ 2 N n=1 u N n (µ)ẑ L 2n - N n,n ,n =1 u N n (µ)u N n (µ)u N n (µ) ẑNL nn n , (94
) where (for example) ẑNL nn n ∈ Y satisfies

(ẑ NL nn n , v) Y = Ω ζ n ζ n ζ n v, ∀ v ∈ Y . (95)
We thus obtain

Yr( • ; µ) 2 Y = • • • + N n,n ,n =1 N m,m ,m =1 u N n (µ) • • • u N m (µ) (ẑ NL nn n , ẑNL mm m ) Y , (96)
where we shall focus only on the highest order (and hence most expensive) summations.

Obviously, naïve treatment of (96) directly yields N 6 operations. However, there are many symmetries that can be exploited. In particular, we note from (95) that ẑNL nn n = ẑNL mm m for any mm m triplet which is a permutation of nn n . We denote by P N 3 the set of unique ordered 3-tuples of integers j ∈ {1, . . . , N } -(j, j , j ) such that j ≤ j ≤ j ; by T the cardinality of P N 3 ; and by

Π k = (Π 1 k , Π 2 k , Π 3 k ), 1 ≤ k ≤ T , the members of P N 3 .
We can thus write our sum of (96) as

T k=1 T k =1 C Π k C Π k u N Π 1 k (µ)u N Π 2 k (µ)u N Π 3 k × u N Π 1 k (µ)u N Π 2 k (µ)u N Π 3 k (ẑ NL Π k , ẑNL Π k ) Y , (97)
where the C Π k are multiplicity constants.

The online complexity of (97) clearly scales as T 2 /2. It is readily shown that the cardinality of P N κ , the set of unique ordered κ-tuples of integers j ∈ {1, . . . , N } is given by

T(N, κ) = ((N -1) + κ)! (N -1)! κ! . (98)
In our particular case, T = T(N, 3) ∼ N 3 /6 for N large; our sum (97) may thus be performed in N 6 /72 operations -a considerable improvement over the naïve estimate of N 6 . In general, T(N, κ) ∼ N 6 /κ!: in relative terms, higher order (e.g., u κ ) nonlinearities thus enjoy greater economies; however, in absolute terms, T(N, κ) will grow very rapidly with N for larger κ. We must be content with relatively low-order (or low-order approximations of) nonlinearities.

Numerical Results

Model Problems

Our model problem has already been specified in Sections 3.1 and 3.2. It remains only to specify the physical domain, Ω ⊂ R 2 -a "tee"-shaped region with ∂Ω N at the root -and the parameter domain -µ min 1 = .1, µ max 1 = 10, µ min 2 = .01, µ max 2 = 1. Note the nonlinearity will be most significant for µ 1 and µ 2 small.

Adaptive Reduced-Basis Approximation

Given the higher powers of N that now appear in our complexity estimates, it is crucial (both as regards online and offline effort) to control N more tightly. To this end, we may gainfully apply our a posteriori error bounds adaptively.

We first construct, offline, an approximation that, over most of the domain, exhibits an error (say, here, in the H 1 -norm) less than ε prior d : we begin with a first sample point µ 1 (S N =1 = {µ 1 }); we next (inexpensively) evaluate ∆ N =1 (µ) over a large test sample of parameter points in D µ , Ξ prior ; we then choose for µ 2 (and hence S N =2 = {µ 1 , µ 2 }) the maximizer of ∆ N =1 (µ) over Ξ prior ; we now repeat this process until the maximum of ∆ N =N prior (µ) over Ξ prior is less than ε prior d . Then, online, given a new value of the parameter, µ, and an error tolerance ε post d (µ), we essentially repeat this adaptive process -but now our sample points are drawn from S N prior , and the test sample is a singleton -µ. We typically choose ε prior d ε post d (µ) since our prior test sample is not exhaustive; and therefore, typically, N post (µ) ≤ N prior .

We present in Table 3 the normalized error e(µ * ) Y / u(µ * ) Y , as a function of N , for the (log) random and adaptive sampling processes (note that, in the results for the random sampling process, the sample S N is different for each N ). We also indicate in Table 3 the online computational time, † which is largely independent of the sampling process but very strongly dependent on N (see Section 3.3.3). Here µ * = µ * (N ) is the point in Ξ prior at which the maximum error bound ∆ N (µ) occurs -note µ * (N ) will be different for the different sampling strategies. We observe that the adaptive sampling procedure yields higher accuracy at lower N ; and that even these modest reductions in N can translate into measurable performance improvements. For purposes of comparison, calculation of s N (µ) requires 50 seconds -and thus, even for an accuracy of 2.0%, the reduced-basis approach is two to three orders of magnitude faster (marginally) than conventional techniques.

N

u(µ) Y ∆ N (µ) u(µ) Y η N (µ)
µ TEST,1 1.45 × 10 -4 5.18 × 10 -3 35.7 µ TEST,2 7.37 × 10 -3 4.33 × 10 -2 5.9 µ TEST,3 6.91 × 10 -3 1.52 × 10 -2 2.2 Of course, in actual practice, the savings indicated in Table 3 can only be realized if our error estimators are true bounds (η N (µ) ≥ 1), and good bounds (η N (µ) ≈ 1). We show in Table 4 the (normalized) error bound ∆ N (µ * )/ u(µ * ) Y , and effectivity, η N (µ * ), as a function of N (for the adaptive case); as before, µ * is the point in Ξ prior at which the maximum error bound occurs. We observe that we do indeed obtain bounds, but that the bounds are not too sharp.

The main cause of the higher effectivities is the relatively small value of α(µ) (and hence α(µ)) for low µ 1 and µ 2 . We present in Table 5 the effectivities for three test points, µ TEST,1 = (0.01, 0.1), µ TEST,2 = (0.1, 1), and µ TEST,3 = (1, 10), for the N = 10 adaptive sample; except near µ 2 = .01, the effectivities are quite close to unity. Numerous remedies exist for low µ 1 , µ 2 ; we thus do not dwell on this further here.

The Burgers Equation

Preliminaries

We consider the domain Ω = ]0, 1[ . We set Y = H 1 0 (Ω), and (

• , • ) Y = ( • , • ) H 1 (Ω) , • Y = • H 1 (Ω) .
The dual space of Y will be denoted Y , with the norm defined as in (3), or, equivalently, (4), (5).

In this case we have a single parameter, µ ∈ D µ ≡ [µ min > 0, µ max ] ⊂ R P =1 + . For any µ ∈ D µ . a( • , • ; µ) : Y × Y → R is given by a(w, v; µ) = a L (w, v; µ) + a NL (w, w, v), ∀ v ∈ Y , (99) where a L (w, v; µ) ≡ µa 0 (w, v) = µ

The development parallels Helmholtz except that we must include our linearization point z in the definitions.

Finally, we assume that we are given a constant ρ such that, for all v ∈ Y (= H 1 0 (Ω)),

v L 4 (Ω) ≤ ρ v Y ; ( 108
)
the existence of a finite ρ (for Ω ⊂ R d=1,2,3 ) is guaranteed by the continuous embedding of Y in L 4 (Ω). [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF] In R 1 we also have Y ⊂ L ∞ (Ω); we thus readily derive, from the Cauchy-Schwarz inequality, ρ = 1 2 . In higher dimensions it is not difficult to develop bounds for ρ in general geometries. Υ N (µ)/ u(µ) Y as a function of N . We observe two difficulties not encountered for µ = 1.0. First, although the reduced-basis approximation in fact converges rather quickly, we can only begin to make an a posteriori error statement for N ≥ 11 -for N < 11, condition (119) is not satisfied. (Recall that these results are for the adaptive sampling procedure; in the case of a random sample, condition (119) is not satisfied for all N .) Second, although the "good" choice, (120), in fact obtains -e(µ) Y ≤ ∆ N (µ) with effectivities O(5-10) for N ≥ 11 -the value of Υ N (µ) is not as large as desirable; it would thus be difficult in practice (when e(µ) Y is not known) to unambiguously rule out the "bad" choice (121). The origin of both these difficulties is the small value of β(µ = 0.01) (and β(µ), µ small, generally). We are hopeful that deflation ideas similar to those successful in the Helmholtz case (see Section 2.4.3) will also prove beneficial here -increasing the effective β(µ), and thereby increasing both our threshold in (119) and the value of Υ N (µ) in ( 118), (121).

  ) It then follows from (60), (64), and (65) that b

Fig. 1

 1 Fig. 1 Membrane with crack of length L o c under harmonic excitation at frequency ω.

  ) for P = 2, µ = (ω 2 , L o c ), D µ = {µ ∈ [25, 50] × [.3, .7] | β(µ) ≥ ε s ∼ = 0.005}, and Q = 8. (Note the relatively large value of Q, relative to P , originates in the (accuracy) requirement that Ω o in and Ω o

Fig. 2

 2 Fig. 2 Contours of β(ω 2 , L o c ) over D µ for the (cracked membrane) model Helmholtz problem.

Fig. 3

 3 Fig.3Polygons R µ j ,τ for cases I, II, III, and V.

Fig. 4 A

 4 Fig. 4 A (more or less) complete "paving" of D µ for case V.

  trilinear forms, respectively. For a given z ∈ Y , we define the bilinear form -associated to the derivative of our operator -d( • , • ; z; µ) : Y × Y → R as d(w, v; z; µ) = a L (w, v; µ) + 2a NL (z, w, v) . (102)It shall prove convenient to introduce the supremizing operator T z;µ : Y → Y such that, for any w ∈ Y ,(T z;µ w, v) Y = d(w, v; z; µ), ∀ v ∈ Y ;(103)it is readily shown thatT z;µ w = arg sup v∈Y β z (µ) = inf w∈Y σ z (w; µ) ≤ sup w∈Y σ z (w; µ) = γ z (µ) ,where σ z (w; µ) ≡ T z,µ w Y w Y .

Fig. 5

 5 Fig.5The inf-sup parameter β(µ), lower bound β(µ), and segments R µ j ,τ =3/4 = [µ j , µ j+1 ] for the model Burgers problem.

Table 2

 2 The effectivity ηN (µ) as a function of N for cases I, II, and V.

		26.7	35.3	9.0
	8	29.3	41.1 10.1
	12 129.4 145.1 12.2
	16 101.3 132.0 11.0
	20	98.2 141.0 11.1
	24	97.9 143.1 11.2
	28 100.6 148.3 11.6
	32 102.8 153.5 12.0
	36 103.2 154.2 12.1
	40 102.3 153.0 12.0

Table 3

 3 Error bound for the cubically nonlinear Poisson problem for random and adaptive samples.

	N	∆ N (µ * ) u(µ * ) Y	η N (µ * )
	2 7.42 × 10 +1	169.3
	4 1.01 × 10 +1	79.3
	6 4.11 × 10 +0	56.6
	8 1.16 × 10 +0	22.8
	10 5.34 × 10 -1	25.6

Table 4

 4 Error bound and effectivity for the cubically nonlinear Poisson problem for an adaptive sample.

µ e(µ) Y

Table 5 Error

 5

, error bound, and effectivity for the cubically nonlinear Poisson problem for the N = 10 adaptive sample; µ TEST,1 = (0.01; 0.1), µ TEST,2 = (0.1; 1), and µ TEST,3 = (1; 10).

Table 6

 6 Error, error bounds, and effectivity as a function of N for the Burgers problem with µ = 1.0.

of 18 American Institute of Aeronautics and Astronautics Paper

† The calculations were performed on a Pentium r 4 2.4GHz processor.

Acknowledgements

We would like to acknowledge our longstanding collaborations with Professor Yvon Maday of University

Problem Formulation

Weak Statement

We introduce an output functional ∈ Y and "data" functional f ∈ Y . Our weak statement of the partial differential equation is then: Given µ, find

where u(µ) ∈ Y satisfies

Equations ( 99), (100), ( 101), (110) are a very good model for the incompressible Navier-Stokes equations (see below), which is our ultimate goal. For sufficiently large µ, (99), (110) -and the incompressible Navier-Stokes equations -have a unique solution; for smaller µ, we can encounter non-uniqueness -multiple solution branches may exist. As for our Helmholtz problem, in actual practice we replace s(µ) and u(µ) with corresponding "truth" Galerkin approximations s N (µ) and u N (µ), respectively (see Section 2.2.1).

Reduced-Basis Approximation

We assume that we are provided with a reducedbasis approximation to u(µ), u N (µ) ∈ W N , where

S N = {µ 1 ∈ D µ , . . . , µ N ∈ D µ }, and u I (µ n ) satisfies (110) for µ = µ n . It follows that u N (µ) may be expressed as

The reduced-basis approximation to the output s(µ), s N (µ), is given by s N (µ) = (u N (µ)). Note u I (µ n ) refers to solutions of (99), (110), which are assumed to reside on a "first" (particular) branch; although we do not dwell here on possible bifurcation structure, other "parametric manifolds" (say, u II (µ)) may, in general, exist.

For the purposes of this paper, we shall consider only standard Galerkin projections:

However, the discrete inf-sup parameter associated with the latter may not be "good," with corresponding detriment to the accuracy of u N (µ) and hence s N (µ). More sophisticated minimumresidual [START_REF] Maday | A blackbox reduced-basis output bound method for noncoercive linear problems[END_REF][START_REF] Rovas | Reduced-Basis Output Bound Methods for Parametrized Partial Differential Equations[END_REF] and in particular Petrov-Galerkin 7, 18 approaches restore (guaranteed) stability, albeit at some additional complexity and cost.

We comment that, for the case in which geometry is fixed and only viscosity varies, our reduced-basis approximation (and associated error estimation) procedure for the Burgers equation directly translates to the full incompressible Navier-Stokes equations -in particular, a divergence-(and hence pressure-) free formulation of the incompressible Navier-Stokes equations.

Error Estimation: Objective

As for the cubically nonlinear Poisson problem, we would like to provide an error bound ∆ N (µ) for e(µ) H 1 (Ω) (and, relatedly, bound ∆ s N (µ) for |s(µ)s N (µ)|) such that the effectivity satisfies (84). However, as we shall see, our error statement will no longer be unqualified -there will be a "choice" that reflects the possible existence of multiple solution branches.

A Posteriori Error Estimation

Preliminaries

We first define, in a slight abuse of notation, T µ ≡ T u N (µ);µ where T z;µ is given by (103); relatedly, we define β(µ) ≡ β u N (µ) (µ), γ(µ) ≡ γ u N (µ) (µ), and σ(w; µ) = σ u N (µ) (w; µ), for β z (µ), γ z (µ), and σ z (w; µ) defined in (105), (106), and (107). We assume that

In what follows, we will explicitly highlight the N -dependence of β(µ), γ(µ), and σ(w; µ) only in those places where this dependence is either not obvious or potentially problematic.

We shall also require operators

where

and

It follows from ( 102), ( 103), ( 113), (114), and (115) that

Note that T 0 and the T n are parameter-independent.

Error Bound

We assume that we are given a β(µ) such that

where

We then define, for Yr(

We note that, as Yr(

thus, ∆ N (µ) → 0 but Υ N (µ) → Const as the residual vanishes.

We can then state Proposition 4 Given µ ∈ D µ , for N sufficiently large such that

where ∆ N (µ) and Υ N (µ) are given by ( 117) and (118), respectively.

Proof It is a simple matter to show that e(µ) = u(µ)u N (µ) satisfies

(122) Since from (101) and (108)

it follows from (122) (with v = T µ e(µ)), (103), (105), and (116) that

The desired result directly follows from solution of this quadratic equation for e(µ) Y .

We note that an alternative proof -which directly places a restriction on e(µ) Y that is subsequently self-consistently determined from the strength of the nonlinearity -is applicable to much more general nonlinearities. However, in the quadratic case, the proof above is simpler and slightly sharper. We do not include here a uniform upper bound for the effectivity, however, it is clear from (123) that, as Yr( • ; µ) Y → 0, η N (µ) ≤ γ(µ)/(1 -τ )ε s (recall that both γ(µ) and ε s (µ) may depend on u N (µ)).

We now turn to an interpretation of the "choice" (120), (121). It is clear that, given W N , (111) (and hence (112)), all evidence would suggest that, as Yr( • ; µ) Y → 0, u N (µ) should well approximate u I (µ); thus (120) -note ∆ N (µ) → 0 as N → ∞ -is the most obvious choice for u I (µ) -u N (µ) Y . However, equally clearly, if a second branch, u II (µ), exists, there is no reason that u N (µ) should -in fact, there is every reason that u N (µ) should not -well approximate u II (µ); thus, for this (possible) second branch, (121) -note Υ N (µ) → Const as N → ∞ -is the most obvious choice. In short, the error bound "sees" only the residual, which in turn "sees" only the branch-independent projection of u I (µ) (or u II (µ)), f (v). Thus, absent other a priori information, the Υ N (µ) option is a nonlinear necessity -a reflection of the potential existence of distinct multiple solutions. † Consistent with these arguments, we note that if d( • , • ; u N (µ); µ) is coercive -certainly the case for sufficiently large µ -then it follows directly from (122), since a NL (e(µ), e(µ); e(µ)) = 0, that we may obtain an unconditional bound on the error in u N (µ) relative to the perforce single branch, u (I) (µ).

However, there is a dark side: we can not rigorously preclude the possibility that u I (µ)-u N (µ) Y ≥ Υ N (µ). Although this is extremely unlikely as N → ∞ -since Υ N (µ) → Const as N → ∞ -it can not be unambiguously ruled out for any fixed N . Clearly, in actual practice, the relative (and absolute) magnitude of Υ N (µ) will directly affect our comfort level in choosing (120). We discuss this again in the context of our numerical results.

Inf-Sup Lower Bound Construction

We assume that we are given a set of J parameter points, L J ≡ {µ 1 ∈ D µ , . . . , µ J ∈ D µ }, and associated segments R µ j ,τ , 1 ≤ j ≤ J, where

and

we further assume that

We then define J : D µ → {1, . . . , J} such that, for a given µ, R µ J (µ) ,τ is that segment (or a segment) which contains µ.

Our lower bound is then:

127) † Note if we include both (say, in the case of two branches) branches, u I (µ n ), u II (µ n ), in W N , then we will typically obtain good reduced-basis approximations to both branches -u I N (µ), u II N (µ). In this case our single bound ∆ N (µ) (→ 0) would apply to both u I (µ) -u I N (µ) Y and u II (µ) -u II N (µ) Y , and Υ N (µ) (→ Const) would apply to both u

for B µ (µ) defined in (125).

We can now state Proposition 5 The construction β(µ) of (127) satisfies the inequality (116).

Proof (Sketch)

The proof is almost identical to the proof of Proposition 2 for the Helmholtz inf-sup lower bound construction. We need only replace the relation ( 44) and (45) with

(128) and

. The continuity result (129) follows from (113)-( 115), from

To derive this last expression we invoke (112), (115), Cauchy-Schwarz, and (108). It thus follows from (105), (106), (128), and (129) that

(130) We now choose, for any given µ, µ = µ J (µ) . We then note that, from (124) and (125), the infimizer of (130) is β(µ J (µ) ); it thus follows that

The desired result, (116), immediately follows from (131), ( 127), (124), and β(µ) ≥ ε s .

Offline/Online Computational Procedure

All the elements of the offline/online procedure for the construction of Burgers a posteriori error bounds have already been introduced in the context of the Helmholtz and cubically nonlinear Poisson problems. We thus restrict ourselves to a few brief comments.

First, in forming the segments R µ j ,τ , 1 ≤ j ≤ J, we do not need to (and could not. . . ) exhaustively verify that for all µ ∈ R µ j ,τ , |µ -µ j | + ρ u N (µ)u N (µ j ) L 4 (Ω) ≤ τ β(µ j ). Rather, we can make plausible continuity assumptions to construct these intervals, and then verify this condition, a posteriori , online. Second, the computationally most intensive online calculation (for large N ) is precisely this u N (µ) -u N (µ j ) L 4 (Ω) evaluation; however, by invoking the symmetry summation techniques developed in Section 3.3.3, we can reduce the relevant operation count to [START_REF] Almroth | Automatic choice of global shape functions in structural analysis[END_REF] 24 N 4 -typically not dominant for the small N realized by our adaptive sampling process. Third, for Burgers equation in R 1 , our reduced-basis approach is not competitive (even as regards marginal cost) with standard techniques, that is, direct computation of s N (µ). However, our complexity estimates also apply to incompressible Navier-Stokes in R 2,3 , in which case we effect very considerable savings relative to finite element calculation of s N (µ).

Numerical Results

Our model problem is given by ( 99)-( 101) and (110); we need only specify f (v) = To begin, we present in Figure 5 β(µ) and β(µ) as a function of µ; we also indicate, on the log(µ)-axis, the segments R µ j ,τ =3/4 , 1 ≤ j ≤ J = 55. There is clearly some deterioration in the length of our segments as µ j decreases -it would appear that J increases more rapidly than ln(µ max /µ min ). For problems with P > 1, this growth would probably not be tolerable. It is possible that deflation techniques -similiar to those introduced in the context of the Helmholtz problem in Section 2.4.3 -could considerably increase the effective inf-sup parameter, and hence considerably decrease J.

The first test case we consider is µ = 1.0. We present in Table 6

and Υ N (µ)/ u(µ) Y as a function of N . We observe that the reduced-basis approximation converges very rapidly; that at least in this particular case, the "good" choice, (120), obtains -e(µ) Y ≤ ∆ N (µ), ∀ N ∈ N; that the effectivities are, as desired, quite close to unity; and that Υ N (µ) is (constant and) very large. From the latter we could very plausibly select (120) over (121) even if -as will be the case in practicewe do not have access to the true error, e(µ) Y .

We now turn to µ = 0.01. We present in Table 7 Table 7 Error, error bounds, and effectivity as a function of N for the Burgers problem with µ = .01.

of Paris VI and Professor Jaime Peraire of MIT. We would also like thank Mr. Nguyen Ngoc Cuong of the National University of Singapore for helpful discussions. This work was supported by DARPA and AFOSR under Grant F49620-01-1-0458 and by the Singapore-MIT Alliance.