多重網格與自調適法於Laplace 方程角奇異解的數值計算

國 立 交 通 大 學

應用數學系

數學建模與科學計算碩士班

碩

士

論

文

多重網格與自調適法於 Laplace 方程角奇異解的

數值計算

Multigrid and Adaptive Methods for Computing

Singular Solutions of Laplace Equation on Corner

Domains

研 究 生：李偉任

指導老師：吳金典 教授

多重網格與自調適法於 Laplace 方程角奇異解的

數值計算

Multigrid and Adaptive Methods for Computing

Singular Solutions of Laplace Equation on Corner

Domains

研 究 生：李偉任 Student：Wei-Jen Lee

指導教授：吳金典 Advisor：Chin-Tien Wu

國 立 交 通 大 學

應用數學系數學建模與科學計算碩士班

碩 士 論 文

A Thesis

Submitted to Department of Applied Mathematics College of Science, Institute of Mathematical Modeling and Scientific Computing

National Chiao Tung University in Partial Fulfillment of the Requirements

for the Degree of Master

In

Applied Mathematics July 2010

Hsinchu, Taiwan, Republic of China

多重網格與自調適法於 Laplace 方程角奇異解的數值計算

學生：李偉任 指導老師：吳金典 教授

國立交通大學應用數學系數學建模與科學計算碩士班

摘 要

橢圓邊界值問題在凹角的地方會有奇異的行為，而這個奇異的行為對於

用有限元素法離散的精確度會受到影響。對於给定 Dirichlet 邊界條件的

Poisson 方程式和在定義域有凹角的情況之下，本論文利用一個奇異解的表

示法

J j=1 u=w+





_js_j

算出較準確的近似值，其中



_j 1{ fs dx₁ u s dx₁ }



    

_

 

_

在工

程上稱之為應力強度因子。這些量的精確計算在許多實際的工程問題上，

是一門很重要的課題。

i

Multigrid and Adaptive Methods for Computing Singular

Solutions of Laplace Equation on Corner Domains

Student: Wei-Jen Lee

Advisor: Chin-Tien Wu

Institute of Mathematical Modeling and Scientific Computing

National Chiao Tung University

Abstract

Elliptic boundary value problems on domain with corners have singular

behavior near the corners. Such singular behavior affect the accuracy of the

finite element method throughout the whole domain. For the Poisson equation

with homogeneous Dirichlet boundary conditions defined on a polygonal

domain with re-entrant corners, it is well known that the solution has the

singular function representation

u=w+



J_j=1



_js_j

,where w is the regular part of the

solution and

sj

are known as singular functions that depend only on the

corresponding re-entrant angles. Coefficients



j

known as the stress intensity

factors in the context of mechanics can be expressed in terms of u by extraction

formula

1 1 { j fs dx u s d1 }





  



_



_

 _ x

, where

sj

are known as dual singular

funciton. Accurate calculation of these quantities is of great importance in many

practical engineering problems. Similar singular function representations hold

for the solutions of interface,biharmonic,elasticity, and evolution problems in [1,

2].

誌 謝

本論文能順利完成，首先要由衷的感謝指導教授吳金典博士在論文的撰寫期間的細 心教導，從整體的架構，資料的蒐集，研究的方法以及邏輯的判斷與思考，甚至於做人 處事及態度上的教誨，使學生在撰寫論文過程中受益良多，對於吾恩師無私的付出，獻 上最誠摯的謝意與感激。 在論文口試中，承蒙口試委員賴明治博士與李國明博士於百忙之中撥冗審閱，提供 了相關多寶貴的意見和指正，使本論文得以更臻嚴謹完整，特此表示由衷的感謝之意。 此外，在兩年的求學生涯中，承蒙學長黃振庭、張裕昇、紀露結、許哲維、方仁洲， 以及同學林威辰、周芳竹、林育賢、蔡玉麟等同學間彼此的支持與勉勵、無論是課業上 或處事上皆能共患難，使我這兩年能夠學習到很多，使我的碩士學位能夠順利完成，我 想這份情感是很難捨的。 最後，絕對必須感謝我的家人，爸爸、媽媽、姊姊及哥哥，長期的教誨和辛苦的養 育，使我在沒有經濟壓力下順利的取得學位。謹此致上我最深的謝意。 國立交通大學應用數學系數學建模與科學計算所碩士班 李偉任 謹誌 中華民國九十九年七月 iii

Contents

Introduction ………1

1 Finite Element Method………..……….………..3

1.1 Introduction of Finite Element Method………..3

1.2 Variational Formulation……….………...….5

1.2.1 Existence and Uniqueness of Solution……….……….5

1.3 Finite Element Discretization……….………6

1.3.1 Linear Interpolation……….………..7

1.3.2 Coordinate Transformations……….……….…8

1.3.3 Linear FEM Discretization………....9

1.3.4 Partial Derivatives……….………...10

1.3.5 Assembling the Element Matrix……….………..12

1.4 Error Estimation ………12

1.4.1 Interpolation Error with Piecewise Linear Functions………..….12

1.4.2 A Priori Error Estimation………...……...18

1.4.3 A Posterior Error Estimation and Adaptive Mesh-Refinement Techniques………..…..…19

1.5 Finite Element Approximation for Singular Functions………..21

1.5.1 Lagrange Interpolation………...…22

v

1.5.2 Singular Element………...………23

2 Multigrid Method………..26

2.1 Introduction of Multigrid Method ……….26

2.2 Relaxation Process……….27

2.2.1 Jacobi Method ………...………...……28

2.2.2 Gauss-Seidel Method………...……….29

2.2.3 Successive Overrelaxation Method………...………30

2.3 Inter-grid Interpolation : Restriction and Prolongation………..32

2.4 Multigrid Algorithm………...34 2.5 Complexity……….37 2.6 Numerical Experiments………..37 3 Research Method………...43 4 Numerical Results……….49 References………55

Introduction

The nite element method has become one of the most popular and effective methods for the numerical solution of partial differential equations, particularly for elliptic equa-tions. In practice many important problems involve polygonal domains. Previously, the geometry of a problem would be restricted so that the triangulation elements t the polyg-onal boundary exactly. Form a theoretical standpoint, under the assumption that solutions were suf ciently smooth, this case has been thoroughly analyzed. Unfortunately, in prac-tice one is not likely to achieve the smoothness required for these previous analyses. It is the basic behavior of elliptic equations that solutions possess singularities at corners. These singularities substantially affect the rate of convergence of numerical approximations.

To handle this problem , here the two main procedures which have been proposed to overcome this dif culty. The rst is based on mesh re nements and has been analysed by Babuska and kellogg [21],Raugel, Schatz and Wahlbin, Thatcher for instance. This method may be applied to most of the practical problems since it requires only a qualitative knowl-edge of the behaviour of the solution near the corners. The second consists in augmenting the space of trail funcitons in which one looks for the approximate solution. This is done by adding some of the singular solutions of the problem to the usual spaces of piecewise polynomial funcitons. For instance, S.C. Brenner and L.Y. SUNG [9], Babuska and Rosen-zweig [22], Kellogg, Lelievre, Djaoua and Ladeveze and Peyret.

In this thesis ,the rst we introduce the singular element to capture the singular point at the corner and see the accuracy reduced.The advantage of the singular element is that

Introduction 2

there are many small tranguler near the singular point , the solution near the singular point can be approximation ef ciency, but the disadvantage is that create large linear systems . To solve the large linear systems , we introduce the Multigrid method.The second,we applied the S.C. Brenner and L. Y. SUNG's method as a standard. The advantage of S.C. Brenner's method is that the stress intensity factors j can be represent by the simple expression and

correct calculation , but the disadvantage is that the lack of accuracy for the whole domain .We will improve the accuracy by introducing the adaptive mesh-re nement and adaptive cut-off function .

Finite element methods and their error estimation are given in section 1. Multigrid method are introduced in chapter 2. The Poisson equation and the singular funciton repre-sentation are given in chapter 3. Numerical results are carried out in chapter 4.

Chapter 1

Finite Element Method

1.1 Introduction of Finite Element Method

The basic idea in any numerical method for solving a differential equation is rst to dis-cretize given continuous problem with in nite degrees of freedom to a discrete problem or with only nite degrees of freedom such that the differential equation is transformed into a system of linear equations which can be solved by using a computer.

Finite element method start from a reformulation of a given differential equation as an equivalent variational problem. In the case of elliptic equations this variational problem in basic case is a minimization problem of the form

F ind u_{2 V such that F (u) 5 F (v) for all v 2 V} (1.1)

where V is a given set of admissible functions and F : V ! R is a functional. F (v) is the total energy associated with v and (1:1) corresponds to an equivalent characterization of the solution of the differential equation as the function in V that minimizes the total energy of the considered system. In general the dimension of V is in nite and thus in general the problem (1:1) can't be solved exactly. To obtain a problem that can be solved on a computer the idea in the nite element method is to replace V by a set Vh consisting

of simple function only depending on nitely many parameters. This leads to a

1.1 Introduction of Finite Element Method 4

dimensional minimization problem of the form:

F ind uh 2 Vh such that F (uh)5 F (v) for all v 2 Vh (1.2)

This problem is equivalent to a linear or nonlinear system of equations. We hope that the solution uhof this problem is suf ciently good approximation of the solution of the original

minimization problem (1:1). Usually one chooses Vh to be a subset of V and in this case

(1:2)corresponds to the classical Ritz-Galerkin method.

The advantage of nite element methods as compared with nite difference meth-ods is that complicated geometry,general boundary conditions and variable or non-linear material properties can be handled relatively easily.In all these cases one meets unnecces-sary arti cial complications with nite difference methodology.Further,the nite element method has a solid theoretical foundation which gives added reliability and in many cases makes it possible to mathematically analyze and estimate the error in the approximate nite element solution.

To solve a given differential or integral equation approximately using the nite ele-ment method, one has to go through basically the following steps:

1. variational formulation of the given problem 2. Mesh Generator

3. discretization using FEM: construction of the nite dimensional space Vh and choose

basis function

1.2 Variational Formulation 5

From step 1~4, one obtain a linear systems we will introduce at section 1.2.

1.2 Variational Formulation

We will now consider the following boundary value problem for the Poisson equation:

u = f in

u = g on @ (1.3)

where _{is a bounded domain in the R}2 _{with boundary @ ,g is a constant, f is a given}

function, where u = @ 2_u @x2 + @2_u @y2 (1.4)

the equivalent variational problem is Z

( u)vdx = Z

f vdx (1.5)

where v is test function in H1

0( ), vj@ = 0.By taking integration by parts ,

Z

OvOudx =Z vf dx + !n Ouvj@ =

Z

vf dx (1.6)

1.2.1 Existence and Uniqueness of Solution

Let a:V V _{! R be a bilinear mapping with following properties:} (1) a(:; :) is symmetric

(2) (Continuity) a(:; :) is continuous,ie,there is a constant > 0 such that ja(v; w)j jjvjjVjjwjjV 8v; w 2 V

1.3 Finite Element Discretization 6

(3) (Coercivity) a(:; :) is V-elliptic,ie,there is a constant > 0such that a(v; v) jjvjj2

V 8v 2 V

(4) L is continuous,ie, there is a constant > 0such that jL(v)j _jjvjjV 8v 2 V

Theorem (Lax-Milgram theorem) Let V be a Hilbert space with scalar product (:; :)V

and corresponding norm jj jjV(the V norm):Suppose that a(:; :) is a bilinear form on

V V and L a linear form on V such that under the assumptions (1)-(4), there exists a unique u 2 V such that

a(u; v) = L(v); for all v 2 V

1.3 Finite Element Discretization

Let Th =fKg is a triangulation of R2, the integral equation can be rewritten as

X K2Th( ) Z K OvOudx = X K2Th( ) Z K vf dx (1.7)

The nite element method is then employed to discretize the termsR OvOudx and R vfdx on element, we rst look the germetry on an element.

The geometry of the 3-node triangle is speci ed by the location of its three corner nodes on the fx; yg plane. The nodes are labeled 1, 2, 3 while traversing the sides in counterclockwise fashion. The location of the corners is de ned by their coordinates:

1.3 Finite Element Discretization 7

the area of triangle is denoted by A and is given by:

2A = det 2 4 x11 x12 x13 y1 y2 y3 3 5 = x21y31 x31y21 where xij = xi xj; yij = yi yj for i; j = 1; 2; 3 i 6= j: Fig1.3.1

1.3.1 Linear Interpolation

One can choose a piecewise ploynomial function to approximate the exact solution u and the test function v .For example , if one choose linear piecewise funciton ,then the function u (x; y)may be expressed as

u (x; y) = a0+ a1x + a2y (1.8)

where a0, a1 and a2 are coef cients to be determined from three conditions. In nite

element work such conditions are often the nodal values taken by u at the corners:

1.3 Finite Element Discretization 8

The expression in triangular coordinates makes direct use of these three values:

u ( ; ) = u1(1 ) + u2 + u3 = u1 u2 u3 2 4 1 3 5 = 1 2 4 uu12 u3 3 5 (1.9)

equation (1:9) is called a linear interpolant for u:

1.3.2 Coordinate Transformations

Consider triangular on regular triangular, points of the triangle may also be located in terms of a parametric coordinate system ;

Fig1.3.2

Cartesian coordinates and triangular coordinates are linked by the relation x y = x1 y1 [1 ] + x2 y2 + x3 y3 x21 x31 y21 y31 + x1 y1 (1.10) These simply apply the linear interpolant formula to the Cartesian coordinates: x = x1(1 ) + x2 + x3 and y = y1(1 ) + y2 + y3 :

1.3 Finite Element Discretization 9 Inversion of (1:10) yields = 1 2A y31 x13 y12 x21 x y 1 2A y31 x13 y12 x21 x1 y1 (1.11)

1.3.3 Linear FEM Discretization

After generating mesh ej;we shall now construct a nite dimensional subspace Vh of the

space V de ned above consisting of piecewise linear function. We now let Vh be the set of

functions v such that v is linear on each subinterval ej, v is continuous on domain and

v = 0on @ :We observe that Vh V: As parameter to describe a function uj = v (xj)

we may choose the values uj = v (xj) at the node points xj; j = 0; :::; m + 1: Let us

introduction the basis function _{j 2 V}h; j = 0; :::; m + 1:de ned by

j(xi) =

1if i = j

0if i 6= j; i; j = 1; :::; M (1.12)

Fig1.3.3

i.e., _j is the continuous piecewise linear function that take the value 1 at node point xj and the value 0 at other node points. A function v 2 Vhthen has the representation

v (x) =

m

X

i=1

1.3 Finite Element Discretization 10

where uj = v (xj), i.e., each uj = v (xj) can be written in a unique way as a linear

combination of the basis function _j: In particular it follow that Vh is a linear space of

dimension m with basis j m i=1:

1.3.4 Partial Derivatives

From equations (1:10) and (1:11) we immediately obtain the following relations between partial derivatives: @x @ = x21; @x @ = x31; @y @ = y21; @y @ = y31: @ @x = y31 2A; @ @x = y12 2A; @ @y = x13 2A; @ @y = x21 2A: The derivatives of function u ( ; ) is

ru = @u @x @u @y = " _@u @ @ @x + @u @ @ @x @u @ @ @y + @u @ @ @y # = 1 2A @u @ y31+ @u @ y12 @u @ x13+ @u @ x21 = 1 2A y31 y12 x13 x21 @u @ @u @ = 1 2A @u @ @u @ (1.14) ,where = y31 y12 x13 x21 :

Similarly,the derivatives of test function v( ; ) is rv = @v @x @v @y = 1 2A y31 y12 x13 x21 @v @ @v @ = 1 2A @v @ @v @ T _(1.15) So,we let

1.3 Finite Element Discretization 11 u = 3 X i=1 ui i ,v = 3 X i=1 vi i (1.16) where 1 = 1 ; 2 = ; 3 = ; @u @ @u @ = " _@ 1 @ @ 2 @ @ 3 @ @ ₁ @ @ ₂ @ @ ₃ @ # 2 4 uu12 u3 3 5 = 1 1 0 1 0 1 2 4 uu12 u3 3 5 = D 2 4 uu12 u3 3 5 (1.17) where D = 1 1 0 1 0 1 : Similarity, @v @ @v @ = v1 v2 v3 2 6 4 @ 1 @ @ 1 @ @ ₂ @ @ ₂ @ @ 3 @ @ 3 @ 3 7 5 = v1 v2 v3 DT:

for one element,we get Z K rvrudx = Z b K v1 v2 v3 DT T 1 4A2 D 2 4 uu12 u3 3 5 jJjd d (L:H:S) (1.18) KK = R b KD T T 1

4A2 DjJjd d is called element stiffness matrix.

and Z K vf dx = Z b K v1 v2 v3 2 4 1₂ 3 3 5 _{1 2 3} 2 4 ff12 f3 3 5 jJjd d (R:H:S) (1.19) MK = R b K 2 4 1₂ 3 3

1.4 Error Estimation 12

1.3.5 Assembling the element matrix

Assembling all the element matrices KK; Mk into global matrices,here !u ; !v ;!fare the

global nodal vector.P_K2T

h( )KK = K;

P

K2Th( )MK = M. The equation (1:7) becomes

!_{v K !u = !v M} !_{f ; for all v 2 Vh} As a result,we obtain a linear system

K !u = M !f (1.20)

1.4 Error Estimation

(Cea Lemma)

Let u 2 V be the solution of a(u; v) = L(v) ; 8v 2 V and uh 2 Vh V:Then

jju uhjjV infjju vhjjV 8vh 2 Vh (1.21)

where is coercivity scalar and is continuity scalar.

For a typical elliptic problem satisfying the conditions (1)-(4) in the section 1.2.1,we choosing vh = hu 2 Vh to be a suitable interpolant of u and estimate the interpolation

error jju hujjV:

1.4.1 Interpolation Error with Piecewise Linear Functions

We rst consider the case V = H1_{( ) and V}

h = fvh 2 V : vhjK 2 P1(K);8K 2 Thg

where Th =fKg is a triangulation of R2,ie,Vh is the standard nite element space of

1.4 Error Estimation 13

hk =the dimeter of K=the longest side of K; k =the dimeter of the circle inscribed in K;

h = max

K2Th

hK:

We shall assume that there is a positive constant independent of h ,such that

K

hK 8K 2 T

h (1.22)

This condition means that the angles of the triangles K are not allowed to be arbitrarily small;the constant is a measure of the smallest angle in any K 2 Th:Let Ni; i = 1; :::; M;

be the nodes of Th:Given u 2 C0( )we de ne the interpolant

hu(Ni) = u(Ni) i = 1; :::; M (1.23)

,thus huis the piecewise linear function agreeing with u at the nodes of Th:

1.4 Error Estimation 14

For j=1,2 and x2 K we have

3 X i=1 i(x) = 1; (1.24) 3 X i=1 pi(x) i(x) = 0; (1.25) 3 X i=1 @ @xj i (x) = @ @xj 3 X i=1 i(x) = 0; (1.26) 3 X i=1 pi(x) @ _i @xj (x) = @v @xj (x): (1.27)

Theorem 1 Let K 2 Th be a triangle with vertices ai,i=1,2,3.Given v 2 C0(K);let the

interpolant v 2 P1(K)be de ned by v(ai) = v(ai); i = 1; 2; 3: (1.28) then jjv v_jjL1_(K) 2h2_Kmax j =2jjjD vjjL 1_(K); (1.29) max j =1jjjD (v v)jjL1(K) 6 h2 K K max j =2jjjD vjjL1(K); (1.30) where jjvjjL1(k) = max x2K jv(x)j: (1.31)

Proof. Let i; i = 1; 2; 3;be the basis functions for P1(K):A general function w 2 P1(K)

then has the representation

w(x) = 3 X i=1 w(ai) _i(x); x_{2 K;} (1.32) so that in particular v(x) = 3 X i=1 v(ai) _i(x); x_{2 K;} (1.33)

1.4 Error Estimation 15

since v(ai_{) = v(a}i_{);we now using the Taylor expansion at x 2 K :}

v(y) = v(x) + 2 X j=1 @v @xj (x)(yj xj) + R(x; y); (1.34) where R(x; y) = 1 2 2 X i;j=1 @2_v @xi@xj ( )(yi xi)(yj xj); (1.35)

is the remainder term of order 2 and is a point on line segment between x and y.In particular by choosing y = ai_{;we have}

v(ai) = v(x) + pi(x) + Ri(x) (1.36) where pi(x) = 2 X j=1 @v @xj (x)(ai_j xj); ai = (ai1; a i 2); Ri(x) = R(x; ai): (1.37) Since jaij xjj hK; i = 1; 2; 3; j = 1; 2; (1.38)

the estimate of the remainder term Ri(x)

Ri(x) 2h2Kmax

j j=2jjD vjjL1(K); i = 1; 2; 3: (1.39)

Now (1:33) and (1:36) combine to give v(x) = v(x) 3 X i=1 i(x) + 3 X i=1 pi(x) i(x) + 3 X i=1 Ri(x) i(x); x2 K: (1.40)

By(1:24) ,(1:25) in Lamma 2 and(1:40) ,we have v(x) = v(x) +

3

X

i=1

1.4 Error Estimation 16

which gives us the following representation of the interpolation error: v(x) v(x) =

3

X

i=1

Ri(x) i(x): (1.42)

Since 0 i(x) 1;if x 2 K; i = 1; 2; 3;we can use the previous estimate (1:39) of the

remainder term Ri to get

jv(x) v(x)_j 3 X i=1 jRi(x)jj i(x)j max i jRi(x)j 3 X i=1 j i(x)j 2h 2 Kmax j j=2jjD vjjL1(K); x2 K:(1.43) we proves (1:29) :

Proof. To prove (1:30) we differentiate (1:33) with respect to x1to get

@ v @x1 (x) = 3 X i=1 v(ai)@ i @x1 (x); (1.44)

which together with (1:36) shows that @ v @x1 (x) = v(x) 3 X i=1 @ _i @x1 (x) + 3 X i=1 pi(x) @ _i @x1 (x) + 3 X i=1 Ri(x) @ _i @x1 (x): (1.45) Hence,by (1:26) and (1:27) we have

@ v @x1 (x) = @v @x1 (x) + 3 X i=1 Ri(x) @ _i @x1 (x); (1.46)

which gives the following representation o the error @v @x1 @ v @x1 : @v @x1 (x) @ v @x1 (x) = 3 X i=1 Ri(x) @ _i @x1 (x); x_{2 K:} (1.47) and max x2K j @ _i @x1 (x)_j 1 K ; (1.48)

which together with (1:39) nally gives j_@x@v 1 (x) @ v @x1 (x)_j 6h 2 K K max j =2jjjD vjjL 1_(K): (1.49)

1.4 Error Estimation 17

In the same way we estimate @v @x2

@ v

@x2 and thus (1:30) follows.The proof of the theorem is

now complete once the lemma is established.

Theorem 2 Under the assumptions of Theorem 1 there is an absolute constant C such that jjv v_jjL2(K) Ch 2 KjvjH2_(K); (1.50) jv v_jH1_(K) C h2 K K jvjH2_(K): (1.51)

Theorem 3 The global interpolation errors

jju hujjH1_{( )} Chjuj_H2_{( )}; (1.52)

jju hujjL2_{( )} Ch2juj_H2_{( )} (1.53)

Proof. We have by summing over K 2 Th;

jju hujj2L2_{( )} = X K2Th jju hujj2L2_(K) X K2Th C2h4_Kjuj2_H2_(K) C2h4 X K2Th juj2H2_(K) = C2h4juj2_H2_{( )}; (1.54)

and similarly using (1:22) ,ie, hK K 1 ; jju hujj2H1_{( )} X K2Th C2h 4 K 2 K juj2H2_(K) X K2Th C2h 2 K 2juj 2 H2_(K) C2h 2 2juj 2 H2_{( )} (1.55) so that jju hujjH1_{( )} Ch jujH2_{( )} = Chjuj_H2_{( )}; (1.56)

1.4 Error Estimation 18

if the constant is included in the constant C;and

jju hujjL2_{( )} Ch2juj_H2_{( )}: (1.57)

1.4.2 A Priori Error Estimation

The bilinear form a(u; v) = R OuOvdx, the nite element space of piecewise linear ele-ment is H1_{( ) = V, the test function space is H}1

0 ( ). The priori error estimate of Possion

equation

jju uhjjH1_{( )} Chjuj_H2_{( )}

Proof. by Céa Lamme, so the inequality is

jju uhjjH1_{( )} infjju v_hjj_H1_{( )} for all v_h 2 V_h

where is coercivity scalar and is continuity scalar ,than we choose vh = hu, where huis interpolant by the piecewise linear basis function

inf_jju vhjjH1_{( )} jju _hujj_H1_{( )}

by interpolation error estimate

1.4 Error Estimation 19

1.4.3 A Posterior Error Estimation and Adaptive Mesh-Re nement

Techniques

On the other way to re ne the mesh ,we introduce a posterior error estimation and the adaptive mesh-re nement techniques. Before presenting the a posteriori error estimators, we introduce some notations.

For any open subset ! of with Lipschitz boundary ,we denote by L2_(!),Hk_(!);

and L2_{( ),k 1, the standard Lebesque and Sobolev spaces ,respectively,equipped with}

the norms jj jj0;! :=jj jjL2_(!), jj jj_k;! :=jj jj_Hk_(!), jj jj_0; :=jj jj_L2_{( )}.

For T 2 Thwe denote by E(T ) and N(T ) the set of its edges and vertices,respectively.Let

Eh := [ T 2Th E(T ); Nh := [ T 2Th N (T );

be the set of all edges and vertices,respectively,in the triangulation. We split Eh and Nh in

the form Eh = Eh; [ Eh;D; Nh = Nh; [ Nh;D; with Eh;D : = fE 2 Eh : E Dg Nh;D : = fx 2 Nh : x Dg

For T 2 Thand E 2 Eh, we denote by hT and hEtheir diamter and length,respectively.

For T 2 Thand E 2 Ehand x 2 Nh ,let

!T := [ E(T )\E(T0₎₆₌₀ T0; !E := [ E2E(T0₎ T0; !x := [ x2N(T0₎ (T0);

1.4 Error Estimation 20

Figure1.4.2: Domain !T; !E and !x

and put VT : =f 2 C(T ) : 2 Y max(k+1;3);2; jE 2 Y k+1;1;8E 2 E(T ) (1.58) (x) = 0;_{8x 2 N(T ); j}E = 0;8E 2 E(T ) \ Eh;Dg; Vx :=f 2 C(!x) : jT0 2 VT0;8T0 !xg: (1.59)

Given E 2 Eh; and 2 L2(!E)with jT0 2 C(T0);8T0 !E;we denote by [ ]E the

jump of across E in an arbitrary,but xed direction.Put RE(uh) := ( [@uh @n]E; 8E 2 Eh; ; 0; _{8E 2 E}h;D: (1.60) and RT := Y T f +4uh; 8T 2 Th; (1.61)

where Q_T is the L2 _{projectors of L}2_{( ) onto the space of piecewise constant functions}

with respect to Th:

The rst error estimator simply is a weighted combination of the residuals RT(uh)

and RE(uh):It is given by

T;R :=fh2TjjRT(uh)jj20;T +

X

E2E(T )nEh;D

hEjjRE(uh)jj20;Eg1=2; 8T 2 Th (1.62)

The second error estimator is based on the solution of local Dirichlet problems. For any x 2 Nh; , it is given by

1.5 Finite Element Approximation for Singular Funcitons 21

For more detialed proof of the lower bound and upper bound for the error estimators is in [17].

According to the error estimator T we obtained above, we want to decide a re

ne-ment of T; there is a strategy that we re ne all T with _T max_{T 2T}

h T:Here 0 1

is a given threshold,typically = 0:5. This strategy is very cheap and often satisfactory results.

Having decided which elements should be re ned, we re ne them by connecting the mid-points of their edges.Triangles may be cut into four new ones by connecting the mid-points of their edges .

Figure1.4.3: triangles with adaptive mesh-re nement

1.5 Finite Element Approximation for Singular Funcitons

The problem on hand is to devise singular elements within which the eld variable would vary as Rt_{(0 t 1)where R is the radial distance from the point where the derivative}

of the function would have a singularity of the type Rt 1_{.For example,let (x}

1; y1)be the

1.5 Finite Element Approximation for Singular Funcitons 22

(y y1)2]1=2:The singularity may be introduced in the geometric transformation between

(x; y)and ( ; ) systems.

1.5.1 Lagrange Interpolation

For one-dimensional domain, consider a straight line along x-axis in a Cartesain reference, such that x0 x xM:Consider the dimensionless coordinate = (x x0)=(xM x0)

such that 0 1:A mth-order polynomial F(m)_{( )can be represented by interpolating}

the discrete values of the function at (m + 1) equidistant points _i[i = 0; :::; m; _i(i=m)] using the Lagrange interpolation function, L(m)

i ( ), as F(m)( ) = m X i=o FiL (m) i ( ) (1.64) where L(m)_i ( ) = m Y j=0;j6=i ( j i j ): (1.65)

It is seen that the functions in (1:65) have the property that

L(m)_i ( _k) = ik (1.66)

where ik = 0 (i6= k) and ik = 1(i = k):

For a simple example, let we consider second-order polynomial F(2)_{( )and the}

trans-formation reduces to

1.5 Finite Element Approximation for Singular Funcitons 23

(x x0) = (xM x0) (1.67)

the Lagrange function can be represent as

L(2)₀ = 2( 1 2)( 1) L(2)₁ = 4 ( 1) L(2)₂ = 2 ( 1 2)

1.5.2 Singular Element

Apart from problems of fracture,at a given point in the domain,say P,the solution funciton F,whould have a singularity in its derivative.If R is the radial length measured from the point P,the solution function F would behave as R (0 < < 1);so that the rst derivative of F at P would be in nite. Such singularities may arise at re-entrant corners, tips of sharp crack,etc.

For One-dimensional cases,consider a straight line along the x-axis,x0 x xM:Let

the geometry of this line be described by an rth-order polynomial as x = r X i=o xiL (r) i ( ) (1.68)

where xi = x( i) and i = (i=r)(i = 0; 1; :::; r) are (r + 1) equidistant points in the

parametric domain 0 1:(xi themselves may not necessarily be equidistant).Also let

the dependent variable,say u,be represented by an mth-order polynomial. u( ) =

m

X

j=0

1.5 Finite Element Approximation for Singular Funcitons 24

where uj = u( j); but now j = (j=m)(j = 0; :::; m) are (m + 1) equidistant points in

0 1:

If the geometric transformation is de ned by

(x x0) = (xM x0) t (1.70)

where t is a positive integer greater than unity, the function u(x) corresponding to u( ) of Eq.(1:69) would be an mth-order polynomial in the variable (x x0)1=t or R1=t,and

hence the derivative (@u=@x) has a term which varies as R(1 t)=t_{,for t > 1,would have a}

singularity at x = x0:

For instance when r = 2 and t = 2; the positioning of nodes at x = x0; x =

x0+ (xM x0)(1₂)2 = x0+ (xM x0)(1₄), and x = xM:it would lead to x 2 and hence

@u=@xwould have a singularity of the type R 1=2:

In general,singularities in (@u=@x) of the type R(1 t)=t_{; t r(t and r are integer),can}

be created by suitably choosing the valus xiin an rth-order geometric transformation of the

type (1:68) :

For two-dimensional domain, consider a standard dimensionless triangle (0 ; 1)as shown in Fig1.5.2. The triangle P1P2P3with Cartesain coordinate (xi; yi)(i = 1; 2; 3)

is then mapped into the standard triangle ,through the relations:

(x x1) = (x2 x1) + (x3 x1)

1.5 Finite Element Approximation for Singular Funcitons 25

Fig1.5.2: Triangular element in (x; y) and ( ; ) coordinate systems

Consider (m + 1) equidistant point along the line = 0 and = 0;respectively,such that _i = (i=m)and _j = (j=m) (i; j = 0; 1; :::; m):A mth-order polynomial (m being an integer) , F(m)_{( ; );can be represented by interpolating the values of the function in the}

standard triangle,as: F(m)( ; ) = m X i=0 m i X j=0 FijL (m) ij ( ; ) (1.72) where L(m)_ij = "_{i 1} Y r=0 r i r # "_{j 1} Y s=0 s j s # 2 4 m (i+j+1) Y t=0 1 ( + ) (t=m) 1 ( _i+ _j) (t=m) 3 5 : (1.73) The function also have the property that

L(m)_ij ( _k; _l) = ik jl: (1.74)

For the detail of singular element geometric transformation in two-dimensional do-main, see the reference [8].

Chapter 2

Multigrid Method

2.1 Introduction of Multigrid Method

Multigrid methods(MG) are among the most ef cient methods of solving the linear sys-tems arising from discretization of eilliptic partial differential equations. There has been intensive research on the convergence of MG since it was introduced by Fedorenko. For symmetric positive-de nite eilliptic problems, thanks to many researchers, such as Bank, Brandt, Braess, Bramble, Dupont, Hackbusch, Mandel and McCormick, etc, the conver-gence theory has matured. The major ingredients for converconver-gence analysis of MG are called the approximation property and the smoothing property. One approach for con-vergence analysis is the so-called compact perturbation technique, which relies on a strong approximation property and treats the lower order terms as a small perturbation of the symmetric positive de ne term. The technique has been successfully applied to diffu-sion problems and Bramble, Pasciak, Wang, Xu have shown robust MG uniform conver-gence. In these studies, uniform convergence of MG can be established with one step of standard Jacobi or Gauss-Seidel smoothing even without regularity assumptions. Another approach requires a strong smoothing property to compensate for poor approximation prop-erty. In this direction, it is very important to nd a robust smoother.Robust smoothers such as the block Jacobi , block Gauss-Seidel method and the incomplete LU factorization (ILU) method are commonly used.

2.2 Relaxation Process 27

The ef ciency of the multigrid algoritthm is achieved from an elegant combination of the smoothing procedure and the coarse grid correction procedure. The smoothing pro-cedure plays the role of reducing highly oscillatory error modes, and the coarse grid is used to correct the remaining smooth error modes. Hackbusch and Braess give the rst rigor-ous proof on the multigrid convergence and identify that the smoothing property and the approximation property are the cornerstones for the convergence analysis of multigrid methods. The smoothing property achieved by stationary iterative method. The approxi-mation property is achieved by choosed intergrid interpolation operator.

In general,the stationary iterative methods can effective eliminating the high-frequency or oscillatory components of the error,but leave the low-frequency or smooth components of the error as we will shown in experiment 1 in section 2.6.This process is called the re-laxation process in MG. Moreover, the smooth error modes on a ne grid appear more oscillatory on the coarse grid. As a result, repeating the relaxation process may conse-quently remove various oscillartory modes of error. To achieve the ef ciency of the MG algorithm, we only have to concern how to restrict smooth error to the coarse grid and how to prolong the errors remained in coarse grid to the ne grid. This process is called the inter-grid interpolation, we will introduce the relaxation and interpolation process in the

llowing sections.

2.2 Relaxation Process

For relaxation process, we introduce the stationary iteration methods.Stationary iterative methods solve a linear system with an operator approximating the original one and based

2.2 Relaxation Process 28

on a measurement of the error in the result (the residual), form a "correction equation" for which this process is repeated. While these methods are simple to derive, implement, and analyze, convergence is only guaranteed for a limited class of matrices. Examples of sta-tionary iterative methods are the Jacobi method, Gauss–Seidel method and the Successive over-relaxation method.

2.2.1 Jacobi Mehtod

One of the methods to solve the linear system is the Jacobi method . The Jacobi method is derived by examining each of the N 1equations in the linear system Av = f in isolation .If in the i thequation

N_X1 j=1

ai;jvj = fi; (2.75)

we solve for the value of viwhile assuming the other entries of v remain xed ,we obtain

vi = (fi

X

j6=i

ai;jvj)=ai;i ; 1 j N 1 (2.76)

this suggests an iterative method de ned by v(k)_i = (fi X j6=i ai;jv (k 1) j )=ai;i ; 1 j N 1 (2.77)

which is the Jacobi method.

In matrix terms, the de nition of the Jacobi method can be expressed as

v(k)= D 1(L + U )v(k 1)+ D 1f; (2.78) where the matrices D; L and U represent the diagonal, the strictly lower-triangular, and the strictly upper-triangular parts of A,respectively.

2.2 Relaxation Process 29

The pseudocode for the Jacobi method is given by following Choose an initial guess v(0) _{to the solution u}

for k = 1; 2; for i = 1; 2; ; N 1 vi = 0 for j = 1; 2; ; i 1; i + 1; ; N 1 vi = vi+ ai;jv (k 1) j end vi = (fi vi)=ai;i end v(k)_{= v}

check convergence;continue if necessary end

There is a simple modi cation which can be made to the Jacobi iteration. v_j = (fi

X

j6=i

ai;jv(k 1)j )=ai;i ; 1 j N 1

However,v_j is now only an intermediate value. The new approximation is given by the weighted average

v(k)_j = (1 !)v(k 1)_j + !v_j ; 1 j N 1

where ! 2 R is a weighting factor which may be chosen.This iteration is called weighted Jacobi method.

2.2.2 Gauss-Seidel Method

Consider again the linear systems Av = f;now assume that the equations are examined one at a time in sequence,and that previously computed results are used as soon as they are available ,we obtain the Gauss-Seidel method:

v_i(k)= (fi X j<i ai;jv(k)j X j>i

2.2 Relaxation Process 30

In matrix form,

Dv(k) Lv(k) U v(k 1) = f

(D L)v(k)= U v(k 1)+ f

v(k)= (D L) 1U v(k 1)+ (D L) 1f

where the matrices D; L and U represent the diagonal, the strictly lower-triangular, and the strictly upper-triangular parts of A,respectively. The pseudocode for the Gauss-Seidel method is given by following

Choose an initial guess v(0) _{to the solution u}

for k = 1; 2; for i = 1; 2; ; N 1 vi = 0 for j = 1; 2; ; i 1 vi = vi+ ai;jv (k) j end for j = i + 1; ; N 1 vi = vi+ ai;jv (k 1) j end vi = (fi vi)=ai;i end v(k) _{= v}

check convergence;continue if necessary end

2.2.3 Successive Overrelaxation Method

The Successive Overrelaxation Method, or SOR ,is devised by applying extrapolation to the Gauss-Seidel method. This extrapolation takes the form of a weighted average between the

2.2 Relaxation Process 31

previous iterate and the computed Gauss-Seidel iterate successively for each component:

v_i(k) = !vi(k)+ (1 !)v (k 1)

i (2.80)

,wherevi(k) denotes a Gauss-Seidel iterate ,and ! is the extrapolation (weighting)

factor.The idea is to choose a value for ! that will accelerate the rate of convergence of the iterates to the solution.

In matrix form, the SOR is written as

Dv(k) = !Dv(k)+ (1 !)Dv(k 1)

= ![f + Lv(k)+ U v(k 1)] + (1 !)Dv(k 1)

=_{) (D} !L)v(k) = [!U + (1 !)D]v(k 1)+ !f

v(k) = (D !L) 1[!U + (1 !)D]v(k 1)+ (D !L) 1!f

where the matrices D; L and U represent the diagonal, the strictly lower-triangular, and the strictly upper-triangular parts of A,respectively.

If the extrapolation factor ! is choosing by one,the SOR method simpli es to the Gauss-Seidel method.If 0 < ! < 1,its called underrelaxation.If 1 < ! < 2, its called overrelaxation and SOR fails to converge when ! is outside the interval (0; 2).

2.3 Inter-grid Interpolation : Restriction and Prolongation 32

The pseudocode for the SOR method is given by following Choose an initial guess v(0) _{to the solution u}

for k = 1; 2; for i = 1; 2; ; N 1 = 0 for j = 1; 2; ; i 1 = + ai;jv (k) j end for j = i + 1; ; N 1 = + ai;jv (k 1) j end = (fi )=ai;i v_i(k)= v_i(k 1)+ !( v(k 1)_i ) end

check convergence;continue if necessary end

2.3 Inter-grid Interpolation : Restriction and Prolongation

First,we consider the linear prolongation, the operator will be denoted Ih

2h.It takes coarse

grid vectors and produces nd grid vectors according to the rule Ih

2hv2h = vh v_2jh = v2h_j vh_2j+1 = 1 2(v 2h j + v 2h j+1), 0 j N 2 1 (2.81)

2.3 Inter-grid Interpolation : Restriction and Prolongation 33

At even-numbered ne grid points ,the values of the vector are transferred directly from 2h_to h_{. At odd-numbered nd grid points,the value of the vector is the average of}

the adjacent coarse grid values.

For the case N = 8; this operator has the form

I_2hh v2h= 1 2 2 6 6 6 6 6 6 6 4 1 2 1 1 2 1 1 2 1 3 7 7 7 7 7 7 7 5 2 4 v1 v2 v3 3 5 2h = 2 6 6 6 6 6 6 6 4 v1 v2 v3 v4 v5 v6 v7 3 7 7 7 7 7 7 7 5 h = vh

The second ,we consider moving vectors from a ne grid to a coarse grid.There are called restriction operators and are denoted by I2h

h :It is de ned by Ih2hvh = v2h,where

v_j2h= 1 4(v h 2j 1+ 2v h 2j+ v h 2j+1); 1 j N 2 1: (2.82)

Figure 2.3.2: Restriction by full weighting of a nd grid vector to the coarse grid. the values of the coarse grid vector are a weighted average of values at neighboring nd grid points.

2.4 Multigrid Algorithm 34

In the case of N = 8,the operator has the form

I_h2hvh = 1 4 2 4 1 2 11 2 1 1 2 1 3 5 2 6 6 6 6 6 6 6 4 v1 v2 v3 v4 v5 v6 v7 3 7 7 7 7 7 7 7 5 h = 2 4 vv12 v3 3 5 2h = v2h:

2.4 Multigrid Algorithm

After kowning the procedure , we compose all the components aobve to make the correction scheme. First , relaxation on the ne grid will eliminate the oscillatory components of the error,leaving a relatively smooth error,then transfer the residual on the ne grid to the coarse grid an solve the residual equation exactly on the 2h_{.Since the error is smooth on}

the 2h_{,we can prolongation the error accurately back to the ne grid.Finally, relaxation on}

the ne grid use the better initial which is correction before. The procedure is given by following

vh

MG(vh_{; f}h₎

Relax 1times on Ahuh = fhon h with initial guess vh:

Compute r2h _{= I}2h

h (fh Ahvh):

Solve A2h_e2h _{= r}2h_on 2h_:

Correct ne grid approximation : vh

vh+ I_2hh e2h: Relax 2times on Ahuh = fhon h with initial guess vh:

Here,in practice, the number of relaxation times 1and 2 are often 1,2 or 3.

Moreover, the correction scheme is not just two level, it can be more level,by re-cursive within itself.The algorithm is called the V cycle.The de nition which is given

2.4 Multigrid Algorithm 35

by

V-cycle Scheme

vh _MVh(vh; fh)

1.Relax 1 times on Ahuh = fh with a given initial guess vh:

2.If h _{=coarsest grid, then go to 4.}

Else f2h I2h h (fh Ahvh): v2h_{= 0} v2h MV2h_(v2h_{; f}2h_): 3.Correct vh vh_{+ I}h 2hv2h:

4.Relax 2 times on Ahuh = fh with initial guess vh:

The V cycleis just one of a family of multigrid cycling schemes.The entire family is called the cyclemethod and is de ned by

-cycle Scheme vh

M h_(vh_{; f}h₎

1.Relax 1 times on Ahuh = fh with a given initial guess vh:

2.If h _{=coarsest grid, then go to 4.}

Else f2h I2h h (fh Ahvh): v2h_{= 0} v2h M 2h_(v2h_{; f}2h_{) times.} 3.Correct vh vh+ I_2hh v2h:

4.Relax 2 times on Ahuh = fh with initial guess vh:

If = 1,which gives the V cycleand = 2is called the W cycle:

Here are other ideas,if it starts on the coarsest discretization with an exact solver, the results are interpolated to the next ner grid with a few cycles (V or W) of the multigrid method are applied .The result is again interpolated to the next ner grid, where again a few cycles of multigrid method are applied.If this is used recursively, the so-called full multigridmethod.The algorithm with V cycleis given by

2.4 Multigrid Algorithm 36

Full Multigrid V-cycle Scheme vh _{F MV}h(vh; fh) Initialize fh_,f2h_,...;vf_,v2h_{,...to zero.}

Solve or relex on coarsest grid ... v4h v4h_{+ I}4h 8hv8h v4h MV4h_(v4h_{; f}4h₎ v2h v2h_{+ I}2h 4hv4h v2h MV2h_(v2h_{; f}2h₎ vh vh_{+ I}h 2hv2h vh MVh_(vh_{; f}h₎

2.6 Numerical Experiments 37

2.5 Complexity

There are various method to solve the linear systems,including stationary iterative methods or nonstationary iterative methods ,the following table gives us the complexity of various method :

Table 3.3

Complexity of Various Methods

Method 2D 3D

JM O(N2_{) O(N}5₃₎

GS O(N2_{) O(N}5₃₎

SOR O(N32) O(N 4 3)

CG O(N32) O(N 4 3)

Multigrid O(N log N) O(N log N)

2.6 Numerical Experiments

Experiment 1

In this experiments we consider the weighted Jacobi method with ! = 2=3 applied to the one-dimensional problem Au = 0 on a grid N=48. We use an initial guess,

vh_j = 1 2 sin( 12j N ) + sin( 30j N )

the results of this calculation are given in Figs 2.6.1(a)-(e). The initial guess is shown in Fig 2.6.1(a) . In Fig 2.6.1(b), the approximation vh_{after one relaxation sweep is superimposed}

on the initial guess. Much of the oscillatory component of the initial guess has already been removed. The maximum norm of the error has decreased signi cation. Fig. 2.6.1(c) shows the approximation after three relaxation sweeps, superimposed on the previous approxima-tions. Further relaxations on the ne grid would provide only a slow improvement at this

2.6 Numerical Experiments 38

point. This signals that it is time to move to the coarse grid.Fig. 2.6.1(d) shows the ne grid error after one relaxation sweep on the coase grid and the error after three coarse grid relaxation sweeps is shown in Fig. 2.6.1(e)

Figure 2.6.1(a): jjejj = 0:8536 Figure 2.6.1(b): jjejj = 0:4300

Figure 2.6.1(c): jjejj = 0:2607 Figure 2.6.1(d): jjejj = 0:1564

2.6 Numerical Experiments 39

Experiments 2

Consider the following two-point boundary value problem. It is given by the second-order ordinary differential equation

( (2 + cos( x))u0(x))0 = f (x); 0 < x < ; and are constant.

u(0) = 0; u( ) = 2 (2.83)

we use nite difference method to discrete the domain.The domain of the problem fx : 0 x _{g is partitioned into N subintervals by introducing the grid points x}j =

jh,where h = =N is the constant width of the subintervals. At each of the N 1 in-terior grid points, the differential equation is replaced by a second-order nite difference approximation. We also introduce vj as an approximation to the exact solution u(xj):

8 > > > > < > > > > : [2 + cos( x_j 1 2)]vj 1 [4 + cos( xj+ 1 2) + cos( xj 1 2)]vj h2 + [2 + cos( x_j+1 2)]vj+1 h2 = f (xj); 1 j N 1: v0 = 0; vN = 2 (2.84)

Let a represent the entry 4 + cos( x_j+1

2) + cos( xj 1

2); brepresent the entry

2 + cos( x_j+1

2); crepresent the entry 2 + cos( xj 1

2);the system can be represent

in matrix form as Av = f: 1 h2 2 6 6 6 6 4 a b c a b . .. ... ... c a b c a 3 7 7 7 7 5 2 6 6 6 6 4 v1 : : : vN 1 3 7 7 7 7 5= 2 6 6 6 6 4 f1 : : : fN 1 3 7 7 7 7 5

The matrix A is tridagonal,symmetric positive de nite and has dimension (N 1) (N 1):

We use multigrid method to solve the linear system.There are several comments in order.First, we choose and the constant one.Second, we use three level v cycleand

2.6 Numerical Experiments 40

the number of relaxation times 1 and 2 is the constant three.Third, we use v = x as an

initial guess.

The solution of this boundary value problem is u(x) = x2_{:The error e}h_{is de ned by}

eh ₌

jjuh _vh

jjL2,where uh is the exact solution and vh is the approximation on h:The

rate of convergence is computed by

= e 2h eh (2.85) Table2.6.1 N e 128 2.7315e-004 256 6.8301e-005 3.9992 512 1.7076e-005 3.9998

The results are tabulated in above .The theoretical convergence rate is 4 for : Experiments 3

In experiment two , we consider the boundary value problem in two-dimension given by

4 u = f in : [0; 1] [0; 1]

u = 0 on boundary (2.86)

and nite element method is used to discretize the domain.To solve the linear sys-tems Ax = b,we applied the Jacobi Method,Gauss-Seidel Method,Successive Overrelax-ation Method, Conjugate Gradient Method and Multigrid Method. Compare the iterOverrelax-ation numbers and the relative residual error between the multigrid method and other methods ,we can nd the advantages of using multigrid method as the linear systems solver.

2.6 Numerical Experiments 41

The exact solution is

u(x; y) = sin(2 x) sin(4 y) (2.87)

The tolerance is 10e-16.The error is de ned by

error =_jjuex bujj1 (2.88)

and the relative residue is de ned by

Relative residue = jjb Axjj

jjbjj (2.89)

2.6 Numerical Experiments 42

Figure3.6.1 show the exact solution and ve kinds of approximation.

Figure2.6.3: Iteration numbers of ve methods

Figure3.6.2 show us the relation between the number of iteration and the log relative residue.With ve kinds of iteration method, Multigrid Method gets the best performance.

Chapter 3

Research Method

Let _{be a bounded polygonal domain in R}2 _{with re-entrant angle.Consider the}

Pois-son equation with homogeneous boundary condition u = f in

u = 0 on @ (3.90)

where f 2 L2_{( ):}

When (3:90) is solved by the P1 nite element method on a quasi-uniform grid ,the

convergence rate in the energy norm is therefore of order O(h( =!) _{);where h is the mesh}

size of the triangulation and ! is the re-entrant angle.

In 1996, S.C.Brenner improves the convergence rate developed in [3].It is based on the full multigrid iteration technique and the following singular function representation of u [4, 5, 6, 7]

u = s + w (3.91)

where w 2 H2_{( )and the s are the singular function associated with the re-entrant angle.}

Note that the coef cient is known as stress intensity factors in elasticity problems. The multigrid method in [3] computes a solution of (3:91) in the form of

uh = hs + wh

where whis a piecewise linear function.It is shown in [3] that

ju uhjH1_{( )} Chjjfjj_L2_{( )}; (3.92)

3 Research Method 44

j hj Ch1+ =! jjfjjL2_{( ):} (3.93)

In 1997 ,S.C.Brenner and L.-Y.Sung [9] extend the results in [3] and [10] to the case of a polygonal domain with crack in R2_{(i.e. the re-entrant angle ! = 2 :)}

Figure3.1: A polygonal domain with cracks.

Let p be the vertice of such that the angle ! associated with p satis es ! > (i.e. the vertice F in the gure1.1).Let polar coordinates (r; ) be chosen at the vertex p so that the angle ! is spanned by the two rays = 0 and = !:

The singular function s is de ned by

s(r; ) = (r)r =!sin(

! ) (3.94)

where (r) is smooth cut-off function which equal 1 identically in a neighborhood of 0,and the supports of the is small enough so that the singular function s vanish identically on @ :

Then the solution u has the representation

3 Research Method 45

The stress intensity fractors can be expressed in terms of u by the following extraction formula = 1 Z f s 1dx + Z u_{4 s} 1dx ; (3.96)

where the dual singular function s 1 is de ned in the polar coordinate system (r; ) as

s 1(r; ) = (r)r =!sin(

! ) ([5, 7, 11, 12, 13, 14, 15, 16]).

There is a idea that we will take advantage of the singular function representation. We substitute the representation (3:95) into (3:90) to obtain the following boundary-value problem for w :

4 w = f + 4 s in

w = 0 on @ (3.97)

If the were known,we could solve (3:97) using piecewise linear nite element method.Unfortunately the is unknown,so we apply the nite element method on the kth level to the following

varational problem:

Find wk 2 H01( )such that

Z

r bwkrvdx =

Z

(f + k4 s) vdx 8v 2 H01( ) (3.98)

where the approximate stress intensity factors k is computed by the extraction formula

(3:96)using the approximate solution uk 1 obtained in the (k 1)stlevel, i.e., k =

1 Z

f s 1dx +

Z

uk 14 s 1dx : (3.99)

We obtain, on the kth level, a piecewise linear approximate solution wkto bwkby applying

3 Research Method 46

to u is de ned to be

uk= ks + wk: (3.100)

In other words we are really computing the regular part w of the solution. The improvement in the convergence rate is possible because w has better regularity than u:

The algorithm of S.C.Brenner's method is given by the following for l = 1:level for k = 1:m if l = 1 , k= 1 l;k = 0 else l;k = 1 R f s 1dx + R ul;k 14 s 1dx end R r bwl;krvdx = R (f + l;k4 s) vdx direct solve bwl;k ul;k= s + bwl;k end end

Although we known that the singular function can replace the solution at the sigular point well, but there still have some problem we can observe.After leaving the singular point, the singular function is now not a correct solution ,so the error occured on this cut-off region.To improve this problem, we introduce two strategys,adaptive mesh-re nement techniques and adaptive cut-off function.

For the rst strategy,adaptive mesh-re nement techniques, there is a relationship be-tween the converge rate and the error estimator,we should nd suitable order of h for the error estimator.The converge rate in the ragular domain is

jju uhjjL2_{( )} Ch2juj_H2_{( )}

3 Research Method 47

the order of error estimator is choose by

E(K) = _jjh(f au)_jjK + ( 1 2 X 2@K h [n (c_ruh)]2)1=2 (3.102)

Now,in the sigular domain ,we also want to nd the same relationship between the converge rate and the error estimator. Since the converge rate in singular domain is

jju uhjjL2_{( )} Ch(3=2) juj_H2_{( )}

jju uhjjH1_{( )} Ch(1=2) juj_H2_{( )}; (3.103)

the order of error estmator is choose by

E(K) = _jjh(3=4)(f au)_jjK+ ( 1 2 X 2@K h(1=4)[n (c_ruh)]2)1=2 (3.104)

we use this order as the rules of adaptive mesh-re nement techniques.

For the second strategy,adaptive cut-off function .Since we realize that the singular function replaced the solution successfully only near the singular point,so the range of cut-off function may affect the error.If we choose the range of cut-cut-off function is too wide, most solution replace by the singular function , that is not a correct .Otherwise, if we choose the range of cut-off functions is too small, the error still big near the singular point.So, this adaptive cut-off function strategy is nd the suitable cut range and move forward slowly .On the one hand, the error will not so big . On the other hand, the solution replace by the singular function more precision with respect to the mesh-re nement steps.

The adaptive cut-off funciton is de ned by

3 Research Method 48

nd the coef cient a0; a1;a2; a3; a4; a5 which is satisfy

8 < : (_3+i1 ) = 1 , (_3+i2 ) = 0 0₍ 1 3+i) = 0 , 0₍ 2 3+i) = 0 00₍ 1 3+i) = 0 , 00₍ 2 3+i) = 0 (3.106) where 1 3+i r 2

3+i; i = 0; 1; 2 is the mesh-re nement steps.

The purpose of this range we choose is that we want the interval of cut range is reduce and close to singular point,also the cut range is move forward slowly.

From the two strategy above , we can hold the stress intensity factorys and also improve the accuracy of the global error.The results shown in the experiment 4 in chapter 4.

Chapter 4

Numerical Results

In this chapter we report the result of some numerical experiments. The rst ex-periment we show the error arise from the singular point associated with the re-entrant angle.We consider the following boundary value problem:

u = f in

u = 0 on @ (4.107)

where is the circle with four angle ! = =2; 3 =2; 7 =4 and =0:51, ! is the maximun re-entrant angle.

4:1a 4:1b

4:1c 4:1d

Figure4.1:Circle with four kinds angle. Let the exact solution u be

u = (1 r2)r sin( ) (4.108)

4 Numerical Results 50

where 0 < r < 1; 0 < < !; =

!; ! is the maximun re-entrant angle.When is solved by piecewise linear nite element method on a quasi-uniform grid, the error is shown in the following:

4:2a 4:2b

4:2c 4:2d

Figure4.2:Error with four kinds of angle

Table 4.1 Results for four kinds of angle ! = 2 ! = 3 2 ! = 7 4 ! = 0:51 jju uhjjH1_{( )} 2.6377E-02 4.1243E-02 7.0759E-02 9.6003E-02

# of points 2205 6469 7791 8323

we can see the Fig4.2a is not a re-entrant angle , so the error of Fig4.2a didn't have peak at the origin,but the others have.We also can nd when the maximun re-entrant angle ! arise, the error also arise with the !:

The second experiment we will concern about the mesh-re nement strategy.There are three different mesh-re nement strategy we used. The rst strategy is uniform mesh

4 Numerical Results 51

with uniformly re nement,the second is singular element mesh with uniformly re nement and the last is uniform mesh with adaptive short cut-region re nement.

The rst mesh-re nement strategy is like as we re ne mesh as usual.The second strategy , we changed the uniform mesh to the singular mesh. According to our mentioned in the singular element ,here the approach is like singular element,we put the exponential grid points (x; y) in the way

(x; y) = (2 i_{; 0)}

(x; y) = (2 i_{cos !; 2} i_{sin !)} ; i = 1; :::; 10 (4.109)

on the two rays.This may couse many grid points located near the singular point.The third strategy is using uniform mesh, but not re ne the whole domain, only re ne the region with radius r = 0:5j _{form the center points (0; 0) , where j represent the times of re nement.}

We consider the maximun re-entrant angle ! =

0:51 case with three different mesh-re nement strategys .The times of mesh-mesh-re nement steps is 2, means that we have the original corasest domain 4h_{; nd domain} 2h_and h_:

Type A Type B Type C

Figure4.3: Error of three types of mesh-re nement we caculate the H1-norm error and L2-norm in the following table

Table4.1

Type A Type B Type C jju uhjjH1_{( )} 9.6003E-02 2.9259E-02 1.3185E-01

jju uhjjL2_{( )} 3.1643E-03 2.3718E-04 4.9345E-03

4 Numerical Results 52

From the Table4.1, we discuss some phenomena. In the second strategy, there are many points around of center.When we re ne the mesh, the points around the center of all will be re ned.Although we reduced the error, but we also pay a high price since too many points cause the large matrix systems.In the third strategy, we only re ne the points inside the cut-region, the error outside the cut-region still not be reduced, therefore, the accuracy performance is not good. The convergence rate for type A and type B in the H1 _normis

therefore of order O(h( =!) _{):(i:e:the theoratical convergence rate is 1.42)}

Table4.2

Type A ratio A Type B ratio B Type C ratio C

4h _{2.3184E-01 1.1065E-01 2.1384E-01}

2h _{1.4147E-01 1.5116 5.6813E-02 1.9476 1.6135E-01 1.3253} h _{9.6003E-02 1.4737 2.9259E-02 1.9417 1.3185E-01 1.2237}

Simultaneously, we apply the multigrid method for solving the linear systems and observed the bene ts for three dif erent mesh-re nement. For the multigrid parameter, we apply V-cycle 2th level iteration and we choose the weight-Jacobi method to relaxation, the max-imum relaxation number is 2.For the restriction operator ,we use full weighting as a restric-tion operator,and for the interpolarestric-tion operator, we consider linear interpolarestric-tion.

Type A Type B Type C

Figure4.4: Multigrid bene ts for three different mesh-re nement

From the Figure4.4 shown in above,we can observe that although the mesh size is dif erent, the slope almost not changed,meens that the converge rate of multigrid method is inden-pentent of mesh size.

4 Numerical Results 53

The third experiment we use the method with S.C.Brenner and L.-Y.Sung.The singu-lar function for is

s(r; ) = (r)r =0:51sin(

0:51 ) (4.110)

where the cut-off function (r) is de ned to be 8 < : 1; 0 r 1₃ 1458r5 _{+ 3645r}4 _3510r3_{+ 1620r}2 _{360r + 32;}1 3 r 2 3 0;2₃ r (4.111) The kare obtained by the extraction formula

k = 1 Z f s 1dx + Z uk 14 s 1dx : (4.112) where s 1(r; ) = (r)r =0:51sin( 0:51 ): (4.113)

Here,for the approximate stress intensity factors k,we choose 3-th level (l = 3) and

itera-tion 5 times (m = 5)

The error is shown in the following

4h 2h h

Figure4.5: Error of three level mesh The theoretical number for the stress intensity factorys is 1.

Table4.3

jju uhjjH1_{( )} # of points

4h _{9.9399E-01 2.3664E-01 556} 2h _{9.9815E-01 1.5423E-01 2128}

h _{9.9944E-01 1.0451E-01 8323}

The fourth experiment we improve the method by including adaptive mesh-re nement techniques and adaptive cut-off function.The following table shown the stress intensity

fac-4 Numerical Results 54

torys and the error in H1 _norm:

Table4.4 jju uhjjH1_{( )} # of points 1 9.9399E-01 2.3664E-01 556 2 1.0017E+01 2.0107E-01 913 3 9.9098E-01 1.7432E-01 1359 4 9.8648E-01 1.4946E-01 1842 5 9.8810E-01 1.3034E-01 2196 6 9.8722E-01 1.2721E-01 2747 7 9.8841E-01 1.1381E-01 3207 8 9.9204E-01 1.0953E-01 3533 9 9.9492E-01 1.0199E-01 4522 10 9.9341E-01 9.8286E-02 5494 11 9.9306E-01 9.7071E-02 6178 12 9.9560E-01 9.5075E-02 7217 13 9.9534E-01 9.5117E-02 7708 14 9.9596E-01 9.3322E-02 9533 15 9.9626E-01 9.3082E-02 10568

we can clearly compare the jju uhjjH1_{( )} with these two results above.Fixed points on

both sides of almost equal,we can clearly see the improvement in error in H1 _norm:

Table4.5

S.C.Brenner Experiment

jju uhjjH1_{( )} # of points jju u_hjj_H1_{( )} # of points

1 2.3664E-01 556 1 2.3664E-01 556 2 1.5423E-01 2128 4 1.4946E-01 1842 3 1.0451E-01 8323 13 9.5117E-02 7708

References

[1]Kellogg RB.Singularities in interface problems .In Symposium on Numerical Solutions of Partial Differential Equations,vol.II. Acadamic Press: New York,1971;351-400 [2]Grisvard P.Singularities in Boundary Value Problems. Masson: Paris, Springer: Berlin,

1992.

[3]S.C.Brenner ,Multigrid methods for the computation of singular solutions and stress intensity factors I: Corner singularities,preprint ,1996.

[4]M.Dauge,Elliptic Boundary Value Problems on Corner Domains,Lecture Notes in Math-ematics 1341,Springer-Verlag,Berlin-Heidelberg,1988.

[5]P.Grisvard,Elliptic Problrms in Non Smooth Domains,Pitman,Boston,1985.

[6]V.Kondratiev,Boundary value problems for elliptic equations in domains with conical or angular points, Tran.Moscow Math.Soc.,16(1967),pp.277-313.

[7]S.A.Nazarov and B.A.Plamenevsky,Elliptic Problems in Domains with Piecewise Smooth Boundaries,de Gruyter,Expositions in Mathematics ,13,Berlin,New York,1994.

[8]Satya N.Atluri, Computational Methods for Plane Problems of Fracture,Georgia Institue of Technology,U.S.A .

[9]S.C.Brenner and L.-Y.Sung ,Multigrid methods for the computation of singular solutions and stress intensity factors II: Crack singularities,1997.

[10]S.C.Brenner,Overcoming corner singularuties by multigrid methods ,preprint,1996. [11]I.Babuska and A. Miller,The post-processing approach in the nite element method

part 2 : The calculation of stress intensity factors,Int.J.Numer.Methods Engrg.,20(1984),pp.1111-1129.

[12]M.S.Birman and G.E.Skvorcov,On the square integrability of the highest derivatives of the solution of the Dirichlet problem on domain with piecewise smooth boundaries(in Russian),Izv.Vyssh.Uchebn.Zaved.Mat., Vol (1962),pp.136-155.

[13]H.Blum and M.Dobrowolski, On nite element methods for elliptic equations on do-mains with corners,Computing,28 (1982),pp.53-63.

References 56

[14]M.Dauge, M.-S. Lubuma, and S. Nicaise, Coef cient des singularities pour le probleme de Dirichlet sur un polygone,C.R.Acad. Sci.Paris Ser. I Math., 304 (1987),pp.483-486. [15]M.Dauge, S. Nicaise, M.Bourlard, and M.-S. Lubuma ,Coef cients des singularities pour des problemes aux limites elliptiques sur un domaine a points coniques I: resultats generaux pour le probleme de Dirichlet, M2_{AN,24 (1990),pp.27-52.}

[16]M.Dauge, S. Nicaise, M.Bourlard, and M.-S. Lubuma ,Coef cients des singularities pour des problemes aux limites elliptiques sur un domaine a points coniques II: quelques operateurs particuliers, M2_{AN,24 (1990),pp.343-367.}

[17]R.Verfurth,A posteriori error estimation and adaptive mesh-re nement techniques, Jour-nal of ComputatioJour-nal and Applied Mathematics 50 (1994) 67-83.

[18]Claes Johnson, ”Numerical solution of partial differential equations by the nite ele-ment method”,Cambridge University Press, 1988.

[19]William L.Briggs,"A Multigrid Tutorial",Department of Mathematics University of Colorado at Denver Denver, Colorado.

[20]Richard Barrett,Michael Berry,Tony F.Chan,James Demmel ,June M. Donato,Jack Don-garra,Victor Eijkhout,Roldan Pozo,Charles Romine,and Henk Van der Vorst"Templates for the Solution of Linear Systems :Building Blocks for Iterative Methods".

[20]I.Babuska, R.B. kellogg, J,Pitkaranta, Direct and Inverse Error Estinates for Finite El-ements with Mesh Re nEl-ements,1972.

[21]I.Babuska, Michael B. Rosenzweig, A nite element scheme for domains with cor-ners,1971.