Alternating-direction collocation for rectangular regions

C O M P U T E R M E T H O D S IN A P P L I E D M E C H A N I C S A N D E N G I N E E R I N G 27 (1981) 265-277 NORTH-HOLLAND PUBLISHING COMPANY

ALTERNATING-DIRECTION COLLOCATION FOR RECTANGULAR REGIONS* L. H A Y E S Department of Aerospace Engineering and Engineering Mechanics. The University of Texas at Austin. USA

and t3. P I N D E R and M. C E L I A Department of Civil Engineering, Princeton Unit'ersity. Princeton. NJ 0~'540. USA Manuscript received 3 March lt~X0

An alternating-direction procedure for tinite elemen! collocation equation,, can be formulated on rectangular regions. Finite element collocation o f f c ~ a great sa,,ings over /inite element techniques in the matrix formation stage of the numerical proces~ Moreo~,er. with an alternating-direction procedure, multidimensional problems can be sol~,.ed as a series of ,me-Jimcn,,ional p r o b l e m s This greatb,' reduces ~he amount of computer time and storage required to ,,oh.t. large problems

I. I n t r o d u c t i o n

The description of many physical problems involves the solution ~f parabolic boundary: value problems. Two numerical techniques ha,~e been widely used in the solution of these problems: finite differences and finite elements. Finite difference cquati,.ms can be formed ,,'cry quickly, but it is often difficult to formulate a finite difference scheme for general geometries or when flux boundary conditions are present. Finite ditterence methods produce a ,,olution only at a discrete set of grid points and not o v e r the entire regi~m. This can be a particular disadvantage in applications where it is the gradient and not the solution itself which is of interest. Finite element methods produce a global solution, they can ea,,ilv handle all types of boundary conditions, and very general geometries can be treated using i,,oparametric elements. However, finite element approximating equations involve integratl()n~, which can be time-consuming, particularly when numerical integration is involved Finite element collocation methods seem to combine the advantages of both linite differences and finite element methods. Collocation equations involve point evaluations. "l'ht ;. as in finite differences, the equations can be formed quickly. Collocation produces a global solution and boundary conditions are readily treated. The vast majority of computati()nal w:~rk with collocation methods has been done for rectangular regions: however, recentb Frind and Pinder [1] have combined the notion of isoparametric finite elements with collocation t(:) generalize the collocation procedure to general domains. *This work was supported in part by a National Science Foundation Grant CME-7915789 at The University t~I Texas.

0045-7825/81/0000-0000/$02.50 © 1981 North-Holland

266

L. Hayes et al., Alternating-direction collocation for rectangular regions

In 1971 Russell and S~ampine investigated the accuracy of the collocation method on a variety of one-d~mensional problems [2]. Over a wide range of problems, they found that the accuracy of the solutions and gradients obtained from collocation were comparable to those obtained from finite element methods using a cubic basis. Houstis, Papathedorou and Rice [3] investigated the efficiency of collocation methods for ten different linear problems in two dimensions. They generally found that collocation was faster and required less storage than finite differences for a solution accuracy of 1.5 digits or more. Pinder, Frind and Celia [1, 7] have applied the collocation method to problems with distorted grids and found that for these problems, collocation produced very good answers, and they observed a 96% decrease in the equation setup time as compared to finite element methods using a bicubic Hermite basis. Classical collocation has been widely used in chemical engineering where the equations often exhibit a strong nonlinearity [4, 5]. The numerical experience with these methods seems to indicate that they produce very good answers for the types of problems in question and they seem to be more efficient than finite element or finite difference methods. While collocation equations can be formed much faster than finite element equations, the time required to solve the resulting system of equations is essentially the same for both methods. Therefore, if the solution time dominates the equation setup time, then the gain in forming collocation equations rather than finite element equations would appear to be negligible. If collocation equations could be solved an order of magnitude faster than standard finite element equations, then collocation methods could have an obvious advantage over both finite element and finite difference methods. One possible candidate for generating an efficient equation solving scheme is the alternating-direction method. This approach, combined with collocation, offers potentially great savings in both the computer storage and computation effort required to solve multidimensional problems. For example, consider a problem with N grid points. In two dimensions, using alternating-direction methods, the solution would require O(N) operations as opposed to O(N log N) operations for standard techniques. In three dimensions the solution using alternating directions again requires O(N) operations as opposed to O(N 4/3) operations for nested dissection. Using an alternating-direction procedure, one solves a series of onedimensional problems so the amount of data required in core at any given time is determined by the one-dimensional size of the multidimensional grid, which can be relatively small. For these reasons, an alternating-direction collocation procedure would appear to ret]uire significantly less computer time and storage than other numerical methods for two- a~d three-dimensional problems. This could greatly extend the range of physical problems that can be treated using a multidimensional modes. We now investigate the development of an alternating-direction collocation formulation. Let us consider the linear parabolic boundary value problem (a)

c Ou(x, t)_ aAu(x, t)+ bVu(x, t)+ du(x, t) = F(x, t) Ot

'

(x, t)~ g~ x (0, T)

,

(b)

u(x, O)= uo(x) ,

x ~ ~,

(c)

u(x, t)= e(x, t),

(x, t)~. Og2, x (0, T),

(d)

Ou Oft (x, t)= oJ(x, t) ,

(x, t)~. Oa2 x (0, T),

(1.1)

L. Hayes et al., Alternating-direction collocation for rectangular regions

267

where the region n is the union of rectangles, aD~ is made up of edges, a£~2 = af~/aJ2~, and 77 is the unit outward normal to the boundary; c, a, b and d are constant. The method extends easily to three or more space dimensions; however, for illustrative purposes, we shall restrict our attention to a single rectangle in R e. The time stepping scheme for this method is a variation of the Crank-Nicolson method. ia section 2 we develop the method. In section 3 we discuss the implementation of the me~.hod when the boundary conditions (1.1.c) and (1.1.d) are non-zero. In section 4 we present estimates of work and storage requirements; and, finally, in section 5, we present numerical results obtained using the alternating-direction collocation method.

2. Simple model problem Consider the simplified problem

(a)

OU c~t

(b)

u(x, 0 ) = 0,

x ~ 1"/.

(c)

u(x, t) = O,

(x. t ) E Of~, × ({). T ) .

(d)

au (x. t)= o.

(x, t ) C c~f~: × (0. T ) .

A u = F ( x , t, u, V u ) ,

(x. t ) ~ f~ × (0, T ) , (2.1)

where/'~ is a rectangle, ,gf~s is made up of edges of/2. and a.f?. = Af?//~f~,. As in the standard finite method, the region 1~ is partitioned into a finite number of subregions or elements. {.fl~}~ t. It is necessary that this be a rectangular grid: however, it does not have to be uniform in either direction. Basis functions '~d,,(x)~,"~~ are con~,tructed ba,,ed on the elements, and the numerical solution U(x. t) is written as a linear combinati,m of these basis functions: N

u(x. t)--- U(x. t ) -

l,,(t)O,(x).

(2.2)

i-!

where ~ j ( t ) } ~ l are undetermined coefficients which depend on time. in the colh,cati,,n method, the basis functions usually are chosen to be C ' ( I ] ) cubic functions We require the basis functions to be tensor products: in other words

,I,,(x, y)= ra(x)#,#(y),

~.2.3)

where t h e / t h basis function corresponds to the ath basis function of x and the /3th basis function of y. A set of points {xi}~=~ is selected in .O and these points are called collocation points. It is well known that the selection of the collocation points can affect the performance of the collocation method. In the standard collocation method, the matrix problem is formed by introducing the numerical solution (2.2) into the diffe,:ential equation (2.1) and evaluating it at each of the collocation points.

268

L. Hayes et al., Alternating-direction collocation for rectangular regions

In order to obtain a matrix problem which will factor into an alternating-direction form, eq. (2.1) is perturbed by a term which is O(At3), that is, eq. (2.1) is replaced by

u , - au + (at)~

0"

40x20y

2 U.

--

F(x, t, u, Vu).

(2.4)

Eq. (2.4) is designed to give an overall accuracy of O((At) 2) in time and is different from the corresponding equation in [6] which was only O(At) in time. We shall use the following notation. Let M be a positive integer and define At = T I M and t" = nAt. In addition, for 0-- n <- M, let U"

=

d,U" = ( U " + ' - U " ) / a t ,

U(x, t"),

U "+'/2 = (U"+! + U")/2 ,

U*.+,/2-- -~(3U" - U"-') ,

(2.5)

where U "+~/2 is an average value for U at the (n + 1/2) time step, and U*+m is an extrapolated value for U at the (n + 1/2) time step. The alternating direction collocation method is defined as follows: For 1 <_ n <--M and {x~}~=~

{

d,U" - A U "+'/2 +

(At) 2 04(d,U " - d,U"-')

4

OX2Oy2

}l = F(x, t n+'/2' U*n+,/2, VU*+t/2)

(2.6)

If U" = X~=mIz,~s(x, y) the matrix problem corresponds to (2.6) has the following form K(/t " + ' - / , " ) = ~ " ,

(2.7)

with the general matrix entry, given by

K . = 4~j(x,. y,)- ~ ) a4,j(x,, y,)+

(aty 4

04O/(Xi,Yi) Ox20y z

and

q,~ = at{F(x,,

04_dtUn-I 1 t "÷''2, U*+,,~, VU*÷l/2)+ au"(x,) + (aty 4 oxeOy 2 J"

Using the tensor product form of the basis functions, the matrix K in (2.7) has a very special form. If Sj(x, y ) = F,,(x)$a(y ), and if the ith collocation point corresponds to the rth point in the x-direction and the sth point in the y-direction, i.e., x~ = (x,, y,), then Ko can be written as rt 1 2 # tt Ku = r~(x,)¢~(y~)- ½at{_r~(x,)¢~(yD + F~(x,)O~(y.)} + ~(at) F~(x,)q,~(y,)

= [r. (x,) - -~a,r-(x,)] [¢,~(y.) - ½atq,~(y,)] =//o(x,)a~(y,),

(2.8)

L. Hayes et aL, Alternating-direction collocation for rectangular regions

269

where

F"(x') =

oW=(x)l Ox 2 I

x,

and ¢,~(y,)= a2q"(Y)] OY 2 It.,'

//=(x,) = [F=(x,)- ½atF:(x,)],

An(y,) =

[q,~,(y,)- ~atg, . ' o(y,)] "

If the collocation points are numbered in a horizontal order, the matrix K factors into

K = KyK+,

(2.9)

with

K~

~, Y= ... , •

"=-

Y~

KX

Ix,, x= O

x33

]

y¥ = [

a,(y,) 0

"".

t+'+,.;+"+i

o 1 , .• Xkl,

IN,N]

,

°1

A,(yq) tr,,pl

x~=

U~(x~)r/,(x~) I H,(x~) r/,(x,) H=(x,) r/,(x,)...... rl,(x:) , % ( x 3 H,(x+,) •.• •

•

•

"

'

"

•

[R~R]

where there are p horizontal grid lines and R vertical grid lines. The matrix K, is block diagonal and each block corresponds to a one-dimensional collocation problem along horizontal collocation lines• The first block corresponds to the first line, the second block corresponds to the second line, etc. This matrix problem can be treated very efficiently since each block is an independent one-dimensional problem. In addition, since the region is a rectangle, each of the diagonal blocks of K= is identical and it corresponds exactly to the standard collocation problem along one line of collocation points. If the collocation points are numbered in a vertical order, then the matrix Ky has exactly the same fo:m as the matrix K,. Again, each of the diagonal blocks of Ky are identical..so only the matrix problem corresponding to one line of collocation points must be calculated. The matrix problem is solved in two steps: The collocation points are taken in a vertical order and one solves

(a)

K,Z = ¢".

The collocation points are then taken in a horizontal order and one solves 0')

Kx(.a"+' = p,") = Z .

Therefore the total matrix problem can be solved as a series of independent one-dimensional problems.

L. Hayes et aL, Alternating-direction collocation for rectangular regions

270

The equation with constant coefficients (1.1) can easily be treated in the same manner. The term which is added to (1.1) so that the matrix will factor is

a2 (A_~3 O4u. c 4 Ox20ye and the coefficient matrix K has the same structure as'(2.9). The entry Kij has the form At22 a

K , = cr~(x,)~,~(y,)- ½aat{r"(x,)¢,~(y,) +

r~(x,)C,~(y~)} + - T

c r~(x,)O~(y,) ,

and the individual factors have the form

l-l~(x,)=[V~cl'o(X,)At a 2 Vc

Aa(y,)= [X/c¢a(y,)

/":(x,)]

At a

'

]

2 V ~ ¢'~(y')

"

Again, in the factored form, each diagonal block of K, and Ky is identical, so only two matrices must be calculated: one f'~r a one-dimensional problem in the x-direction and one for a one-dimensional problem in the y-direction.

3. Implementation of boundary conditions We shall consider three types of boundary conditions for eq. (2.1). The extension to eq. (1.1) is trivial. (1) Zero Dirichlet or Neumann boundary conditions

u(x, t)= o,

(x, t)E 012, x (13,T),

O.u_(x, t) = o

(x, t) ~. 002 x (0, T).

brl

(2) Space variable Dirlchlet or Neumann boundary conditions

u(x, t)= e(x),

(x, t)E 0/]1 × (0, T),

8U (x, t)= w(x) Oft

(x, t)E 002 x (0, T).

(3) Space and time variable Dirichlet of Neumann boundary conditions

u(x, t)= e(x, t),

(x, t ) ~ 012, x (0, T ) ,

Ou O~l (x, t) = w(s, t) ,

(x, t) ~. 0a2 x (0, T) .

L. Hayes et ai., Alternating-direction collocation for rectangular regions

271

Table 1 Correspondence between one-dimensional and two-dimensional boundary conditioas

(A) (B) (C) (D)

Type of B.C. along vert.ical edge"

Type of B.C. along horizontal edge"

I~.C. at corner

One-dimensional basis function to be ignored at the corner

Dirichlet

u(u,)

Dirichlet

u

F(x)

u(u~)

u~, u~

4,(y)

Neumann ux (u~).) Dirichlet u(uy ) Neumann

Neumann uy(u,y ) Neumann uy(u~y ) Dirichlet

u~ u),, u~. u

u

aF(x)/c~x /~b(y)/0y F(x) ,'~(y)/c~y OF(x)/,Ox

u,(u,)

u(u~)

u.. u,,

~(y)

u . uxy

"Explicit Boundary Condition (Implicit Boundary Condition).

The strategy will be to first form the matrix problem ignoring the boundary conditions and then to apply the boundary conditions and adjust the right hand side, which in effect removes the boundary nodes from the matrix problem. Case I: Zero Boundary Conditions. Since the prescribed value at these boundary nodes is zero, there is no correction to be made to the right hand side. One simply ignores the boundary nodes and the corresponding basis function.~ in the creation of the matrix problem. Table I shows the correspondence of the one-dimensional basis lunctions to the applica).ion of the two-dimensional boundary conditions. For instance, on a rectangle, if one applies a Dirichlet boundary condition along a horizontal and a vertical edge. the tangential derivative is then fixed along these edges, and at the corner one is in fact prescribing u. u, and u,. If these boundary conditions are set to zero. this corresponds to setting l'(x)= ~,(y)= () along each edge. Consequently, the one-dimensional basis functions F(x) and ¢(y), which correspond to function values at the corner, will be ignored in the formulation of the matrix problem. The matrix problem is then solved exactly as described previously. Case 2: Space Variable Boundary Conditions. One selects U" to satisfy the boundaD' conditions. At each time step, we are solving for differences. U " " - - U " and since the boundary conditions are time independent, the difference at boundary n(~es is zero. Thus. there will not be a correction to the right hand side because of differences in bounda~ values. Therefore, as in Case 1, at each time step one simpfiy ignores the boundary n,.~es and corresponding basis functions in forming the matrices K, and K~. However, in forming the right hand side, one must recall that N

NBC

j=,

I:!

where {4,j(x)}~Y=, are basis functions corresponding te nodes without boundary conditions. NBC is the number of boundary conditions, {Bt(x)}~ c, are basis functions corresponding to n(×les with boundary conditions, and {p~}~c are the known boundary values. The right hand side, ~ " has the form

L. Hayes et al., Alternating-direction collocation for rectangular regions

272

NBC

N

O i"= A t { F ( x i , .tn+l/2, U*,,.+,:,, v u , "* , + , : , ) + (j~-i p.TAq,~(x,) + ~'~ PtdB,(x,)) l~l

+-T ~_

(~

a'q,j(x,)~l

At

/ ax'ay ~]J"

The term ~_-~ P~dB~(xi) is added to the right hand side because U" satisfies the boundary conditions. Case 3" Space and Time Variable Boundary Conditions. In this case, the boundary conditions are changing in time. As before, U" satisfies boundary conditions so N

NBC

U" - ~, p,';~j(x) + ~, PTB,(x), jffil

lffil

where P7 is the /th boundary condition at time t", and B: is a boundary basis function. In order to remove t h e / t h node from the matrix problem, we must correct the right-hand side. Let g[d,] be the operator At

~:[q,] = ¢ , ( x ) - -~ aq,(x) +

A t 2 O40(x)

4 8x20y 2"

It is easy to see that in (2.7) Kij = g[C~/(xi)]. The ith entry of • ~ is corrected LSy~ c

(p~+~_

PT)g[Ba(x,)] and d~j- = At { F(xi, t "+I/2, U.+I/2, *

v.U.+i/2) +

N

NBC

N

+-T ~mc

at p~-i

+ ~ P7 t=l At

/Ox~Oy•

a'8,(x,)lx 1 oxe~y ~-]].[

NBC

(P7 +I- P't)g[Bz(x,

)l .

/ffil

4. Work and storage requirements Consider the case where the region is a rectangle and there are R vertical lines of collocation points and p horizontal lines of collocation points. If the equation has constant coefficients then the diagonal blocks of Kx and Ky in (2.9) are identical and only one decomposition in the x-direction and one decomposition in the y-direction will be required. This requires only O(R + p) operations. This decomposition is done only once per problem. The decomposition of the standard collocation matrix requires O(N 3/2) operations in two dimensions, and requires O(N 2) operations in three dimensions. At each time step, the

L. Hayes et al., Alternating-direc.ion collocation for rectangular regions

273

solution of the alternating-direction matrix problem requires only O(N) operations; whereas in two dimensions, standard techniques require O(N log N ) operations and in three dimensions George's nested dissection requires O ( N 4/3) operations. The total solution time for alternating directions in both two and three dimensions is O(N + ,~fV) where M is the number of time steps. Standard collocation requires O(N3/2+ M N log N) in two dimensions and O(N 2 + M N 4/3) in three dimensions. The alternating-direction method deals with only one line of data at a time, so the storage requirements are quite low. This technique offers a great savings in both storage and execution time for large problems and may greatly extend the range of three-dimensional problems which can be modeled using numerical techniques.

5. Numerical results

Let us now consider the results of several numerical experiments to demonstrate some properties of the solution scheme. The first experiment addresses the question of convergence, while the second illustrates an application to a field equation. For each problem cubic Hermite basis functions are used. The model problem for the convergence study is

c3u V 2 u + u = f at .

0
.

0
t>0

(5.1)

subject to

u(x, t)= O,

x E cgl'l,

u(x, 0)= 0,

x~O.

The homo£,,-,neous boundary condition forces all tangential gradients along the boundaries, and the cross-derivatives at the corners to be zero in the Hermite approximation (as shown in table 1). The forcing function f is chosen so that the analytical solution is u(x, y . t ) = ( x - x2)(y - y2)e"t e'. The problem is solved using square elements of equal size. Convergence in both space and time is examined, and the results are summarized in fig. I. Fig. la illustrates spatial convergence; the mesh is successively refined while holding the time step constant at a very small value. The error norms displayed in fig. 1 are defined a~ follow~:

Ilell '° = fn e= da,

Ilell ' = f,,

[e 2 + (Oe/ax) 2 + (Oe/OY) 2] da

where the error e is defined as

e-u-O,

%

,:2_

'0'

¢0

'0.

q.

!0

I . . . . .

AX

0.125

! . . . . . . . . . .

0.25

0.0625

,'~ X

YI~

bJ

=

to..

od

. . . . . . . . . . . . .

1"--" . . . . . .

1.0

T . . . . . . . .

Q5

1" . . . . . . . .

0 2

~

,'~t

0 1

.....

"1"-'-

0.05

. . . . .

•

0.02

Fig. 1. (a) Convergence in space for model problem (5.1), (b) convergence in time for model oroblem (5.1).

0.5

//ellH t

(a)

0.01

AX= 1/16 T = 1.00

(b)

w

L. Hayes et al., Alternating-direction collocation for rectangular regions

275

with u and 0 being the analytical and numerical solutions, respectively. Fig. l a clearly shows fourth order convergence in the H0 norm and third order convergence in the/-/1 norm; these results are as expected, since our basis functions are Hermite cubics. , Fig. lb shows the results of the study of convergence in time. These values are computed using a very fine spatial mesh and successively decreasing time steps. The Crank-Nic~lson finite difference approximation used here is second order correct in time; fig. lb shows a convergence rate of slightly greater than two. Thus, convergence in time in this particular problem is slightly better than the theoretical minimum rate. The second problem considered is a typical equation describing the transport of a contaminant with decay,

Ou + 20u 02U 02U 0-7 ~-20-~x-~-10-~y2+0.005u=0,

0
0
t>0,

(5.2)

subject to a 'patch' source boundary condition

O<-y
b
Ou/Oy = 0 along y = O, y = L, u,

Oulax -.,0 as x -., oo,

u(x,O)= o,

x

n.

All other necessary boundary conditions are implemented as per table 1. The analytical solution to this problem is

u(x, y, t)=

x

'

1

2

4x/T6× l - erfc(:.~,/2"~) ÷ e r f c ( ~ )

} d'r.

The problem is solved using rectangular elements of unequal size (see fig. 2). The domain is truncated at x equal to 90. The boundary conditions at infinity are then applied at this finite boundary. The error measure used is the maximum norm. E - max

le j I,

l

where e, = u , - O~ and i is the node number. Note that for this problem, the necessary nodal approximations are calculated using harmonic mean values from adjacent elements. Several time steps are used and the results are given in table 2, which shows that this method behaves well for this type of problem. The results converge in time to a limiting value, which is nothing more than t h e spatial error for this grid configuration. This spatial error arises from two sources. The first is the fact that the boundary conditions at infinity are applied to a finitt.

L. Hayes et al., Alternating-direction collocation for rectangular regions

276

0

0

/

04"

/i /i /I i

0 ¢D ,,,,i m

0

~m o

/

/

01

20

:50

5

60

90

X

Fig. 2. Mesh for field problem (5.2).

Table 2 Results of second test problem

At

E = max levi

5.0 2.5 1.0 0.5 0.2

0.264 0.077 0.012 0.007 0.007

l

domain, and the second is th.~.t the finite element discretization results in a spatial error. This error could be decreased by taking a larger domain and refining the spatial grid. 6. Summary An alternating-direction technique for finite element collocation has been presented. Convergence studies in time and space on a model problem show that this method produces solutions which are comparable to standard finite element methods using cubic basis and using

L, Hayes et aL, Altematit~g-direetion collocation for rectangular regions

277

a Crank-Nicolson time-stepping scheme. Numerical results show that this method works well on a field equation which is typical of those describing transport of a contaminant with decays The numerical method discussed here has a great advantage over standard finite element and finite difference methods in terms of computation time and storage. Using this alternatingdirection collocation method, multidimensional problems can be solved as a series of onedimensional collocation problems.

References [1] E.O, Frind and G.F. Pinder, A collocation finite element method for potential problems in irregular domains, J. Numer. Meths. Engrg. 14 (1979) 681-701. [2] R.D. Russell and L.F Shampine, A collocation method for boundary value problems, Numer. Math. 19 (1972) 1-28, [3] E.N. Houstis, T.S. Papatheodorou and J.R. Rice, Algorithms for elliptic partial differential equations: metalogorithms and selection, Preliminary Report, 1974. [4] B.A. Finlayson, The Method of Weighted Residu~is and Variational Principais (Academic Press, New York, 1972). [5] J.V. Villadsen and F. Michelsen, Solution of Differential Equation Models ~y Polynomial Approximation (Prentice-Hall, Englewood Cliffs. N J, 1978j. [6] L.J. Hayes, An alternating-direction collocation method for f'~ite element approximations on rectangles. Comput. Math. Appl, 6 (1980) 45-50, [7] E.O, Frind, G.F. Pinder and M. Celia, Ground water flow simulation using collocation finite elements, m: Finite Elements in Wator Resources, Vol. !! (Pentech Press. 1978) pp. 171-186