- Email: [email protected]
Optimal error analysis of spectral methods with emphasis on non-constant coefficients and deformed geometries
COMPD’IER METHODS NORM-HOLLAND
IN APPLIED
MECHANICS
AND ENGINEERING
80 (1990) 91-115
OPTIMAL ERROR ANALYSIS OF SPECTRAL METHODS WITH EMPHASIS ON NON-CONSTANT COEFFICIENTS AND DEFORMED GEOMETRIES Yvon MADAY ~b~ra~oire C. N. R. S. d’ Arouse ~~~r~q~ de I’ Uni~e~i~~Pierre et Marie Curie, Place Jussieu, Par&, France and Departmend of Mechanical Engineering, M.I. T., Cambridge, MA, U.S. A.
Einar M. R@NQUIST Department of Mec~~ca~ Engineering, M. I. T., ~a~ridge~ Received 26 June 1989 Revised manuscript received December
MA, U.S. A.
1989
In this paper we present the mimer&I anafysis of spectral methods when non-constant coefficients appear in the equation, either due to the original statement of the equations or to take into account the deformed geometry. A particular attention is devoted to the optimality of the discretization even for low values of the discretization parameter. The effect of some ‘overintegration’ is also addressed, in order to possibly improve the accuracy of the discretization.
1. Introduction Spectral methods were introduced around 20 years ago by S.A. Orszag in order to provide a high order accuracy for the numerical simulation of partial differential equations. The original idea was to use truncated Fourier series to approximate the (smooth) solution when the problem was supposed to be provided with periodic boundary conditions. In order to tackle problems with more general boundary con~tions (Dirichlet or Neumann type), the set of (algebraic) polynomials replaced the set of truncated series, but the characterization of the unique discrete function that would provide the numerical solution was still achieved following the original strategy. This enters in the context of the collocation method where the numerical solution is chosen so as to satisfy the original partial differential equation at some suitably chosen collocation points. Of course the choice of the set of collocation points is of fundamental importance for the accuracy of the method and the first remark is to notice that the number of collocation points must be equal to the dimension of the space of approximation. Otherwise, the problem could, in general, be overspecified. When other authors arrived to the analysis of the problem, ([l-4] among others)- they realized that the collocation method could be interpreted as a variational problem with numerical integration. In fact, it is quite common now to realize that the spectral methods are very close to the finite element method in its p or h.~p version, where convergence is achieved by increasing the order of the ~01~0~~ degree and not by diminishing the size of the elements. It is in this framework that the domain of application of the spectral methods has been generalized. The plain spectral method suffers from being constrained to very simple 004%7825/90/$03.50
@ 1990 - Elsevier Science Publishers B.V. (North-Holland)
92
Y. Maday, E.M.
R)nquist,
Optimal error analysis of spectral methodr
domains: they are limited to be slightly deformed squares (in 2-D) or cubes (in 3-D). The idea to couple domain decomposition techniques to the spectral discretization developed rapidly in order to cope with this initial drawback [5,6]. However, by starting from the strong formulation of the equations, this can only produce the Schwarz algorithm [6,7] or a strong coupling between the elements [5], where the solution (in the case of second order elliptic problems) is searched as a global %’ function that is piecewise polynomial. This results in some drawbacks. First, in the Schwarz algorithm, there is an increase of the work due to the double computation over the overlapping region (recall that this one has to be large enough in order to achieve a good convergence rate of the algorithm [7,8]). Second there is a lack of optimality of the approximation in the strong coupling formulation-both from the numerical analysis point of view and from the algorithmic point of view-a consequence of an overconstrained problem due to the %’ matching. In this context, understanding of the similarity between the collocation method and a variational formulation used with consistent numerical quadrature brings a lot of flexibility. Indeed, the coupling required in the nonoverlapping decomposition of the domain is weaker (only %“) and allows for constructing an optimal method. This has lead to the spectral element method [9-111 and more recently to the mortar element method [12-141. The variational method involves integrals that can be computed (with or without numerical quadratures) separately over each subdomain. The spectral element method uses consistent quadrature formulas and in this respect (but not only in this one) conserves the spirit of former spectral discretization (the other aspects being the use of tensorial basis and tensorial evaluations of the residuals). It has been shown [lo, 151 that the method can still be interpreted as a collocation method within each element whereas a suitable equation is satisfied at the interface of each subdomain. There are still fundamental differences between the spectral and the finite element methods that are important as regards the numerical implementation of the method, however the general philosophy is the same. The point that we want to address in this paper is related to one of the differences between the p or h-p version of the finite element method and the spectral methods, more precisely to one possible drawback of the spectral method. Indeed, as we have said, the spectral method can be interpreted in a finite element framework when a particular numerical integration formula is used. Derived from the collocation method, the numerical integration formula, used to compute the integrals in the variational formulation of this discrete problem, is constrained to be based on a certain number of points. This number is not related to the fact that some quantity must be well (or exactly) computed, but has to be equal to the dimension of the discrete space. Using a vocabulary that is standard in the finite element context, the points of the numerical quadrature formula have to be unisolvent with respect to the discrete space. This is a restriction with respect to the fact that in finite element methods, people use more general integration formulas based on accuracy considerations. In this sense, there is often an over-integration in the finite element statement. The current spectral methods use these ‘consistent’ (with respect to the numerical discretization) quadrature rules with a priori no lack of accuracy. As shown in the appendix, the possibility to use an overintegration in the context of spectral methods is however possible, but is not of importance in most interesting cases as demonstrated in this paper. The situations where overintegration could be of importance are those where non-constant coefficients are involved, either because they are present in the original formulation of the equations or because they are introduced when a deformed geometry is considered and a
Y. Maday, E.M. Rbnquist,
Optimal error analysis of spectral methods
93
mapping to a domain of reference is used. Indeed, it is in these cases that the numerical integration based on N + 1 Gauss-Lobatto points is not accurate enough. We analyze in details these two situations and point out the cases where overintegration is required. The analysis is illustrated with many numerical experiments. This allows to strengthen our theoretical analysis since the exact value of some constants could be of importance and are not exactly evaluated by the theory. The analysis of the effect of overintegration has never been addressed in the spectral context. Previous analysis of Legendre spectral approximations of problems with non-constant coefficients [lo] or with deformed geometries [7] could in this sense be misleading as the results presented there were proving a spectral type convergence (i.e., faster than any algebraic rate if the solution is analytic) but were not optimal. The basic ingredient to get optimality is a new result of [16] and deals with the fact that the interpolation at Gauss-type points is optimal. The paper is organized as follows. In Section 2, we consider the case of original non-constant coefficients equations and we analyze in which cases the overintegration may improve the accuracy of the approximation to allow for optimal results. An appendix presents some numerical considerations on the implementation that shows that the tensor product evaluations can still be preserved even in the case of overintegrations. In Section 3, we analyze the case where the non-constant coefficients are induced by the treatment of a deformed geometry. We prove that in this situation the overintegration is unnecessary to provide the optimality of the approximation. In Sections 2 and 3 we also present some remarks related to the impact of our analysis as regards the design of the best schemes in the case of deformed geometries. 2. Analysis of spectral approximation with non-constant coefficients 2.1. The one-dimensional case Let us consider the following problem: -(pu,),
+ qu =f
Find u E Hi(-1,
in (-1,l).
l), such that
(1)
This problem is supposed to be well posed in Ht(-1, 1) in the sense that the bilinear form a, defined over (Ht(-1, 1))’ as follows: 1
1
a(u, u) =
I
-1
pu,u,
dx +
I
-1
quu
dx ,
(2)
is elliptic and continuous over (Hi(-1, l))*. In fact, for sake of simplicity, we shall assume that there exist three constants pl, p2 and q2, such that O
In addition, approximate
VxE(-l,l), VxE(-1,l).
(3)
(4)
the forcing function f is supposed to be at least in L*(- 1,l). We want to this problem by a spectral element method. Our first work will be to design the
94
Y. Maday, E.M.
Rlnquist,
Optimal error analysis of spectral methodr
discrete space that will approximate Hi(-1, 1). Given a fixed integer K, we shall first consider a partition of (-1,l) in K subintervals Zk, where Zk = (a,, %+*) 7
(5)
and the ak are (K + 1) points in (-1,l)
such that
-1 = a, < a, < - - - < aK_l < aK = I
(and the length ak+l - ak is supposed to be O(K-‘)). The spaces of discretization are imbedded in YN = (9 E WI7
1) I ‘pllk E y&H
7
(6)
where N is some integer and lPN(A) denotes the set of all polynomials of degree SN over A. This space will discretize L “(- 1,l). In order to discretize Hi(- 1, 1), let us introduce
x, = r, n HA(-1, 1)
(7)
and cp vanishes at +l). (which means that qr,(a,+, ) = qIk+l(ak+l) The second step in the discretization process is to define the discrete problem. As indicated in the introduction, the point of departure is the following continuous variational formulation of the problem (1): Find u E Ht(- 1, l), such that VvEHA(-l,l),
a(u,v)=(f,v)
(8)
where the notation (e, .) stands for the L2(-1, l)-scalar product. the discrete problem is: Find U, E X,,,, such that
Given an integer M, M 3 N,
The definition of aM and (e, s)~ from a and (s, *) will use a composed Gauss-Lobatto quadrature formula. More precisely, this one is based on the data, over each segment Zk of a set of M + 1 points 53
and a set of associated weights c CD= 2 GL
z
p:
that are such that
@(&p:=I,@(x)dx
VC?EY,,_,
(11)
.
k=O i=O
The definition of aM and (e, *)M is now as follows: %(!W
4
= c PN*u, GL
wcp,
+ cGL W#N @) E
y;
*
W,,
UN)E Xi
3
(12) (13)
Y. Maday, E.M. R@nquist, Optimal error analysk of spectral methods
95
REMARK 2.1. The standard spectral element method utilizes a consistent Gauss-Lobatto quadrature rule, i.e. M = N. The case M > N corresponds to an ove~ntegration with respect to the standard spectral element method. Such an overintegration might be of interest in the case of non-constant coefficients. The analysis that follows indicates the situations for which this is the case. The interesting feature about (a, .)M is that it defines over Y, a discrete scalar product that is uniformly equivalent to the L2(-1, 1) one. Indeed, it has been proven in [4] that
The following Lemma is then completely
standard.
LEMMA 2.1. There exists one and only one solution u, to problem (9), and there exists a co~tant C, such that
PROOF.The
Lemma follows easily from the constatation that a, is uniformly and elliptic over X,, i.e., there exist two positive constants LYand y such that a
IIh4i+1,1~Cad+, 4 6 YI141&-1,1j7
and can be deduced
from (3)) (4) and (14).
continuous
(15)
q
Our further analysis is devoted to the derivation of optimal error bounds for the numerical solution. The following Lemma is also standard and follows from (15). LEMMA 2.2. There exists a ~u~t~~t C independent of N, such that
The following Lemma proves that the last term in (16) is of the same order as the best L2(-1, 1)tit off by polynomials of degree 2M - N. Let us first denote by 92:(-l, 1) the space Q-1,1)
= {$ E L2(-1,
and similarly, for any CL,by Hs(-1,
Then we have the following lemma.
1) I Wk E ~“o(f,,): $,I, = $k> , 1)
the space
(17)
Y. Maday, E.M.
96
Rbnquist,
Optimal error analysis of spectral methods
LEMMA 2.3. Let p be a real number >i.
For any 4 E Hc(-1,
l), we have (19)
(20) PROOF. Let us first consider the case where M = N. Let us introduce, for any integer L, the L2(-1, 1)-projection $L = II,(+) of + over P,(-1,l). Then we have
It is an easy matter to note that
(22)
I<+,4 - WA0 %)I Gw - ~~lIL2~-l.~~llVNllL2~-1,I) ’ Besides, let us denote by ,%Nthe operator such that $,I,+$)=+(&;)
Vi=O,.
of interpolation
. . ,N Vk=l,..
from the set %i(-1,
. , K.
1) onto
Y,
(23)
which, thanks to (14), gives
It remains to estimate the middle term in (21). Let us first assume that K = 1, corresponding to a plain spectral method with no domain decomposition. The term in question can then be written as follows: W/7 UN) - (V&Y%A = ($h&)[(LN,
LN) - (LN,
LJNI 7
where, for any i E N, ,. +j =
($7
tvN7
Lj)
(Lj, Li) ’
+ =
Lj)
(Lj, Lj)
(25)
Y. Maday, E.M. R@nquist, Optimal error analysis of spectral methods
97
are the coefficients of $ ( and thus of JIN) and uN, respectively, in the basis of the Legendre It is then an easy matter to note, from (14) that polynomials (Lj)jERI’
K&v, u,) - MN9 %),I Sz ~l~~~~lll~~ll2L~~-~.~)
442Nl ll~~ll~2~-~,~~l~~l IIb4L~e1.1)* Also, it follows easily that
So we deduce
The multidomain case where K > 1 is treated similarly over each subdomain Ik and the error is the same (the notations are in this case just more comp~cated). Let us recall that we have for any pa0 and any L .
We derive by sugaring
(21), (22), (24) and (26) that
‘Ibis completes Lemma 2.3 in the case M = N, after we recall that for any p > $ and any L we have [16]
II40- %4lL2(-1,1)6 CL-pIIdfP(-I,I)
WPE ~Y--~, 1) ,
(28)
so that the contribution of N-‘II$ - elN_*IILzt_l,l) is much smaller than the other terms. Let US now consider the ca:e where M > N and let us use the same strategy as before, with $zM_Nin place of +!+, where JI,,_, is some good approximation of JI of degree s2M - N that will not be necessarily its L2(-1, 1)projection and will be precised in the following. In order to prove Lemma 2.3, we must now estimate
The first term and the second term are treated exactly in the same way as previously and we
98
Y. Maday, E.M.
R)nquist,
Optimal error analysis of spectral methods
obtain
From the basic properties
of the Legendre
polynomials
we recall that
and that (n + l)&+,(X)
= (2n + l)&(X)
- nL,_,(x) .
A reiterate use of these two formulas yields
tL2M-N,LN)M=
4M-2N-1 244
_
N+l 2N
N
(&u-,-P
L,+&
= **
Now using (14), we derive that (L,,_,,
~3llLMllt+I.
L,),
1)
~3ll~ZM-NIIL~~-1.I)ll~NllL~~-l,l~
so that, as before,
J(&M-N, u,)
-
<&,_,,u,),I s CN-‘IIICI,,-, - &M-N-I IIG-IA.
(32)
The treatment of the last term in (29) has to be different. Indeed the same procedure as before would result in a bound of this term by the quantity ))I(,- $~+l)~2~__~,~)and thus a loss with respect to the optimality we could expect from the other terms. Besides, the numerical experiment (see Fig. 2.1) clearly shows a behaviour like II+ - $zM_NIIL2C_1,1).In order to prove this type of convergence, let us go back to the case K = 1 and note that (28) yields, for any CL> 1,
114,(PllL2(-1.1) 6 lI~lIL~(-*.1) + M-%lIH~(-1,1) v9 E fm-l~ 1).
(33)
From (14) it now follows that ($2,~,,
%)M - ($7 UN)M
= (G&j-N - $7 GM ~3lI$M(~Wv
- ~)lIL2~-1,1)ll~NllL2~-1.1,
+ M-‘lh,-rv - ~‘~H’(-~,I)I~~uN(IL~~-~,~) . <3[ll ICTZM-N - ~lIL~~-I,l, At this point the correct &,_, will be chosen. We see that it has to be a good approximation norm as well. This will be the case of the projection of $ with respect to of + in the H"(-1, l)-
Y. Maday, E.M. Rgnquist, Optimal error analysis of spectral methods
::I0 IL&____ 0
2
4
6
8
10
12
14
16
18
99
20
N
Fig. 2.1 A plot of various numerical results relevant to the approximation of the one-dimensional problem (l), in order to analyze the effect of overintegration. We have here considered the case where q = 0 and different values of p and f are simulated and the effect of overintegration is addressed. The relation between p and f is always chosen (of course, we do not use the knowledge of the exact in such a way that the exact solution is u = sin(fIrx)e-” solution in the simulations!!). The discretization is based on K = 4 equal spectral elements for different values of the polynomial degree, N. First of all, since this is the aim of the numerical scheme, we have computed the best fit of u in Ht(O, 1) by piecewise polynomials of degree d N (plot A); the plot proves a convergence better than exponential that is consistent with the fact that the solution is entire. Next we consider the case p = 1 with no overintegration (plot 0) and we can check that the approximation is optimal and very close to the best fit. We then choose p = 1 l(x + a) for various values of a. The case where a = 0.25, treated with no overintegration, is represented on the plot @ and we can remark that the approximation is still very good. The convergence rate however does not follow the best fit of u, but the best fit of p (consistently with Corollary 2.5) that is a straight line in this scale since p has a pole outside (0,l). Let us now use a more singular p, corresponding to a = 0.05. The plot 0 represents the results when no overintegration is used. It is very close to three other plots representing the results when overintegration (M = 3N/2) is used to compute only the right-hand side (the contribution off) (0) or only the left-hand side (corresponding to p) (plot A), while the third one represents the best fit of p (plot l) in the L*(O, 1)-norm by the discrete functions. In order to recover optimal results, overintegration must be used on both the left and the right-hand side of the equation (plot 0).
the W(-
1, 1)-norm over YN. It is proven in [17] that this element
satisfies
(34) and from (29), (30) and (32) the proof of Lemma 2.3 is complete M>N.
also in the case where
100
Y. Maday, E.M.
The second inequality vanish. 0
Ranquist,
Optimal errcv analysis
of spectral
methodv
(20) is quite trivial since the middle terms in (21) and (29) simply
We are now in a position to state the following Theorem. TUTORED 2.4. Let us suppose that the suluti~n u to problem (1) belongs to H”,( - 1,l) , that p belongs to Hk(-1, l), q belongs to Hcf-l,l) and that f belongs to Hi(-1, l), and in addition that the four real numbers CT- 1, CL,Yand p are larger than 1. Then the following error estimate holds :-
PROOF. It is a direct application
of Lemma 2.3 to derive
Let us choose w, as being the best fit of u in the H’(- 1, l)-norm.
It follows from (17) that
II@ - %?llH’(-l,I)si c~*-“ll~ll~~~-~,~~* If we now let G, be any polynomial
(37)
in X,, it follows from (3) and (4) that
so that sup
a@, UN)- %A%
qv)
uNEXh
The interest of this decomposition is that we can now choose GM, still close to u but such that the higher norm (i.e. H”,(-1, 1)) is unifo~ly bounded. This will now be the case if we choose GN as being the projection of u with respect to the H-(-l, 1)-norm over X,. It is proven in [17] that this element satisfies
101
Y. Maday, E.M. R)nquist, Optimal error analysis ofspectral methods
From the definition of a and uM we have
Using Lemma 2.3, first with 4 = pGNx and then with $ = qfiN, we derive that
l4L UN)- %t~;N,GI S C[(2M - N)-
min'u'o-l'~~P~Nxll~~in(v.o-l)(_~,~)
(2M - N)-
+
6 C[(2M
mi”(p,u)ll ~~,lln~“cr.~~(-~,~)lll~,IlH’(-1,I)
- N)- mi”(“.rr-l)llpUxIIHmi”c”,~-l,c_l,l,
+(2M - N)-
+ C[(2M
mi”(p9u)II 4~II~~~~r.~~(-~,l)lll~NIIH~~-~,1)
- N)- min(“‘ul-l)))p(U-
+(2M - N)-
~N)xlJHgi"(Y.~-')(_l,l)
mi”(p~rr)llq(~ - ~~~Il~~“c~.~~(-1,1)111~~ll~‘(-l,l) *
Let us recall that min(v, p, u - 1) > 4 so that H~in(vyo-l)(-l,l) algebras. This proves that
and I~IKrni”(~*~)(-l,1) are
146 NT4 - QMGL 4 SC[(2M
- N) -min(v,o-l)llpUxIla~“(u.~-l)(_l,l)
+(2M - N)-
+ C[(2M
- N)- min~“~u-l~llpll~~(_~,~~ll~~ - ~N)x(lHpin(u.o-l)(_l,l)
+(2M - N)which together
mi”~a~o~ll~~ll~~~~~~~~~~-~,~~lll~~II~~~-~,~~ min~~~u~ll~llHf;(-l,l)II~~ - ~,)o,~~~r.~~(-,,,)lll~NIIH~~-*,1) f
with (39) and the fact that M 3 N (so that 2M - N 2 N) gives
l4w”w %I - %fNv~ QJI
- N) -mi”(u*o-l)IJPUxJIHKmi”(“.a-l,(_l,l) -
N)-mi”‘~~o’~~q~~~H~.cr~u,~_l,l~]~~~N~~H~~_I,I~
+~~~l-(rll~ll~;;(-*,~)IIUll~~(-*,l) +~-“l14118~~-1,1~ll~IlHg(-1,1)lll~NllH~~-~,*~ * Theorem
2.4 now follows from the fact that H’(- 1,l)
is an algebra when r > 1.
Cl
Y. Maday, E.M.
102
R#nquist,
Optimal error analysis of spectral methods
COROLLARY 2.5. Let us suppose that the hypotheses of Theorem 2.4 are fullfilled and in addition let us suppose that M = N. Then the following error estimate holds:
REMARK
2.2. Interpretation of the results. In the practical computations, the more common situation is the one where the regularity of the solution is limited by the regularity of the non-constant data. More precisely, and especially if we consider the extension to the 2-D case that will be treated in the forthcoming section, we shall have (T - 1 s min(v, F, p). In this case, the estimates given in Corollary 2.5 prove that the optimality of the scheme is achieved with no use of over-integration. This statement is confirmed by the numerical experiments shown in Fig. 2.1. In the case where the regularity of the solution is better than the regularity of the non-constant data, we expect, from Theorem 2.4 that some amount of overintegration is necessary to recover the optimality of the discretization. The numerical experiments of Fig. 2.1 also confirm this statement with exactly the rate theoretically proven. It is important to note that the possibility to include some amount of overintegration is only true in the case where the variational formulation is the starting point of the discretization. The pure collocation method is incompatible with the concept of overintegration.
2.2. The two-dimensional case Let us suppose here that the problem is set on a square R = (-1, l)* and that no domain decomposition is used. The extension of our results to the spectral element case is very simple and it is only for sake of limitation of the notations that we restrict our analysis to this particular case. Let us consider now the following test problem: Find u E Hi(O), such that
V-A-Vu=f,
(42)
where f is a given force (supposed to be in L*(0)), A = (ajj) is a 2 X 2 matrix with non-constant coefficients and we assume, for sake of simplicity, that it is bounded, symmetric, positive and non-degenerate, i .e . , such that there exists 2 real numbers p1 > 0 and p2 3 0 such that AV-Vsp,(V;+V;) laij( Sp2
The variational
VV=(V,,V,)E[W*,
Vi, j=
formulation
a(4 u) = ((f,
v))
(43)
1,2.
of this problem:
(44 Find u E HA(O), such that
(45)
Vu E Hi(R)
will be used as a starting point to define the discrete problem.
Here
Y. Maday, E.M. R@nquist, Optimal error analysis of spectral methodr
a(u, u) = ((AVU, Vu))
103
V(U, u) E (H’(O))* .
The space of discretization, X,, consists of all polynomials (~~(~)~ that vanish at the boundary
of degree SN in each variable
x;, = UqL?) n H;(R). The quadrature formula is based on (M + 1)” points derived by a tensor product of the one-dimensional Gauss-Lobatto formula, and consists of the evaluation CGL,i CGL,2, where the indices 1 and 2 refer to the first and second spatial direction, respectively. The discrete problem is then: Find uN E X, such that
%f(%v~ %v)= (If7 MM where we have introduced
and
(46)
the following discrete scaIar product over PN(fi):
((R JI)), = o?, o;, @ I
WV E xiv>
,
V(cp, $1 f P&W2
%(%v, Qiv)= ((AV% V%)),
(47)
*
It is standard to note that, from (43) and (44), this problem has a unique solution and that in the case where M = N it is equivalent to a collocation problem based on the (N + 1)2 Gauss-Lobatto points. The analysis of the error will be based on Lemma 2.3 and on the follo~ng complement of this Lemma. LEMMA 2.6. Let p be a real number > 1. For any $ E Hc(- 1, l), we have
(48) PROOF. This proof is derived from (29) by taking $2M_N_1instead of I&~_~, which has the effect to cancel the middle difference in the right-hand side of (29). REMARK 2.3. It is an easy matter to note that (48) cannot be improved to yield a bound by the best fit by polynomials of degree 2M - N as it was in Lemma 2.3. We are now in a position to give the following error estimate. THEOREM 2.7. Let us suppose that the solution u to problem (42) belongs to H”(O), that the coefficients (aij) of the matrix A belong to H’(a) and that f belongs to HP(a) and let us suppose in addition that the real numbers p, u - 1, p are larger than 1. Then the following error
Y. Maday, E.M.
104
R@quist,
Optimal error analysis of spectral methods
estimate holds:
(4%
+w4 - WPllfllHP(~)] *
PROOF. The proof is derived from Lemma 2.2 (with obvious changes in the notations) following exactly the same lines as in the previous subsection. The only difference relies on the fact that Lemma 2.6 has to be used in the analysis of the quantity sup
4% 4 - %fh
4
vNExN
Here wN is a approximation bidimensional As in (38) we
sup
polynomial of degree N (and not N - 1) suitably chosen to be an optimal of u in the norm H”(R), say, and in any other lower norm from the equivalent to (39) (see [ 161). can write 44
UN) - 5&h,
uNEXN
IIhAil(*)
UN)
a(w,, 4 - a&h
uNEXN
where the first term on the right-hand
which can be bounded
s sup
II%hi~(*)
vd
+ clb - %&f’(R)7
side is a sum of four terms of the type
from Lemma 2.6 as
The reason why (20) is not sufficient to analyze this term is that in the second direction, of degree sN. It follows easily from the analysis in the one-dimensional case that we can prove
sup
(f’ ‘A’)- (f7vN)MC(2M Q
- N)-PllflI
HP(~)
’
uNXis
(51)
vNEXN
which allows us to derive Theorem
2.7.
Cl
COROLLARY 2.8. Let us suppose that the hypotheses of Theorem 2.7 are fulfilled and let us suppose in addition that M = N. Then the following error estimate holds:
IIu- u,IIHltn) s C(N-~~~min~Y~~~l~p~~ll~ll~~~~~II~ll~~~~~ + IIf IIwyn)) .
Y. Maday, E.M. R@nquist, Optimal error analysis of spectral methods
105
REMARK 2.4. Interpretation of the results. As pointed out in the previous subsection, there is no need to over-integrate in the general case where the solution has the same regularity as the non-constant coefficients. In two dimensions this conclusion still holds, however more weakly. Indeed, we have only been able to prove that the plain spectral or spectral element method provides the same accuracy as the best fit by the polynomials of degree one less. An overintegration with just one more point (M = N + 1) is sufficient to recover the optimal accuracy. Although the results are asymptotically the same, this remark could suggest the use of an overintegration. However, the associated increase in work makes it preferable to use a standard spectral element method with polynomial degree of N + 1 (instead of N) since this is less expensive than an over-integration, with M = N + 1 (see Appendix). It is only in the case where the spectral method would require much higher values of N that over-integration, coupled with the preservation of the tensor product formula, could be of interest. Finally, note that the numerical experiments again confirm the theoretical statements with an even less important difference between the cases M = N + 1 and M = N; see Fig. 2.2.
10-l
.
c
* lo-*
.
8
c
2
.
0
1o-3
0
E
. IO“ t LL
& 10-s p kz lo4
c
lo-’
p
lo+
c
lo+
r
.
Cl
8
10-I
-I
r
1
’
’
2
’
’
3
’
’
4
8
’
’
5
’
6
8
’
7
’
’
8
’
’ 9
’
1
10
N
Fig. 2.2. This figure represents the equivalent of the previous case in a two-dimensional domain. The equation is here (42) with A = pId and different values of p and f, computed so that the solution is always u = sin(nx)y(leye’). The domain of computation 0 is (0, 1)2 and decomposed in four squares. The plot 0 represents the case p = 1 with M = N, M = N + 1, M = 3N (there is no difference between the two last experiments, and the first case is a bit less precise but cannot really be distinguished from the two last cases). The plot 0 represents the experiment with p = ex+q that is, a very smooth function, and a treatment with no overintegration gives the same accuracy as the previous plot (that certainly also corresponds to the best fit). Now comes the treatment of the case p = 1 + cos(27rx + 37ry), which is a smooth function but is worse than u. The plot A corresponds to no over-integration and we clearly see that some degree of overintegration is necessary to recover the optimal accuracy, as shown in the case where M = N + 1 (plot 0) and M = 3N (plot A).
106
Y. Maday, E.M.
[email protected], Optimal error analysis of spectral methods
3. Analysis of spectral approximation in deformed geometries Our purpose is now to extend the previous analysis to the case where the non-constant coefficients are induced from the treatment of a non-rectangular geometry. Since the point here is not to analyze the behaviour of the approximation when original non-constant coefficients are present in the equation, we shall restrict our analysis to the simple problem of a Poisson equation in a deformed domain 0 which reads as follows: Find u E HA(a), such that
II
n
VuVv dx dy =
fodxdy
VVEH@).
(52)
In order to explain an important issue related to the error analysis, let us first consider the simple case where the geometry is rectilinear but not rectangular. 3.1. Rectilinear geometries Let us suppose that the domain 0 is the trapezoid with vertices A = (0,0), B = (4,0), C = (0,l) and D = (4 - 4/a, l), where a is a real number >l. It is an easy matter to check that the transformation F : (r, s) I+ (x, y) with 4 a
(
x=r4--s
> ,
(53)
y=s,
is a one-to-one mapping from the square fi = (0, 1)2 onto the domain a. We can now define a natural correspondence between functions 9 defined over 0 and functions C$defined over fi as follows:
G(r, 4 = cp(F(r,4) . This will allow to set the discrete problem. Indeed as we have indicated in the introduction, the spectral method has to be used on square domains. The other domains have to be mapped on such a reference square in order to allow for a spectral computation (the spectral element method allows to deal with more complex geometries, however, the elements must still be mapped to such a reference square). The variational formulation transferred to the domain fi is the following: Find li E Hi(b), such that
+$
2
o,(r,s)]drds=//fi
@J(r,s)drds
Vi,EHi(fi).
The coefficients J and oi, i = 1, 2, 3, are respectively the Jacobian of the transformation some smooth geometric factor given by
(54) F and
Y. Maday,E.M. R&quid,
J=Z(a-S),
US=(4-;)2.
(55)
and consists of the following:
Find
U2=(4-f)(2),
The discretization of the problem Li, E lPN(a) n HA(a), such that
107
Optimal error analysis of spectral methods
can now proceed
vo E P,(A)
)
(56)
where ((s, .))M has been defined in (47) and a,(ii,,
ti,) =
(( -
-
w2
-J
aii,-ai3, ar ' as M
The coefficients in the right-hand side of (56) are extremely regular. We can deduce from the previous analysis that no overintegration is required to compute accurately this term in order to achieve optimal approximation. On the contrary, the coefficients that arise in aM on the left-hand side, although analytic, may give rise to ,a problem since the factor q/J has a pole where J vanishes, i.e. at s = a. This pole is not in 0 but can be very close to it. The similarity of this problem with the one treated in the second section (Fig. 2.2) seems to imply that, when a is sufficiently close to 1, an overintegration is required to compute correctly aM in order to get an optimal approximation. However, surprisingly, the experiments performed in Fig. 3.1 show that no over-integration is required, even with those values of a that in the onedimensional case induced a large error associated with the consistency error (50). The interpretation of this fact can be easily understood, a posteriori, as the pole of l/J is compensated by the fact that the other part of the integrand vanishes and makes the term
more regular than it appears. In order to be more precise, let us state the following lemma. LEMMA 3.1. Let us assume that the solution u to problem real number u > 2. Then we have
(52) belongs to H”(a)
for some
(57)
108
Y. Maday, E.M.
Renquist,
Optimal error analysis of spectral methods
B 0
t lo-‘O e
lo-”
a 0
F
1o-L2
j ’
2
’ 3
I
’ 4
’
’ 5
’
’ 6
’
’ 7
’
’ 6
’
’ 9
’
’ 10
’
’ 11
’
12
N
Fig. 3.1. On this plot we prove that the singularity of the function l/Jacobian is not relevant for the need of overintegration. We consider problem (52) for two different domain decompositions A and B. They are
A
B
The three plots represent the discretization error in the semi-norm with no overintegration. The plot 0 corresponds to the solution u = sin(2y) for the two cases A and B. They both coincide since the deformation is in the x-direction and the solution only depends on y. Next, we repeat the experiments with u = sin( )x and we arrive to the (a priori) surprising result: plot a for A and q for B. We refer to Section 3.1 for the explanation. The fact that B is even better in this case is due to the fact that the discretization on the two domains is more equilibrated and that the left domain in case A is too long.
f
PROOF.
Let us go back to the original equation
(52) and remark that
and that a simple change of variable formula yields
Y. Maday, E.M. R)nqutit, Optimal error analysis of spectral methods
109
so that Q(k UN) - %.#i),,
u,) = A, + B, + A, + B,,, + B,,, + A, + B,,, + B,,, ,
(58)
where
A,=
G n + (r, s) 2
II
(6 s)(4 - z s)drds-(((4-zs):,?)),,
B~,1=(((~)~,~)),-(((4)~,~)),, B~,~=(((4)~,~))M-(((g)~,~))M. In the decomposition (58), the terms A_ are small from the accuracy of the quadrature formula; the terms B_ are small since GN is a polynomial that is close ,fo ii. More precisely, let us denote by C the supremum of 11/JI, I(1 - ol) lJ, (021Jl, 1t.03/Jlon 0. It is an easy matter to derive that, for any index cr, where (Y= 1, (2,1), (2,2), (3, l), (3,2), we have
where t in the derivation
stands for s or r. Using (33), with p = 1, we derive
(59)
110
Y. Maday, E.M. IQnquist,
Let us denote now by C?the supremum 2.2 that
Optimal error analysis of spectral methods
of 1, /(4 - 4sfa)l, 14rlal on fi. We derive as in Section
and the lemma follows from (SS), (59) and (60).
•!
In order to analyse the error in the approximation of problem (52) with (56), we first note that (16) (with obvious changes of notations) also holds in this case. This leads to the following theorem. THEOREM 3.2. Let us suppose that the solution u of problem (52) belongs to H”(0) and that f belongs to HP(R) and let us suppose in addition that the real numbers (+ - 1, p are larger than 1. Then the following inequality holds: II; - unrllul(~)s C~l-“llu~ln~~n~ + C(2M - N-- l)l-%IIn~~~~ + G(2M - WP
IIfllHPcn,*
(61)
PROOF. This is an easy consequence of (48) and of Lemma 3.1, if we choose tii), as being the best fit of ii in H”(b) by elements of PN(&). 3.2. General ~forma~ons Let us suppose in this subsection that the domain 0 is the image of & by some mapping F : (r, s) I-+ (x = F,(r, s), y = _F2(r,s)). We no longer assume tha! F is bilinear, nevertheless we assume it has some regularity, more precisely, that it is in Hail) for some p > 1, and that its Jacobian J is larger than some constant p1 > 0. The problem (52) is transferred onto fi as before and reads:
Find ti E HA(h), such that (54) is satisfied, where the geometric
factors oi, i = 1,2,3
are defined as follows:
The discrete problem is now: Find tiN E ~~(~)
n Hi(h)
Using here the same technique Lemma 3.3.
, such that (56) is satisfied.
as in the proof of Lemma 3.1 and Theorem
2.4, we obtain
LEMMA 3.3. Let us suppose that the solution u to problem (52) belongs to H”(a), for some real number cr > 2, that the deformation mapping F belongs to H’(h) for some real number
Y. Maday, E.M. Renquist, Optimal error analysk of spectral methods
111
u > 1. Then the following inequality holds :
REMARK 3.2. It is important, at this stage, to note that the regularity of G and the regularity of am where t is either s or t is bounded by the regularity of both u and F, which means that
is the best result we can expect. Note also (but this is natural after the analysis of Section 3.1) that it is only the regularity of F that is important and not at all the regularity of its inverse (as one could previously fear from the factor 1 lJ). From Lemma 3.3 and (65), we derive in a (now) standard manner the following theorem. THEOREM 3.4. Let us suppose that the solution u belongs to problem (52), belongs to H”(O), for some real number u > 2, that the deformation mapping F belongs to HF(@ for some real number p > 1 and that the forcing function f belongs to HP(a) for some real number p > 1. Then the following inequality holds:
+(2M -
N)-min(C’p)I(FII,~(n)lI f (It+n(r.p+a)} .
W)
REMARK 3.3. In the light of the previous theorem, it appears that there is no need to overintegrate in the case where the deformation is responsible for non-constant coefficients in the formulation of the discrete problem, even if the geometry is very distorted. This is readily seen by comparing the results obtained in the theorem by choosing M >> N and M = N (the gain of just one degree is not, as explained in Remark 2.4, sufficient to start the overintegration machinery). We can explain it a posteriori now as overintegration could allow to improve part of the consistency error term (57) in the case where the geometry mapping F is less regular then the solution u itself (contrary to the case in Section 2 this can happen, especially if domain decomposition is used). However, as pointed out in the previous remark, loss of regularity of F induces also a loss of regularity of & which is the function that is of importance since it is this one that is approximated by polynomials. REMARK 3.4. We have not considered here the case where the geometry mapping F itself is approximated by isoparametric polynomials. This effect can easily be analysed by following the same lines as previously and including the standard arguments adapted to the finite element discretization [ 181.
112
Y. Maday, E.M. loo
Rtinquist,
Optimal error analysis of spectral methods
k
I .
10-l
c
$2 A
sa . .
10-2 c
. IO”
A
0 0
c
A r
(L
w”
.
0
A
0
A
.
10-7
. A0
1o-5 c lo4
.
0 A
8
3 .
A A0
.
104
.
s
c
.
E loa
c
.
10-O c . 1
lo-
10' lo2 N Fig. 3.2. This plot is to prove that the overintegration can even have a negative effect, Again we consider (52) for p = 1, with a decomposition of a= [0,7] x [0, l] into two subdomains: 100
The interface is described by the equation x = $ - 3 ]y13. This equation is not regular and its interpolant is more rough when the degree is higher. Since this interpolant is the only geometric factor for the discrete problem we can understand the plots A and A that correspond to M = N and M = 2N, respectively, for the approximation of the solution u = sin( 4x). As before, if the deformation and the solution are in different directions (this is a very particular case!!!), for example u = sin 2y, the behaviour of the approximation is particular also and this is the only case where overintegration can help (in the case of deformed geometries) as is illustrated in the plots 0 and 0 corresponding to M = N and M = 2N, respectively.
Acknowledgment The question of the importance of overintegration has been raised during a meeting with I. BabuSka. The many discussions with A.T. Patera have been very enlightening. We thank them both for their help. This research has been supported in part by NSF under ASC-8806925, ONR and DARPA under contract NOOO14-85-K0208. Part of this work has been done while the first author was in residence at ICASE, NASA Langley Research Center, Hampton, VA, U.S.A.
Y. Maday, E.M. R)nquist, Optimal error analysis of spectral methods 10-Z
,
I
113
I
!a 1o-3
r
z
OA
A
lo+ :
A
-j
A
0
A
10-5 c
A
A
0
A
8
oz lo6
7
25
0 1o-'O loo
I
1
8
I-I I
h
lo*
10' N Fig. 3.3. A plot of the dkcretization error in the semi-norm 0 = [0,2] x [0, l] 1s * d ecomposed into two subdomains:
when solving
(52) with p = 0. The domain
I
1
and this time the geometry is tortured in the following way. The subdomain on the left is mapped onto the reference square by the mapping x = r + OS(lrl - 1) r”, while the mapping is affine in the other subdomain. The plot A corresponds to a solution u = sin( 4x) that only depends on x and is treated both with M = N and M = 2N. The plot 0 corresponds again to a particular case as it is related to u = sin 2y in both cases M = N and M = 2N.
Appendix (in collaboration
with A.T. Patera)
In this appendix we consider the computational cost of using overintegration. We show that the matrix-vector products can still be evaluated efficiently using tensor product forms, and that for a fixed discretization (polynomial degree N), the increase in operation count is at least a factor of d in R’. Overintegration is attractive to use only if the convergence rate (error as a function of polynomial degree) improves substantially. In this case we can achieve a fixed accuracy either with a low polynomial degree and overintegration, or with a higher polynomial degree with consistent quadrature. With no loss of generality we consider here the evaluation of the left-hand side of (46),
%M(%Q %J = ((PV%&J), *
Y. Maday, E.M.
114
R)nquist,
Optimal error analysis of spectral methoak
Here p represents a non-constant coefficient which either comes from the original strong formulation (as in (42)), or represents a geometric factor in the case of deformed geometry, or both. The evaluation of (67) represents the computationally most expensive part in an iterative solution of (46). In order to implement (46) we require a basis for our high-order polynomial space X,. The choice of basis does not effect the error estimates, however it greatly effects the conditioning and sparsity of the resulting set of algebraic equations, and is critical for the efficiency of parallel iterative solution procedures. We choose an interpolant basis to represent wN E X,:
p=o q=o
w,qq~)~q(4
(68)
9
where (x, y) E 0 -+ (r, S) E l-1, l[‘. Here the hp(.z) are the one-dimensional Nth-order Lagrangian interpolants through the Gauss-Lobatto Legendre points 5, (hp E P,(]-1, l[), hp( 5p) = $,,,), and wPq is the value of wN at the local node ( 5, , 6,). In addition to (68) we requn-e w, C X,,, to honor the homogeneous Dirichlet boundary conditions. The bases (68) are now inserted into (67), and the test functions are systematically chosen to be unity at one global node and zero at all the other Gauss-Lobatto Legendre points. We then arrive at the following discrete statement:
(69) where we for convenience only have included the derivatives with respect to the first (local) spatial direction. Here pP are the one-dimensional quadrature weights, and D and I are the one-dimensional derivative and interpolation operators defined as
Dpq=
2
(70)
c&J,
zpq =h,G,) Using tensor-product as follows:
(71)
* sum-factorization
techniques
where the expression in the innermost partenthesis order results in the following operation count:
[5] we can now evaluate (69) efficiently
is evaluated
first. An evaluation
in this
MN2+M2N+M2+M2+M2N+MN2=2(M2N+MN2+M2). In the case where M = N (no overintegration), and (69) reduces to the following expression:
the interpolation
operator
becomes IPq = aPp4
Y. Maday, E.M. R)nquist, Optimal errar analysis of spectral methodr
115
which can be efficiently evaluated in 2(N3 + N2) operations. We are now in a position to compare the computational cost with and without overintegration. In the case A4 = N + 1 the operation count is approximately twice the cost of using consistent quadrature (A4 = N) (in three dimensions the cost would increase with a factor of three). If A4 = 2iV, say, the cost to evaluate (67) is four times the cost with A4 = N. The factors which determine whether over-integration is economical are the difference in convergence rate with and without overintegration, and the specified accuracy of the discrete solution. References D. Gottlieb, The stability of pseudospectral Chebyshev methods, Math. Comp. 36 (1981) 107-118. ;;; B. Mercier, Stabilite des methodes spectrales polynamiales application B l’bquation d’advection, RAIRO Anal. NumQ. 16 (1982) 67-100. Y. Maday and A. Quarteroni, Approximation of Burgers equations by pseudospectral methods, RAIRO [31 Anal. Numtr. 16 (1982) 375-404. 141C. Canuto and A. Quarteroni, Approximation results for orthogonal polynomials in Sobolev spaces, Math. Comp. 38 (1982) 67-86. S.A. Orszag, Spectral methods for problems in complex geometries, J. Comput. Phys. 37 (1980) 70-92. ;i; Y. Morchoisne, Resolution des equations de Navier Stokes par une methode spectrale de sous domaines, in: Proc. 3rd Intemat. Conf. Numer. Methods Sci. Engrg. (Gamni, Paris, 1983). [71 B. MCtivet, Resolution spectrales des equations de Navier-Stokes par une methode de sous-domaines courbes, These de Doctorat, Universite Pierre et Marie Curie, 1987. PI C. Canuto and D. Funaro, The Schwarz algorithm for spectral methods, SIAM J. Numer. Anal. 25 (1) (1988) 24-40. 191A.T. Patera, A spectral element method for fluid dynamics: Laminar flow in a channel expension, J. Comput. Phys. 54 (1984) 468-488. [101 Y. Maday and A.T. Patera, Spectral element methods for the incompressible Navier Stokes equations, in: A.K. Noor, J.T. Oden, eds., State of the Art Survey in Computational Mechanics (1989) 71-143. incompressible Pll E.M. Ronquist, Optimal spectral element methods for the unsteady three-dimensional Navier-Stokes equations, Ph.D. Thesis, MIT, 1988. P21 C. Bemardi, Y. Maday and A.T. Patera, A new nonconforming approach to domain decomposition: The mortar element method, in: H. Brezis and J.-L. Lions, eds., Nonlinear Partial Differential Equations and their Applications (Pitman and Wiley, 1989). 1131Y. Maday, C. Mavriplis and A.T. Patera, Nonconforming mortar element methods: Application to spectral discretizations, in: T. Chan, ed., Proceedings of the 7th International Conference on Domain Decomposition Techniques for P.D.E. (SIAM, Philadelphia, 1988). 1141C. Mavriplis, Nonconforming discretizations and a posteriori error estimators for adaptative spectral element techniques, Ph.D. Thesis, MIT, 1989. 1151D. Funaro, A multidomain spectral approximation of elliptic equations, Numer. Methods in P.D.E. 2 (1986) 187-202. Y. Maday, Resultats d’approximations des operateurs d’interpolation aux points de type Gauss, Notes aux WI Comptes rendus de 1’Academie des Sciences de Paris, submitted. 1171C. Bemardi and Y. Maday, Approximation results for spectral methods with domain decomposition, J. Appl. Numer. Math. 6 (1989-90) 33-52. WI P.G. Ciarlet and P.A. Raviart, The combined effect of curved boundaries and numerical integration in isoparametric finite element methods, in: A.K. Aziz, ed., The Mathematical Foundation of the Finite Element Method with Applications to Partial Differential Equations (Academic Press, New York, 1970) 409-474.
IN APPLIED
MECHANICS
AND ENGINEERING
80 (1990) 91-115
OPTIMAL ERROR ANALYSIS OF SPECTRAL METHODS WITH EMPHASIS ON NON-CONSTANT COEFFICIENTS AND DEFORMED GEOMETRIES Yvon MADAY ~b~ra~oire C. N. R. S. d’ Arouse ~~~r~q~ de I’ Uni~e~i~~Pierre et Marie Curie, Place Jussieu, Par&, France and Departmend of Mechanical Engineering, M.I. T., Cambridge, MA, U.S. A.
Einar M. R@NQUIST Department of Mec~~ca~ Engineering, M. I. T., ~a~ridge~ Received 26 June 1989 Revised manuscript received December
MA, U.S. A.
1989
In this paper we present the mimer&I anafysis of spectral methods when non-constant coefficients appear in the equation, either due to the original statement of the equations or to take into account the deformed geometry. A particular attention is devoted to the optimality of the discretization even for low values of the discretization parameter. The effect of some ‘overintegration’ is also addressed, in order to possibly improve the accuracy of the discretization.
1. Introduction Spectral methods were introduced around 20 years ago by S.A. Orszag in order to provide a high order accuracy for the numerical simulation of partial differential equations. The original idea was to use truncated Fourier series to approximate the (smooth) solution when the problem was supposed to be provided with periodic boundary conditions. In order to tackle problems with more general boundary con~tions (Dirichlet or Neumann type), the set of (algebraic) polynomials replaced the set of truncated series, but the characterization of the unique discrete function that would provide the numerical solution was still achieved following the original strategy. This enters in the context of the collocation method where the numerical solution is chosen so as to satisfy the original partial differential equation at some suitably chosen collocation points. Of course the choice of the set of collocation points is of fundamental importance for the accuracy of the method and the first remark is to notice that the number of collocation points must be equal to the dimension of the space of approximation. Otherwise, the problem could, in general, be overspecified. When other authors arrived to the analysis of the problem, ([l-4] among others)- they realized that the collocation method could be interpreted as a variational problem with numerical integration. In fact, it is quite common now to realize that the spectral methods are very close to the finite element method in its p or h.~p version, where convergence is achieved by increasing the order of the ~01~0~~ degree and not by diminishing the size of the elements. It is in this framework that the domain of application of the spectral methods has been generalized. The plain spectral method suffers from being constrained to very simple 004%7825/90/$03.50
@ 1990 - Elsevier Science Publishers B.V. (North-Holland)
92
Y. Maday, E.M.
R)nquist,
Optimal error analysis of spectral methodr
domains: they are limited to be slightly deformed squares (in 2-D) or cubes (in 3-D). The idea to couple domain decomposition techniques to the spectral discretization developed rapidly in order to cope with this initial drawback [5,6]. However, by starting from the strong formulation of the equations, this can only produce the Schwarz algorithm [6,7] or a strong coupling between the elements [5], where the solution (in the case of second order elliptic problems) is searched as a global %’ function that is piecewise polynomial. This results in some drawbacks. First, in the Schwarz algorithm, there is an increase of the work due to the double computation over the overlapping region (recall that this one has to be large enough in order to achieve a good convergence rate of the algorithm [7,8]). Second there is a lack of optimality of the approximation in the strong coupling formulation-both from the numerical analysis point of view and from the algorithmic point of view-a consequence of an overconstrained problem due to the %’ matching. In this context, understanding of the similarity between the collocation method and a variational formulation used with consistent numerical quadrature brings a lot of flexibility. Indeed, the coupling required in the nonoverlapping decomposition of the domain is weaker (only %“) and allows for constructing an optimal method. This has lead to the spectral element method [9-111 and more recently to the mortar element method [12-141. The variational method involves integrals that can be computed (with or without numerical quadratures) separately over each subdomain. The spectral element method uses consistent quadrature formulas and in this respect (but not only in this one) conserves the spirit of former spectral discretization (the other aspects being the use of tensorial basis and tensorial evaluations of the residuals). It has been shown [lo, 151 that the method can still be interpreted as a collocation method within each element whereas a suitable equation is satisfied at the interface of each subdomain. There are still fundamental differences between the spectral and the finite element methods that are important as regards the numerical implementation of the method, however the general philosophy is the same. The point that we want to address in this paper is related to one of the differences between the p or h-p version of the finite element method and the spectral methods, more precisely to one possible drawback of the spectral method. Indeed, as we have said, the spectral method can be interpreted in a finite element framework when a particular numerical integration formula is used. Derived from the collocation method, the numerical integration formula, used to compute the integrals in the variational formulation of this discrete problem, is constrained to be based on a certain number of points. This number is not related to the fact that some quantity must be well (or exactly) computed, but has to be equal to the dimension of the discrete space. Using a vocabulary that is standard in the finite element context, the points of the numerical quadrature formula have to be unisolvent with respect to the discrete space. This is a restriction with respect to the fact that in finite element methods, people use more general integration formulas based on accuracy considerations. In this sense, there is often an over-integration in the finite element statement. The current spectral methods use these ‘consistent’ (with respect to the numerical discretization) quadrature rules with a priori no lack of accuracy. As shown in the appendix, the possibility to use an overintegration in the context of spectral methods is however possible, but is not of importance in most interesting cases as demonstrated in this paper. The situations where overintegration could be of importance are those where non-constant coefficients are involved, either because they are present in the original formulation of the equations or because they are introduced when a deformed geometry is considered and a
Y. Maday, E.M. Rbnquist,
Optimal error analysis of spectral methods
93
mapping to a domain of reference is used. Indeed, it is in these cases that the numerical integration based on N + 1 Gauss-Lobatto points is not accurate enough. We analyze in details these two situations and point out the cases where overintegration is required. The analysis is illustrated with many numerical experiments. This allows to strengthen our theoretical analysis since the exact value of some constants could be of importance and are not exactly evaluated by the theory. The analysis of the effect of overintegration has never been addressed in the spectral context. Previous analysis of Legendre spectral approximations of problems with non-constant coefficients [lo] or with deformed geometries [7] could in this sense be misleading as the results presented there were proving a spectral type convergence (i.e., faster than any algebraic rate if the solution is analytic) but were not optimal. The basic ingredient to get optimality is a new result of [16] and deals with the fact that the interpolation at Gauss-type points is optimal. The paper is organized as follows. In Section 2, we consider the case of original non-constant coefficients equations and we analyze in which cases the overintegration may improve the accuracy of the approximation to allow for optimal results. An appendix presents some numerical considerations on the implementation that shows that the tensor product evaluations can still be preserved even in the case of overintegrations. In Section 3, we analyze the case where the non-constant coefficients are induced by the treatment of a deformed geometry. We prove that in this situation the overintegration is unnecessary to provide the optimality of the approximation. In Sections 2 and 3 we also present some remarks related to the impact of our analysis as regards the design of the best schemes in the case of deformed geometries. 2. Analysis of spectral approximation with non-constant coefficients 2.1. The one-dimensional case Let us consider the following problem: -(pu,),
+ qu =f
Find u E Hi(-1,
in (-1,l).
l), such that
(1)
This problem is supposed to be well posed in Ht(-1, 1) in the sense that the bilinear form a, defined over (Ht(-1, 1))’ as follows: 1
1
a(u, u) =
I
-1
pu,u,
dx +
I
-1
quu
dx ,
(2)
is elliptic and continuous over (Hi(-1, l))*. In fact, for sake of simplicity, we shall assume that there exist three constants pl, p2 and q2, such that O
In addition, approximate
VxE(-l,l), VxE(-1,l).
(3)
(4)
the forcing function f is supposed to be at least in L*(- 1,l). We want to this problem by a spectral element method. Our first work will be to design the
94
Y. Maday, E.M.
Rlnquist,
Optimal error analysis of spectral methodr
discrete space that will approximate Hi(-1, 1). Given a fixed integer K, we shall first consider a partition of (-1,l) in K subintervals Zk, where Zk = (a,, %+*) 7
(5)
and the ak are (K + 1) points in (-1,l)
such that
-1 = a, < a, < - - - < aK_l < aK = I
(and the length ak+l - ak is supposed to be O(K-‘)). The spaces of discretization are imbedded in YN = (9 E WI7
1) I ‘pllk E y&H
7
(6)
where N is some integer and lPN(A) denotes the set of all polynomials of degree SN over A. This space will discretize L “(- 1,l). In order to discretize Hi(- 1, 1), let us introduce
x, = r, n HA(-1, 1)
(7)
and cp vanishes at +l). (which means that qr,(a,+, ) = qIk+l(ak+l) The second step in the discretization process is to define the discrete problem. As indicated in the introduction, the point of departure is the following continuous variational formulation of the problem (1): Find u E Ht(- 1, l), such that VvEHA(-l,l),
a(u,v)=(f,v)
(8)
where the notation (e, .) stands for the L2(-1, l)-scalar product. the discrete problem is: Find U, E X,,,, such that
Given an integer M, M 3 N,
The definition of aM and (e, s)~ from a and (s, *) will use a composed Gauss-Lobatto quadrature formula. More precisely, this one is based on the data, over each segment Zk of a set of M + 1 points 53
and a set of associated weights c CD= 2 GL
z
p:
that are such that
@(&p:=I,@(x)dx
VC?EY,,_,
(11)
.
k=O i=O
The definition of aM and (e, *)M is now as follows: %(!W
4
= c PN*u, GL
wcp,
+ cGL W#N @) E
y;
*
W,,
UN)E Xi
3
(12) (13)
Y. Maday, E.M. R@nquist, Optimal error analysk of spectral methods
95
REMARK 2.1. The standard spectral element method utilizes a consistent Gauss-Lobatto quadrature rule, i.e. M = N. The case M > N corresponds to an ove~ntegration with respect to the standard spectral element method. Such an overintegration might be of interest in the case of non-constant coefficients. The analysis that follows indicates the situations for which this is the case. The interesting feature about (a, .)M is that it defines over Y, a discrete scalar product that is uniformly equivalent to the L2(-1, 1) one. Indeed, it has been proven in [4] that
The following Lemma is then completely
standard.
LEMMA 2.1. There exists one and only one solution u, to problem (9), and there exists a co~tant C, such that
PROOF.The
Lemma follows easily from the constatation that a, is uniformly and elliptic over X,, i.e., there exist two positive constants LYand y such that a
IIh4i+1,1~Cad+, 4 6 YI141&-1,1j7
and can be deduced
from (3)) (4) and (14).
continuous
(15)
q
Our further analysis is devoted to the derivation of optimal error bounds for the numerical solution. The following Lemma is also standard and follows from (15). LEMMA 2.2. There exists a ~u~t~~t C independent of N, such that
The following Lemma proves that the last term in (16) is of the same order as the best L2(-1, 1)tit off by polynomials of degree 2M - N. Let us first denote by 92:(-l, 1) the space Q-1,1)
= {$ E L2(-1,
and similarly, for any CL,by Hs(-1,
Then we have the following lemma.
1) I Wk E ~“o(f,,): $,I, = $k> , 1)
the space
(17)
Y. Maday, E.M.
96
Rbnquist,
Optimal error analysis of spectral methods
LEMMA 2.3. Let p be a real number >i.
For any 4 E Hc(-1,
l), we have (19)
(20) PROOF. Let us first consider the case where M = N. Let us introduce, for any integer L, the L2(-1, 1)-projection $L = II,(+) of + over P,(-1,l). Then we have
It is an easy matter to note that
(22)
I<+,4 - WA0 %)I Gw - ~~lIL2~-l.~~llVNllL2~-1,I) ’ Besides, let us denote by ,%Nthe operator such that $,I,+$)=+(&;)
Vi=O,.
of interpolation
. . ,N Vk=l,..
from the set %i(-1,
. , K.
1) onto
Y,
(23)
which, thanks to (14), gives
It remains to estimate the middle term in (21). Let us first assume that K = 1, corresponding to a plain spectral method with no domain decomposition. The term in question can then be written as follows: W/7 UN) - (V&Y%A = ($h&)[(LN,
LN) - (LN,
LJNI 7
where, for any i E N, ,. +j =
($7
tvN7
Lj)
(Lj, Li) ’
+ =
Lj)
(Lj, Lj)
(25)
Y. Maday, E.M. R@nquist, Optimal error analysis of spectral methods
97
are the coefficients of $ ( and thus of JIN) and uN, respectively, in the basis of the Legendre It is then an easy matter to note, from (14) that polynomials (Lj)jERI’
K&v, u,) - MN9 %),I Sz ~l~~~~lll~~ll2L~~-~.~)
442Nl ll~~ll~2~-~,~~l~~l IIb4L~e1.1)* Also, it follows easily that
So we deduce
The multidomain case where K > 1 is treated similarly over each subdomain Ik and the error is the same (the notations are in this case just more comp~cated). Let us recall that we have for any pa0 and any L .
We derive by sugaring
(21), (22), (24) and (26) that
‘Ibis completes Lemma 2.3 in the case M = N, after we recall that for any p > $ and any L we have [16]
II40- %4lL2(-1,1)6 CL-pIIdfP(-I,I)
WPE ~Y--~, 1) ,
(28)
so that the contribution of N-‘II$ - elN_*IILzt_l,l) is much smaller than the other terms. Let US now consider the ca:e where M > N and let us use the same strategy as before, with $zM_Nin place of +!+, where JI,,_, is some good approximation of JI of degree s2M - N that will not be necessarily its L2(-1, 1)projection and will be precised in the following. In order to prove Lemma 2.3, we must now estimate
The first term and the second term are treated exactly in the same way as previously and we
98
Y. Maday, E.M.
R)nquist,
Optimal error analysis of spectral methods
obtain
From the basic properties
of the Legendre
polynomials
we recall that
and that (n + l)&+,(X)
= (2n + l)&(X)
- nL,_,(x) .
A reiterate use of these two formulas yields
tL2M-N,LN)M=
4M-2N-1 244
_
N+l 2N
N
(&u-,-P
L,+&
= **
Now using (14), we derive that (L,,_,,
~3llLMllt+I.
L,),
1)
~3ll~ZM-NIIL~~-1.I)ll~NllL~~-l,l~
so that, as before,
J(&M-N, u,)
-
<&,_,,u,),I s CN-‘IIICI,,-, - &M-N-I IIG-IA.
(32)
The treatment of the last term in (29) has to be different. Indeed the same procedure as before would result in a bound of this term by the quantity ))I(,- $~+l)~2~__~,~)and thus a loss with respect to the optimality we could expect from the other terms. Besides, the numerical experiment (see Fig. 2.1) clearly shows a behaviour like II+ - $zM_NIIL2C_1,1).In order to prove this type of convergence, let us go back to the case K = 1 and note that (28) yields, for any CL> 1,
114,(PllL2(-1.1) 6 lI~lIL~(-*.1) + M-%lIH~(-1,1) v9 E fm-l~ 1).
(33)
From (14) it now follows that ($2,~,,
%)M - ($7 UN)M
= (G&j-N - $7 GM ~3lI$M(~Wv
- ~)lIL2~-1,1)ll~NllL2~-1.1,
+ M-‘lh,-rv - ~‘~H’(-~,I)I~~uN(IL~~-~,~) . <3[ll ICTZM-N - ~lIL~~-I,l, At this point the correct &,_, will be chosen. We see that it has to be a good approximation norm as well. This will be the case of the projection of $ with respect to of + in the H"(-1, l)-
Y. Maday, E.M. Rgnquist, Optimal error analysis of spectral methods
::I0 IL&____ 0
2
4
6
8
10
12
14
16
18
99
20
N
Fig. 2.1 A plot of various numerical results relevant to the approximation of the one-dimensional problem (l), in order to analyze the effect of overintegration. We have here considered the case where q = 0 and different values of p and f are simulated and the effect of overintegration is addressed. The relation between p and f is always chosen (of course, we do not use the knowledge of the exact in such a way that the exact solution is u = sin(fIrx)e-” solution in the simulations!!). The discretization is based on K = 4 equal spectral elements for different values of the polynomial degree, N. First of all, since this is the aim of the numerical scheme, we have computed the best fit of u in Ht(O, 1) by piecewise polynomials of degree d N (plot A); the plot proves a convergence better than exponential that is consistent with the fact that the solution is entire. Next we consider the case p = 1 with no overintegration (plot 0) and we can check that the approximation is optimal and very close to the best fit. We then choose p = 1 l(x + a) for various values of a. The case where a = 0.25, treated with no overintegration, is represented on the plot @ and we can remark that the approximation is still very good. The convergence rate however does not follow the best fit of u, but the best fit of p (consistently with Corollary 2.5) that is a straight line in this scale since p has a pole outside (0,l). Let us now use a more singular p, corresponding to a = 0.05. The plot 0 represents the results when no overintegration is used. It is very close to three other plots representing the results when overintegration (M = 3N/2) is used to compute only the right-hand side (the contribution off) (0) or only the left-hand side (corresponding to p) (plot A), while the third one represents the best fit of p (plot l) in the L*(O, 1)-norm by the discrete functions. In order to recover optimal results, overintegration must be used on both the left and the right-hand side of the equation (plot 0).
the W(-
1, 1)-norm over YN. It is proven in [17] that this element
satisfies
(34) and from (29), (30) and (32) the proof of Lemma 2.3 is complete M>N.
also in the case where
100
Y. Maday, E.M.
The second inequality vanish. 0
Ranquist,
Optimal errcv analysis
of spectral
methodv
(20) is quite trivial since the middle terms in (21) and (29) simply
We are now in a position to state the following Theorem. TUTORED 2.4. Let us suppose that the suluti~n u to problem (1) belongs to H”,( - 1,l) , that p belongs to Hk(-1, l), q belongs to Hcf-l,l) and that f belongs to Hi(-1, l), and in addition that the four real numbers CT- 1, CL,Yand p are larger than 1. Then the following error estimate holds :-
PROOF. It is a direct application
of Lemma 2.3 to derive
Let us choose w, as being the best fit of u in the H’(- 1, l)-norm.
It follows from (17) that
II@ - %?llH’(-l,I)si c~*-“ll~ll~~~-~,~~* If we now let G, be any polynomial
(37)
in X,, it follows from (3) and (4) that
so that sup
a@, UN)- %A%
qv)
uNEXh
The interest of this decomposition is that we can now choose GM, still close to u but such that the higher norm (i.e. H”,(-1, 1)) is unifo~ly bounded. This will now be the case if we choose GN as being the projection of u with respect to the H-(-l, 1)-norm over X,. It is proven in [17] that this element satisfies
101
Y. Maday, E.M. R)nquist, Optimal error analysis ofspectral methods
From the definition of a and uM we have
Using Lemma 2.3, first with 4 = pGNx and then with $ = qfiN, we derive that
l4L UN)- %t~;N,GI S C[(2M - N)-
min'u'o-l'~~P~Nxll~~in(v.o-l)(_~,~)
(2M - N)-
+
6 C[(2M
mi”(p,u)ll ~~,lln~“cr.~~(-~,~)lll~,IlH’(-1,I)
- N)- mi”(“.rr-l)llpUxIIHmi”c”,~-l,c_l,l,
+(2M - N)-
+ C[(2M
mi”(p9u)II 4~II~~~~r.~~(-~,l)lll~NIIH~~-~,1)
- N)- min(“‘ul-l)))p(U-
+(2M - N)-
~N)xlJHgi"(Y.~-')(_l,l)
mi”(p~rr)llq(~ - ~~~Il~~“c~.~~(-1,1)111~~ll~‘(-l,l) *
Let us recall that min(v, p, u - 1) > 4 so that H~in(vyo-l)(-l,l) algebras. This proves that
and I~IKrni”(~*~)(-l,1) are
146 NT4 - QMGL 4 SC[(2M
- N) -min(v,o-l)llpUxIla~“(u.~-l)(_l,l)
+(2M - N)-
+ C[(2M
- N)- min~“~u-l~llpll~~(_~,~~ll~~ - ~N)x(lHpin(u.o-l)(_l,l)
+(2M - N)which together
mi”~a~o~ll~~ll~~~~~~~~~~-~,~~lll~~II~~~-~,~~ min~~~u~ll~llHf;(-l,l)II~~ - ~,)o,~~~r.~~(-,,,)lll~NIIH~~-*,1) f
with (39) and the fact that M 3 N (so that 2M - N 2 N) gives
l4w”w %I - %fNv~ QJI
- N) -mi”(u*o-l)IJPUxJIHKmi”(“.a-l,(_l,l) -
N)-mi”‘~~o’~~q~~~H~.cr~u,~_l,l~]~~~N~~H~~_I,I~
+~~~l-(rll~ll~;;(-*,~)IIUll~~(-*,l) +~-“l14118~~-1,1~ll~IlHg(-1,1)lll~NllH~~-~,*~ * Theorem
2.4 now follows from the fact that H’(- 1,l)
is an algebra when r > 1.
Cl
Y. Maday, E.M.
102
R#nquist,
Optimal error analysis of spectral methods
COROLLARY 2.5. Let us suppose that the hypotheses of Theorem 2.4 are fullfilled and in addition let us suppose that M = N. Then the following error estimate holds:
REMARK
2.2. Interpretation of the results. In the practical computations, the more common situation is the one where the regularity of the solution is limited by the regularity of the non-constant data. More precisely, and especially if we consider the extension to the 2-D case that will be treated in the forthcoming section, we shall have (T - 1 s min(v, F, p). In this case, the estimates given in Corollary 2.5 prove that the optimality of the scheme is achieved with no use of over-integration. This statement is confirmed by the numerical experiments shown in Fig. 2.1. In the case where the regularity of the solution is better than the regularity of the non-constant data, we expect, from Theorem 2.4 that some amount of overintegration is necessary to recover the optimality of the discretization. The numerical experiments of Fig. 2.1 also confirm this statement with exactly the rate theoretically proven. It is important to note that the possibility to include some amount of overintegration is only true in the case where the variational formulation is the starting point of the discretization. The pure collocation method is incompatible with the concept of overintegration.
2.2. The two-dimensional case Let us suppose here that the problem is set on a square R = (-1, l)* and that no domain decomposition is used. The extension of our results to the spectral element case is very simple and it is only for sake of limitation of the notations that we restrict our analysis to this particular case. Let us consider now the following test problem: Find u E Hi(O), such that
V-A-Vu=f,
(42)
where f is a given force (supposed to be in L*(0)), A = (ajj) is a 2 X 2 matrix with non-constant coefficients and we assume, for sake of simplicity, that it is bounded, symmetric, positive and non-degenerate, i .e . , such that there exists 2 real numbers p1 > 0 and p2 3 0 such that AV-Vsp,(V;+V;) laij( Sp2
The variational
VV=(V,,V,)E[W*,
Vi, j=
formulation
a(4 u) = ((f,
v))
(43)
1,2.
of this problem:
(44 Find u E HA(O), such that
(45)
Vu E Hi(R)
will be used as a starting point to define the discrete problem.
Here
Y. Maday, E.M. R@nquist, Optimal error analysis of spectral methodr
a(u, u) = ((AVU, Vu))
103
V(U, u) E (H’(O))* .
The space of discretization, X,, consists of all polynomials (~~(~)~ that vanish at the boundary
of degree SN in each variable
x;, = UqL?) n H;(R). The quadrature formula is based on (M + 1)” points derived by a tensor product of the one-dimensional Gauss-Lobatto formula, and consists of the evaluation CGL,i CGL,2, where the indices 1 and 2 refer to the first and second spatial direction, respectively. The discrete problem is then: Find uN E X, such that
%f(%v~ %v)= (If7 MM where we have introduced
and
(46)
the following discrete scaIar product over PN(fi):
((R JI)), = o?, o;, @ I
WV E xiv>
,
V(cp, $1 f P&W2
%(%v, Qiv)= ((AV% V%)),
(47)
*
It is standard to note that, from (43) and (44), this problem has a unique solution and that in the case where M = N it is equivalent to a collocation problem based on the (N + 1)2 Gauss-Lobatto points. The analysis of the error will be based on Lemma 2.3 and on the follo~ng complement of this Lemma. LEMMA 2.6. Let p be a real number > 1. For any $ E Hc(- 1, l), we have
(48) PROOF. This proof is derived from (29) by taking $2M_N_1instead of I&~_~, which has the effect to cancel the middle difference in the right-hand side of (29). REMARK 2.3. It is an easy matter to note that (48) cannot be improved to yield a bound by the best fit by polynomials of degree 2M - N as it was in Lemma 2.3. We are now in a position to give the following error estimate. THEOREM 2.7. Let us suppose that the solution u to problem (42) belongs to H”(O), that the coefficients (aij) of the matrix A belong to H’(a) and that f belongs to HP(a) and let us suppose in addition that the real numbers p, u - 1, p are larger than 1. Then the following error
Y. Maday, E.M.
104
R@quist,
Optimal error analysis of spectral methods
estimate holds:
(4%
+w4 - WPllfllHP(~)] *
PROOF. The proof is derived from Lemma 2.2 (with obvious changes in the notations) following exactly the same lines as in the previous subsection. The only difference relies on the fact that Lemma 2.6 has to be used in the analysis of the quantity sup
4% 4 - %fh
4
vNExN
Here wN is a approximation bidimensional As in (38) we
sup
polynomial of degree N (and not N - 1) suitably chosen to be an optimal of u in the norm H”(R), say, and in any other lower norm from the equivalent to (39) (see [ 161). can write 44
UN) - 5&h,
uNEXN
IIhAil(*)
UN)
a(w,, 4 - a&h
uNEXN
where the first term on the right-hand
which can be bounded
s sup
II%hi~(*)
vd
+ clb - %&f’(R)7
side is a sum of four terms of the type
from Lemma 2.6 as
The reason why (20) is not sufficient to analyze this term is that in the second direction, of degree sN. It follows easily from the analysis in the one-dimensional case that we can prove
sup
(f’ ‘A’)- (f7vN)MC(2M Q
- N)-PllflI
HP(~)
’
uNXis
(51)
vNEXN
which allows us to derive Theorem
2.7.
Cl
COROLLARY 2.8. Let us suppose that the hypotheses of Theorem 2.7 are fulfilled and let us suppose in addition that M = N. Then the following error estimate holds:
IIu- u,IIHltn) s C(N-~~~min~Y~~~l~p~~ll~ll~~~~~II~ll~~~~~ + IIf IIwyn)) .
Y. Maday, E.M. R@nquist, Optimal error analysis of spectral methods
105
REMARK 2.4. Interpretation of the results. As pointed out in the previous subsection, there is no need to over-integrate in the general case where the solution has the same regularity as the non-constant coefficients. In two dimensions this conclusion still holds, however more weakly. Indeed, we have only been able to prove that the plain spectral or spectral element method provides the same accuracy as the best fit by the polynomials of degree one less. An overintegration with just one more point (M = N + 1) is sufficient to recover the optimal accuracy. Although the results are asymptotically the same, this remark could suggest the use of an overintegration. However, the associated increase in work makes it preferable to use a standard spectral element method with polynomial degree of N + 1 (instead of N) since this is less expensive than an over-integration, with M = N + 1 (see Appendix). It is only in the case where the spectral method would require much higher values of N that over-integration, coupled with the preservation of the tensor product formula, could be of interest. Finally, note that the numerical experiments again confirm the theoretical statements with an even less important difference between the cases M = N + 1 and M = N; see Fig. 2.2.
10-l
.
c
* lo-*
.
8
c
2
.
0
1o-3
0
E
. IO“ t LL
& 10-s p kz lo4
c
lo-’
p
lo+
c
lo+
r
.
Cl
8
10-I
-I
r
1
’
’
2
’
’
3
’
’
4
8
’
’
5
’
6
8
’
7
’
’
8
’
’ 9
’
1
10
N
Fig. 2.2. This figure represents the equivalent of the previous case in a two-dimensional domain. The equation is here (42) with A = pId and different values of p and f, computed so that the solution is always u = sin(nx)y(leye’). The domain of computation 0 is (0, 1)2 and decomposed in four squares. The plot 0 represents the case p = 1 with M = N, M = N + 1, M = 3N (there is no difference between the two last experiments, and the first case is a bit less precise but cannot really be distinguished from the two last cases). The plot 0 represents the experiment with p = ex+q that is, a very smooth function, and a treatment with no overintegration gives the same accuracy as the previous plot (that certainly also corresponds to the best fit). Now comes the treatment of the case p = 1 + cos(27rx + 37ry), which is a smooth function but is worse than u. The plot A corresponds to no over-integration and we clearly see that some degree of overintegration is necessary to recover the optimal accuracy, as shown in the case where M = N + 1 (plot 0) and M = 3N (plot A).
106
Y. Maday, E.M.
[email protected], Optimal error analysis of spectral methods
3. Analysis of spectral approximation in deformed geometries Our purpose is now to extend the previous analysis to the case where the non-constant coefficients are induced from the treatment of a non-rectangular geometry. Since the point here is not to analyze the behaviour of the approximation when original non-constant coefficients are present in the equation, we shall restrict our analysis to the simple problem of a Poisson equation in a deformed domain 0 which reads as follows: Find u E HA(a), such that
II
n
VuVv dx dy =
fodxdy
VVEH@).
(52)
In order to explain an important issue related to the error analysis, let us first consider the simple case where the geometry is rectilinear but not rectangular. 3.1. Rectilinear geometries Let us suppose that the domain 0 is the trapezoid with vertices A = (0,0), B = (4,0), C = (0,l) and D = (4 - 4/a, l), where a is a real number >l. It is an easy matter to check that the transformation F : (r, s) I+ (x, y) with 4 a
(
x=r4--s
> ,
(53)
y=s,
is a one-to-one mapping from the square fi = (0, 1)2 onto the domain a. We can now define a natural correspondence between functions 9 defined over 0 and functions C$defined over fi as follows:
G(r, 4 = cp(F(r,4) . This will allow to set the discrete problem. Indeed as we have indicated in the introduction, the spectral method has to be used on square domains. The other domains have to be mapped on such a reference square in order to allow for a spectral computation (the spectral element method allows to deal with more complex geometries, however, the elements must still be mapped to such a reference square). The variational formulation transferred to the domain fi is the following: Find li E Hi(b), such that
+$
2
o,(r,s)]drds=//fi
@J(r,s)drds
Vi,EHi(fi).
The coefficients J and oi, i = 1, 2, 3, are respectively the Jacobian of the transformation some smooth geometric factor given by
(54) F and
Y. Maday,E.M. R&quid,
J=Z(a-S),
US=(4-;)2.
(55)
and consists of the following:
Find
U2=(4-f)(2),
The discretization of the problem Li, E lPN(a) n HA(a), such that
107
Optimal error analysis of spectral methods
can now proceed
vo E P,(A)
)
(56)
where ((s, .))M has been defined in (47) and a,(ii,,
ti,) =
(( -
-
w2
-J
aii,-ai3, ar ' as M
The coefficients in the right-hand side of (56) are extremely regular. We can deduce from the previous analysis that no overintegration is required to compute accurately this term in order to achieve optimal approximation. On the contrary, the coefficients that arise in aM on the left-hand side, although analytic, may give rise to ,a problem since the factor q/J has a pole where J vanishes, i.e. at s = a. This pole is not in 0 but can be very close to it. The similarity of this problem with the one treated in the second section (Fig. 2.2) seems to imply that, when a is sufficiently close to 1, an overintegration is required to compute correctly aM in order to get an optimal approximation. However, surprisingly, the experiments performed in Fig. 3.1 show that no over-integration is required, even with those values of a that in the onedimensional case induced a large error associated with the consistency error (50). The interpretation of this fact can be easily understood, a posteriori, as the pole of l/J is compensated by the fact that the other part of the integrand vanishes and makes the term
more regular than it appears. In order to be more precise, let us state the following lemma. LEMMA 3.1. Let us assume that the solution u to problem real number u > 2. Then we have
(52) belongs to H”(a)
for some
(57)
108
Y. Maday, E.M.
Renquist,
Optimal error analysis of spectral methods
B 0
t lo-‘O e
lo-”
a 0
F
1o-L2
j ’
2
’ 3
I
’ 4
’
’ 5
’
’ 6
’
’ 7
’
’ 6
’
’ 9
’
’ 10
’
’ 11
’
12
N
Fig. 3.1. On this plot we prove that the singularity of the function l/Jacobian is not relevant for the need of overintegration. We consider problem (52) for two different domain decompositions A and B. They are
A
B
The three plots represent the discretization error in the semi-norm with no overintegration. The plot 0 corresponds to the solution u = sin(2y) for the two cases A and B. They both coincide since the deformation is in the x-direction and the solution only depends on y. Next, we repeat the experiments with u = sin( )x and we arrive to the (a priori) surprising result: plot a for A and q for B. We refer to Section 3.1 for the explanation. The fact that B is even better in this case is due to the fact that the discretization on the two domains is more equilibrated and that the left domain in case A is too long.
f
PROOF.
Let us go back to the original equation
(52) and remark that
and that a simple change of variable formula yields
Y. Maday, E.M. R)nqutit, Optimal error analysis of spectral methods
109
so that Q(k UN) - %.#i),,
u,) = A, + B, + A, + B,,, + B,,, + A, + B,,, + B,,, ,
(58)
where
A,=
G n + (r, s) 2
II
(6 s)(4 - z s)drds-(((4-zs):,?)),,
B~,1=(((~)~,~)),-(((4)~,~)),, B~,~=(((4)~,~))M-(((g)~,~))M. In the decomposition (58), the terms A_ are small from the accuracy of the quadrature formula; the terms B_ are small since GN is a polynomial that is close ,fo ii. More precisely, let us denote by C the supremum of 11/JI, I(1 - ol) lJ, (021Jl, 1t.03/Jlon 0. It is an easy matter to derive that, for any index cr, where (Y= 1, (2,1), (2,2), (3, l), (3,2), we have
where t in the derivation
stands for s or r. Using (33), with p = 1, we derive
(59)
110
Y. Maday, E.M. IQnquist,
Let us denote now by C?the supremum 2.2 that
Optimal error analysis of spectral methods
of 1, /(4 - 4sfa)l, 14rlal on fi. We derive as in Section
and the lemma follows from (SS), (59) and (60).
•!
In order to analyse the error in the approximation of problem (52) with (56), we first note that (16) (with obvious changes of notations) also holds in this case. This leads to the following theorem. THEOREM 3.2. Let us suppose that the solution u of problem (52) belongs to H”(0) and that f belongs to HP(R) and let us suppose in addition that the real numbers (+ - 1, p are larger than 1. Then the following inequality holds: II; - unrllul(~)s C~l-“llu~ln~~n~ + C(2M - N-- l)l-%IIn~~~~ + G(2M - WP
IIfllHPcn,*
(61)
PROOF. This is an easy consequence of (48) and of Lemma 3.1, if we choose tii), as being the best fit of ii in H”(b) by elements of PN(&). 3.2. General ~forma~ons Let us suppose in this subsection that the domain 0 is the image of & by some mapping F : (r, s) I-+ (x = F,(r, s), y = _F2(r,s)). We no longer assume tha! F is bilinear, nevertheless we assume it has some regularity, more precisely, that it is in Hail) for some p > 1, and that its Jacobian J is larger than some constant p1 > 0. The problem (52) is transferred onto fi as before and reads:
Find ti E HA(h), such that (54) is satisfied, where the geometric
factors oi, i = 1,2,3
are defined as follows:
The discrete problem is now: Find tiN E ~~(~)
n Hi(h)
Using here the same technique Lemma 3.3.
, such that (56) is satisfied.
as in the proof of Lemma 3.1 and Theorem
2.4, we obtain
LEMMA 3.3. Let us suppose that the solution u to problem (52) belongs to H”(a), for some real number cr > 2, that the deformation mapping F belongs to H’(h) for some real number
Y. Maday, E.M. Renquist, Optimal error analysk of spectral methods
111
u > 1. Then the following inequality holds :
REMARK 3.2. It is important, at this stage, to note that the regularity of G and the regularity of am where t is either s or t is bounded by the regularity of both u and F, which means that
is the best result we can expect. Note also (but this is natural after the analysis of Section 3.1) that it is only the regularity of F that is important and not at all the regularity of its inverse (as one could previously fear from the factor 1 lJ). From Lemma 3.3 and (65), we derive in a (now) standard manner the following theorem. THEOREM 3.4. Let us suppose that the solution u belongs to problem (52), belongs to H”(O), for some real number u > 2, that the deformation mapping F belongs to HF(@ for some real number p > 1 and that the forcing function f belongs to HP(a) for some real number p > 1. Then the following inequality holds:
+(2M -
N)-min(C’p)I(FII,~(n)lI f (It+n(r.p+a)} .
W)
REMARK 3.3. In the light of the previous theorem, it appears that there is no need to overintegrate in the case where the deformation is responsible for non-constant coefficients in the formulation of the discrete problem, even if the geometry is very distorted. This is readily seen by comparing the results obtained in the theorem by choosing M >> N and M = N (the gain of just one degree is not, as explained in Remark 2.4, sufficient to start the overintegration machinery). We can explain it a posteriori now as overintegration could allow to improve part of the consistency error term (57) in the case where the geometry mapping F is less regular then the solution u itself (contrary to the case in Section 2 this can happen, especially if domain decomposition is used). However, as pointed out in the previous remark, loss of regularity of F induces also a loss of regularity of & which is the function that is of importance since it is this one that is approximated by polynomials. REMARK 3.4. We have not considered here the case where the geometry mapping F itself is approximated by isoparametric polynomials. This effect can easily be analysed by following the same lines as previously and including the standard arguments adapted to the finite element discretization [ 181.
112
Y. Maday, E.M. loo
Rtinquist,
Optimal error analysis of spectral methods
k
I .
10-l
c
$2 A
sa . .
10-2 c
. IO”
A
0 0
c
A r
(L
w”
.
0
A
0
A
.
10-7
. A0
1o-5 c lo4
.
0 A
8
3 .
A A0
.
104
.
s
c
.
E loa
c
.
10-O c . 1
lo-
10' lo2 N Fig. 3.2. This plot is to prove that the overintegration can even have a negative effect, Again we consider (52) for p = 1, with a decomposition of a= [0,7] x [0, l] into two subdomains: 100
The interface is described by the equation x = $ - 3 ]y13. This equation is not regular and its interpolant is more rough when the degree is higher. Since this interpolant is the only geometric factor for the discrete problem we can understand the plots A and A that correspond to M = N and M = 2N, respectively, for the approximation of the solution u = sin( 4x). As before, if the deformation and the solution are in different directions (this is a very particular case!!!), for example u = sin 2y, the behaviour of the approximation is particular also and this is the only case where overintegration can help (in the case of deformed geometries) as is illustrated in the plots 0 and 0 corresponding to M = N and M = 2N, respectively.
Acknowledgment The question of the importance of overintegration has been raised during a meeting with I. BabuSka. The many discussions with A.T. Patera have been very enlightening. We thank them both for their help. This research has been supported in part by NSF under ASC-8806925, ONR and DARPA under contract NOOO14-85-K0208. Part of this work has been done while the first author was in residence at ICASE, NASA Langley Research Center, Hampton, VA, U.S.A.
Y. Maday, E.M. R)nquist, Optimal error analysis of spectral methods 10-Z
,
I
113
I
!a 1o-3
r
z
OA
A
lo+ :
A
-j
A
0
A
10-5 c
A
A
0
A
8
oz lo6
7
25
0 1o-'O loo
I
1
8
I-I I
h
lo*
10' N Fig. 3.3. A plot of the dkcretization error in the semi-norm 0 = [0,2] x [0, l] 1s * d ecomposed into two subdomains:
when solving
(52) with p = 0. The domain
I
1
and this time the geometry is tortured in the following way. The subdomain on the left is mapped onto the reference square by the mapping x = r + OS(lrl - 1) r”, while the mapping is affine in the other subdomain. The plot A corresponds to a solution u = sin( 4x) that only depends on x and is treated both with M = N and M = 2N. The plot 0 corresponds again to a particular case as it is related to u = sin 2y in both cases M = N and M = 2N.
Appendix (in collaboration
with A.T. Patera)
In this appendix we consider the computational cost of using overintegration. We show that the matrix-vector products can still be evaluated efficiently using tensor product forms, and that for a fixed discretization (polynomial degree N), the increase in operation count is at least a factor of d in R’. Overintegration is attractive to use only if the convergence rate (error as a function of polynomial degree) improves substantially. In this case we can achieve a fixed accuracy either with a low polynomial degree and overintegration, or with a higher polynomial degree with consistent quadrature. With no loss of generality we consider here the evaluation of the left-hand side of (46),
%M(%Q %J = ((PV%&J), *
Y. Maday, E.M.
114
R)nquist,
Optimal error analysis of spectral methoak
Here p represents a non-constant coefficient which either comes from the original strong formulation (as in (42)), or represents a geometric factor in the case of deformed geometry, or both. The evaluation of (67) represents the computationally most expensive part in an iterative solution of (46). In order to implement (46) we require a basis for our high-order polynomial space X,. The choice of basis does not effect the error estimates, however it greatly effects the conditioning and sparsity of the resulting set of algebraic equations, and is critical for the efficiency of parallel iterative solution procedures. We choose an interpolant basis to represent wN E X,:
p=o q=o
w,qq~)~q(4
(68)
9
where (x, y) E 0 -+ (r, S) E l-1, l[‘. Here the hp(.z) are the one-dimensional Nth-order Lagrangian interpolants through the Gauss-Lobatto Legendre points 5, (hp E P,(]-1, l[), hp( 5p) = $,,,), and wPq is the value of wN at the local node ( 5, , 6,). In addition to (68) we requn-e w, C X,,, to honor the homogeneous Dirichlet boundary conditions. The bases (68) are now inserted into (67), and the test functions are systematically chosen to be unity at one global node and zero at all the other Gauss-Lobatto Legendre points. We then arrive at the following discrete statement:
(69) where we for convenience only have included the derivatives with respect to the first (local) spatial direction. Here pP are the one-dimensional quadrature weights, and D and I are the one-dimensional derivative and interpolation operators defined as
Dpq=
2
(70)
c&J,
zpq =h,G,) Using tensor-product as follows:
(71)
* sum-factorization
techniques
where the expression in the innermost partenthesis order results in the following operation count:
[5] we can now evaluate (69) efficiently
is evaluated
first. An evaluation
in this
MN2+M2N+M2+M2+M2N+MN2=2(M2N+MN2+M2). In the case where M = N (no overintegration), and (69) reduces to the following expression:
the interpolation
operator
becomes IPq = aPp4
Y. Maday, E.M. R)nquist, Optimal errar analysis of spectral methodr
115
which can be efficiently evaluated in 2(N3 + N2) operations. We are now in a position to compare the computational cost with and without overintegration. In the case A4 = N + 1 the operation count is approximately twice the cost of using consistent quadrature (A4 = N) (in three dimensions the cost would increase with a factor of three). If A4 = 2iV, say, the cost to evaluate (67) is four times the cost with A4 = N. The factors which determine whether over-integration is economical are the difference in convergence rate with and without overintegration, and the specified accuracy of the discrete solution. References D. Gottlieb, The stability of pseudospectral Chebyshev methods, Math. Comp. 36 (1981) 107-118. ;;; B. Mercier, Stabilite des methodes spectrales polynamiales application B l’bquation d’advection, RAIRO Anal. NumQ. 16 (1982) 67-100. Y. Maday and A. Quarteroni, Approximation of Burgers equations by pseudospectral methods, RAIRO [31 Anal. Numtr. 16 (1982) 375-404. 141C. Canuto and A. Quarteroni, Approximation results for orthogonal polynomials in Sobolev spaces, Math. Comp. 38 (1982) 67-86. S.A. Orszag, Spectral methods for problems in complex geometries, J. Comput. Phys. 37 (1980) 70-92. ;i; Y. Morchoisne, Resolution des equations de Navier Stokes par une methode spectrale de sous domaines, in: Proc. 3rd Intemat. Conf. Numer. Methods Sci. Engrg. (Gamni, Paris, 1983). [71 B. MCtivet, Resolution spectrales des equations de Navier-Stokes par une methode de sous-domaines courbes, These de Doctorat, Universite Pierre et Marie Curie, 1987. PI C. Canuto and D. Funaro, The Schwarz algorithm for spectral methods, SIAM J. Numer. Anal. 25 (1) (1988) 24-40. 191A.T. Patera, A spectral element method for fluid dynamics: Laminar flow in a channel expension, J. Comput. Phys. 54 (1984) 468-488. [101 Y. Maday and A.T. Patera, Spectral element methods for the incompressible Navier Stokes equations, in: A.K. Noor, J.T. Oden, eds., State of the Art Survey in Computational Mechanics (1989) 71-143. incompressible Pll E.M. Ronquist, Optimal spectral element methods for the unsteady three-dimensional Navier-Stokes equations, Ph.D. Thesis, MIT, 1988. P21 C. Bemardi, Y. Maday and A.T. Patera, A new nonconforming approach to domain decomposition: The mortar element method, in: H. Brezis and J.-L. Lions, eds., Nonlinear Partial Differential Equations and their Applications (Pitman and Wiley, 1989). 1131Y. Maday, C. Mavriplis and A.T. Patera, Nonconforming mortar element methods: Application to spectral discretizations, in: T. Chan, ed., Proceedings of the 7th International Conference on Domain Decomposition Techniques for P.D.E. (SIAM, Philadelphia, 1988). 1141C. Mavriplis, Nonconforming discretizations and a posteriori error estimators for adaptative spectral element techniques, Ph.D. Thesis, MIT, 1989. 1151D. Funaro, A multidomain spectral approximation of elliptic equations, Numer. Methods in P.D.E. 2 (1986) 187-202. Y. Maday, Resultats d’approximations des operateurs d’interpolation aux points de type Gauss, Notes aux WI Comptes rendus de 1’Academie des Sciences de Paris, submitted. 1171C. Bemardi and Y. Maday, Approximation results for spectral methods with domain decomposition, J. Appl. Numer. Math. 6 (1989-90) 33-52. WI P.G. Ciarlet and P.A. Raviart, The combined effect of curved boundaries and numerical integration in isoparametric finite element methods, in: A.K. Aziz, ed., The Mathematical Foundation of the Finite Element Method with Applications to Partial Differential Equations (Academic Press, New York, 1970) 409-474.