Incremental unknowns for convection–diffusion equations

Applied Numerical North-Holland

APNUM

Mathematics

365

11 (1993) 365-383

378

Incremental unknowns for convection-diffusion

equations

Min Chen The Institute for Scientific Computing Bloomington, IN 47405, USA

and Applied Mathematics,

Indiana

Unillersity,

Roger Temam The Institute for Scientific Computing and Applied Mathematics, Indiana University, Bloomington, IN 47405, USA; and Lahoratoire d’ilnalyse Nume’rique, UnirlersitP Paris-&d, Bat. 425, 91405 Orsay Cedex, France

Abstract Chen, M. and R. Temam, Incremental Mathematics 11 (1993) 365-383.

unknowns

for convection-diffusion

equations,

Applied

Numerical

When iterative methods (e.g. MR, GCR, Orthomin(kI, see [3]) arc used to approximate the solution of a nonsymmetric linear system AU = b, where A is an N X N matrix with positive symmetric part, i.e. M = (A + A’)/2 is positive-definite, the convergence rates depend on the number v(A) = A,,,(A’A)‘/*/h,,,(M) which is the condition number of A when A is symmetric positive-definite. In this paper, we use incremental unknowns (IU) see [1,2,5]) in conjunction with the above iterative methods to approximate the solution of the nonsymmetric linear system generated from the discretization of a convection-diffusion equation. We show theoretically that, with the use of IU, I, is of order O((log(h))‘) instead of 0(1/h”) which is the order of Y with the USC of the usual nodal unknowns, where h is the mesh size for the finite differences. With the use of (2.3) in Theorem 2.1, we can show that af most O((log(h)J4) iterations arc needed to attain an acceptable solution. Numerical results are also presented. They show that actually, O( Ilog(h)l) iterations are needed with the use of IU and 0(1/h) iterations are needed with the USC of nodal unknowns to obtain the approximate solution. The algorithms arc efficient and easy to implement.

1. Introduction For a class of convection-diffusion equations, we use finite differences discretization. The problem of finding the solutions of these equations becomes approximating the solutions of large sparse systems of linear equations

for the spatial the problem of

AU=b,

(1.1) where U,b E RN and A is a nonsymmetric matrix of order N. The dimension N is the number of nodal points and U is the vector corresponding to the nodal values of the unknown function. Correspondence Cedex, France.

to: R. Temam, Laboratoire d’Analyse Numerique, Telephone: (33-l) 69417177. Fax: (33-l) 69416718.

0168-9274/93/$06.00

0 1993 - Elsevier

Science

Publishers

Universite

Paris-Sud,

B.V. All rights reserved

Bat. 425, 91405 Orsay

M. Chen, R. Temam / Incremental unknowns for comection-diffusion

366

equations

Under some assumptions, the discretization matrix A is regular and has positive symmetric part, namely M = +(At +A) is positive-definite. In order to approximate the solution of (l.l), several iterative methods (e.g. Minimum Residual (MR), Generalized Conjugate Residual (GCR), and Orthomin(k) methods) can be used (see [3]). These methods are in particular suitable for the above linear system since all references to A are in the form of a matrix-vector product Au. There is no need to store the matrix A. The convergence rates of these methods depend on v(A) (cf. [3]). When these methods are used directly to solve the system (l.l), v(A) is of order 0(1/h*), which is the order of the condition number of A when A is the discretized matrix of the Laplace operator. Numerical experiments are used to examine the convergence rate. Results presented in Section 6 show that in order to reduce the error, which is the difference between the computational solution and the exact solution, by a certain factor, at least 0(1/h) iterations are needed. In this paper, we introduce incremental unknowns (IU) as we did for the symmetric linear systems (see [1,2]). If we let U E [WN be the IU and let S be the transfer matrix from 0 to U, U=KJ,

(I.21

then (1.1) becomes -AU=&,

(1.3) with A= StAS and b = S’b. We show theoretically that v(x) is of order O((log(h))‘) (which is smaller than v(A)). Furthermore, with the use of (2.3) in Section 2 (which is not optimal), we can show that at most O((log(h))4) iterations are needed for reducing the error by a certain factor if we perform the above iteration methods on (1.3). This number is not optimal either. Our numerical results show that only 0( ]log(h) I> iterations (which is smaller than 0(1/h), when h is small enough) are needed. Due to the construction of S (see [l]), the computational flops needed for the matrix-vector product Au is 0(1/h) which is the order of the flops required for the matrix-vector product Au. The computation time needed to perform one iteration on (1.3) is comparable to that on (1.1). Therefore, the total CPU time needed to approximate the solution of (1.1) can be reduced by using incremental unknowns in conjunction with iterative methods, i.e. performing these methods on (1.3). Our numerical results agree with these results. Furthermore, our algorithm is easy to implement.

2. MR, GCR, GCR( k), and Orthomint k) methods For approximating

the solution

of a linear system

Ax=b,

(2.1)

where A is a large sparse nonsymmetric matrix of order N with positive-definite symmetric part, we consider the following four algorithms: MR, GCR, GCR(k), and Orthomin(k). All these algorithms have the following general form: Choose x(0). Compute r(0) = b -Ax(O). Set p(0) = r(O).

M. Chen, R. Temam / Incremental unknowns for conuection-diffusion

For i = 0, 1,. . . , compute: a(i) = (r(i), Ap(i)l/(Ap(i), x(i + 1) =x(i) + a(i)p(i), r(i + 1) = r(i) - a(i)Ap(i), compute p(i + 11.

367

equutions

A&)),

The four algorithms differ in the way they compute p(i + 1): (9 MR algorithm: p(i + 1) = r(i + 1).

This algorithm resembles the steepest descent method in the symmetric case. (ii) GCR algorithm: p(i + 1) = r(i + 1) + h bj”p(j), j=O

where bj’) = -(Ar(i

+ 11, Apt j)l/(Ap(j),

Ap( j)>,

j < i.

(2.2)

x(i + 11 minimizes

the residual IIr(i + 1) II2 over the space n(O) + Therefore, with the absence of computer error, we can obtain the exact solution of (2.1) in N steps. (iii) GCR(k) algorithm: The algorithm restarts the GCR algorithm every k + 1 iterations (i.e. after computing x(i), we set x(O) =x(i) and start GCR using this x(0) as the initial guess). This algorithm requires less storage space than the GCR algorithm. (iv) Orthomin(k) algorithm: In this algorithm

span1 p(O), . . . ,p(i)}.

p(i + 1) = r(i + 1) +

bj’)p( j).

i j=i-k+

1

In this algorithm x(i + 1) minimizes , . . ,p(i)} for i 2 k.

IIr(i + 1) II2 over the space x(i - k) + span(p(i -k),

These four algorithms have some similarities. They are all equivalent to the conjugate residual method, which is a variant of the conjugate gradient method, when A is symmetric and positive-definite. We observe that the MR algorithm is a special case of the GCR(k) and Orthomimk) algorithms with k = 0. For the convergence of these methods, we have the following result which is proved in [3]: Theorem 2.1. If {r(i)} is the sequence of residuals generated by the MR, GCR, GCR(k) Orthomin( k) algorithms, then

II r(i) II 2 G

{1 -

~)“‘ll~(~)ll2,

where v(A) = A,,(A’A)‘/2/Amin(M) and A4 = +(A +A’) is the symmetric part ofA. Hence in order to reduce the residual by a factor 6, O( I log(S) I) operations are needed.

or

(2.3)

M. Chen, R. Temam / Incremental unknowns for corwection-diffusion equations

368

3. Properties

of incremental

We consider

a Dirichlet

problem

in 0;

u=O

Lu=f

With finite differences, L&l

unknowns

0naR.

the discretized

(3.1)

equation

can be written

as

=fh>

(3.2)

where fh and uh are the approximate values of f and u at grid points respectively. Incremental unknowns can be defined when multi-level discretizations are used. When two levels of discretization are used and R is the square (0, 1>2, we let h = 1/(2N), NE N, and 2h be the fine mesh size and coarse mesh size, respectively. If uj,j = u(ih, jh> is the approximate value of the unknown function u at the mesh point (ih, jh), then the incremental unknowns consist of the nodal values y2i,2j = u2i,2j at the points of the coarse grid (2ih, 2jh), i,j = 0,. . . , N, and of appropriate incremental quantities z,,~ at the other points: ‘2i,2j+

i=l

1 =U2i,2j+l_3(‘2i,2j+U*i,2,+2),

=2i+

1,zj =~2i+l,2j-i(~2i,2j+~2i+2,2j),

X2i+

1,2j+ 1 =

‘2i+l,Zj+l

-

y...>N-l,

i=O

i(“2i,2j

+

‘21+2,2j

+

‘2i,2;:2

j=o,***>N-l;

..,N-l, +

i,j=O ,..., Readers are referred to [5] for the initial motivations the motivations for linear equations. We now consider d + 1 levels of meshes in [w2: h0 = (hL”Y

h2,lA

h, =

h2,,),

(h,,,>

hi,r =

hi,o/‘2’,

j=l,...,N-1;

‘2i+2,2j+2),

N-l.

for nonlinear

equations

and to [1,2] for

O
Hence h, is the coarsest mesh and h, the finest mesh. Different meshes are allowed in both directions x1 and x2. To these meshes we associate the grids 3, =LZ?~) consisting of points (j,h, I, j2h2,J, where jl,j2 E Z. We denote by %‘[ the set of nodal points %[ =9?[ n a. For j = CJ,, j,> E .Z2 we denote also by ~j,I the rectangle q,/

= (j1h,,lY (jl + I)h,,J

X (j2h2,17 (j2 + IP2,$

and by Y1 the set of all rectangles Xj,[, 0 < 1~ d. For simplicity we shall emphasize the case where R is a rectangle (0, a,) X (0, a2> and hi,o = ai/No, i = 1,2, N, E N. More general open sets 0 and more general meshes can be considered similarly. Although we are interested in finite differences, a space of finite element functions on R will be used. More precisely for 0 G 1 G d we denote by V, the space of continuous real functions on n that are Q, (i.e. affine with respect to x, and x2 respectively) on each rectangle q.,[ c 0. Since the rectangles 3 E 7, are obtained by dividing the rectangles of 7/_, into four equal rectangles, we observe that V, c V, c . . . c V,. W, of V,_, in V, so that VI = V,_, @ I+‘,, Since V,_ 1 c V,, it is useful to define a supplement 1 G 1 G d. By reiteration we then have

M. Chen, R. Ternam / Incremental unknowns for conuection-diffusion

We can also define the decomposition

of u E V, corresponding

equations

369

to (3.3); it reads:

d u =u,+

c

Ll/,

U[E

U,El/o,

w,.

(3.4)

I=1

A detailed analysis shows that the u[(x), x E ZYl\Z!_l, I= 1,. . . , d, are exactly the incremental values of u at different levels (see [1,2]). We then endow I$ with a semi-norm corresponding to the incremental values: [u]‘,=

;

c

I=1

XErz/\Z,_,

Iu[(q”.

To a function u E V, we associate the step function fi such that G(x)

=U(jlhl,dT

j2h2,d),

vx

Ezj,dT

jJ, 0 < ji
j = (jl,

region ‘Rh*,,where

Lemma 3.1. For every u E V,, u = u0 + C;‘_lul as in (3.4), there exist constants cl, c2, c3, and c4 depending only on 8 = h,,,/h,,,, such that

and therefore

Proof. Formula (3.5) and (3.6) are proved as [2, Lemma 4.31 and [2, Lemma 4.41. Formulas (3.7) can be obtained by simply combining (3.5) and (3.6). 0 We then define the finite difference operators V,,h,= (Vl,hd, Vz2,Jt: %J$$(~> = +--{+(”

+

hi,dei)

-

+b>},

i =

1,2,

where e, = (1, 0) and e2 = (0, 1). We observe again that V1,h,fi is defined in the extended region q, = (0, a,) x (0, a2 + h,,,) and V2,hdGis defined in the extended region a&d = (0, al + hl,hd) X (0, a2). Now we assume that a multi-level discretization is used with meshes h,, . . . , h, as before, we then have several systems (3.2), L,&,,

= j-/,,,

0 < 1
In order to make the notations simple, we also write: L,U,=f,,

O
(3.8)

M. Chen, R. Temam / Incremental unknowns for comection-diffusion

370

At each level I, the set of incremental

unknowns

0, consists

equations

of the the following:

‘-V Y’ U, = 2’ i i’

z’=

:

)

\ -& )

where l Y’ = Y, = the properly ordered set of the approximate nodal values of u at the coarsest grid; l 2, = the properly ordered set of the incremental unknowns at the level j. Hence using the notations above, if u E V, is the approximate function, then Y’ = YO consists of the values of U& x), x E ZO, while Zj consists of the values of uj(x), for x E Zj\ Zj_ r. We are in fact interested in the finest mesh, i.e. when 1= d. For system (3.81, we can pass from U, to u, by using a transfer matrix S,. Hence we find as in (1.2) and (1.31, u, = s,u,, _LJJd =fd,

(3.9)

with z, = SiL,S, and fd = Stfd. Besides associating a step function G to a function u E V, as above, we can also associate a step function ii to a vector U, in the following way. Since U, can be considered as a vector consisting of the nodal values at the grid points, there is a unique function u E V, which takes the nodal values from U,. We then associate a step function fi with this u. Proposition

3.2. There exists a constant c5 depending only on 6’ = h,,,/h,,,, IIs, II2

(3.10)

G c5 * zd.

Proof. From the definition

of matrix norm, (sdu,)

II sd

such that

11; =

cud,

‘dT7,>

sup c!?~#O

(u,,

Od)

=

where ( * , . ) denotes the usual Euclidean to ud, we use (3.7) and obtain (U,, U,) = h

:, l,d

/ 2,d

;;:o

(u,, ud)

scalar product.

lfil* dx<

’

Let ii be the step function

(I~~l~~~~~+h~,~h2,~[~1~)~

h ‘; 1,d

ah*,

ud)

2,d

With the help of (3.5), we get

I uo I&2)

G

I fi, I * dx c:/ ah*,

and

[u]“, = (.zd, zd>.

G

+y0h2,&‘o, Yo),

associated

M. Chen, R. Temam / Incremental unknowns for convection-diffusion

371

equations

llSdh/2d * 1.00E+Oa

-

9.00E-01

_

8.00E-01 7.00E-01

-

6.OOE-01 5.00t-01

-

4.OOE-01

-

*p*

XOOE-01 Z.OOE-01

-

l.OOE-01

-

o.ooE+OO

Fig. 1. (1S, 112/zd

*d

0.

1.

d is plotted

against

2.

3.

4.

5.

6.

7.

8.

here with 0 = (0, 1)2. This figure confirms

the result in Proposition

3.2.

Therefore (U,, u,> < &

* 4d{c;h,,,h2,0(Y0,

l,o

yo> + h,,,h,,,(E

2”)).

2,0

By taking c5 = (c,’ max(cz, 1)>‘/2, we obtain cud, &>
ad),

and (3.10) follows (see also Fig. 1). 4. Incremental

unknowns

0

for convection-diffusion

Consider the convection-diffusion Lu=

equations

equation

-(U,,+Uyy)+bUx=f,

b>O

(4.1)

on the square (0, a,> X (0, u2) with homogeneous Dirichlet boundary conditions. five-point discretization scheme with mesh size h,: ( Ld”d)O

=

&%,p

-

Ua-1,p

-

4+1.p>

l,d

+ =

u a,P = 0,

b(* 2h,,,

a+l,P

-

Ua-l,6

f a,p’ cu,p=1,...,

if cy or p = 0 or 2dNo.

+

&a$ 2,d

-

UcrJ-l

-

We use the

%p+1>

>

2dNo-l,

(4.2) (4.3)

M. Chen, R. Temam / Incremental unknowns for convection-diffusion

372

equations

are the approximate values of S and u at (ah, @z) respectively. Let A,, Here fa,@ andua,@

Bi,d, Ci,d, i = 1,2, the matrix forms of the difference

(AdUdLP= ~~(ZU,$

operators, be defined by:

- Ua_r$ - uol+r,p) + >(2%,,

(Bl,dUd)a,P

=

(Bz,dU,)(%

ua+l,p P)

(c,,dud)a,,

=

(c2,dud)(a,

=

-

-

Ua,P+l)’

Ua,i3~

Ua,p+~

‘a,/3 P>

- u,,p-I

24

l,d

ua,p,

Un-l,P’

-

=

-

‘a,@

ucx,p-l,

-

for cr,p = 1,. . . , 2dN - 1. NOW we can write (4.2) and (4.3) in matrix form: 1

LdUd =

AdUd

$-cBl,d l,d

+

h,,dh2,d

When b f 0, L, calculation:

is nonsymmetric.

(LdUd, u,> = h

(4.4)

=fd*

Using boundary

condition

(4.3) we obtain

by simple

(4.5)

2,d

Hence L, is positive-definite @(Ed +z;)ud,

‘I,,)

(A&&, ud>*

:, 1,d

+

because A, is. Since

od) = (Ed&, &) = (LdUd, u,),

(4.6)

where Ed and ud are as in (3.9), Ed has a positive symmetric part. Therefore using the iterative methods in Section 2 on equation -LdUd =fd, where fd = Stfd, the methods are convergent and (2.3) holds. Theorem 4.1. There exist constants c6, c7, and c8 which depend on L,, but not on h,, such that ‘6 Amin

(4.7)

a hl,dh2,d(d

hmax( z:zd)1’2

+

’

+ b)

4 ’

‘1’

h,

dh2 ,

d 3

(4.8)

’

and hence +d)

=Gc,(l + b)(d + l)‘,

where M = i(zd + zd). Therefore in order to reduce the residual of the approximate solution by a factor E, we need at most k iterations, where k = 2c;(l

+ b)‘(d

+ 1)” Ilog

I.

(4.9)

M. Chen, R. Temam / Incremental unknowns for convection-diffusion

Proof. In [2], we have shown that there h,, such that for any I!?~

exist two constants

equations

373

c; and c; which do not depend

on

(4.10) and also

(AdUd, u,>+ (U,,

4 (d+ 1)2 G where

(&,

il is the step function

a

hl,dh2,dV

With cg = c;, we obtain Since

II2

to U,.Using (4.51, (4.61, and (4.10), we obtain

*

the first inequality

(4.7).

we can find the upper bound of A,,(~,~,)‘/2 of the above inequality. Using (4.10) we find,

With the help of Proposition

(4.11)


CT)

associated

+

u,> ,

by b ounding

each term on the right-hand

3.2,

Because @%,dsdu,, iI Bl,dSd

b2 =

Bl,dSd’d)

( ‘l,dUd,

SUP U&O

(q,

ud)

=

;;:o

‘l,dUd)

(od,

od)

side

374

M. Chen, R. Ternam / Incremental unknowns for convection-diffusion

equations

we obtain, by using (4.10) again, IIB,,,S, II2 G

(c#“,

11El,,

112 =G (c;y2cs

* 2d.

Similarly, we can obtain IIc&

< (c;)1’2cs * zd.

Therefore

A,,,( L;z,y2< h

“1 +g 1,d

2,d

(c;)1/2cs

1,d

1 c;

’

+ bh,,,(

c;y2c,).

hl,dh2,d

Now with c7 = max{c;, h2,0(c~)‘~2cs}, inequality (4.8) follows, and therefore

with cs = c7/c6. Now using (2.3), we see that for reducing the residual of the solution by a factor F, we need at most k steps, where k satisfies: k/2

1

1-p

<&. ‘(‘d)’

i

I

Then with (4.7) and (4.81, it is sufficient to have k/2

1

l-

c&I + 1)4(1 + b)2 I

i

<&,

or equivalently, 2 Ilog(E) I

ka

1 log lI

i

Since (log(1 -n) I ax,

Ci(d + 1)4(1 + h)2 iI. it is sufficient to have

k>2c,2(d+1)“(1+b2)llog(&)I.

KI

4.1 that O((log( I h, 1>)4> iterations are needed for the solution to be acceptable, namely for the residual of the solution to be smaller than a certain number, where I h, I =
Remark 4.2. We can see from Theorem

Remark 4.3. Since (2.3) is not optimal, (4.9) is not optimal either. Our numerical tests show that

O(log( I h, I)) instead of O((log( I h, l>j4) iterations are needed for an acceptable solution, in the case of the examples we have tested.

M. Chen, R. Temam / Incremental unknowns for convection-diffusion

375

equations

b = 1, h,,, = h,,, usually makes the constant c8 small. But for b sufficiently have to be chosen properly (depending on b) to make c&b t l)* small. large, h,,, and h,,,

Remark 4.4. When

5. More general elliptic operators

More general elliptic Dirichlet problem

operators

can be handled

in the same way. For instance

consider

the

P-1) u=O

ondo,

a,,, aI2 = LZ*~,u2*, b,, b,, and c are given in B?fi>

where 0 = (0, a,) x (0, u2). The functions and satisfy:

V[

=

(tl,

C$[f<

t

i=l

i,j=l

t2)

E I$*,

uij(x)[i[j

(5.2)

i=l

for x E 0 a.e., and with 0 < (y GE,

-ii~~~b,(r)+C(x)rO;

c(x)~y’O.

(5.3)

I

Using the usual five-point &fUdLP

-

= i =

U a,P =

0,

finite difference 5

0, &zijV&~)

scheme, + +t

’

i,j=l

i=l

we find

hi@ )hd + q 9/,d)u + cu I a#

(54

f ff>P’

if (Y or /3 = 0 or 2dN,

where

We then have similar results

as in Theorem

4.1.

Theorem 5.1. There exist constants cg, c,~, and cl1 which depend on L,, g, Z, b,, b,, and c(x),

but not on h, and d, such that Amin

c9

a h,,dh2,d(d

and hence $-,)

< cll(d + I)‘,

(5.5) +

‘1”

’

M. Chen, R. Temam / Incremental unknowns for convection-diffusion

376

equations

where M = i( z, + Et,). And therefore in order to reduce the residual of the approximate solution by a factor E, we need at most k iterations, where k = 2c;,(d + I)” [log(&) I.

(5 J)

Proof. As in [2], the quadratic form corresponding

+

From the assumption

to (5.4) is

d

n C’(X) - +

2 (y,h,+ c,,,)6i(x>lS2dx

i 1

i=l

*

P-8)

(5.3), we know that for h small

C”(x) - + z? (ol,,Z<, + &J&(x))

2 0,

i=l

Using (5.2) we obtain that dx>O

for U,#O.

Therefore, using (4.61, we see that $(z, + z:,) is positive-definite. GCR, GCR( k), or Orthomin( k) method on equation

Then we can use the MR,

--

LdUd=L and we know that these methods converge and (2.3) holds. Now, using (4.10), we see that

With cg = c;g, (5.5) holds. On the other hand, we assume that Ic(x) I G c;,

lb,(x)1 =Gci,

Let Md and Nd be such that

i= 1,2.

(5.9)

M. Chen, R. Temam / Incremental unknowns for convection-diffusion

equations

377

then

= IIEd II2GII@i

~,,,( wd)1’2

II2+ II@d112.

We again can find the upper bound of A,,,(z>z,)‘/2 by bounding IIad II2 and II#d II2. From the definition of Md, it is clear that Md is symmetric and positive-definite. Therefore, using (5.8), (5.9), and (5.2) and then using (4.11) we obtain --

(M,u,, Ud)

IIii& II2= sup

ii&O

(L

0,)

1 G

Qih2,d

max( Cu, c;)c; .

Because

1 t

=

hI,dh2,di=l

-Di(Bi,d l 2hi,d

+

ci,d)ud,

‘d

) )

where A,, Bi,d, Ci,d, i = 1,2, are defined as in Section 4 and the Di are diagonal i = 1,2. It is easy to see by using (5.9) that Di = diag(bi(x)),

i= 1,2.

II Di II 2a;, Now,

As in Section

4, we can prove that

II Bi,dSd11 2 < Therefore,

(C;)1’2yII Ci,dsdII 2< (C;)“2>

using Proposition

3.2, we get

2c,c;(c;)1’2 IIN,/?

G

min(hI,O,

h2,“).

1 hl,dh2,d

’

i= 1,2.

matrices,

378

M. Chen, R. Temam / Incremental unknowns for convection-diffusion

equations

Hence, 2c5c;(c;)“’ c; max(Z,

c;) + mw,,m

i

hz,,)

I

.

Choosing cl0 properly, we obtain (5.6). Choosing c1 I = cJc9, we find

V(Ld)< C&f +

1)2.

Exactly as in the proof of Theorem reduce the residual of the approximate

4.1, we can prove that we need at most k iterations q solution by a factor I, where k is as in (5.7).

to

Example 5.1. Lu = -(u,, u=O

+ uyy) + bu, + au =f

in 0,

onafl,

satisfy conditions where a and b are constants, a 2 0. With a > 0 the coefficients equation discussed in Section (5.3). With a = 0, it becomes the convection-diffusion

(5.2) and 4.

Example 5.2. Lu = -u,,

+ u, - (1 +y)(u,,

+ uy) =f

in a,

on NJ.

u=O

It is easy to verify that the coefficients

of this problem

satisfy conditions

(5.2) and (5.3).

Example 5.3. Lu = -(b(x, Y)

=g(x, u=O

y)uJx

- (c(x,

Y)UJy

+ (d(x,

Y>U)x + (e(-G

Y)U)Y +f(x,

Y)

in fi,

ondo,

where b(x,

y) = ePXy,

e(x,

Y> = Y(X +Y),

c(x,

Y> = P(x

f(x,

d(x,

+Y + 1)

1 Y) = (1 +xy)

y) = P(.x +Y)eeXY,

*

It is again easy to verify that the coefficients satisfy conditions (5.2) and (5.3) when e(r + 1) - p > 0. This is a standard elliptic test problem. Controlling p and y changes the degree of nonsymmetry.

6. Numerical

results

The comparison the MR algorithm

between the efficiencies of the MR algorithm and the MRIU algorithm (i.e. on (1.3)) is done numerically. We compare first the numbers of iterations

379

M. Chen, R. Temam / Incremental unknowns for conrlection-diffusion equations Table 1 Comparison d 2d MR MRIU

of the numbers

of iterations

for equation

(4.1) with h = 1

3 8

4 16

5 32

6 64

6 128

8 256

9 512

28 15

60 17

120 20

245 24

495 27

995 30

1993 33

-d -l/

lh,l

-O(l,‘lh,l) -O(d)

needed by using both algorithms to approximate the solution. We then compare the CPU time used to get the approximate solution. The results are very clear. With the same requirement of accuracy, by using the MRIU algorithm, fewer iterations are needed and less CPU time is used to obtain the approximate solution than by using the MR algorithm. For our numerical test, 0 = (0, l>* and h, o = h,,, = 1 are taken. We choose ug = 0 as our initial guess. The discrete values at nodal points of u(x, y) = 10 sin(xy*(l -x)(1 - y>‘12) are taken as the exact solution. The right-hand side fd is obtained by computing fd = AdUd, with U, being the exact solution.

Test 1. The convection-diffusion which is the difference between 10p6, the numbers of iterations column shows the approximate

equation (4.1) with b = 1 is solved. For the relative error, the exact solution and the approximate solution, to be less than needed by both algorithms are compared in Table 1. The last relations with respect to the grid level d.

(Ermr(t)(/(Error(O)(

1.OOE-07 1 0.00

(Log1 0 scale)

1co.00

200.00

300.00

400.00

9 500.00 No. of Cycles

Fig. 2. Relative Z* norm of the error is plotted against the iteration numbers. For problem Lu = -(u,, + ales)+ u, = f, the efficiencies of the MRIU algorithm (line 1) and the MR algorithm (line 2) arc compared with the initial error, Error(O), being the value of 10 sin((xy(1 - x)(1 - y>>‘/*) at grid points. Here WC use the uniform 128 x 128 mesh.

M. Chen, R. Temam / Incremental unknowns for conoection-diffision

380

(Ermr(t)l/lError(O)l

equations

(Log1 0 stole)

l.ooE+oa f

1 .OOE-01

I .oaE-02 1.OOE-0.J

l.ooE-04

l.oclE-05

l.ooE-06

0.00

9.95

19.90

2985

39.80

49.75

CPU (second)

Fig. 3. Relative 1’ norm of the error is plotted against the CPU time used. For problem Lu = -(u,, + u,,)+ u, = f, the efficiencies of the MRIU algorithm (line 1) and the MR algorithm (line 2) are compared with the initial error, Error(O), being the value of 10 sin((xy(l - x)(1 - y))1/21 at grid points. Here we use the uniform 128 x 128 mesh.

Based on Table 1, we can see that the number of iterations needed by using MRIU is of order O(d) = O(log 1h, I> which is less than O(l/ 1h, I>, the number of iterations needed by using the MR algorithm. For d = 9, only 33 iterations are needed by using MRIU instead of 1993 iterations required by using MR algorithm. The histories of error decrease against the number of iterations for both algorithms are plotted in Fig. 2 with 128 X 128 meshes. The advantage of MRIU can be clearly seen. We then plot the error decreases with respect to the CPU time used on Fig. 3. By drawing a horizontal line at lo-“, the line intersects the two curves at points with coordinates approximately (10, 10p5) and (3.5, 10p5). This means that 10 seconds are needed for the error to be less 10m5 by using the MRIU method instead of 35 seconds needed by using the MR method. The difference gets bigger when the number of grids gets larger. Therefore, the MRIU algorithm is very efficient and can save CPU time when we compare it with MR algorithm. Test

2. This test is still dealing with the convection-diffusion equation (4.1). The only difference with Test 1 is that here we take b = 10 instead of b = 1. For the error to be less than 10p6, the numbers of iterations needed by both algorithms are given in Table 2.

Table 2 Comparison d

2d MR MRIU

of the numbers

of iterations

for equation

(4.1) with b = 10

3 8

4 16

5 32

6 64

7 128

8 256

9 512

47 33

85 42

162 48

358 62

1053 77

2931 101

> 4000 109

-d -l/lb,/ -O(l/ lh,l> - O(d)

M. Chen, R. Temam / incremental unknowns for comection-diffusion

lError(t)l/lError(a)l

l.ooE-04

(Log10

equations

381

scale)

?

’

0.00

2ooJ30

u)o.oa

600.00

aoom

1om.M) No. of Cycles

Fig. 4. Relative I* norm of the error is plotted against the iteration numbers. For problem ,!,I* = -(II,, + IA,,,)+ lOu, = f, the efficiencies of the MRIU algorithm (line 1) and the MR algorithm (line 2) are compared with the initial error, Error(O), being the value of 10 sin((xy(lx)(1 - y))“‘) at grid points. Here we use the uniform 128 X 128 mesh.

IError(t)l/lError(O)I ) 1.ooE+oa

-

l.ooE-05

-

(Lag10

scale)

9

I.OOE-OB 0.00

20.60

41.21

61.61

62.41

103.02

CPU (second)

Fig. 5. Relative I* norm of the error is plotted against the CPU time used. For problem LA = -(u,, + Use)+ lOu, = f, the efficiencies of the MRIU algorithm (line 1) and the MR algorithm (line 2) are compared with the initial error, Error(O), being the value of 10 sin((xy(1 - xX1 - y))r/‘) at grid points. Here we use the uniform 128X 128 mesh.

M. Chen, R. Ternam / Incremental unknowns for convection-diffusion

382 Table 3 Comparison

of the numbers

of iterations

for Example

5.1 with a = 1 and b = 1

d

3

4

5

6

7

8

9

2d

8

16

32

64

128

256

512

24 16

61 17

108 20

218 24

437 27

875 30

1745 33

MR MRIU

equations

-d -l/l&l -o(l/ lh,l) - O(d)

The, conclusions for Test 1 are valid here. Figures 4 and 5 shows the error decrease with respect to the number of iterations and the CPU time used with 128 X 128 meshes respectively. We can see, in addition, from Fig. 5 that the convergence rate decreases when the error gets smaller when the MR algorithm is used. We think that this is due to the instability of the scheme and the accumulation of the computer error. The MRIU algorithm is better in this aspect too (cf. Fig. 5).

Test 3. This time we take Example 5.1 with a = 1 and b = 1 as our test problem. Again for the error to be less than 10p6, the numbers of iterations needed by both algorithms are given in Table 3. The advantages of the MRIU algorithm can be easily seen as in Test 1. Figures 6 and 7 show the error decrease against number of iterations and the CPU time used.

IErmr(t)l/IErmr(O)(

(Log10 scale)

1.OaE+oa T

3.96f-02 l.saE-03

l.oOC-07 3.98E-Ohl

.

2 l.xIC-10 6JlE-12 291E-13

1 l.OOE-14

P

’ 0.00

1801x)

360.00

540.00

720.00

wO.00 No.

al Cycles

Fig. 6. Relative 1* norm of the error is plotted against the iteration numbers. For problem Lu = -(u,, + u,,)+ uY + u = f, the efficiencies of the MRIU algorithm (line 1) and the MR algorithm (line 2) are compared with the initial error, Error(O), being the value of 10 sin((xy(lx)(1 - y>)‘/*) at grid points. Here we use the uniform 128 X 128 mesh.

M. Chen, R. Temam / Incremental unknowns for conuection-diffusion iError(t)l/lError(O)I

(Log10

equnrions

383

scale)

+ 1.ooE+ocl

2.51E-06 1.Oot-07

&3lE-12

.

2.!ilE-13

B

l.ooE-14 0.00

18.77

37.53

5630

75.07

93.83

CPU (second)

Fig. 7. Relative 1’ norm of the error is plotted against the CPU time used. For problem Lu = -(u,, + uY,,) + uY + u = f, the efficients of the MRIU algorithm (line 1) and the MR algorithm (line 2) are compared with the initial error, Error(O), being the value of 10 sin(xy(1 - xx1 - y>)‘/‘> at grid points. Here we use the uniform 128 x 128 mesh.

Acknowledgement

This work was supported Indiana University.

in part by NSF grant DMS-8802596 and by the Research Fund of

References [l] M. Chen and R. Temam, Incremental unknowns for solving partial differential equations, Numer. Math. 59 (1991) 255-271. 121M. Chen and R. Temam, Incremental unknowns in finite differences: condition number of the matrix, SIAMJ. Matrix Anal. Appl. 14 (2) (1993). iterative methods for nonsymmetric system of linear [31 S.C. Eisenstat, H.C. Elman and M.H. Schultz, Variational equations, SZAM J. Numer. Anal. 20 (1983) 345-357. 141 G.H. Golub and C.F. Van Loan, Matrix Computations (Johns Hopkins, University Press, Baltimore, MD, 1989). 151 R. Temam, Inertial manifolds and multigrid methods, SZAMJ. Math. Anal. 21 (1990) 154-178.