- "i 

NASA Con&oi Report 181620 

ICASE REPORT NO. 88-11 

ICASE 

DOMAIN DECOMPOS ITION PRECONDITIONERS 

FOR THE SPECTRAL COLLOCATION METHOD ~ ~ 

(#iISA-CR-l816Zo) DGflAIIH &ECCCfCS.I!XION 3a8-18345 

IbECCIDIIICIEBS EC6 TLL SFECSEBL CCLLOCATIOH 

itErHeC Einal Regcrt (sass) TE F CSCL 12s 

Unclas 

G3/64 0126002 

Alfio Quarteroni 

G iov anni . Sacc h i-Land r i ani 

Contract No. NAS1-18107 

January 1988 

INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING 

NASA Langley Research Center, Hampton, Virginia 23665 

Operated by the Universities Space Research Association 
DOMAIN DECOMPOSITION PRECONDITIONERS 

FOR THE SPECTRAL COLLOCATION METHOD 

ABSTRACT 

We propose and analyze several block iteration preconditioners for the solution of elliptic problema by spectral collocation methods in a region partitioned into several rectangles. It is shown that 

convergence is achieved with a rate which does not depend on the polynomial degree of the spectral 

solution. The iterative methods here presented can be effectively implemented on multiprocessor 

system due to their high degree of parallelism. 

Research was supported under the National Aeronautics and Space Administration under NASA Contract No. 

NAS1-18107 while the first author was in residence at the Institute for Computer Applications in Science and 

Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 

P) Dipanimento di Matematica. Universita Cattolica di Brescia. 25121 BRESCIA 

(**) Istit~to di Analisi Numerica del C.N.R.. C.so Carlo Albem 5, 27100 PAVIA 

i 
ODUCTION 

We consider the numerical solution of elliptic problems by spectral collocation methods in regions 

partitioned into several rectangles (also called subdomains). The interaction between the 

subdomains is handed by propa iterative methods. 

Such iterative domain decomposition methods am particulary attractive in the context of spectral 

approximations. Actually, they allow the reduction of a problem set in a complicated geomeuy to a 

sequence of problems of similar type but with smaller size in every subdomain. The latter can 

therefore be faced by standard monodomain spectral solvers. Besides, due to their (hopefully high) 

degree of parallelism, these methods can be advantageously implemented on multiprocessor 

systems. 

To be effective, the rate of convergence of the iterative method should not deteriorate as the 

polynomial degree N of the numerical solution in each subdomain increases. This is what we prove 

for our method which alternates the solution of Dirichlet problems on the odd subdomains to that of 

mixed Neumann-Dirichlet problems on the even ones. Our proof applies to any domain R 

decomposed into rectangles without internal crosspoints, i.e., whose vertices lie on the boundary 

of R (see fig.0.1). In the last section we consider a decomposition with internal cross points; in 

this case it is possible to prove the convergence of the iterative method for fmed N, although the 

rate of convergence decreases if N increases. This convergence behavior with respect to the 

decomposition is shared by other multidomain approaches based on finite elements or finite 

differences. We refer, e.g., to [41, [51. 

RG.O.l a) Decomposition without internal cross points b) Decomposition with an internal cross 

point 

The results we obtain generalize those of [15] which were established with a different proof 

technique and were limited to the case of a rectangle divided into two subdomains only. 

The outline of this paper is the following. In section 1 we write the differential problem in a 

variational form which is set up on the muladomain decomposition of the physical domain. 
~~ 

In section 2 we introduce th 

-2 - 

multidomain ollocation problem. This amounts to collocating the differential problem at the Chebyshev (or Legendre) collocation points internal 

to each subdomain, and to enforcing, at the interface-boundary points, the continuity of 

the solution as well as of a suitable combination of the residual of the equation and the normal flux. This method has been proposed by Funaro [13] and retains the same convergence 

properties of the more classical patching method introduced by Orszag [23]. In particular, as shown in [13] and [19], the numerical solution converges with spectral accuracy to 

the physical solution as the polynomial degree N grows on each subdomain. In section 2 

we also define an “iteration-by-subdomains” method for solving the discrete multidomain 

decomposition problem. Each step of this iterative method is based on the solution of 

Dirichlet problems in the odd subdomains, and mixed problems in the even ones. In section 3 we give a variational formulation of both the multidomain decomposition problem 

and the iteration-by-subdomain method. These variational formulations allow us to prove, 

in section 4, the convergence of the iterative method, the rate of convergence being independent of the degree of the polynomial solution in each subdomain. In section 5 we show 

that the iteration-by-subdomain method here analyzed amounts to solving the capacitance 

system governing the interface unknowns by a linear iteration procedure using a proper 

block diagonal preconditioner. Going further along this equivalence, we propose in section 

6 other preconditioners for the capacitance system, which in turn gives rise to some new 

iteration-by-subdomains algorithms. Finally in section 7 we consider a decomposition of 

the domain with internal cross points and we state how our method can be formulated 

and analyzed in this case. 

Numerical experiments on this itation-by-subdomain method are presented in [ 131, [28]. 

Earlier works on algorithms of this nature for finite element and finite difference domain 

decomposition methods are reviewed and analyzed in (261. 
-3 - 

1. ~.~Im.mFO RMULA "ION OF THE Dm RENTIAL PROBLEM 

Let R be an open two dimensional domain with boundary 20. We consider the 

boundary-value-problem 

(1.1) Lu=f inR, u=O on2R 

where f is a given function of L*(Q and 

In variational form (1.1) reads as 

UE& (Q): a(u,v): = (VU vv + %UV)~X = fvdx v v E H: (Q) . i I n 

(1.3) 

The Lax-Milgram theorem ensures that there exists a unique solution to (1.3); moreover, 

ue@(R), hence (1.1) holds almost everywhere in R (see, e.g., [18]). 

Let us consider a partition of R by non intersecting open subdomains R, i = 1, ..., M. as in fig.l.1, 

n 

aqd denote by ri the common boundary between R, and Qi+*, i I M - 1. Then set r = UT;. and, 

for each i, we denote with 0 := n 

i=l 

M- 1 

(ri) the space of functions defined on r, which are the 

traces on r of functions of HI@) (see [18]). Note that these functions vanish at the endpoints of 

ri for i=l, ..., M-1. 

r i- 1 

Q. 

I 

Fig.l.l: The partition of R 
-4 - 

For any cpe0, let cp E H'o(f2) denote any extension of 9 to 

"harmonic" extension Rcp, wbosc ddon \cp to Oi satisfes the Dirichlet problem: 

An example is furnished by the 

wherc we have set 

% (u,v) = (vu vv + Qo uv) dx . I ni 

The following equivalence statement is dy proven 

problem (1.3) is eau ivalent to fi nd UE Hlo(R) so tha& u;= ula. 1 5 i I M, 

3 (up) = fvdx VvcHi(Ri) , lSi5M , 

4 

on ri , ISiSM-I, 

(1.5) I 

ui = ui+l 

ui-0 onaRinaR , IliSM , ( 1.7) 

Clearly, (1.5) amounts to require that 

(1.9) Lui=f inRi , 1SiSM 

while (1.8) is equivalent to the nansmisswn codnon 

(1.10) 

aui aui+l -=- onri , 15iSM-1 

hdu 

(u is the outward unit vector to Ti ). 
-5 - 

In this section we consider a spectral collocation app-roximation of problem (1 .S)-( 1.8) 

which consists of collocating the differential equations (1.9) at Gaussian collocation points internal 

to Ri and to enforcing the interface conditions at some selected points of r. Either 

Gaussian-Legendre and Gaussian-Chebyshev points will be considered. We also give a variational 

formulation of the multidomain problem, then we detine an "iteration-by-subdomains" procedure to 

solve effectively this problem. 

We confbe ourselves to the case of a rectangular domain R partitioned into M non intersecting 

rectangles ni without internal cross points (see fig.2.11, but what we are proving can be extended 

to more general situations as well (cf. remarks 3.1 and 4.2). 

I 

FIG.2.1 Decomposition by aligned subdomains 

Let D be the reference domain (- 1, 1)2 and let N,, N,, be two given natural numbers. The Legendre 

Gauss Lobatto collocation points in D are the roots (Ckm, 0 I k I N,, 0 I m I NY) of the 

polpornid 

where h(t) is the Legendre polynomial of degree k in [-1,1]. In the Chebyshev case, the Gauss 

Lobatto collocation points in D arc given by 

, forOSkSN, , OSmSNy . 
-6 - 

We remind that, in both Legendre and Chebyshev cases, lies on the boundary of D if either k 

equals 0 or N, or m equals 0 or N,,. Moreover 

Ci := [(xi, ym) , 0 S k 5 N: , 0 S m 5 N,,) , for i=l, ..., M , 

will denote the set of collocation points in the domain Ri. Given M natural numbers Nix, i = 

1, ..., M, we take the points of Ci as the images of the Legendre (or Chebyshev) points in D, with 

Nx = Nix. through the affme transformation which maps D into Ri. For the convenience of 

exposition we set 

Cy := Ci n Ri , C: := Ci n aR , Cri := CiYr u e) for i=l, ..., M . 

See fig.2.2 for an example with M=2. 

0 pointsofcb~ 

0 pointsofCb2 

+ points of Clint 

e points of $1 

0 points n cb2 

points of ~2int 

Fig.2.2 Interior and Boundary Collocation Points 

For i=l, ..., M we denote by PN(Ri) the space of polynodals of degree N', with respect to the x 

variable and degree Ny with respect to the y one. Morcover we set 

For the convenience of the reader we recall the Legendre and Chebyshev Gauss-Lobano formulae 

that will be extensively used in the sequel (see also [lo]). 
-7 - 

e Gaius-Lobatto fonnulg 

where 

(2.2) 

We have set = xi - xi-l, where xi denotes the abscissa of ri. 

(2.3) 

where 

'ni 

2 -lR [ai (x,y) = I(X - XiJXi - x)]-'R (1 - y ) , 

x x for m=l, ..., N, - 1 , W, = - for m=O and m=Ny . 2NY 

We can now introduce the specnaf coffocurion upproximarion to the problem (1.1), which 

fits well the multidomain formulation (1.94 1.8). The spectral solution uN verifies: 

(2.5) U* := uN~ E PN(S2,) , i=l, ..., M , 

(2.6) LUiN = f in q' , i=l, ..., M , 

(2.7) ua - u~+~,~ in C , i=l, .... M-1. 

(2.8) 'ifi =O in C: , i=l, .... M , 

- Ti 

ri &i,N aui+l.N + -ax - 7 = - (LU~ - f) wi - (LU~+~,,, - f) wi in C , i=1, ..., M-1 
-8 - 

-i + i+l 

Hcre we have set ai = o . and ai = a,, for i-1, ..,h4- 1. We are assuming f to be a continuous 

function in der that (2.6) and (2.9) make sense. Note that at the interface collocation points, the 

jump of the normal derivative is not asked to vanish (as in (1 .lo)), but rather to balance a suitable 

linear combination of the residual of the equation from both sides. We observe,however, that the 

coefficients of such combination tend to zero if N', tends to infinity. This follows from the 

defintions of the weights dk of the Gauss-Lobatto famula given in (2.2) and (2.4). 

From a computational point of view, dealing with condition (2.9) rather than with the pure flux 

condition a;; uiN = m; ui+lJ does not bring any appreciable extra work. 

considerable advantage is achieved for the theoretical properties of the scheme, as is 

* 

a a 

On the other hand, a 

will be pointed out in the fonhcm 

The collocation scheme (2.5)-(2.9) has been first proposed by Funaro [13], who 

sassed its equivalence with a suitable variational approximation to (1 3-( 1.8) (see below).This 

equivalence allows one primanIy to investigate the convergence properties, for Nix. Ny + -, of 

themethod. 

The most relevant result (for the Legendre collocation) is the following (see [ 131 and [19]). 

Semng Ni = min('N; , Ny) and assuming that ui E H '(QJ with bi 2 1, for i=l, ..., M, the following 

MOT estimate holds 

a. 

M 

(2:lO) 

where C is a positive constant independent of Nj and u-, for j=l, ..., M. In particular this estimate 

shows that if f = 0 then u~,~= 0 for i = 1, ..., M, hence the discrete problem has a unique 

solution. 

I 

Another implication of (2.10) is that as WX,Ny + - the jumps in the normal derivatives 

a;; us - bx u~+~~ at the interface points vanish with spectral rate, hence so does the right hand side 

of (2.9). In 1131 it is also shown by numerical experiments that using (2.9) rather than the pure 

flux condition does not reduce the accuracy of the results. 0 

a a 

We now introduce an iteration-by-subdomains procedure for the solution of the above 

collocation problem. Let us suppose, here and in the sequel, that M is an odd number. the case of 

M even can be studied analogously. We set r = u Ti, and define 
-9 - 

PON 0 := (p : r + R : pi := pln is a polynomial of degree S Ny on ri. 

vanishing at the endpoints of ri) . 

00 For a given g E PN (Q we look for a sequence &,n 2 1, satisfying: u>EP~(RJ, i=l, ..., M,and: 

for: 

(2.11) Luna = f in Cih' , 

(2.13) u"a-0 inCbi ; 

(2.14) LuniN=f inqint , 

- 

ri inc , 

(2.17) una = 0 in Cbi 

Then set: 

r .- I en UY.1.N + (1 - en) gn" on C' i-l for i=3, ..., M , i odd 

(2.18) gn = .- 

[e, u:+lS + (1 - en) gn-' on C' for i=1, ..., M-2 , i odd 

and rcstaxt from (2.1 1) with n+l instead of n. In (2.18). 8, is a positive relaxation parameter, how 

to choose it will be dealt with in the remark 5.4. 

Following this method, at each step one solves the M/2+1 independent Dirichler problems 

(2.1 1)-(2:13); then, after computing the (independent) right-hand sides of (2.1342.16). one has 

to solve the M12 independent mixed problems (2.14)-(2.17). The degree of parallelism of this 

algorithm is therefore W. 

In the simple situation of a &comDosiaon bv t wo SUMO mains, the above method simplifies .. 
- 10- 

to: 

(2.19) LIP18 = f in Clint , 

(2.20) un18 ,g"-l inC,' , 

(2.21) un18 =o , 

(2.22) Lunm =f in cz'" 

(2.24) unm = 0 9 

and gn = en unW + (1 - 0") g"'' on Cr. Here r is the interface between R, and "2. 

B€lwu2 (Mo=eeneraldiff-torS) The muhidomain approach (2.5)-(2.9). and the 

corresponding iteration-by-subdomain procedure (2.1 1)-(2.18). apply also to the more general 

differential opemm. 

with a. symmetric. uniformly negative definite, bounded, and smooth on R. In such a case, the 

operator L in Ri should be replaced by its pseudo-sgecrrd realisation 

U 

(2.26) 

Here pN v denotes the polynomial of PN(Ri) which interpolates v at the collocation nodes of Ci (see 

[a], Ch. 2). Moreover, the conormal derivative allah + a,,d/ay should be used instead of the 

nOnnal derivative ah at the subdomain interfaces. 

In particular, we remind that this situation occurs whenever we consider the Laplace equation 

within a domain R subdivided into subdomains with curved boundaries. The mapping of the 

physical subdomain R, into the computational one D=(- 1 ,1)2 introduces a transformed operator 

which has the form (2.25).0 

Remark 2.3 The iteration procedure for two subdomains (2.19)-(2.24) has been introduced in 

[lS] for Chebyshev collocation methods, and subsequently applied in [21] to finite element 

approximations. The convergence analysis was carried out in [ 151 for the Helmholtz operator using 

a separation of variable argument. In section 4, we apply a proof technique based on a suitable 

extension theorem, similar to the one presented in [21] for finite elements. The higher generality of 
- 11- 

this technique allows us to deal with the case of several subdomains.0 

3 VARIATIONAL FORM OF THE CO LLOCATION MULTIDO MAIN PROBLEM 

In this section we give a variational formulation of the collocation multidomain problem 

(2.5)-(2.9), and, accordingly, of the iterative procedure (2.1 1-)-(2.18). The new formulation is 

more suited for the convergence analysis that will be canid out in section 4. 

We confine ourselves to the case of Legendre collocation points, where we will make use of the 

following notation. For functions continuous in Ri we define a discrete inner product, 

approximating that of L2(Ri), as follows: 

(3.1) (w,z)* := (wz) (XI. y,) dk om , i=l, ..., M 

k=O mlO 

Due to the exacmess of the quadrature formula (2.1) for polynomials g E Pm-l, we have 

(w,z)iN= wzdx ZWZE P*N-* (Ri) , i=l. ..., M I (3.2) 

n, 

We also define the discrete bilinear f&m 

(3.3) 

For any cp E PoN (r), we denote with E PN (R) any piecewise polynomial extension of cp to R 

which is determined by the values of cpi-l and ‘pi solely. can be, e.g., either the discrete 

harmonic extension, or the interpolant extension of 9. We can now prove the following 

equivalence result. 

The coll-n multldomain problem (2.5) - (2.9) keauivalent to loo foE ‘N 

SUMUS u~N := uN~ E PN(Ri), i=l, .... Mo 1 
mf This proof generalizes the proof of an analogous equivalence statement given in [ 131 for a 

decomposition by two subdomains. 

We start proving the following discrete integration by parts formula: 

for i=l, ..., M , V W, z E P&) such that w=z=O on aRin aR . 

Indeed, by (3.3) and (3.1) we obtain, using the propexty (3.2): 

Now we get (3.8) using again (3.2). Let now uN be the solution to (2.5)-(2.9). Equations (2.6) 

give: 

This yields (3.4) by virtue of (3.8). 

We verify now (3.7). For this, let i be an odd integer between 1 and M, and let cp be any function 

of Po, 0.h Ri one has = @-l) + 4') , where @) is a polynomial extension of the function 
- 13- 

r. r1 ti) which coincides with cp on C and vanishes on C , 1 + j. In particular, we can set w = 0 on 

if 1 # j, j + 1. Using (3.8) we have 

In a similar way we find that: 

By summation of the last two equations it follows: 

It is now an easy task to see that (3.7) can be obtained by summing the relations (3.10) on all odd i 

between 1 and M. 

We shall now prove the converse, namely that the solution of (3.4)-(3.7) is also solution of 

(2.5)-(2.9). To this end, for each i-1, ..., M let P denote a point of Ck The discrete characteristic 

function associated with P is the polynomial of PN(Ri) which equals 1 at P and vanishes at the 

remaining points of Ci. Equations (2.6) arc obtained fnt~ (3.4) &rig 9 to be either of the discrete 

characteristic functions associated with either point of Clint. Similarly, (2.9) follows from (3.7) 
-14- 

taking as cp a function of PNy cri, which equals 1 at a point of cri, and vanishes at all the other 

points of Cri, then working as done for proving (3.10). Clearly, the relations (2.7) and (2.8) 

follow from (3.5) and (3.6).0 

In view of the Proposition 3.1, it is now clear that the iterative procedure (2.11)-(2. IS) admits 

the following equivalent variational formulation: 

For all iadd (i=l, ... N) solve 

(3.1 1) 

(3.12) 

(3.13) 

aiN (una, v) = (f. vIiSJ V v E PoN(Qi) , i=l, ..., M , 

en-1 u:-& + (1 - en-l) ut; ul&= 

on , i=3 ,.... M-1 

ui+lsJ + (1 - u$ on Ti , i=l, ..., M-2 , 

then for all ieven. (i=2, ..M- I) solve 

(3.15) Una = 0 

For each i odd, (3.11)-(3.13) is a Dirichlet problem, while, for i even, (3.14)-(3.17) is a mixed 

type problem. The above variational formulation of the multidomain collocation iterations 

(2.11)-(2.18) will be used in the next section to investigate the convergence properties 

of this method. 

Pernark 3.1 If R is a plurirectangle, partitioned into subrectangles L$ without internal cross 

points (as in fig.0.l.a). the above multidomain problem and the relative iteration procedure 

(2.11)-(2.18) can still be formulated. In such case each venex of the decomposition lies on the 

boundary of R, hence the value of each iterate u"i~ is set to zero (the prescribed boundary data) 

there. 0 

Remark 3.2 For one dimensional multidomain methods using the Chebyshev collocation points an 

equivalent variational formulation is still available (see [ 141). However, the bilinear form a*(-,.) 

that holds such an equivalence is much more involved than the Legendre one (3.3). (Actually it 

must take into account the Chebyshev weight functions which blow up at each subd+ 

main interface). The generalization of this equivalence to two dimensional problem is 

not yet available. Thus the convergence analysis we carry out in the next section will 

deal with Legendre collocation points only.O 
-15- 

4 CONE RGENCE ANN.YS 15 

t 

In this section we confine ourselves to the case of the Legendre collocation. We pve that the 

sequence (u",&), solution to (2.1 1)-(2.18) (or, equivalently to (3.1 1)-(3.17)) converges, as n -+ 

QO, to u, and that the rate of convergence is independent of N. 

For any (PE PN(T) we recall that cponT,isdenoted by cp,, i=l, ..., M-1; then we define 

R, 9 E PN (ni) through the equations: 

(4.4) R,(P=O on a R n a R,. 

If i=l or i=M the conditions (4.3) and (4.2) respectivey must be dropped. In view of the analogy 

between this problem and (1.4). R, cp is referred as the discrere harmonic exrenrion of cp to Ri. 

We now defme the following norms. 

fM 

(4.5) 

(4.6) 

, 112 

In the next lemma we state that the norms (4.5) and (4.6) are uniformly equivalent on P", 0. This 

is an extension type theorem for spectral collocation approximations similar to the one stated in 131 

and [21] for finite elements. Its proof is essentialy based on a result due to Bemardi and Maday 

PI. 

Theorem 4.1 ylereex ist two pos itive u 'L independe nt ofN s!Kuia 
(4.9) 

-1 

Far a general rectangle R, it is easy to set that C, is proportional to 1 + meas <Ri) 

Setting v = RiN 9 - QiN cp in (4.1) we obtain 

On the other hand we wall that (see [ 121 and [ 6 * Ch. 1 11) 

(4.1 1) 

and 

(4.12) 

where C, and C, arc positive constants independent of N. Hence, from (4.9) and (4.10) we 

deduce 

(4.13) 
-17- 

In order IO estimate the right hand side of (4.13) we usc l~nw thc mtinuity of the b'tce operator 

from H1 (Qi) to H'" (a Qi) (see [181) 

where C, is proportional to p (M) 

Using again (4.1 1) wc also get 

For i=l or i-M working in the saine way we have 

Finally summing up (4.14) and (4.15) for i odd and for i even we deduce the desired result (4.7).0 

In order to prove that the solution un, of (3.11)-(3.18) converges, as n -+ to the one 

ug of (3.4)-(3.7)- we look for the operator T, such that 

(4.16) qaN n+l - uiN = Te (uiJ n - us) on r * i even. 

The operator T, can k defined as follows. For any cp E Po, 0 let wih' E P, (Q,), i=1 ,... *M. i 

m, be defined through: 

(4.17) a. bh . (wiN . v) = 0 v v E pc, (R,) 

M M 

i even i odd 

(4.19) w,=O on a R n a R, 

Then we define 

(4.20) 

urd 

(4.21) 

Tq := wiN for i=l, ... *M-l * iri 

Tgq:=eT9+(1-e)cp V~EPO~~ . 
-18- 

Due to (4.1) and (4.17) one could substitute RLN yf in (4.18) with any other extension 

of y in Qp without affecting the definition of w,.O 

Remark 4.2 The operator T can be interpreted as follows (suppose M > 2). For a given cp E Po, 

0 and an even number i with 2 S i 5 M-1 one can find Ri-lJ cp and Ri+l,N cp by solving the 

Dixichlet problems (4.1)-(4.4). Then compute the left flux at the boundary side ri-] and the right 

flux and the boundary side ri as follows 

for m=l, ..., N-1. The functions ws arc now obtained by solving the mixed pmblem 

Finally Tcp on Ti-1 and Ti is obtained by computing wiN on Tiel and Ti. 

The operator which associates to the values of Di-l, and Di+lm the values of wiN on Ti-1 and Ti 

is a discrete Poincd-Steklov operator (see [ 11). 

R 

i- 1 

R 1 

R 

i+ 1 

FIG.4.1 Computation of T cp 
-19- 

In the following theorem we prove that for suitable values of 8, Te is a contraction. This, together 

with (4.16). will imply the convergence of una to uiN for n + -. 

rem42 ForMM22ereex &B*>Osuchthatfo r any N > 0 the followine holds 

(4.23) v e E [e, e**] k(e)SK<l . 

From (4.5) and (4.21) we have, using the symmeny of the forms aa (e,.) 

Moreover, using (4.17) and (4.18). we have 

hence, from theorem 4.1 and (4.24) we deduce 

From (4.25). the Cauchy-Schwarz inequality and (4.7) we also have 

(4.27) I1 Tc~ l,, S Jo I1 ~p I,, , 

In order to find a lower bound for I1 Tc~ I,, we remark that from (4.17)-(4.18) we have the 

fundamental relation 
-20- 

M M M 

iodd i even i even 

Using the Cauchy-Schwarz inequality, the property (4.7) and the definitions (4.5) and (4.6) we 

deduct: 

From (4.26),(4.27),(4.28) it follows 

If we &fine 

We can readily sce that (4.22) holds and that 

k(8)<1 if€ 04<8*=mh 

Then (4.23) follows from the continuity of k(8).0 

We can now derive the following convergence result. 

4.1 & ui3, i=l, ... , M, be the so lution of (3.4)-(3.7) and let uni,~, for n 2 1 and 

i=l, ... , M, be the so lution of (3.11)-(3.18). Lfene [e, e"] we havt 

moreme r 

(4.30) 

whm C is a consta nt deDe ndent on go but indebendent of bot h N and M. 

Proof 

4.1 and to (4.11). the estimate (4.30) can be deduced from (4.29).0 

The estimate (4.29) follows using repeatidly (4.16) and theorem 4.1. Due to the theorem 

Remark 49 In view of (4.8), if H is much smaller than one, then ts and I: grow like H-'. If 
-21 - 

t+ 1 

0 r+7+2 

is the value of 8 that minimizts K (e), one has 8- = 

(1 - CH2)ln, where C is a positive constant. 0 

The proof of theorem 3.1 and corollary 3.1 can also be applied to the case of a 

plurircctangle R partitioned into rectangles i2, without internal cross points (see fig.O.1 and remark 

3.1).0 

and K (eqJ behaves like 

To solve the multidomain problem (2.5)-(2.9) we have proposed, in section 2, an iterative 

method based on a sequence of differential probl&s to be solved in each subdomain. An (a priori) 

different point of view consists of applying the influence (or capacitance) mamx method for 

solving (2.5)-(2.9). This approach, in the framework of spectral approximations, has been 

extensively pursued by Peyrct and his Coworkers (see e.g. [ 171, [24]) and by other authors (see 

[20],[22]). The influence matrix coincides with the Schur complement of the matrix of 

the system (2.5)-(2.9) with respect to the interface variables (see e.g. [9]). It is precisely 

the matrix of the system of the interface unknowns and it is derived from the global 

system by block Gaussian elimination (see e.g. [3]). In this section we show that the 

iteration-by-subdomain method introduced in section 2 is equivalent to a preconditioned 

Rlchardson iterative method for the resolution of the influence system. In particular we 

will give the precise structure of the preconditioner for the influence matrix. 

We proceed now to construct the influence matrix associated with the problem (2.342.5 j. 

Let (qie Po, (r.), i =1, ..., M-1) be M-1 given polynomials. Correspondingly, for each 

i=2, ...$l, denote by va the solution to the problem 

ri 

b 

Vi& = 0 inC 

ViJq = 0 in ci 

I 

Moreova, for each i=l, ... ,M-1 we let za be such that 
-22- 

0 0 

FIGS. 1 Boundary conditions for vW arid zD 

Let us denote with qi the vector of the values attained by the polynomial qi at the collocation points 

Cri of Ti. Then we define Sil, Siz, Si3. Si4 to be the square mamces of order N-1 such that 

for m=l, ..., N-1. Si1 is defined for i=l, ..., M-1, S'z for i=2 ,..., M and Si3, Sig for i=2 ,..., M-I. 

We shall call for brevity interfaceflux on ri of a function w defined on Ri the quantity 

a. 

a;; w (XL. Y,) + Lw (xh 9 Y,) ai for m=l, ..., N-1 . 

Then the mamx Sil represents the operator which associates to qi on Ti the interface flux on Ti of 

the left discrete harmonic extension of qi on Ri. Si2, Si3, Si4 can be interpreted analogously. 
-23 - 

i Si a 

f 9 

of si*, si,, si3, si.+ 

We now define the block nidiagonal matrix S, of order (M-1) (N-1). as follows. 

We claim that the matrix S in (5.7) is actually the influence (or capacitance) mamx of the system 

(2.5)-(2.9). As a matter of fact, let us solve the homogeneous collocation problems 
-24- 

far i=l, ... ,M and compute the vcctor b = [b, ,..., bi ,..., bM-I]Tdcfincd by 

Then define the vector 6 = [5~,...,6i,...,5~.1]T as the solution of the linear system - 

(5.10) S$=b . 

It is now readily seen (using (5.1)-(5.10)) that the solution of (2.5)-(2.9) is given by 

where wiN and ziJ are obtained from (5.1) and (5.2) taking qi = $. Thus, the matrix S defined 

in (5.7) is actually the influence matrix associated with the problem (2.5)-(2.9). 

It can be proved, using repeatedlythe integration by parts, that in the case of the 

Llgendre collocation the matrix S is symmetric. On the other hand S is not symmemc if the 

Chebyshev collocation points an used. S is again non symmetric in the case of both Legendre or 

Chebyshev collocation if the pure flux condition ax = -T on ri is used instead of (2.9).0 %.N 

The relations (5.8)-(5.11) suggest a possible algorithm, based on the influence mamx 

system, for the solution of the multidomain problem (2.5)-(2.9). The crucial step of this algorithm 

is the use of an efficient method for the solution of the system (5.10). To this end, preconditioned 

iterative methods (which do not require the explicit construction of the maaix S but simply the 

calculation of the product of S time a given vector) are currently used (see e.g. 131, [4], [SI, 171, 

[8],[11],[16]). In this section we show that our iteration by subdomain method amounts 

to the use of the Richardson iterative method for (5.10) preconditioned with a proper 

matrix P (see (5.15) below). 

Let us mall that if P is any non singular manix of order (M-1) (N-1). the Rirbrdson iterative 

method for solving (5.10) preconditioned with P is: 

(5.12) 

The iteration maaix at the step n is 

. 

P (r' - 53, = 8,(b - Sen) - . 

(5.13) (1 - e,,) I +e, (I - P-1 S) . 
-25 - 

Let us define the matrix P as follows (we assume, e.g., M odd). 

(5.14) P= 

S2 

s2 

2 

4 

0 

S2 

S2 

3 

1 

s4 s4 

s 4 s 4 

2 3 

4 1 

0 

sM-l sMsM-l sM2 3 

4 1 

We denote with Te and T the matrices associated With the optm defined in (4.20) and (4.21) 

andbearing the samc name. 

.. msmon 5.1 The followine identitv holds 

Etpnf In view of (4.20) and (4.21) it is sufficient to prove that T = I - P1 S or equivalently that 

PT = P-S. Recalling (5.3)-(5.6) the action of P consists of solving a sequence of Dirichlet 

problems in Rifor i even, and then to compute the interface fluxes on r. Recalling the 

interpretation of the operator T given in the remark 4.2, it is readily seen that applying the matrix 

PT to the values of a polynomial on Cr amounts to solve a sequence of Dirichlet problems in Ri 

for i odd with such boundary conditions, and then to compute the interface fluxes on r of the 

obtained solution. 

On the other hand the matrix P-S is given by 
P-s = 

1 

1 

-S ii 3 

2 

-S 

s3 

s3 

3 

3 

1 

-S 

-26- 

M O4 -S 

2 

ROm (5.3)-(5.6) we deduce that the matrix P-S is associated with the solution of Dirichlet 

problems in Qifm i odd and to the computation of the interface fluxes on r. This shows that 

PT=P-s. 0 

lbnaku rrhee of two subdomams) ' If M is even it is possible to give an analogous 

interpretation of the itcrafion-by-suMd method. In the particular case in which M=2 we have 

(5.16) S =S'1 +S22 and P=S22 , 

and-therefore (5.15) becomes T, = (1 - 0) I - 0 (S2$' S',). If R, and-%have the same measure 

and L = -A + aI, with a E R+, then S22 = Sl1, hence taking e,, = 1/2, Qn, yields exact 

convergence in two iterations (this result was already found in [ 151). 0 

.. he nrecond inoned mam 'x) Let us consider the case of the 

Legendre collocation with M=2. Recalling the definitions of S1, and SZ2 (see (5.3) (5.4)) and 

using the discrete integration by parts formula (3.8) we deduce from (5.16) that 

.. Fl.maU2 (Condltlon number of t 

where cp is the vector such that Q~ = cp(ym) m=O, ..., N,,, and 

(5.18) 

Similarly the preconditioner P = s2, satisfies 
From (5.18),(5.19) and theorem 4.1 we obtain the fundamental relation 

t 

where r is a constant independenr ofN. This ensures that the eigenvalues of the matrix P-lS are 

real and positive. Moreover the ratio Laxlhmin between the maximum and the minimum 

eigenvalue of PIS is bounded by 7.0 

RcmiKu& choice of the - en). The interpretation of the iterative 

method(2.11)-(2.18) as a Richardson scheme suggests some convenient choices of the parameter 

9, appearing in (2.18). When information is available on the exmme eigenvalues of the ponditioned influence ma& the optimal pararneter 8 = 

Otherwise a dynamical choice inspired bya"minimal residual'' strategy can be used (see [ 11 for the 

case of multidomain finite element approximation). 

Another strategy is suggested by the proof thmm 4.2. The minimal value of the conmcuon constant k(0) is achieved for 9, = (r + l/a % + .e + 2). In [21] the values of e,, are choosen so that 

can be used for all n (see e.g.[7]). 1 

1- + A& 

-2 

sup (On] = 8, and no extra computational work is needed to obtain 6,. The same strategy 

can be applied to the case of spectral collocation approximations here discussed.0 

5 OTHERPRECONDITIOh% R S 

Based on the interpretation of the preconditioned Richardson method for solving the 

influence system (5. lo), one can formulate other iteration-by-subdomain procedures by simply 

taking a preconditioner P in (5.12) different than the one in (5.14). For instance if we take 

S2 0 

(6.1) P= [ s' 

2 
-28 - 

then the ireration (5.12) yields the following iteration by subdomain pmccdurc. (We identify, for 

simplicity of notation, the vector 5" with the poho~al r E PN (T') whose values at the 

collocation points arc the components of 5"). 

Given eo, look for ulli E PN (Q,) for i=l, ...,M- I such that for any n 2 1 

LU+f inq , 

(6.2) in c'i , Up=O incand,ifi22, 

Then, for i=2, ...,M look for V"i E PN (QJ such that 

Finally cnis defined by 

(6.4) 

Once cn-' is available, one can compute Vi. i=2, ...&I by solving the M-1 independent problems 

(6.2). Then the M-1 independent systems (6.3) can be solved simultaneously.+At each step of this 

algorithm one has to solve 2(M-1) differential problems, instead of M as prescribed by the method 

(2.1 1)-(2.18). however in this case the degree of parallelism of the algorithm is M-1, instead of 

M/2. 

The choice of the preconditioner 

- 

n-1 

(=e,, V:+l + (1 - en) ti on Cri for i=l, ..., M-I 

L 
(6.5) P= 

I 

s: 

gives rise to an algorithm analogous to (6.2)-(6.4). In this case the Dirichlet problems are solved in 

Ri for i=2 ,..., M, and the mixed Dirichlet-Neumann problems in Ri for i=l, ..., M-1. 

Another prtconditiona which comsponds to an iteration-by-suMomain technique is 

(6.6) P = 

.I 

1 

S 

s2 s2 

41 

0- 

sM-l s'- 

4 1- 

The differential interpretation of the preconditioned Richardson iteration is now the following. 

First one solves M-1 independent Dirichlet problems in R,, i=2, ..., M, with 6:-1 as boundary data. 

Then computes the interface flux Qi-l on ri-, of the solution in Ri. Finally one solves sequentially 

M-1 mixed problems going from R, to RM-*. For the mixed problem in Ri one enforces the flux Q, 

in ri, while on ri-l the value of the solution coming from Ri-* is enforced as Dirichlet data. We 

note that these mixed problems arc not solvable in parallel. 

An analogous algorithm is obtained by using the upper triangular prcconditioner 

t 

I 
P= 

'2 s s2 

s3 

2 3 

2 

- 

-30- 

- 

s'-l s'-l 

SM 

2 3 

2 

The convergence analysis for these prccondiaoners is under investigation. 

Notice that in the iteration by subdomain methods we have considered so far the subdomains can 

be handled in parallel whenever the preconditioncr P has a block diagonal structure. The degree of 

parallelism is precisely the number of diagonal blocks. 

In this section we consider the case of decompositions of R with internal cross points (see 

fig.O.l b). 

For simplicity we consider the case of a square R partitioned into four squares Ri. i=1, ...., 4 (see 

fig.6.1). The general case can be studied combining the results of this case with those of the 

privious sections. 

FIG.7.1 The decomposition of R. 
t 

-31 - 

4 

We denote by @ the space of the functions which are the traces on r := u rj of the functions of 

H1o(R). The variational multidomain formulation of problem (1.1) is given by (1.5)-(1.8) 

provided u1 is identified with u5. 

The collocation points in Ri , i=l, ..., 4, are defined in section 2. We recall that here Cri is the set 

j=l 

of the internal collocation points of the interval Ti, i=1, ..., A The cross point 0 is the unique 

collocation point common to all the subdomains. For simplicity we look for a spectral solution 

which is of equal degree (N) in the four subdomains. Hence the weights defined in (2.2) and (2.4) 

are such hat aik= y , V k = 0, ...a, moreover we-have = oo. and we denote this value with 

0. 

Let us consider the following multidomain spectral approximation of problem (1.34 1.8). For 

i=l ,..., 4 we look for uiN E PN(Ri) such that for i-1, ... ,4 (we identify us& with u18) 

aUiN - - aui+l,N - = -&u* - 00 - (Lui+ls - 00 in c 'i , (74) hi hi 

4 

=-c(L"a-f)(o)o . i= 1 

The vis are the outward normal directions to ri, ordered clockwise as indicated in fig.6.1. 

The conditions (7.4) are analogous to the ones in (2.9). The "cross point condition" (7.5) gives a 

relation between the four polynomials uiH, i=l, - ,4. In the spirit of the method (2.1 1)-(2.18) it is 

possible to formulate an iterative method for solving (7.1)-(7.5). The idea is to solve, at each 

iteration, two Dirichlet problems in R, and R,, and two mixed problems in C+ and R4 as follows. 

For i=1, ..., 4 let goi E PN (Ti> be such that goi = 0 on aR and go1 = go2 = go3 = g04 in 0 we look 

for a sequence un, n 2 1, satisfying: una E PN (Qi). i=l, ..., 4 and 

inQp!@: 

(7.6) Luys = f in cff. uys = grl in C rl , u:,~ = d" in C r4 , uyeN = gy-' in 0, u:,~ = 0 in C: 
The polynomial fi arc defined recursively as follows: 

Note that problems in $2, and & are independent, while those in 

the two cross point conditions (7.14) and (7.15). 

If the Legendre collocation points arc used, the problem (7.1)-(7.5) can be written in a variational 

and S2, are coupled through 
-33 - 

i 

form which coincide with (3.4)-(3.7) provided M4, usy is identified with u,~ , and PO, 0 is 

substituted with the space 

4 

QN := (Q : v ri + R : cp E PN cg , i=l, ..., 4, cp is continuous in 0, q4 in aR n (v r.)). 

i=l 'i 

We skip the proof which is analogous to the one of proposition 3.1. With this variational 

foxmulation we can now follow the guideline of section 4 to prove the convergence, as n + 00, of 

the iterative scheme (7.6)-(7.16). In this case however the extension theorem 4.1 takes a different 

form. Rcisely (4.7) holds now with two constants o and r which depend on N as follows 

o=a0logN, r= blogN. 

A consequence of this is that the constant 8* of theorem 4.2, as well as the mr reduction factor 

k(B), far 8 E (0, e*), depend now on N. The convergence result (4.29) is still me, but K tends to 

one as N + - with a logarithmical speed. The precise form of this result and a detailed proof of it 

can be found in [25]. 

REFERENCES 

V.I.Agoshkov, Poincd-Steklov's operators and domain decompostiion methods in finite 

dimensional spaces, to appear in [26]. 

C.Bernardi, Y.Maday-Ropcrties of some weighted Sobolev spaces and applicamon to spectral 

approximations, 1987 - Submitted to SIAM J. Numer. Anal. 

PEBjorstad, 0.B.Widlund - Itcxative methods for the solution of elliptic problems on regions 

paxtitioned into subsrmcturcs. SIAM J. Numer. Anal., 23 (1986), pp.1097-1121. 

J.H.BrambIe, J.E.Pascia!., A.H.Schatz - The construction of preconditionen for elliptic 

problems by substrumring I, Math. of Comp. 47 (1986), pp.103-134. 
(71 T.F.Chon - Analysis of preconditionen for domain dwmposition, SUM J. Numn. hal. 

2 4 (1987) PP. 382 - 3 3 0 

(81 T.F.Chan, D.C.Refasco - A suntey of preconditioners for domain decomposition - Research 

Repon Yale /DCS/RR-414, (1985). 

[ 9) R.Cottle - Manifestations of the Schur complement. Linear Algebra Appl.8, (1974). 

~p.189-211. 

(IO] PlDavis, P.Rabinowitz - Methods of numerical integration. Academic Press (New York), 

(1975). 

(1 11 MDxyja - A finite elementcapacitance method far clliptic problems on regions partitioned in 

subregions. Numcr.Math. 44, (1984). pp.153-168. 

[ 121 D.Fum - Approssimazione numenca di problcmi parabolici e ipcrbolici con wtd spemali. 

Thesis, University of Pavia, 198 1. 

[ 13) D.Funaro - A rnuludornain spectral approximation of elliptic problems. Numer. Meth. for 

PDE, 2. (1986). ~p.187-205. 

[I41 D.Funaro - Domain decomposition methods for pseudo spectral approximations. (1987) 

Submitted to Numer. Math. 

(IS] D.Funaro, A.Quamroni, P.Zanolli - An iterative procedure with interface relaxation for 

damain decomposition methods, SIAM J.Numcr. Anal. (1987) To appear. 

(161 G.HGolub, D.Maycn - The use of preconditioning over irregular regions. at Sixth 

IntConf. an Camputing Methods in Applied Sciences md Engjndng, VtrsaiUer, Dec.1983. 

(17) J.M.Lacroix, RPeyret, J.P.Pulicuti - A pseudo-spectral multidomain method for the 

N8Via-Stdres quation with rpplicatims to doublediffusion problan~. 7th GAMM Cod. on 

 NU^ Mahods in Fluid MCCMCS. LaWrin-h-Nt~e Stpt 9-1 I, 1987. 

I181 JLlAons, EMagenes - -am and aoD1,Fations. 
-35 - 

Vol.I.Springer-Vcrlag, Berlin 1972. 

[ 191 Y .Maday, A.Patera, and E.Ronquist - A well-posed optimal spectral element approximation 

for the Stokes problem. NASA CR-178343,1987, submitted to SIAM J. Numer. Anal. 

(201 HMacaraeg, L.Strcett - Improvements in spectral collocation through a multiple domain 

technique. Appl. Numer. Math. To appear. 

[21] D.Ma.~ini, A.Quarteroni - A relaxation procedure for domain decomposition methods using 

finite elements. (1987), submitted to Numer. Math. 

[22] Y.Marion,. B.Gay - Resolution des equation de Navier-Stokes par une methode 

pseudo-spectrale via une technique de coordination. Sixth Int. Symp. on Finite Elements 

Methods in Flow Problems. Anabes, June 1620,1986. 

[23] S.A.Orszag - Spectral methods for problems in complex geomemes. J. Comput. Phys., 37, 

(1980), pp.70-92. 

[24] J.P.Pulicani - A spectral multimethod for the solution of ID Helmotz and Stokes type 

equation. Preprint 11.122, Universitd de Nice, (1986). To appear in Computers and Fluids. 

[25] A.Quarteroni, G.Sacchi-Landriani, Parallel Algorithms for the Capacitance Mamx Method in 

Domain Decompositions, to appear in calcolo. 

[26] SIAM Proceedings of the First Int. Symp. on Domain Decomposition Methods, Paris, 1987. 

R.Glowinski, G.H.Golub, Eds. 

[27] 0.B.Widlund - Iterative substructuring methods: algorithms and theory for elliptic problems in 

the plane. In [26]. 

[28] Zampieri. P.Zanolli - Domain decomposition preconditioners for elliptic equations using 

spectral methods. 1987, to be submitted to CALCOLO. 

I 
Report Documentation Page 

1. Report No. 2. Government Accession No. 

NASA CR-181620 

ICASE Report No. 88-11 

4. Title and Subtitle 

3. Recipient's Catalog No. 

5. Report Date 

DOMAIN DECOMPOS ITION PRECONDITIONERS FOR 

THE SPECTRAL COLLOCATION METHOD 

7. Authoris) 

Alfio Ouarteroni and Giovanni Sacchi-Land riani 

9. Performing Organization Name and Address 

Institute for Computer Applications in Science 

Mail Stop 132C, NASA Langley Research Center 

and Engineering 

HanDton, VA 33665 - 5735 

12. Sporkoring Agency Name and Address 

National Aeronautics and Space Administration 

Langley Research Center 

Hampton, VA 23665-5225 

January 1988 

8. Performing Organization Report No. 

88-1 1 

10. Work Unit No. 

505-90-21-01 

11. Contract or Grant No. 

NAS 1-1 8 1 0 7 

13. Type of Report and Period Covered 

Contractor Report 

14. Sponsoring 4gency Code 

17. Key Words (Suggested by Authork)) 

iomain decomposition, spectral 

collocation methods, parallel 

calculations, approximation of elliptic 

?quat ions 

Final Report 

16. Abstract 

We propose and analyze several block iteration preconditioners for the solution of elliptic problems by spectral collocation methods in a region partitioned 

into several rectangles. It is shown that convergence is achieved with a rate 

which does not depend on the polynomial degree of the spectral solution. The 

iterative methods here presented can be effectively implemented on multiprocessor 

systems due to their high degree of parallelism. 

18. Distribution Statement 

64 - Numerical Analysis 

Unclassified - unlimited 

19. Security Classif. (of this report) 20. Security Classif. (of this page) 21. No. of pages 

Unclassified Un c 1 as si f i ed 37 

NASA-Langley, 1988 

22. Price 

A0 3 

! 

i 

1 

I 

I 

i 

1 

I 

1'