MOVING MESH FINITE ELEMENT APPROXIMATIONS ... - CiteSeerX

0 downloads 0 Views 1MB Size Report
HKBU FRG: (Grant numbers: 98-99/II-14 and I-32). ... z Department of Mathematics, The Hong Kong Baptist University, Kowloon Tong, ..... ((div(Ar ) + f)(1 ? h u h.
MOVING MESH FINITE ELEMENT APPROXIMATIONS FOR VARIATIONAL INEQUALITY I: STATIC OBSTACLE PROBLEM RUO LI y , WENBIN LIU

x 

, AND TAO TANG

z

Abstract. Finite element schemes loose accuracy when approximating a class of vari-

ational inequalities, elliptic obstacle problems, due to the existence of free boundaries. In this paper, moving mesh nite element method is applied to solve the elliptic obstacle problems. Computational meshes are constructed by combining harmonic mapping and sharper a posteriori error estimators which are normally used in h-version adaptive nite element approximation. Some important issues such as the selection of monitor functions and monitor smoothing are addressed. The numerical schemes are applied to a number of test problems in two dimensions. It is shown that the moving mesh nite element methods with appropriate monitor functions eliminate (to leading order) the errors arising from the free boundaries.

Funded by:  Special Funds for Major State Basic Research Projects of China.  EPSRC (Grant number: GR/L67387)  Hong Kong RGC (Grant number: 2033/99P)  HKBU FRG: (Grant numbers: 98-99/II-14 and I-32). Subject Classi cation: 65N15, 65N30, 41A05, 41A29, 41A36. Keywords: Finite element approximation, adaptivity, moving mesh method, a posteriori

error estimators, obstacle problems.

School of Mathematical Sciences, Peking University, Beijing 100871, P. R. China. CBS & Institute of Mathematics and Statistics, The University of Kent, Canterbury, CT2 7NF England, Email: [email protected]. z Department of Mathematics, The Hong Kong Baptist University, Kowloon Tong, Hong Kong. Email: [email protected] . 1 y

x

2

R. LI, W.-B. LIU, AND T. TANG

1. Introduction Variational inequalities represent an important class of nonlinear physical models. In particular, variational inequalities have been used to study and compute various free or moving boundary problems, and have been an important subject of studies since 60s. Obstacle problem is among the most studied variational inequalities due to its importance in the theory and application of the variational inequalities. Some of the progress in this area can be found in, for example, [1, 2, 3, 4, 5], and the references cited therein. Finite element approximation of variational inequalities is among important topics in variational inequality theory, and has been extensively studied. The literature is too large to give an even very brief review here, though some of earlier work has been summarized in [6] and [7]. In recent years adaptive algorithms for the nite element approximation have been extensively investigated. They ensure a higher density of nodes in certain area of the given domain, where the solution is more dicult to approximate. At the heart of any adaptive nite element method is an error monitor or indicator. The decision of whether further adjustment of meshes is necessary is based on the estimate of the discretization error. If further adjustment is to be performed then the error monitor or indicator is used as a guide as to show the adjustment might be accomplished most eciently. The literature in this area is also large. Some of techniques directly relevant to our work can be found in [8, 9, 10, 11, 12] and [13, 14, 15]. Adaptive nite element approximation is among the most important means to boost accuracy and eciency of the nite element discretization. There has been some work on h-version adaptive nite element computation of the variational inequalities (particularly of static obstacle problems) in the last decade, see, e.g., [16] using local solution, and [17] using gradient recovery. Recently there is a renewed interest in deriving a posteriori error estimates of residual type and in developing ecient adaptive schemes for nite element approximation of the variational inequalities, see [18, 19, 20] where various residual type estimates are derived. In particular, equivalent a posteriori error estimates of residual type are given in [20] for static obstacle problem, and implemented in numerical tests. Based on these estimates, various h-version adaptive nite element schemes have been implemented quite successfully in tracking the coincidence set of static obstacle problems. One of the important conclusions from the existing theoretical analysis and numerical experiment is that in general the maximum approximation error is concentrated around free boundaries (i.e., the boundary of the coincidence set), as to be seen later. Thus it is essential to use sharper a posteriori error estimators, which can identify the free boundaries, as error indicators in h-version adaptive schemes, in order to re ne the meshes correctly around the free boundary. It is still a challenging problem to implement h-version adaptive nite element scheme eciently for evolutional obstacle problems which frequently arise in practice such as the American Option pricing. For evolutional problems, the eciency can be very poor even with sharper a posteriori error estimators or indicators. The essential diculty seems

MOVING MESH APPROXIMATIONS FOR OBSTACLE PROBLEM

3

from the existing implementation techniques of h-version mesh construction. Matured techniques have been developed to re ne the meshes according to a given error indicator aiming to equi-distribute the approximation errors, see, [21, 22], but it is dicult to coarsen the meshes except in one-dimension. More precisely, once an area is re ned, then one can only either leave it or further re ne it unless a completely new mesh is used. This defect has little impact in static problems. However, in the time dependent problems, once a computational diculty area (e.g. a singular point of the solution) is re ned at a time, the mesh will be kept at least so dense in that area in all computational time afterwards, even after the singularity is moved away. This may seriously reduce any bene ts brought by h-version adaptive nite element schemes for evolutional obstacle problems: Since a computational ecient mesh is very ne around the moving boundary, a large stripe of very ne mesh may be developed, which may be completely unnecessary for the approximation. One of the possible cures is to use the moving mesh method in the adaptive nite element schemes. Moving mesh nite element methods have been widely studied. In Miller [23, 24], methods based on Galerkin formulation were given. There are many applications and extensions of Miller's moving nite element methods, see e.g. Baines [25], Russell and co-authors [26, 27] and Flaherty and co-authors [28, 29]. Numerical error analysis and convergence theory have been investigated by Dupont, Banks, Liu etc, see e.g. [30, 31, 32]. On the practical side, several moving mesh techniques have been introduced in the past, in which one of the most advocated methods is the one based on solving elliptic partial di erential equations (PDEs) rst proposed by Winslow [33]. There are also many applications and extensions of Winslow's method, see e.g. Godunov and Prokopov [34], Brackbill et al. [35, 36], Thompson et al. [37], and Ren and Wang [38]. In most of these works, the moving mesh methods have been applied to partial di erential equations with large solution gradients. To our knowledge, there seems no existing work on implementing the moving mesh method to approximate variational inequalities either static or evolutional. Therefore, as the rst step we will apply the moving mesh methods to the static obstacle problem before we attack the more complicated case of the evolutional problem. The purpose of this paper is to investigate the moving mesh adaptive nite element approximation for the following elliptic obstacle problem:

(1.1)

8 > > > > < > > > > :

?div(Aru)  f; u  ; (?div(Aru) ? f )(u ? ) = 0 ;

where A is a given matrix function, and f and  are given functions. The present work should pave the way to further investigation of ecient adaptive numerical methods for more general evolutional obstacle problems and variational inequalities.

4

R. LI, W.-B. LIU, AND T. TANG

The moving mesh methods to be used in this research will be based on harmonic mapping, see e.g. [35, 39, 40]. In order to obtain the harmonic maps eciently, we construct the harmonic map between the physical mesh and the logical mesh by an iteration procedure. Each iteration step is to move the mesh closer to the harmonic map. The idea of iteration mapping was proposed by Tang and Trummer [41] and Liu and Tang [42]. Selecting appropriate error monitors is one of the most important issues that have to be addressed in implementing the moving mesh methods. It seems that the commonly used gradient based monitors are not suitable in this case. In this work, we will show that the monitor based on sharper a posteriori error estimates can handle the free boundaries very well. The plan of the paper is as follows: In Section 2, we formulate nite element approximation of elliptic obstacle problems. In Section 3, the moving mesh method which is based on harmonic mapping is introduced. In Section 4, we discuss several possible selections of the error monitors to be used in computations. Four supporting numerical examples are given in Section 5. 2. Finite Element Approximation Let be a bounded open set in R2 with a Lipschitz boundary @ . In this paper we adopt the standard notation W m;q ( ) for Sobolev spaces on with norm k  km;q; and semi-norm j  jm;q; . We set W0m;q ( )  fw 2 W m;q ( ) : wj@ = 0g. We denote W m;2( ) by H m( ) with norm kkm; and semi-norm jjm; . In addition c or C denotes a general positive constant independent of h. We rst review some known results about variational inequalities and in particular, the obstacle problem (1.1). Let f 2 L2( ); g 2 H 1( ) and Vg = fv 2 H 1( ) : v ? g 2 H01( )g: Let Z

a(u; v) = (Aru; rv)Rn ; (f; v) = where

Z



fv;

8u; v 2 H 1( );

8v 2 H 1( );

A() = (ai;j ())nn 2 (L1( ))nn: Assume that there is a constant c > 0 such that for any vector X 2 Rn X tAX  ckX k2Rn : Let K be a closed convex set in Vg . Then the following inequality is commonly referred to as variational inequality: For given f; g, nd u 2 K such that (2.1) a(u; v ? u)  (f; v ? u); 8v 2 K:

MOVING MESH APPROXIMATIONS FOR OBSTACLE PROBLEM

5

When the matrix A is symmetric, the above variational inequality problem is equivalent to the following minimization problem (see e.g. [5])): 1 a(u; u) ? (f; u): min (2.2) u2K 2 It is well known that the variational inequality (2.1) has a unique solution in K , see [5] for example. In particular, let  2 H 1( ) such that (g ? )j@  0. Then the following variational inequality is a weak formulation of the obstacle problem (1.1) with the boundary condition u = g on @ : For given f; ; g, nd u 2 Kg () such that (2.3) a(u; v ? u)  (f; v ? u); 8v 2 Kg (); where Kg () = fv 2 Vg : v  g: It is well known that the variational inequality (2.3) has a unique solution u 2 Kg (), and furthermore that the solution is regular, say in H 2( ), if the data are so. For instance, u 2 W 2;p( ) for any p > 1 if A 2 (C 1 (  ))nn;  2 C 2(  ); f 2 C (  ); g 2 C 2+ (  ); @ 2 C 2+ for an > 0 (see e.g. [4]). If u 2 H 2( ), then (2.3) is equivalent to (1.1) (see e.g. [4]). In the rest of the paper we assume that the solution of (2.1) is continuous, which implies that the noncoincidence set of the obstacle problem (1.1)

+ = fx 2 : u(x) > (x)g is open, and the coincidence set

? = fx 2 : u(x) = (x)g is closed. The boundary of the coincidence set in is referred to as free boundary. In this section we study the nite element approximation of the elliptic obstacle problem. For ease of exposition we only consider triangular elements in this paper. Let h be a polygonal approximation to with boundary @ h . Let T h be a partitioning of h S into disjoint regular triangular  , so that  h =  2T h  . Each element has at most one face on @ h , and  and 0 have either only one common vertex or a whole face if  and  0 2 T h. We further require that Pi 2 @ h ) Pi 2 @ where fPig(i = 1::::J ) is the vertex set associated with the triangulation T h. Let h denote the maximum diameter of the element  in T h. We shall only discuss the piecewise linear element due to the limited higher order regularity for the solution of a variational inequality. Associated with T h is a nite dimensional subspace V h of C 0(  h), such that j are ane functions for all  2 V h and  2 T h. For ease of exposition we will assume that h = , though all the results can be extended to more general cases. Let V0h = f 2 V h : (Pi) = 0; 8Pi 2 @ g; and Vgh = f 2 V h : (Pi) = g(Pi); 8Pi 2 @ g:

6

R. LI, W.-B. LIU, AND T. TANG

Let h : C 0 (  h) to V h, denote the standard Lagrange interpolation operator such that for any v 2 C 0 (  h), h v(Pi) = v(Pi) for all Pi (i = 1;    ; J ). We now consider the nite element approximation for (2.1). We only consider conforming linear Lagrange elements here. Let K h be a closed convex set in Vgh. Then a possible nite element approximation of (2.1) reads: For given f and g, nd uh 2 K h such that (2.4) a(uh; vh ? uh)  (f; vh ? uh); 8vh 2 K h: It follows from the standard theory that there exists a unique solution to this approximation scheme. We shall devote the rest of this section to the nite element approximation of the obstacle problem (2.3). We rst consider a special obstacle problem where  = 0: For given f; g, nd u 2 Kg such that (2.5) a(u; v ? u)  (f; v ? u); 8v 2 Kg ; where Kg = fv 2 Vg : v  0g: Note that here we have gj@  0. More general obstacle problems may be transferred into this form by using the transformation u = u ? . Let Kgh be a closed convex set in Vgh. Then the nite element approximation of (2.5) reads: Find uh 2 Kgh such that (2.6) a(uh; vh ? uh)  (f; vh ? uh) 8vh 2 Kgh: It is known that there is a unique solution to this approximation scheme which converges to the solution of (2.5) as h ! 0, if Kgh converges to Kg in the sense given in, say, [7]. In this paper we shall take Kgh = fvh 2 Vgh : vh  0g: It should be noted that Kgh may not be contained in Kg if g is not an ane function. Finite element approximation of the obstacle problem has been extensively studied in the past, see, e.g, [7]. The main purpose of this paper is to approximate the obstacle problems using the moving mesh nite element method based on harmonic mapping, which will be discussed in the next section. 3. Moving Mesh Finite Element Method Based on Harmonic Mapping In this section, we brie y review the moving mesh methods introduced in Li, Tang and Zhang [40], where the moving mesh methods were proposed to solve time-dependent PDEs. In this work, the moving mesh methods will be modi ed to solve the timeindependent variational inequalities in Sections 4 and 5. Consider the following PDEs in a domain  R n (3.1) ~ut = L(~u) in  [0; T ]

MOVING MESH APPROXIMATIONS FOR OBSTACLE PROBLEM

7

with boundary and initial conditions B~uj@ = ~ub in @  [0; T ] ~ujt=0 = ~u0 in : The basic idea of the moving mesh method (also called r-method) is to equidistribute some functional of the numerical solution on the mesh. In 1-D case, it is called equidistribution principle as proposed by de Boor [43]. In the moving mesh method, the functional to be equidistributed is named monitor. The natural extension of the equidistribution principle is to construct the harmonic mappings from physical domain to the computational domain c:  : x 7! 

! c by minimizing the energy functional (3.2)

E ( ) =

XZ

k



k @ k Gij @ @xi @xj dx

in which the standard summation convention is assumed. In the functional (3.2), Gij is the inverse matrix of the monitor and can be regarded as the metric matrix on . The existence and uniqueness of harmonic maps are proved by Hamilton [44], Schoen and Yau [45]. A framework for adaptive grid generation based on the Hamilton-SchoenYau theory was proposed by Dvinsky [39]. In Li-Tang-Zhang [40], Dvinsky's method was extended to provide an ecient solver for mesh-redistribution. The key idea is to construct the harmonic map between the physical space and a parameter space by an iteration procedure. Each iteration step is to move the mesh closer to the harmonic map. In this work, we will apply the numerical procedure proposed in [40] to solve the elliptic obstacle problem (1.1). We only outline the basic idea of the moving mesh method { the detailed procedure, derivations and formulas for the moving mesh strategy based on harmonic mapping can be found in [40]. The Euler-Lagrange equation of the functional (3.2) is of the form 



@ Gij @ k = 0 : @xi @xj

(3.3)

In principle, the desired harmonic mapping can be obtained numerically. However, in the numerical implementation what we need to compute is the inverse of the harmonic mapping instead of the harmonic mapping itself. The numerical approach to sway o the spot is to deform the given mapping to the harmonic one by solving the heat equation (3.4)





@ k = @ Gij @ k : @ @xi @xj

Transforming back to the physical domain according to (3.5) will lead to the map from ~ to ~x.

@xl = ? @xl @ k ; @ @ k @

8

R. LI, W.-B. LIU, AND T. TANG

In the rst step, we prepare the mesh on the logical domain c. Let X = (Xi) be the nodes, and  = (j ) be basis functions such that j (Xi) = ij , where ij is the Kronecker delta. We also choose a convex domain c as the logical domain. By solving the Poisson equation (3.6) ~ = 0 with Dirichlet boundary condition: (3.7) ~j@ = ~b; we obtain a mesh Tc in the logical domain. The above harmonic map is the same as the one originally proposed by Winslow [33]. By solving equations derived from (3.5), we will obtain the mesh-moving vector ~x. With this information, we move the mesh to a new location by updating the location of nodes and updating the values of those functions de ned on the mesh by interpolation. By repeating the procedure iteratively, we can obtain the inverse of the harmonic mapping within the desired accuracy. The complete mesh-moving process can be summarized as the following:

Algorithm 1

(i) compute the mesh-moving vector ~x; (ii) move the mesh to a new location; (iii) interpolate those functions from the old mesh to the new mesh; (iv) judge if the mesh is close to the harmonic map. if not, go to (i).

In part (iv), in order to ensure the quality of the harmonic map, we repeat the operation of the mesh-moving until the L2 -norm for the distance between the solution of equation (3.3) and the logical mesh Tc is smaller than a preassigned tolerance TOL. One of the important issues requiring some attention is the monitor function Gij used in the equation (3.2). There are several excellent papers in this direction, including Brackbill [35] and Cao et al. [27]. Most of the works in this direction are concerned with the solutions with large gradients (or thin layers), so it is natural to use the solution gradients in the monitor function (see gradient monitor in Sect. 4). However, it is quite complicated to choose a proper monitor in solving the obstacle problem (1.1). The standard monitor function involving gradients can not produce a satisfactory mesh. It is found that a good choice of the monitor function should be based on sharp a posteriori error estimates. We will devote the next section to the study of the monitors suitable to the obstacle problem. 4. Construction of Error Monitors It is vital to select a suitable error monitor in the moving mesh method. Traditionally, some weighted gradient or second order derivative of the solutions are used as error

MOVING MESH APPROXIMATIONS FOR OBSTACLE PROBLEM

9

monitors in the moving mesh methods. However, it can be clearly seen form Figure 2 that larger errors occur on the free boundary, and the approximation errors around the free boundary obtained from the moving mesh methods with the gradient-based monitors can be even worse than those from the nite element method with the uniform mesh. The reason is that the gradients are not necessarily large around the free boundary. It is clear that a suitable monitor should be able to identify the free boundary. Our approach is to combine existing a posteriori error estimation of residual type with the moving mesh method. We rst introduce two a posteriori error estimates for the nite element approximation of the obstacle problem with zero obstacle. The rst estimate can be found in Chen and Nochetto [18] and Liu and Yan [20].

Theorem 4.1. Let Kg = fv 2 Vg : v  0g and Kgh = fvh 2 Vgh : vh  0g. Assume that g 2 C 2(  ). Let u and uh be the solutions of (2.5) and (2.6), respectively. Then ju ? uhj21;  C (12 + 22 + 32 ); (4.1) where

12 22

= =

X



h2

X

lZ\@ =;

Z



(f + div(Aruh))2 ; Z

hl [(Aruh)  n]2 ; l

j(Aruh)  njjg ? hgj; where hl is the size of an edge l, and [(Aruh)  n] denotes the jump of A-normal derivative 32 =

@

on an interior edge l, de ned by

[(Aruh)  n]l = (Aruhjl1 ? Aruhjl2 )  n; l = l1 \ l2 ;

where n is the outer normal vector of l1 .

The estimate (4.1) is actually held for more general variational inequalities, see [20]. The key ingredient in deriving the estimate is to construct suitable positive-keeping and H 1 stable interpolator, see [18]. The main disadvantage of this estimate is that it only provides a upper bound for the error. It is clear that in general the estimate cannot be a lower bound since one could change the value of f in the coincidence set provided f  0 without changing that of the solution. Thus the value of f inside the coincidence set should take little e ect in sharper error estimates. By exploiting the fact that the obstacle problem may be equivalent to the following nonlinear equation: For given f; g, nd u 2 Kg such that (4.2)

a(u; v) + (fu; v) = (f; v)

8v 2 H01( );

10

R. LI, W.-B. LIU, AND T. TANG

where



0 u(x) = 10 uu((xx)) = >0

sharper a posteriori error estimates have indeed been derived in [20]. The equivalence certainly holds if the solution u 2 H 2( ) and the free boundary has Lebesgue measure zero, see [20]. Theorem 4.2. Assume that A is a constant matrix. Let u, uh be the solutions of (2.5) and (2.6) with f  0 on ? . Let Kgh = fv 2 Vgh : vh  0g, and g 2 H 2( ). Assume that u 2 W 1;p( ) with p > n and that there is an 2 (0; 2] such that ku ? uhk0;1;  Ch jln(h)j: Suppose that u is also the solution of (4.2). Then (4.3) c(^12 + ^22 + 32 + 1 )  ju ? uhj21;  C (^12 + ^22 + 32 + 2 ); where

^12

=

X



h2

Z



(f (1 ? huh ))2;

huh = h h+ u ;

^22

=

X

l\@ =;

Z

hl [(Aruh)  n]2 ; l

0< < ;

h

1 and 2 are some higher order terms when the solution is in H 2( ) and the free boundary

is consist of nite pieces of smooth curves.

It follows from the proof of Theorem 5.1 in [20] that the upper bound hold for the case where A is a variable matrix if one replaces ^1 with X



h2

Z



(f (1 ? uh ) + div(Aruh))2 :

The case  6= 0 can be easily dealt with by using the transformation u = u ?  provided  2 H 2( ). In this case one has (see [20]) (4.4) c(^12 + ^22 + 32 + 1 )  ju ? uhj21;  C (^12 + ^22 + 32 + 2 ); where

^12 ^22

= =

X



h2

X

l\@ =;

Z



((div(Ar) + f )(1 ? huh ))2; Z

hl [(Aruh)  n]2 : l

It is seen that the estimates in Theorem 4.2 are sharper than that given by Theorem 4.1, in the sense that the inequalities (4.3) provide an equivalent a posteriori error estimate.

MOVING MESH APPROXIMATIONS FOR OBSTACLE PROBLEM

11

Remark 4.1. When  is not smooth, then it is much more dicult to obtain an equiv-

alent a posteriori error estimator if the solution is not smooth. If the soution is still smooth, then the nonsmooth part of the obstacle is normally outside the coincidence set. In that case, we have the following a posteriori error estimate: (4.5) ju ? uhj21;  C (^12 + ^22 + 32 + ); where  is a higher order term,

^12 ^22 where

= =

X



h2

X

l\@ =;

Z



(f (1 ? huh ) ? div(Ar)huh )2; Z

hl [(Aruh)  n]2 ; l

huh = h + hu ?  : h

h

In this paper, we construct the error monitors in the following way: Assume that P 2 =   is an a posteriori error estimator, where  may contain integral terms around its edges. Then the error monitor M is de ned by p M = c + 2 ; where c  0 is a constant to be decided in computation. It can be shown that 3 above is normally of higher order than the others. Thus it will be omitted in our computations. One of the diculties in applying the moving mesh techniques to the obstacle problem is that on the coincidence set the error estimators may be zero. This is a good feature for h-version re nement as it increases computational eciency. However, in the moving mesh methods the monitor will be used as metric matrix, which should be positive de ned, and no-zero at any point. In case the error estimator is zero at some point, we need to introduce a small positive number c within the estimator. The magnitude of c will decide the density of the nodes on the coincidence set. If c is too large, then there will be too many elements inside the set, which are not needed, as main approximation error is concentrated around the free boundary. On the other hand, if c is small, there will have a large jump of density of elements over the free boundary as the monitor is relatively large just outside the coincidence set and near zero inside the set. This, as observed in our numerical experiments, leads to very singular meshes and higher approximation error around the free boundary. We overcome this diculty by smoothing the monitors around free boundary. We have tried two strategies:  Monitor smoothing strategy I: To make the monitor more smooth to cope with the free boundary, we rst interpolate M from L2( ) into H1;h( ), from piecewise constant to piecewise linear:

2

(hM )jat P =

P

P :P is vertex of  M jon  j j :  :P is vertex of  j j

12

R. LI, W.-B. LIU, AND T. TANG

We then project it back into L2 ( ): X M~ jon  = 31 (h M )at P : P :P is vertex of   Monitor smoothing strategy II: An alternative approach for smoothing the monitor is the following: Given d > 0. If M  = 0, let M = max 2Th M =2, where Th = f : dist(;  )  dg. Then we replace M  by M . 5. Numerical Experiments Let h be a polygonal approximation to with boundary @S h . Let T h be a partitioning of h into disjoint regular triangular  , so that  h =  2T h  . We shall only consider the piecewise linear element due to the limited higher order regularity for the solution of a variational inequality. Denote Xi as nodes and j as basis functions, such that j (Xi) = ij . In all our experiments, A in the variational inequality (2.1) is chosen as the unit matrix. Thus it follows from the equivalence between (2.1) and the minimization problem (2.2) that the variational inequality (2.3) can be approximated by the following optimization problem: (5.1) min 1 U T HU + qT U 2 s:t: U  where Z i ij uh = Ui  ; H = ri  rj dx ;

qi =



Z



f idx ; i = (Xi) :

The functions f and  used above are given by the elliptic obstacle problem (1.1). 5.1. The algorithm. In order to solve the problem (5.1) with adaptive mesh, we use the following algorithm:

Algorithm 2 (i) solve the optimization problem (5.1) on the current mesh; (ii) compute the mesh-moving vector ~x; (iii) move the mesh to a new location, as described in [40]; (iv) judge if the mesh is close to the harmonic map. if not, go to (i).

In part (i), we use a projection gradient method to solve the optimization problem (5.1). Here we brie y describe the method { the detailed description of the method can be found in [46]. For given H and q, if the projection operator is denoted by proj (), then  (1) g = HU + q;

MOVING MESH APPROXIMATIONS FOR OBSTACLE PROBLEM

13

 (2) e = U ? proj (U ? g);  (3) error = kek22, if e <  then exit;  (4)  = kek2k+eke2T He ; 2  (5) U = U ? e, go to (1). In part (iv), in order to ensure the quality of the harmonic map, we repeat the operation of mesh-moving until the L2 -norm for the distance between the solution of equation (3.3) and the logical mesh Tc is smaller than a preassigned tolerance TOL. The essential di erence between Algorithm 1 described in Section 3 and Algorithm 2 is that we do not update the value of u on the new mesh by interpolation. Instead we obtain the value of u on the new mesh by solving the optimization problem (5.1) which is an approximation to the variational inequality (2.3). 5.2. The monitor function. In this subsection, we choose three di erent types of monitor functions.

 Monitor I: Gradient Monitor

p

M = 1 + jruj2

where is a positive constant. This monitor has been used by many authors, see e.g. [27, 35, 38, 40]. However, as the numerical results shown later, it is not a good choice for the obstacle problem.

 Monitor II: Monitor based on a non-equivalent a posterior error estimator M jon  =

s

h2

Z



f2 +

Z

X

l: edge of 

hl [(Auh)  n]2 : l

As indicated in Theorem 4.1 the above monitor function provides an upper bound for the nite element approximation to the variational inequality (2.3).

 Monitor III: Monitor based on an equivalent a posterior error estimator s

M jon  = c + 12; +

X

l: edge of 

22;l ;

where c is a positive constant, 1; is the integration in every element

12;

= h2

Z



(f (1 ? huh ))2;

and 2 is the integration on every side

22;l

Z

= hl [(Auh)  n]2 : l

The reason for adding a positive constant c in the de nition for Monitor III is that the monitor is to be used as metric matrix so it can not be zero at any point.

14

R. LI, W.-B. LIU, AND T. TANG

5.3. Numerical Examples. In order to illustrate the foregoing numerical schemes, we describe the results for four particular problems.

Example 5.1. In the rst example, we choose a square physical domain and a constant f : = [?1:5; 1:5]  [?1:5; 1:5], f = ?2. The obstacle function  is zero. The exact solution is

 r2

1 u = 02 ? ln(r) + ln(1) ? 2

if r  1 if r < 1

p

where r = x2 + y2. The free boundary is u(x; y) = 0 along the circle r = 1. The exact solution is plotted in Figure 1.

1

0.8

0.6

0.4

0.2

0 −1.5

−1.5 −1

−1 −0.5

−0.5 0

0 0.5

0.5 1

1 1.5

1.5

Figure 1. The exact solution for Example 5.1.

We solve Example 5.1 by using the three monitors described in the last subsection. In using Monitor I, namely the gradient monitor, the constant is set to be 10. In using Monitor III, the constant c is set to be 10?3 . In Figure 2, we display the numerical results obtained by using 30  30 elements. For comparison, we also solve this problem using a uniform grid. It is observed from Figure 2 that with the uniform grid and Monitors I and II the large errors occur on the free boundary. In particular, it is seen that with the gradient monitor (Monitor I) the large errors occur only on a set of one order lower dimensionality, i.e. the free boundary. This implies that the gradient monitor can not catch the free boundary and therefore is not appropriate for the obstacle problem. On the other hand, the monitor based on the equivalent a posterior error estimator, namely Monitor III, removes the errors from the free boundary so that the overall error level is decreased signi cantly. More importantly, we observe from the plots on the left hand side of Figure 2 that the moving mesh method with Monitor III provides a satisfactory approximation for the free boundary, while the other three approaches either give a wrong indication of the free boundary or give no information on it.

MOVING MESH APPROXIMATIONS FOR OBSTACLE PROBLEM

15

(a) Uniform mesh: −3

1.5 2

x 10

1.8

1 1.6

1.4

0.5

1.2

0

1

0.8

−0.5 0.6

0.4

−1 0.2

−1.5 −1.5

−1

−0.5

0

0.5

1

1.5

0 −1.5

−1

−0.5

0

0.5

1

1.5

= 10:

(b) Monitor I with −3

1.5 7

1

x 10

6

5

0.5

4

0 3

−0.5 2

−1 1

−1.5 −1.5

−1

−0.5

0

0.5

1

1.5

0 −1.5

−1

−0.5

0

0.5

1

1.5

(c) Monitor II: −3

1.5 2

x 10

1.8

1 1.6

1.4

0.5

1.2

0

1

0.8

−0.5 0.6

0.4

−1 0.2

−1.5 −1.5

−1

−0.5

0

0.5

1

1.5

0 −1.5

−1

(d) Monitor III with

−0.5

0

0.5

1

1.5

0.5

1

1.5

c = 10?3:

−3

1.5 2

x 10

1.8

1 1.6

1.4

0.5

1.2

0

1

0.8

−0.5 0.6

0.4

−1 0.2

−1.5 −1.5

−1

−0.5

0

0.5

1

1.5

0 −1.5

−1

−0.5

0

Figure 2. Adaptive meshes and the solution errors for Example 5.1. The

maximum errors are approximately 1:9  10?3; 6:5  10?3; 1:7  10?3 and 0:7  10?3 for the above four cases, respectively. Example 5.2. This example(is as same as the Example 5.1, except that f is given by r ?4 0:5 2 f = 40:95 ? ? 2 0:95 ifif rr  < 0:5: r

16

R. LI, W.-B. LIU, AND T. TANG (a) Uniform mesh: −3

1.5

x 10

2

1

0.5

1.5

0 1

−0.5

0.5

−1

−1.5 −1.5

−1

−0.5

0

0.5

1

1.5

0 −1.5

−1

−0.5

0

0.5

1

1.5

= 10:

(b) Monitor I with −3

1.5 9

x 10

8

1 7

0.5

6

5

0 4

−0.5

3

2

−1 1

−1.5 −1.5

−1

−0.5

0

0.5

1

1.5

0 −1.5

−1

−0.5

0

0.5

1

1.5

(c) Monitor II: −3

1.5

x 10

2

1

0.5

1.5

0 1

−0.5

0.5

−1

−1.5 −1.5

−1

−0.5

0

0.5

1

1.5

0 −1.5

−1

(d) Monitor III with

−0.5

0

0.5

1

1.5

0.5

1

1.5

c = 10?3:

−3

1.5

x 10

2

1

0.5

1.5

0 1

−0.5

0.5

−1

−1.5 −1.5

−1

−0.5

0

0.5

1

1.5

0 −1.5

−1

−0.5

0

Figure 3. Adaptive meshes and the solution errors for Example 5.2. The

maximum errors are approximately 2:0  10?3; 9:5  10?3; 1:5  10?3 and 1:0  10?3 for the above four cases, respectively.

MOVING MESH APPROXIMATIONS FOR OBSTACLE PROBLEM

The exact solution is p

17



2 1 u = (0r ? 1) rr < 1

where r = x2 + y2. Again, we compare the performance of the di erent monitor functions in Figure 3, with 30  30 elements used. Observations similar to Example 5.1 are made from Figure 3. Furthermore, it is clear from Figures 2 and 3 that the mesh re nement scheme based on the non-equivalent a posterior error estimates, namely Monitor II, re nes the meshes almost uniformly everywhere in the computational domain. It can not recognize the coincident set. However, the moving mesh scheme based on the equivalent a posterior error estimates, namely Monitor III, only re nes the domain outside the coincidence set and therefore is more ecient. To demonstrate the convergence of the moving mesh method with Monitor III, we compute Example 5.2 with 20  20, 30  30 and 40  40 elements, and display the adaptive meshes and the solution errors in Figure 4. Indeed we observe from Figure 4 that the moving mesh methods with Monitor III have increased the accuracy signi cantly. The leading errors at the free boundary are eliminated with the moving mesh method. Example 5.3. In this example, we choose = [?1; 1]  [?1; 1], f = 1. Let K = fu 2 H01( ) : u  g with (x) = dist(x; @ ). Note that the above problem is a variation of type (2.3) and the obstacle function is piecewise linear. It was studies by Kornhuber [17]. The exact solution is unknown. For this nonzero-constraint problem, we modify the characteristic function huh by letting

huh = u ? h + h ; h h

where is some positive constant, as pointed out in Remark 4.1. To demonstrate the convergence of our moving mesh methods with Monitor III, we compute Example 5.3 with 20  20, 30  30 and 40  40 elements and display the adaptive meshes and the numerical solutions in Figure 5. It indicates that the moving mesh method with Monitor III works satisfactorily for this test problem. Example 5.4. In the last example, the obstacle  = 0; f = ?1 and 8 0:5h21 ? qx 0  x  l; y = 0 > > < 2 h1 ? y ) + x = l; 0 < y  h1 uj@ = > 00::5( 2 5(h2 ? y)+ x = 0; 0 < y  h1 > : 0 0 < x < l; y = h1 where h1 = 10, h2 = 2, l = 5, q = (h21 ? h22 )=(2l), and (h ? y)+ = max(h ? y; 0). The problem is de ned in the domain = [0; l]  [0; h1]. This problem arises in the modeling of the seepage of an incompressible, inviscid uid through an unsaturated rectangular dam. The true solution of this problem is not readily accessible. This problem was studied numerically by Ainsworth, Oden and Lee [16]

18

R. LI, W.-B. LIU, AND T. TANG −3

1.5 3

1

2.5

0.5

2

0

1.5

−0.5

1

−1

0.5

−1.5 −1.5

−1

−0.5

0

0.5

1

1.5

0 −1.5

3

1

−0.5

0

0.5

1

1.5

−1

−0.5

0

0.5

1

1.5

−1

−0.5

0

0.5

1

1.5

x 10

2.5

0.5

2

0

1.5

−0.5

1

−1

0.5

−1

−0.5

0

0.5

1

1.5

0 −1.5

−3

1.5 3

1

x 10

2.5

0.5

2

0

1.5

−0.5

1

−1

0.5

−1.5 −1.5

−1

−3

1.5

−1.5 −1.5

x 10

−1

−0.5

0

0.5

1

1.5

0 −1.5

Figure 4. Adaptive meshes and solution errors with 20  20, 30  30 and

40  40 grid points (from top to bottom), respectively, for Example 5.2.

who derived local a posteriori error estimator for assessing the error in a nite element approximation. Again, we compute Example 5.4 with di erent number of elements. The adaptive meshes and the numerical solutions are shown in Figure 6. The numerical solutions as well as the curve of free boundary seem convergent as the number of elements is suciently large. This indicates that the moving mesh method with Monitor III work well for this test problem. 6. Conclusions The main idea of this work is to solve the obstacle problem by combining the moving mesh technique and the equivalenet a posteriori error estimators developed by the authors recently [40, 20]. Some technical parts on how to implement the moving mesh

MOVING MESH APPROXIMATIONS FOR OBSTACLE PROBLEM

19

1

0.8

0.6 1

0.4 0.8

0.2 0.6

0 0.4

−0.2 0.2

−0.4 0 −1

−1

−0.6 −0.5

−0.5 0

−0.8

0 0.5

−1 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

0.5 1

1

1

1

0.8

0.6 1

0.4 0.8

0.2 0.6

0 0.4

−0.2 0.2

−0.4 0 −1

−1

−0.6 −0.5

−0.5 0

−0.8

0 0.5

−1 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

0.5 1

1

1

1

0.8

0.6 1

0.4 0.8

0.2 0.6

0 0.4

−0.2 0.2

−0.4 0 −1

−1

−0.6 −0.5

−0.5 0

−0.8

0 0.5

−1 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

0.5 1

1

Figure 5. Adaptive meshes and numerical solution with 20  20, 30  30

and 40  40 elements, respectively, for Example 5.3. The moving mesh method with Monitor III is used.

methods, such as the selection of monitor functions and monitor smoothing, are addressed. The numerical examples indicate that the numerical scheme suggested in this work, i.e. moving nite element method with Monitor III, eliminates (to leading order) the errors arising from the free boundaries, while the moving nite element method with some commonly used monitors either re nes the meshes almost uniformly everywhere in the computational domain (such as Monitor II) or give wrong information on the free boundaries (such as Monitor I). In the second part of this work, we will extend the numerical techniques developed in this work to solve evolutional obstacle problems. To this end, we need to obtain an equivalent a posteriori error estimator, similar to that in Theorem 4.2, which will be used to provide an appropriate monitor function.

20

R. LI, W.-B. LIU, AND T. TANG 10

9

8 50

7 40

6 30

5 20

4 10

3 0 0

2

0 2

1 4

1

2 6

3 8

0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

4 10

5

5

10

9

8 50

7 40

6 30

5 20

4 10

3 0 0

2

0 2

1 4

1

2 6

3 8

0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

4 10

5

5

10

9

8 50

7 40

6 30

5 20

4 10

3 0 0

2

0 2

1 4

1

2 6

3 8

0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

4 10

5

Figure 6. Adaptive meshes and numerical solutions with 10  20, 20  40

and 30  60 elements (from top to bottom), respectively, for Example 5.4. The solution domain is [0; 5]  [0; 10]. The moving mesh method with Monitor III is used.

Acknowledgment. We thank Dr. B. S. He of Nanjing University for introducing us

the optimization code used in Section 5.

MOVING MESH APPROXIMATIONS FOR OBSTACLE PROBLEM

21

References [1] Brezis H. Problemes unilateraux. Journal de Mathematiques Pures et Appliquees 1972; 51:1-168. [2] Duvaut G, Lions JL. The inequalities in mechanics and physics. Springer-Verlag, 1973. [3] Elliott CM, Ockendon JR. Weak and variational methods for moving boundary problems. Research Notes in Mathematics 59. Pitman: Boston, 1982. [4] Friedman A. Variational principles and free-boundary problems, Academic Press, New York, 1982. [5] Kinderlehrer D, Stampacchia, G. An introduction to variational inequalities and their applications, Academic Press, New York, 1980. [6] Ciarlet PG, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [7] Glowinski, R., Lions, J.L., and Tremolieres, R.: Numerical analysis of variational inequalities, North-Holland, Amsterdam, 1976. [8] Ainsworth M, Oden JT. A uni ed approach to a posteriori error estimation using element residual methods. Numerische Mathematik 1993; 65:23-50. [9] Ainsworth M, Oden JT. A posteriori error estimators in nite element analysis. Computer Methods in Applied Mechanics and Engineering 1997; 142:1-88. [10] Bank RE, Weiser A. Some a posteriori error estimators for elliptic partial di erential equations. Mathematics of Computation 1985; 44:283-301. [11] Baranger J, Amri HE. A posteriori error estimators in nite element approximation of quasiNewtonian ows. Mathematical Modelling and Numerical Analysis 1991; 25:31-48. [12] Duran RD, etc. On the asymptotic exactness of error estimators for linear triangular nite elements. Numerische Mathematik 1991; 59:107-127. [13] Verfurth R. A posteriori error estimators for the Stokes equations. Numerische Mathematik 1989. 55:309-325. [14] Verfurth R. A posteriori error estimates for non-linear problems. Mathematics of Computation 1994; 62:445-475. [15] Zienkiewicz OC, Zhu JZ. A simple error estimator and adaptive procedure for practical engineering analysis. International Journal for Numerical Methods in Engineering, 1987; 24:337-35. [16] Ainsworth M, Oden JT, Lee CY. Local a posteriori error estimators for variational inequalities. Numerical Methods for Partial Di erential Equations 1993; 9:22-33. [17] Kornhuber R. A posteriori error estimates for elliptic variational inequalities. Computers and Mathematics with Applications 1996; 31:49-60. [18] Chen ZM, Nochetto RH. Residual type a posteriori error estimates for elliptic obstacle problems. To appear. [19] French DA, Larsson S, Nochetto RH. Pointwise a posteriori error analysis for an adaptive penalty nite element method for the obstacle problem. In preparation. [20] Liu WB, Yan N. A posteriori error estimates for a class of variational inequalities. To appera. [21] Verfurth R. A posteriori error estimators for the Stokes equations. Numerische Mathematik 1989; 55:309-325. [22] Verfurth R. A Review of a Posteriori Error Estimation and Adaptive Mesh-Re nement Techniques. Wiley-Teubner, 1996. [23] Miller K, Miller RN. Moving nite element methods I. SIAM Journal on Numerical Analysis 1981; 18:1019-1032. [24] Miller K. Moving nite element methods II. SIAM Journal on Numerical Analysis 1981; 18:10331057. [25] Baines MJ. Moving Finite Elements. Oxford University Press, 1994. [26] Cao WM, Huang WZ, Russell RD. An r-adaptive nite element method based upon moving mesh PDEs. Journal of Computational Physics 1999; 149:221-244. [27] Cao WM, Huang WZ, Russell RD. A study of monitor functions for two dimensional adaptive mesh generation, SIAM Journal on Scienti c Computing 1999; 20:1978-1994.

22

R. LI, W.-B. LIU, AND T. TANG

[28] Davis SF, Flaherty JE. An adaptive nite element method for initial-boundary value problems for partial di erential equations. SIAM Journal on Scienti c and Statistical Computing 1982; 3:6-27. [29] Moore PK, Flaherty JE. Adaptive local overlapping grid methods for parabolic systems in two space dimensions. Journal of Computational Physics 1992; 98:54-63. [30] Dupont TF. Mesh modi cation for evolution equations. Mathematics of Computation 1985; 44:3952. [31] Bank RE, Santos RF. Analysis of some moving space-time nite element methods. SIAM Journal of Numerical Analysis 1993; 30(1):1-18. [32] Liu Y, Bank RE, Dupont TF, Garcia S, Santos RF. Symmetric error estimates for moving mesh mixed methods for advection-di usion equations, Preprint, 1999. [33] Winslow A. Numerical solution of the quasi-linear Possion equation in a nonuniform triangle mesh. Journal of Computational Physics 1967; 1:149-172. [34] Godunov SK, Prokopov. Computational Mathematics and Mathematical Physics 1971; 12:182. [35] Brackbill JU, Saltzman JS. Adaptive zoning for singular problems in two dimensions. Journal of Computational Physics 1982; 46:342-368. [36] Brackbill JU. An adaptive grid with direction control. Journal of Computational Physics 1993; 108:38-50. [37] Thompson JF, Warsi ZUA, Mastin CW. Numerical Grid Generation. North Holland: New York, 1985. [38] Ren W, Wang XP. An iterative grid redistribution method for singular problems in multiple dimensions. To appear in Journal of Computational Physics. [39] Dvinsky AS. Adaptive grid generation from harmonic maps on Riemannian manifolds. Journal of Computational Physics 1991; 95:450-476. [40] Li R, Tang T, Zhang PW. Moving mesh methods in multiple dimensions based on harmonic maps. Submitted, 2000. [41] Tang T, Trummer MR. Boundary layer resolving pseudospectral methods for singular perturbation problems. SIAM Journal on Scienti c Computing 1996; 17:430-438. [42] Liu WB, Tang T. Spectral methods for singular perturbation problems. Proceedings of Symposia in Applied Mathematics, AMS 1994; 48:323-326. [43] de Boor C. Good approximation by splines with variable knots, II, in Conf. On the Numerical Solution of Di erential Equations. Lecture Notes in Math. 363. Springer-Verlag:Berlin, 1973; 1220. [44] Hamilton R. Harmonic maps of manifolds with boundary. Lecture Notes in Math. 471. SpringerVerlag:New-York, 1975. [45] Schoen R, Yau ST. Inventiones Mathematicae 1978. 44:265. [46] He BS. Solving a class of linear projection equations. Numerische Mathematik 1994; 68:71-80.