- Email: [email protected]
Effective condition number for the finite element method using local mesh refinements
Applied Numerical Mathematics 59 (2009) 1779–1795
Contents lists available at ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
Effective condition number for the finite element method using local mesh refinements Zi-Cai Li a,b , Hung-Tsai Huang c,∗ a b c
Department of Applied Mathematics and Department of Computer Science and Engineering, National Sun Yat-sen University, Kaohsiung, Taiwan 80424 Department of Applied Mathematics, Chung-Hua University, HsinChu, Taiwan Department of Applied Mathematics, I-Shou University, Taiwan 840
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 12 May 2008 Received in revised form 30 December 2008 Accepted 23 January 2009 Available online 14 February 2009 MSC: 65N10 65N30 Keywords: Condition number Effective condition number Finite element method Local refinements Singularity problem
This is a continued study but at advanced levels of effective condition number in [Z.C. Li, C.S. Chien, H.T. Huang, Effective condition number for finite difference method, J. Comput. Appl. Math. 198 (2007) 208–235; Z.C. Li, H.T. Huang, Effective condition number for numerical partial differential equations, Numer. Linear Algebra Appl. 15 (2008) 575–594] for stability analysis. To approximate Poisson’s equation with singularity by the finite element method (FEM), the adaptive mesh refinements are an important and popular technique, by which, the FEM solutions with optimal convergence rates can be obtained. The local mesh refinements are essential to FEM for solving complicated problems with singularities, and they have been used for three decades. However, the traditional condition number is given 2 ) in Strang and Fix [G. Strang, G.J. Fix, An Analysis of the Finite Element by Cond = O(h− min Method, Prentice-Hall, Englewood Cliffs, NJ, 1973], where hmin is the minimal length of elements. Since hmin is infinitesimal near the singular points, Cond is huge. Such a dilemma can be bypassed by small effective condition number. In this paper, the bounds of the simplified effective condition number Cond_EE are derived as O(1), O(h−1.5 ) or O(h−0.5 ), where h(< 1) is the maximal length of elements. Evidently, Cond_EE is much smaller than Cond. The numerical experiments are carried out, to verify the stability analysis. Small effective condition numbers explain well the satisfactory FEM solutions obtained. This paper provides a stability justification for the adaptive mesh refinements used in FEM. Compared with [Z.C. Li, C.S. Chien, H.T. Huang, Effective condition number for finite difference method, J. Comput. Appl. Math. 198 (2007) 208–235; Z.C. Li, H.T. Huang, Effective condition number for numerical partial differential equations, Numer. Linear Algebra Appl. 15 (2008) 575–594], the analysis in this paper is more difficult and challenging, its proof techniques are new and intriguing, and the results are more important and useful. © 2009 IMACS. Published by Elsevier B.V. All rights reserved.
1. Introduction Consider the linear algebraic equations Ax = b,
(1.1) n×n
where x ∈ R and b ∈ R are the unknown and known vectors, respectively, and the matrix A ∈ R positive definite. The traditional condition number in the 2-norm is defined in [3,16] by n
*
n
Corresponding author. E-mail address: [email protected] (H.-T. Huang).
0168-9274/$30.00 © 2009 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2009.01.006
is symmetric and
1780
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
Cond =
λmax (A) , λmin (A)
(1.2)
where λmax and λmin are the maximal and the minimal eigenvalues of matrix A, respectively. For the perturbed system A(x + x) = b + b, there exists the error bound in [3,16],
b x Cond × , x b
(1.3)
where x is the 2-norm. The condition number is used to measure the solution errors resulting from data disturbance, in particular from rounding errors. Hence, the Cond plays a role in measuring numerical stability. Recently in Li et al. [8], we propose a better criterion to measure numerical stability. Let the eigenvalues of A be arranged in a descending order: λ1 λ2 · · · λn > 0. The eigenvectors ui satisfy Aui = λi ui , where {ui } are orthogonal, with uiT u j = δi j , where δi j = 1 if i = j and δi j = 0 if i = j. The new effective condition number is defined in [8] by Cond_eff =
λn
b
n 2 2 i =1 (βi /λi )
=
b , λn x
(1.4)
where the projection coefficients βi = uiT b. In [8], for A(x + x) = b + b,1 we have derived the following bound,
x b Cond_eff × . x b
(1.5)
Since there always exists the inequality: Cond_eff Cond, the effective condition number may provide a lower bound of the solution errors in data disturbance. In particular, when βn = 0, the simplest effective condition number is defined in [8] by Cond_EE =
b . |βn |
(1.6)
There also exists the inequality [8]: Cond_eff Cond_EE. Hence we may derive the bounds of Cond_EE, instead of those of Cond_eff for some problems. In this paper, we apply the effective condition number for the finite element method for Poisson’s equation. For the singular solutions of PDEs, when the local refinements of grids are used, the optimal convergence rates can be obtained. The local mesh refinements are essential to FEM for solving complicated problems with singularities, 2 and they have been used for three decades. However, the traditional condition number is given by Cond = O(h − ) in min Strang and Fix [13], where hmin is the minimal boundary length of elements. Since hmin is infinitesimal near the singular points, Cond is huge. Although a renovated bound of Cond is given in (6.17), Cond is still large. Such a dilemma can be bypassed by effective condition numbers. In this paper, the bounds of the simplified effective condition number Cond_EE are derived as O(1), O(h−1.5 ) or O(h−0.5 ), where h is the maximal boundary length of elements. Evidently, Cond_EE are much smaller than Cond, and so are Cond_eff. The numerical experiments are carried out, to verify the stability analysis. Small effective condition numbers explain well the satisfactory solutions obtained. This paper provides a stability justification for the important and popular adaptive mesh refinements used in FEM. Here we briefly review the condition number and the effective condition number. A rather detailed review is provided in [8]. The definition of the traditional condition number was given in Wilkinson [16], and then used in many books and papers. The condition number is used to provide the bounds of relative errors from the perturbation of both A and all b. However, in practical applications, we only deal with a certain vectors b, and the true relative errors may be smaller, or even much smaller than the worst Cond. Such a case was studied in Chan and Faulser [2], and called the effective condition number. However, this case was first proposed in Rice [12] in 1981, but called the natural condition number. We choose the terminology, the effective condition number used in [2]. This paper is organized as follows. In Section 2, the optimal convergence rates are briefly introduced for FEM using local mesh refinements. In Sections 3 and 4, the bounds of the simplest effective condition number are derived for the homogeneous and the non-homogeneous boundary conditions, respectively. In Section 5, the relations of assumptions are discussed, and the stability analysis is improved. In Section 6, numerical experiments are carried out for Motz’s problem. Below let us describe the key results of the following sections, and their linkage. We also compare the results in this paper with those in our previous study [8,9]. To explore the local mesh refinements used in FEM, the next section is for error analysis, and Sections 3–5 are for stability analysis. Based on (1.3) and (1.5), both Cond and Cond_eff (or Cond_EE) provide the bounds of amplified factors of the solution errors, with respect to the errors in disturbance data. For the same matrix A in (1.1), the bounds of effective condition numbers may be different for different b by (1.4) and (1.6). This is distinct from Cond in (1.3), which reaches the upper bounds of amplified factors for all b. For Poisson’s equation by the finite difference method (FDM) in [8,9] and the finite elements method in this paper, the bounds of Cond_eff and Cond_EE are different for different boundary conditions. Interestingly, for the homogeneous boundary conditions, the constant bounds, Cond_eff = O(1) and Cond_EE = O(1), 1
Errors are discussed in [9] for the general perturbed system: (A + A)(x + x) = b + b.
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
1781
remain even when the minimal boundary length hmin → 0. Such an extra-ordinary feature of effective condition number is remarkably beneficial to local mesh refinements in FEM, since hmin is infinitesimal in application. Also since the case of homogeneous boundary conditions is very interesting and important, they are first studied in Section 3. In general, the homogeneous conditions cannot be satisfied for Poisson’s equation. By a transformation, however, partial homogeneous boundary conditions may be satisfied. In Section 4, we assume that the non-homogeneous boundary conditions are subjected to the domain boundary far from the corner singular point. Under some simple partition, the bounds Cond_EE = O(h−3/2 ) and Cond_EE = O(h−1/2 ) are derived in (4.19) later for the non-homogeneous Dirichlet and the nonhomogeneous Neumann boundary conditions, respectively, where h is the maximal boundary length of elements. From (1.5) and Cond_eff Cond_EE, we have
x b Cond_EE × , x b
(1.7)
to give
x b Ch−q × , x b
(1.8)
where q = 0, 32 , 12 are for the homogeneous, the non-homogeneous Dirichlet, and the non-homogeneous Neumann boundary conditions, respectively. The amplified factors O(h−q ) in (1.8) are not large, since h is not small in application. Evidently, the stability of local mesh refinements in FEM is not as severe as we thought before. Contrarily, its stability is quite well. This paper enhances the local mesh refinements which have been used in FEM for three decades. In Section 5, the analysis for assumption A2 in (4.17) is twofold. First, we may answer the question: why does the good stability happen for local mesh refinements? The reason comes from the equivalent assumption A2* in (5.6) (i.e. (5.8)) that for the PDE solutions, their portion of low frequency is not insignificant. The real amplified factors in (1.8) are much smaller 2 ), where only the high frequent portion is dominant in the PDE solutions. This fact leads than the traditional Cond = O(h− min
to the miracle of good stability analysis, if using effective condition number. Second, by assumption A2* , when the simple partition is replaced by the common quasiuniform partition in S 1 , the bounds in (4.19) are also valid. For the key results in (4.19), the analysis in Section 4 is simple and straightforward, and the analysis via Section 5 is complete and sophisticated. Sections 3–5 provide a comprehensive study for the new stability analysis of local mesh refinements used in FEM. Moreover, let us also compare this paper with our previous papers [8,9]. In papers [8,9] and this paper, the effective condition numbers are applied to the finite difference method (FDM) and the finite element method (FEM), respectively. In [8,9] the solutions of PDE are supposed to be highly smooth without singularity. In this paper, the solutions with corner singularity are discussed. Compared with [8,9], the analysis in this paper is more difficult and challenging, and its proof techniques are new and intriguing. Moreover, since the local mesh refinements are essential to FEM for singularity solutions and complicated problems, the new stability analysis in this paper is more significant and helpful to modern engineering computation. Hence, the stability analysis in this paper is the study of effective condition numbers at advanced levels. 2. Optimal convergence rates of FEM using local mesh refinements Local mesh refinements were discussed early in Gregory et al. [4] in 1978, and error analysis is provided in Wahlbin [14] in 1991. A review of local mesh refinements are provided in [14]. Hence, the local mesh refinements in finite elements have been used for three decades. This paper considers the benchmark Motz’s problem of singularity problems of Laplace’s equation, and focuses on stability analysis by the effective condition number. The stability analysis will be carried out for the local mesh refinements, to maintain the optimal convergence rates. Since the optimal convergence rates are closely related to our stability analysis, their results are introduced before an exploration of effective condition numbers is unfolded. In this section, we consider the finite element method (FEM) for solving the homogeneous boundary conditions of Poisson’s equation,
−u = f u=0
in S ,
on Γ D ,
(2.1)
∂u = 0 on Γ N , ∂n
(2.2)
where ∂∂ nu is the outward derivatives on ∂ S, and S is a polygon, having only a corner O with the interior angle Θ ∈ ( π2 , 2π ], see Fig. 1. Eqs. (2.2) are the mixed type of the Dirichlet and the Neumann boundary conditions. From Lehman [5] and Wigley [15] (also see [11], Chapter 11), the solution on the sector S 0 = {(r , θ) | 0 r r0 , 0 θ Θ} has the expansion u=
∞
di r μi cos(μi θ) (r , θ) ∈ S 0 ,
(2.3)
i =0
where (r , θ) are the polar coordinates with the origin O , di are the true coefficients, and
μi =
(i + 12 )π , Θ
i = 0, 1, . . . .
(2.4)
1782
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
(a)
(b)
Fig. 1. (a) A concave polygon, (b) the local refinements near O .
Hence when r → 0 we have
u (r , θ) = O r σ ,
(2.5)
where
σ = μ0 =
π 2Θ
1 4
(2.6)
.
Then in this paper, we consider only Θ ∈ ( π2 , 2π ] and
σ ∈ [ 14 , 1) for singularity problems.
Let S be divided into small triangles i j such that S = ments of triangulation near O are chosen as follows: r i = (ih) p 1 , where p 1 1, h =
i, j
i j . In S 0 the nodes (r i , θ j ) of triangles i j in local refine-
θ j = j θ, r0 Na
and θ =
(2.7) Θ N
. The parameter p 1 ( 1) should be chosen to reach the optimal convergence rates as
p 1 > σ1 from Theorem 2.1 given below. In the rest of the region S 1 = S \ S 0 , since the solution is smooth enough, we may still use the quasiuniform triangles i j with max(i , j )∈ S 1 h i j min(i , j )∈ S 1 h i j
C,
(2.8)
where h i j are the maximal boundary length of i j , and C is a constant independent of h. Also denote h = max{h, θ, max(i , j )∈ S 1 h i j } and assume h < 1 in application. Denote by V h0 ⊂ H 1 ( S ) the finite dimensional collection of all continuous piecewise linear interpolation functions on S h satisfying v |Γ D = 0, where H 1 ( S ) is the Sobolev space. Then the finite element method for (2.1) and (2.2) can be expressed as: To seek u h ∈ V h0 such that
∇ uh ∇ v ds = S
f v,
∀ v ∈ V h0 .
(2.9)
S
near the origin O as shown in Fig. 2, which are formulated from the four nodes A (r i , θ j ), Denote the pair of ± ij
B (r i +1 , θ j ), C (r i , θ j +1 ) and D (r i +1 , θ j +1 ). Denote u I the piecewise interpolant function on ± based on the true values ij u A , u B , u C and u D . Below, we give a lemma and an important theorem, whose simple proof is given in [10], and is omitted, since these results are well known in Wahlbin [14]. be given with (2.7) in Fig. 2. Then there exists the bound, Lemma 2.1. Let σ ∈ [ 14 , 1) and the triangles ± ij
|u − u I |21,+ ∪− C θ (r i )3 r i2σ −3 + (r i )(θ)2 r i2σ −1 , ij
(2.10)
ij
where u I is the piecewise linear interpolant of u on ± . ij Theorem 2.1. Let σ ∈ [ 14 , 1), and the local mesh refinements as (2.7) be given in S 0 and the quasiuniform triangulation in S 1 . Then the solution u h from (2.9) has the following convergence rates,
u − u h 1, S = O h τ ,
1 p1 <
1
σ
,
(2.11)
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
1783
Fig. 2. The triangles near O .
where τ = p 1 σ , and
u − uh 1, S = O h | ln h| , u − uh 1, S = O(h),
p1 >
p1 = 1
σ
1
σ
(2.12)
,
(2.13)
,
where v 1, S is the Sobolev norm. From Theorem 2.1, when p 1 σ1 and p 1 > σ1 , the reduced and the optimal convergence rates are obtained, respectively. Hence in the next sections for stability, we will study the local mesh refinements with p 1 > σ1 only. 3. Homogeneous boundary conditions Consider the corresponding eigenvalue problem with the homogeneous boundary conditions,
− w = λ w in S , w =0
∂w = 0 on Γ N . ∂n
on Γ D ,
(3.1)
From [7,11], the lowest frequent wave as the eigenfunction near O is given by
w min = w min (r , θ) = O J σ (r ) = O r σ
as r → 0,
(3.2)
where J σ (r ) is the Bessel function defined by J σ (r ) =
∞ i =0
2i +σ (−1)i r .
(i + 1) (i + σ + 1) 2
(3.3)
Hence we may have the non-homogeneous function in (2.1)
f = f (r , θ) = O r σ −2 ,
σ
1 4
as r → 0.
(3.4)
We have the following lemma. Lemma 3.1. Let (3.4) be given. Then
β¯n =
f w min ds C , S
where C is a bounded constant independent of h.
(3.5)
1784
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
Proof. Since S = S 0 ∪ S 1 , we have
β¯n =
f w min ds =
f w min ds +
S
S0
f w min ds.
(3.6)
S1
For the first term in the above equation,
Θr0 f w min ds C
S0
r0
r σ −2 r σ r dr dθ = C Θ r 2σ 0 = O(1).
(3.7)
0 0
It is easily to see
f w min ds C .
(3.8)
S1
The desired result (3.4) follows. This completes the proof of Lemma 3.1.
2
Denote the vector
un = yhmin = c 0 h . . . , w min (xi , y j ), . . .
T
(3.9)
,
where the coefficient c 0 h results from the orthogonal vector un = 1, and c 0 is a constant independent of h. From (2.9), the linear algebraic equations are obtained, Ax = b,
(3.10)
where A ∈ R is the symmetric and positive definite matrix, b results from S ˆf ds, and ˆf is the piecewise linear interT is an approximation of the eigenfunction w min corresponding to λmin of (3.1). Such an polant of f . Suppose that un = ymin n×n
assumption is supported by our numerical experiments, and more exploration is given in Remark 3.1. Hence we have
T βn = yhmin b = c 0 h
ˆf w h
min ds ≈ c 0 h
S
f w min ds.
(3.11)
S
From Lemma 3.1, we may give the following assumption A1,
β¯ =
f w min ds O(1).
(3.12)
S
The notation a b or a O(b) implies that there exists two positive constants C 1 and C 2 such that 0 < C 1 b |a| C 1 b, where b > 0. Under A1 we have
βn ≈ c 0 β¯ h O(h).
(3.13)
Besides, we have the following lemma. Lemma 3.2. Let FEM,
σ ∈ [ 14 , 1), p 1 > σ1 , and the local refinements with (2.7) and the quasiuniform partition in S 1 . Then for the linear
b = O(h).
(3.14)
Proof. Denote the linear interpolation function on
ˆf h = f 1 φ1 + f 2 φ2 + f 3 φ3 on ,
(3.15)
where 1, 2 and 3 are three vertices of , and the hat-like basic functions φi satisfy |φi | 1. Then we have
3 3 3 ˆ ( f h ) ds = f i φi ds | f i | ds = | f i | × ||.
i =1
i =1
(3.16)
i =1
Let (i , j ) = (r i , θ j ) be an interior node with a few neighboring i , j , see Fig. 3, and their control regions S i j are formulated by connecting the midpoints of edges of i , j . Hence, we obtain the discrete vector b = {b i , j } with the components,
|bi , j | c i j S i , j | f i , j |,
(3.17)
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
1785
Fig. 3. The control region S i j .
where c i j are the constants independent of h. Then we have
b 2 C
f i2j S 2i j C
r i r0
ij
f i2, j S 2i , j +
r i >r0
j
f i2, j S 2i , j .
(3.18)
j
For r ∈ (0, r0 ), there exist the bounds,
| f i , j | Cr iσ −2 , S i , j Cr i r i θ,
p r i = r i +1 − r i = (i + 1)h 1 − (ih) p 1 Ch(ih) p 1 −1 .
(3.19)
Then for the first term in (3.18) we have
r i r0
f i2, j S 2i , j C
r i2σ −4 r i2 h2 (ih)2p 1 −2
r i r0
j
= C θ
(θ)2
j
(ih)(2σ −2) p 1 h2 (ih)2p 1 −2
r i r0
= C θ
(ih)2σ p 1 r i r0
i2
(3.20)
.
Since p 1 > σ1 we obtain
(ih)2σ p 1 i2
r i r0
Ch2σ p 1
(i )2σ p 1 r i r0
i2
Ch2σ p 1 ( N a )2σ p 1 −1 = C ( N a h)2σ p 1 C where N a =
r0 h
r i r0
1 Na
1 Na
Ch,
(3.21)
= O( h1 ). Combining (3.20) and (3.21) gives f i2, j S 2i , j Ch2 .
(3.22)
j
Since the solution is smooth for r > r0 , we can show
r i >r0
f i2, j S 2i , j Ch2 .
(3.23)
j
We obtain the desired result (3.14) from (3.18), (3.22) and (3.23). This completes the proof of Lemma 3.2. From (1.6), (3.13), and Lemma 3.2, we obtain the following theorem, immediately.
2
1786
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
Theorem 3.1. Let σ ∈ [ 14 , 1), p 1 > σ1 , and the local refinements with (2.7) and the quasiuniform partition in S 1 . Under assumption A1, for the homogeneous boundary conditions in (3.1), there exists the bound, Cond_EE = O(1).
(3.24)
First, the constant in (3.24) is extra-ordinary, against our past intuition for the case of local mesh refinements before. Second, the condition p 1 > 1 in Theorem 3.1 grants the optimal convergence rate given in Theorem 2.1. From our past knowledge, the p 1 should be chosen as small as possible, but to still retain the optimal convergence rates. A natural choice is given by p 1 = σ1 + δ , where 0 < δ 1. Since the worry about instability is now removed, in computation, we may choose rather lager p 1 ,2 to provide more a flexibility in application. In theory, we may even choose p 1 → ∞, and obtain (2.13) and (3.24) simultaneously. Such a result is very interesting and important, if noting the traditional condition number
2 = O h−2p 1 → ∞ Cond = O h− min
when p 1 → ∞.
(3.25)
The renovated bounds of Cond are also provided in Remark 6.1 later. Remark 3.1. The corresponding matrix eigenvalues of (3.1) may be expressed by3 Ay = μy,
(3.26)
or Az = λBz,
(3.27) n×n
where the matrix B ∈ R is symmetric and positive. The approximation of w min in (3.1) is zmin in (3.27), but not ymin in (3.26). Strictly speaking, the equality ymin = zmin holds if and if only two matrices A and B are exchangeable. In fact, we may consider the perturbation of x and f in the following linear algebraic equations, instead of (3.10) Ax = Bf,
(3.28)
where f = {. . . , f i j , . . .} and Bf = b ∈ R . Denote the weight norm T
x B =
n
(Bx, x) = B1/2 x,
(3.29)
where matrix B is symmetric and positive definite. The eigenvectors zi of (3.27) are orthogonal in the norm · B ,
(Bzi , z j ) = δi j ,
(3.30)
where δi j = 0 if i = j and δi j = 1 if i = j. By following [8,9], we have
f B x B CondB_eff × , x B f B
(3.31)
where the effective condition number is defined by CondB_eff =
f B λhmin x B
(3.32)
,
and λhmin is the minimal eigenvalue of (3.27), to satisfy
λhmin ≈ λmin O(1),
(3.33)
and λmin is the minimal eigenvalue of (3.1). Moreover, the simplest effective condition number can also be defined by CondB_EE =
f B , βnB
where
βnB = Bf, zhmin ≈
(3.34)
f zhmin ds ≈ S
f w min ds.
(3.35)
S
By following the arguments in Lemma 3.1, we may derive In Section 6, we choose p 1 = 4 > σ1 , where σ = 12 , and good numerical results are obtained. For simplicity, we only consider the homogeneous boundary conditions in this section; the conclusions for other problems discussed in the next section are similar. 2 3
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
f B = O(1) when p 1 > Also we have βnB ≈ we obtain
1
σ
(3.36)
.
f w min ds O(1), based on A2* :
S
CondB_EE = O(1)
1787
when p 1 >
1
σ
S
u w min ds O(1), shown in (5.6) later. Hence from (3.34)–(3.36)
(3.37)
.
The same constant effective condition number as in (3.24) is achieved. Although the analysis for CondB_EE is more strict for the equality: ymin = zmin , we prefer Cond_EE to CondB_EE, because Eq. (3.10) is easier and simpler in application than (3.28). On the other hand, for the goal β¯ O(1) in assumption A1, the approximation between ymin and zmin may be greatly relaxed. Denote the real approximate function w ∗min (= w min ) for ymin , and also let
β¯ =
f w ∗min ds.
(3.38)
S
As long as S f w ∗min ds O(1), we have β¯ O(1) and βn O(h), and then Theorem 3.1 holds. Hence, the approximation: ymin ≈ zmin is not a necessary condition for β¯ O(1). In fact, we have carried some computational comparisons, to show that there does exist an approximation: ymin ≈ zmin , where their maximal relative errors are less than 4%. 4. Non-homogeneous boundary conditions Below, let us consider the non-homogeneous boundary condition far from the corner O ,
−u = f , in S , ∂ u = 0, u = u |θ=Θ = 0, ∂ n θ=0 ∂ u ∗ u = u |Γ = g , = g∗, D ∂ n Γ ∗
(4.1)
N
where Γ D∗ = Γ D ∩ ∂ S¯ 1 , Γ N∗ = Γ N ∩ ∂ S¯ 1 and ∂ S = Γ D ∪ Γ N . When ∂∂ nu |θ=Θ = g ∗ and u = u |Γ ∗ = g, we can find the harmonic D functions u¯ from [11] to satisfy
−u¯ = 0, in S 0 , ∂ u¯ = g∗, u = u¯ |θ=Θ = g . ∂ n θ=0
(4.2)
¯ we obtain the problem (4.1). Hence by a transformation u − u, The FEM solution is expressed as: To seek u h ∈ V h such that
∇ uh ∇ v ds = S
f v ds +
g ∗ v d,
∀ v ∈ V h0 ,
(4.3)
Γ∗
S
N
where V h and V h0 denote the finite collections of all continuous piecewise linear functions on S h = i j i j satisfying the non-homogeneous and the homogeneous Dirichlet boundary conditions, respectively. In S 1 , the quasiuniform i j at Γ D∗ are first chosen to be right triangles,4 as shown in Fig. 4(a). In Fig. 4(b), denote A ∈ Γ D∗ and B ∈ (Γ D∗ )h , where (Γ D∗ )h is the nearest element edge to (Γ D∗ ). For the piecewise linear function v, we have
|∇ v |2 dx dy = 1 ∪2
1 ∪2
∂v ∂n
2 +
∂v ∂s
2 dx dy ,
(4.4)
where
1 ∪2
∂v ∂n
2 dx dy =
vB − vA d
2 |1 ∪ 2 |.
In (4.5), d is the distance between Γ D∗ and (Γ D∗ )h in Fig. 4, and 4
Such a limitation is removed in Section 5.
(4.5)
1788
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
(a)
(b)
Fig. 4. Right triangles near the non-homogeneous boundary condition.
si −1 + si
|1 ∪ 2 | = d
(4.6)
,
2 with si = si +1 − si . Hence we have
1 ∪2
∂v ∂n
2
dx dy =
( v B − v A )2 d
si −1 + si
×
2
(4.7)
.
∂v 2 1 ∪2 ( ∂ s ) dx dy has no effects on b, the vector
Since the integral
b = b f + bg + bg∗ ,
(4.8)
where the bound of b f is given in (3.18), and
bg = . . . ,
si −1 + si
2d
gi, j , . . . ,
Also we have from Γ ∗ g ∗ v
bg∗ = . . . ,
si −1 + si 2
(i , j ) ∈ Γ D∗ .
g i∗, j , . . . ,
(4.9)
(i , j ) ∈ Γ N∗ .
(4.10)
Below we consider
yhmin
T
b g = c0h
gi, j
si −1 + si 2d
(i , j )∈Γ D∗
w hmin ¯i , ¯j ,
(4.11)
where (¯i , ¯j ) ∈ (Γ D∗ )h . Since w min (Γ D∗ ) = 0 we have w min ((Γ D∗ )h ) d
≈−
∂ w min , ∂n
as d → 0.
(4.12)
Then Eq. (4.11) leads to
yhmin
T
b g ≈ c0h
w min ((Γ D∗ )h )
g
d
Γ D∗
d ≈ −c 0 h Γ D∗
g
∂ w min d. ∂n
(4.13)
Also since (si −1 + si )/d2 O( h1 ), we have from (4.9)
2 si −1 + si b g = g i2, j 2d
(i , j )∈Γ D∗
−0.5 Ch
si −1 + si
(i , j )∈Γ ∗
2
g i2 Ch−0.5 g 0,Γ ∗ . D
(4.14)
D
Next, we have
yhmin
T
b g∗ = c0h
(i , j )∈Γ N∗
≈ c0h Γ N∗
si −1 + si 2
g ∗ w min d.
× g i∗, j w hmin i , j (4.15)
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
1789
Also from (4.10)
b g ∗ =
si −1 + si 2 2
(i , j )∈Γ ∗
( g ∗ )2i , j Ch0.5 g ∗ 0,ΓN∗ .
(4.16)
N
Assume the following assumption A2,
β¯ ∗ =
f w min ds −
g
w min
∂n
ΓD ∗
S
d +
g ∗ w min d O(1).
(4.17)
Γ∗ N
Then we have from (3.11), (4.13) and (4.15)
T βn = yhmin b c 0 hβ¯ ∗ O(h),
(4.18)
and from (4.8), (3.14), (4.14) and (4.16) Cond_EE =
b C h−1.5 g 0,Γ D∗ + h−0.5 g ∗ 0,ΓN∗ . |βn |
(4.19)
We have the following theorem. Theorem 4.1. Let σ ∈ [ 14 , 1), p 1 > σ1 , and the local refinements with (2.7) hold. Suppose that the triangles i j in S 1 are quasiuniform, but the i j near Γ D∗ are right triangles. Under assumption A2, for the non-homogeneous boundary conditions in (4.1), the simplest effective condition number has the bound (4.19). From Theorem 4.1 for the non-homogeneous Dirichlet boundary condition far from O , we have
Cond_EE = O h−1.5 ,
(4.20)
and for the non-homogeneous Neumann boundary condition (i.e. g = 0 and g ∗
= 0),
Cond_EE = O h−0.5 .
(4.21)
Note that the bounds of Cond_EE in (4.20) and (4.21) do not depend on p 1 (i.e. the minimal boundary length hmin ). In theory, even when p 1 → ∞, the values of Cond_EE remain to be moderate, for non-homogeneous boundary conditions. Hence, based on the simplest effective condition number Cond_EE, the instability is not severe for the FEM using the local mesh refinements near the singularity at O . 5. Intrinsic view of assumption A2 and improvements of Theorem 4.1 5.1. Intrinsic view of assumption A2 The eigenvalue problem for (4.1) is given by
− w = λ w in S , ∂ w w |Γ D = 0, ∂n
= 0.
(5.1)
ΓN
We have the following lemma. Lemma 5.1. Let u and w be the solutions of (4.1) and (5.1), respectively. Then there exists an equality
f w ds − S
g
∂w d + ∂n
ΓD
g ∗ w d = λ
u w ds.
(5.2)
S
ΓN
Proof. For the Green formulas we have
( w u − u w ) ds = S
w
(5.3)
∂S
There exist the equalities,
( w u − u w ) ds = − S
∂u ∂w −u d. ∂n ∂n w f ds + λ
S
u w ds, S
(5.4)
1790
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
and
w
∂u ∂w −u ∂n ∂n
g ∗ w d −
d =
∂S
ΓN
g
∂w d . ∂n
(5.5)
ΓD
Combining (5.3)–(5.5) gives the desired result (5.1). This completes the proof of Lemma 5.1.
2
Based on Lemma 5.1, assumption A2 in (4.17) is rewritten as A2* :
u w min O(1),
(5.6)
S
due to λmin O(1). Denote x is the solution vector,
x = . . . , u h (xi , y j ), . . .
T
(5.7)
,
assumption A2* leads to unT x O(1),
(5.8)
where (σn , un ) is the minimal eigenpair of the discrete matrix A, which is symmetric and positive definite. Note that Eq. (5.6) resembles (5.8). Eq. (5.6) implies that the solution u is rich in the eigenfunction w min , and Eq. (5.8) also implies that the solution x is rich in the eigenvector un . For the real application, the PDE solution u can only vary in a certain region. Therefore, the opposite case of (5.6) is very seldom, and the analysis of the effective condition number in this paper is useful. 5.2. Improvements of Theorem 4.1 In Section 4, the triangles i j near Γ D are to be right triangles. In the rest of this section, such a limitation can be removed. Denote the vector b g = {. . . , αi j g i j , . . .}, where
(i , j ) ∈ Γ D∗ ,
(5.9)
αi j are constants. We have the following lemma.
Lemma 5.2. For the quasiuniform i j near Γ D∗ , the constants αi j in (5.9) have the bound,
|α i j | C ,
(5.10)
where C is a constant independent of h. Proof. Define the energy I (v ) =
1
|∇ v |2 ds −
2 S
f v ds −
g ∗ v d.
Γ N∗
S
Also denote the finite dimensional collection V h∗ (⊂ H 1 ( S )) of all continuous piecewise linear functions on S h1 = satisfying v |Γ ∗ = g. Then the FEM can be rewritten as: To seek u h ∈ V h∗ such that
(5.11)
ij
i j
D
I (u h ) = min I ( v ), v ∈ V h∗
i.e.,
∂ I (v ) = 0. ∂ v i, j
(5.12)
From Zienkiewicz [17], the linear interpolant functions are given by v = ϕi v i + ϕ j v j + ϕm v m , where i, j, and m are three nodes of , and
ϕi =
ai + b i x + c i y 2| |
,
ai = x j ym − xm y j , b i = y j − ym , c i = xm − x j ,
1 1 | | = area( ) = xi 2 y
i
1 xj yj
. ym
1 xm
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
1791
Similarly, the coefficients a j , b j , c j , am , bm and cm are obtained using the counter-clock rotation, see Zienkiewicz [16]. The discrete energy of (5.11) is given by I (v ) =
1 2
∂v ∂x
2
+
∂v ∂y
2 dx dy −
f v ds −
g ∗ v d.
(5.13)
∩Γ N∗
In (5.13), only the first term on the right hand is related to the vector b g , denoted by I1(v ) =
=
1 2
1
8| |
∂v ∂x
2 +
∂v ∂y
2 dx dy
(bi v i + b j v j + bm v m )2 + (c i v i + c j v j + cm v m )2 ,
and from (5.12)
∂ I 1 ( v ) 1 2 = b i + c 2i v i + (b i b j + c i c j ) v j + (b i bm + c i cm ) v m . ∂ vi 4| |
(5.14)
Since are quasiuniform, there exist the bounds: |b i | = | y i − ym | Ch, |c i | = |xm − y i | Ch and O(h2 ). Hence the coefficients in (5.14) have the bounds
|bi b j + c i c j | C,
|bi bm + c i cm | C.
(5.15)
In (5.14), when i ∈ S 1 and j ∈ Γ D (or m ∈ Γ D ) then v i = g i (or v m = gm ) from the Dirichlet boundary conditions. Then the coefficients of g i obtained form (5.14) are O(1), based on (5.15). Since the boundary nodes (i , j ) ∈ Γ D are connected to finite quasiuniform . Hence the coefficients of g i j in (5.9) satisfy (5.10). This completes the proof of Lemma 5.2. 2 Below, let us renovate Theorem 4.1. First, assumption A2 is replaced by A2∗ in (5.6). From Lemma 5.2, we can easily show
b g Ch−0.5 g 0,Γ D∗ . To remove the limitation of right triangles i j near Γ D∗ by the common assumption: quasiuniform elements, we need only prove
yhmin
T
bg ≈ − Γ∗
∂ w min g d, ∂n
(5.16)
D
to replace (4.13), since other bounds (4.14)–(4.16) also hold. First, consider the Neumann boundary condition on Γ N∗ , which can be obtained approximately by (4.3). Denote the element node A (i , j ) ∈ Γ N∗ . The integration of elements in (4.3) related to node A are given by 3 i =1
∇ uh ∇ v ds =
3 i =1
i
g ∗ v d,
i ∩Γ N∗
which gives
Lu i j +
S i +1 + S i 2
g ∗A v A = 0,
where i +1 and i are not necessarily to be the right triangles as 1 and 2 as shown in Fig. 4(b). Since v A is arbitrary, we obtained the boundary equation u A ∈ Γ N∗ ,
S i +1 + S i ∗ S i +1 + S i ∂ v Lu i j = − gA = − . 2 2 ∂n A
Next, consider the node A ∈ Γ D∗ . We obtain from (5.17) 3 i =1
i
S i +1 + S i ∂ v ∇ uh ∇ v ds = ( Lu i j · v ) A = ( Lu i j ) · g A ≈ − · gA. 2 ∂n A
(5.17)
1792
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
Hence we obtain
yhmin
T
b g = −c 0 h
S i +1 + S i ∂ w min ∂ w min g g d. ≈ −c 0 h · 2 ∂ n ∗ ∂n (i j )∈Γ ∗
(i , j )∈Γ D
D
Γ D∗
This is the desired result (5.16). We summarize this new result as a theorem. Theorem 5.1. Let σ ∈ [ 14 , 1), p 1 > σ1 and the local mesh refinements with (2.7) hold. Suppose that the triangles i j in S 1 are quasiuniform. Under assumption A2∗ (i.e. A2), the simplified effective condition number has the bound (4.19). 6. Numerical experiments In this section, first we consider Motz’s problem (see Fig. 5)
u = 0 in S , u = 0 on O D ,
u = 500 on A B , ∂u = 0 on BC ∪ C D ∪ O A , ∂n
(6.1)
where S = {(x, y ) | −1 < x < 1, 0 < y < 1}. The true functions are found as in [6] u=
∞
1
di r i + 2 cos i +
i =0
1 2
θ
in S ,
(6.2)
with the computed expansion coefficients di . The origin O is a singular point due to the intersection between the Neumann and the Dirichlet boundary conditions. To deal with the singularity, we choose the local refinements oftriangulation near O as shown in Fig. 6. Let the solution domain S be divided into small triangles i j such that S = i , j i j . The nodes
(r i , θ j ) of the small triangles i j near O are chosen as (2.7). In the our computation, we use the larger p 1 = 4(> σ1 = 2). √ where N a and N are the division numbers along the radial direction and the semi-circle, Denote h = 2/ N a and θ = π N respectively. The numerical results are listed in Table 1. Since Cond_eff Cond, and since Eq. (1.4) is easier than (1.6) in computation, we provide only the numerical data of Cond_eff. In the table, λmax (A) and λmin (A) are the maximal and the minimal eigenvalues of A, respectively, and Ratio(Cond_eff) = Cond_eff|2N /Cond_eff|2N −1 . From Table 1 we can see
Fig. 5. Motz’s problem.
Table 1 Errors and condition numbers for the non-homogeneous Dirichlet boundary condition with p 1 = 4. N (= N a )
4
8
16
32
64
128
maxi , j |εi , j | u − uh 0, S u − uh 1, S λmax (A) λmin (A) b x Cond(A) Cond_eff Ratio(Cond_eff)
3.03(1) 2.12(1) 1.10(2) 3.97(1) 5.44(−1) 5.04(2) 4.38(2) 7.30(1) 2.12 –
8.47 6.77 5.95(1) 8.27(1) 1.50(−1) 6.88(2) 1.03(3) 5.51(2) 4.45 2.099
4.67 2.01 30.9(1) 1.68(2) 3.76(−2) 9.92(2) 2.16(3) 4.47(3) 1.22(1) 2.742
2.44 5.58(−1) 15.7(1) 3.36(2) 9.42(−3) 1.42(3) 4.35(3) 3.57(4) 3.47(1) 2.844
1.04 1.47(−1) 7.88 6.65(2) 2.39(−3) 2.03(3) 8.72(3) 2.78(5) 9.74(1) 2.807
3.87(−1) 3.72(−2) 3.95 1.32(3) 6.03(−4) 2.88(3) 1.74(4) 2.19(6) 2.74(2) 2.813
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
1793
Fig. 6. Partition for Motz’s problem with N = 32.
Table 2 Errors and condition numbers for the non-homogeneous Neumann boundary condition with p 1 = 4. N (= N a )
4
8
16
32
64
128
maxi , j |εi , j | u − uh 0, S u − uh 1, S λmax (A) λmin (A) b x Cond(A) Cond_eff Ratio(Cond_eff)
4.72(1) 4.34(1) 1.08(2) 3.97(1) 2.31(−1) 2.41(2) 7.62(2) 1.72(2) 1.37 –
1.54(1) 1.28(1) 6.00(1) 8.27(1) 7.70(−2) 2.09(2) 1.32(3) 1.07(3) 2.06 1.5036
4.52 3.62 3.12(1) 1.68(2) 1.91(−2) 1.60(2) 2.41(3) 8.79(3) 3.47 1.6845
2.26 9.66(−1) 1.58(1) 3.36(2) 5.04(−3) 1.17(2) 4.60(3) 6.67(4) 5.05 1.4553
9.93(−1) 2.49(−1) 7.90 6.65(2) 1.30(−3) 8.46(1) 8.96(3) 5.12(5) 7.26 1.4376
3.74(−1) 6.28(−2) 3.96 1.32(3) 3.30(−4) 6.05(1) 1.77(4) 4.00(6) 1.03(1) 1.4187
max |εi , j | = O h2 | ln h| , i, j
u − uh 0, S = O h2 , u − uh 1, S = O(h), −1.5 Cond_eff = O h , −1 λmax (A) = O h , λmin (A) = O h2 , Cond(A) = O h−3 ,
(6.3) (6.4) (6.5) (6.6) (6.7)
where εi , j = u (i , j ) − u h (i , j ) and u (i , j ) = u (r i , θ j ). The asymptotic rates of errors in (6.4) are the optimal, to coincide with Theorem 2.1. The effective condition number in (6.5) also verifies Theorem 5.1 for the non-homogeneous Dirichlet boundary condition. Next, we consider the non-homogeneous Neumann boundary condition (see Fig. 5),
u = 0 in S , u=0
on O D ,
∂u = g ∗ on A B , ∂n
∂u = 0 on BC ∪ C D ∪ O A , ∂n
(6.8)
where g ∗ (= 0) is the normal derivative of (6.2) on A B, to replace the non-homogeneous Dirichlet condition u = 500 on A B in (6.1). The numerical results are listed in Table 2. From Table 2 we can see clearly
max |εi , j | = O h2 | ln h| , i, j
(6.9)
1794
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
u − uh 0, S = O h2 , u − uh 1, S = O(h), −0.5 , Cond_eff = O h −1 λmax (A) = O h , λmin (A) = O h2 , Cond(A) = O h−3 .
(6.10) (6.11) (6.12) (6.13) O(h −0.5 ) again
The asymptotic rates of numerical errors are also optimal, and the effective condition number Cond_eff = verifies Theorem 5.1 for the non-homogeneous Neumann boundary condition. From the analysis and the computation in this paper, the instability is not severe for the FEM using local refinements of elements near the singularity, based on 2 ) in [13]. Moreover, we obtain the ratio from Table 2 for Cond_eff = O(1), O(h−1.5 ) or O(h−0.5 ), in contrast to Cond = O(h− min N = 128, Cond(A) Cond_eff
=
4.01(6) 1.04(1)
= 3.85(5).
(6.14)
The above facts also confirm that the effective condition number is a better criterion for numerical PDE than the traditional condition number. Remark 6.1. In Bank and Scott [1], the local mesh refinements are limited under certain restrictions, to give an improved bound of Cond,
Cond(A) Ch−2 1 + ln
hmax
hmin
,
(6.15)
where C is a constant independent of h. Evidently, Eq. (6.15) does not coincide with the numerical results in (6.7) and (6.13):
Cond(A) = O h−3 .
(6.16)
Since the mesh family based on (2.7) is degenerate, thus the assumption of non-degenerate triangulations in [1] is not complied with. Note that the local mesh refinements (2.7) with the simplest node location is much simpler than those in [1,14]. For (2.7), a further exploration displays the new bound (details appear elsewhere): Cond(A) Ch−3 | ln hmin |,
(6.17)
where C is a constant independent of h. Eq. (6.17) is consistent perfectly with (6.16). Cond in (6.17) is large, although it is acceptable in application. In summary, based on effective condition numbers, but not on the traditional condition number, the instability is not severe. This paper provides a stability justification, to support the important adaptive mesh refinements in FEM existing for three decades. Remark 6.2. Motz’s problem is a typical corner singularity with mixed boundary conditions. The other kind of corner singularity is for the Dirichlet condition, where S 0 a concave polygon with Θ ∈ (π , 2π ]. The harmonic functions in S 0 with u = 0 at θ = 0 and θ = Θ have the expansion u=
∞
∗
c i r μi sin(μ∗i θ),
(r , θ) ∈ S 0 ,
(6.18)
i =0
where c i are the true coefficients, and
μ∗i =
iπ
Θ
,
i = 1, 2, . . . .
(6.19)
Hence when r → 0 we have
∗ u (r , θ) = O r σ ,
(6.20)
where
σ ∗ = μ∗0 =
π Θ
1 2
,
(6.21)
to give σ ∗ ∈ [ 12 , 1) for π < Θ 2π of the concave polygon. Compared with (2.5) for Motz’s problem, the corner singularity of Dirichlet conditions is weaker. Moreover, the Motz’s singularity is, indeed, of the singularity of concave corners with Dirichlet conditions. Let us consider an interior crack at O D = (−1 < x < 0) ∩ ( y = 0) within S ∗ = {(x, y ) | −1 < x < 1, −1 < y < 1} with the Dirichlet condition on ∂ S ∗ where u = 0 at O D (see Fig. 7). Suppose the solution in S ∗ is symmetric with respect to the axis X , i.e., u (x, y ) =
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
1795
Fig. 7. The interior crack problem.
u (x, − y ). By using this symmetry, the interior crack problem may be simplified to the problem in its half domain S (see Fig. 5) with the following mixed type of Dirichlet and Neumann boundary conditions at y = 0
∂u = 0, on O A , ∂n u = 0, on O D .
(6.22) (6.23)
Eqs. (6.22) and (6.23) are exactly the same corner boundary conditions (6.1) for Motz’s problem. In this paper, we choose Motz’s problem as the corner singularity for easy comparisons with our previous study for singularity problems in [6,11]. Acknowledgements We are grateful to Professors Robert Beauwens, Peter Jimack and the referees for their valuable comments and suggestions. We are also indebted to Mr. D.T.C. Wu for the computation in Section 6. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
R.E. Bank, R. Scott, On the conditioning of finite element equations with highly refined meshes, SIAM J. Numer. Anal. 26 (1989) 1383–1394. F.C. Chan, D.E. Foulser, Effectively well-conditioned linear systems, SIAM J. Stat. Comput. 9 (1988) 963–969. G.H. Golub, C.F. van Loan, Matrix Computations, second ed., The Johns Hopkins, Baltimore and London, 1989. J.A. Gregory, D. Frishelov, B. Schiff, J.R. Whiteman, Local mesh refinement with finite elements for elliptic problems, J. Comput. Phys. 29 (1978) 133–140. R.S. Lehman, Developments near an analytic corner of solution of elliptic differential equation, J. Math. Meth. 8 (1959) 727–760. Z.C. Li, Combined Methods for Elliptic Equations with Singularities, Interfaces and Infinities, Kluwer Academic Publishers, Dordrecht, Boston, 1998. Z.C. Li, Error analysis for boundary approximation methods solving eigenvalue problems, J. Comput. Appl. Math. 200 (2007) 231–254. Z.C. Li, C.S. Chien, H.T. Huang, Effective condition number for finite difference method, J. Comput. Appl. Math. 198 (2007) 208–235. Z.C. Li, H.T. Huang, Effective condition number for numerical partial differential equations, Numer. Linear Algebra Appl. 15 (2008) 575–594. Z.C. Li, H.T. Huang, Effective condition number for the finite element method using local mesh refinements, Technical report, Department of Applied Mathematics, National Sun Yet-sen University, Kaohsiung, Taiwan, 2007. Z.C. Li, T.T. Lu, H.Y. Hu, A.H.D. Cheng, Trefftz and Collocation Methods, WIT Press, Southampton, Boston, 2008. J.R. Rice, Matrix Computations and Mathematical Software, McGraw-Hill Book Company, New York, 1981. G. Strang, G.J. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cliffs, NJ, 1973. L.B. Wahlbin, Local behavior in finite element method, in: P.G. Ciarlet, J.L. Lions (Eds.), in: Finite Element Methods, vol. I, North-Holland, 1991, pp. 352– 522. N.M. Wigley, Asymptotic expansions at a corner of solutions of mixed boundary value problems, J. Math. Meth. 4 (1964) 549–579. J.H. Wilkinson, The Algebraic Eigenvalue Problem, Clarendon Press, Oxford, 1965, 191 p. O.C. Zienkiewicz, The Finite Element Method, third ed., McGraw-Hill, London, 1977.
Contents lists available at ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
Effective condition number for the finite element method using local mesh refinements Zi-Cai Li a,b , Hung-Tsai Huang c,∗ a b c
Department of Applied Mathematics and Department of Computer Science and Engineering, National Sun Yat-sen University, Kaohsiung, Taiwan 80424 Department of Applied Mathematics, Chung-Hua University, HsinChu, Taiwan Department of Applied Mathematics, I-Shou University, Taiwan 840
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 12 May 2008 Received in revised form 30 December 2008 Accepted 23 January 2009 Available online 14 February 2009 MSC: 65N10 65N30 Keywords: Condition number Effective condition number Finite element method Local refinements Singularity problem
This is a continued study but at advanced levels of effective condition number in [Z.C. Li, C.S. Chien, H.T. Huang, Effective condition number for finite difference method, J. Comput. Appl. Math. 198 (2007) 208–235; Z.C. Li, H.T. Huang, Effective condition number for numerical partial differential equations, Numer. Linear Algebra Appl. 15 (2008) 575–594] for stability analysis. To approximate Poisson’s equation with singularity by the finite element method (FEM), the adaptive mesh refinements are an important and popular technique, by which, the FEM solutions with optimal convergence rates can be obtained. The local mesh refinements are essential to FEM for solving complicated problems with singularities, and they have been used for three decades. However, the traditional condition number is given 2 ) in Strang and Fix [G. Strang, G.J. Fix, An Analysis of the Finite Element by Cond = O(h− min Method, Prentice-Hall, Englewood Cliffs, NJ, 1973], where hmin is the minimal length of elements. Since hmin is infinitesimal near the singular points, Cond is huge. Such a dilemma can be bypassed by small effective condition number. In this paper, the bounds of the simplified effective condition number Cond_EE are derived as O(1), O(h−1.5 ) or O(h−0.5 ), where h(< 1) is the maximal length of elements. Evidently, Cond_EE is much smaller than Cond. The numerical experiments are carried out, to verify the stability analysis. Small effective condition numbers explain well the satisfactory FEM solutions obtained. This paper provides a stability justification for the adaptive mesh refinements used in FEM. Compared with [Z.C. Li, C.S. Chien, H.T. Huang, Effective condition number for finite difference method, J. Comput. Appl. Math. 198 (2007) 208–235; Z.C. Li, H.T. Huang, Effective condition number for numerical partial differential equations, Numer. Linear Algebra Appl. 15 (2008) 575–594], the analysis in this paper is more difficult and challenging, its proof techniques are new and intriguing, and the results are more important and useful. © 2009 IMACS. Published by Elsevier B.V. All rights reserved.
1. Introduction Consider the linear algebraic equations Ax = b,
(1.1) n×n
where x ∈ R and b ∈ R are the unknown and known vectors, respectively, and the matrix A ∈ R positive definite. The traditional condition number in the 2-norm is defined in [3,16] by n
*
n
Corresponding author. E-mail address: [email protected] (H.-T. Huang).
0168-9274/$30.00 © 2009 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2009.01.006
is symmetric and
1780
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
Cond =
λmax (A) , λmin (A)
(1.2)
where λmax and λmin are the maximal and the minimal eigenvalues of matrix A, respectively. For the perturbed system A(x + x) = b + b, there exists the error bound in [3,16],
b x Cond × , x b
(1.3)
where x is the 2-norm. The condition number is used to measure the solution errors resulting from data disturbance, in particular from rounding errors. Hence, the Cond plays a role in measuring numerical stability. Recently in Li et al. [8], we propose a better criterion to measure numerical stability. Let the eigenvalues of A be arranged in a descending order: λ1 λ2 · · · λn > 0. The eigenvectors ui satisfy Aui = λi ui , where {ui } are orthogonal, with uiT u j = δi j , where δi j = 1 if i = j and δi j = 0 if i = j. The new effective condition number is defined in [8] by Cond_eff =
λn
b
n 2 2 i =1 (βi /λi )
=
b , λn x
(1.4)
where the projection coefficients βi = uiT b. In [8], for A(x + x) = b + b,1 we have derived the following bound,
x b Cond_eff × . x b
(1.5)
Since there always exists the inequality: Cond_eff Cond, the effective condition number may provide a lower bound of the solution errors in data disturbance. In particular, when βn = 0, the simplest effective condition number is defined in [8] by Cond_EE =
b . |βn |
(1.6)
There also exists the inequality [8]: Cond_eff Cond_EE. Hence we may derive the bounds of Cond_EE, instead of those of Cond_eff for some problems. In this paper, we apply the effective condition number for the finite element method for Poisson’s equation. For the singular solutions of PDEs, when the local refinements of grids are used, the optimal convergence rates can be obtained. The local mesh refinements are essential to FEM for solving complicated problems with singularities, 2 and they have been used for three decades. However, the traditional condition number is given by Cond = O(h − ) in min Strang and Fix [13], where hmin is the minimal boundary length of elements. Since hmin is infinitesimal near the singular points, Cond is huge. Although a renovated bound of Cond is given in (6.17), Cond is still large. Such a dilemma can be bypassed by effective condition numbers. In this paper, the bounds of the simplified effective condition number Cond_EE are derived as O(1), O(h−1.5 ) or O(h−0.5 ), where h is the maximal boundary length of elements. Evidently, Cond_EE are much smaller than Cond, and so are Cond_eff. The numerical experiments are carried out, to verify the stability analysis. Small effective condition numbers explain well the satisfactory solutions obtained. This paper provides a stability justification for the important and popular adaptive mesh refinements used in FEM. Here we briefly review the condition number and the effective condition number. A rather detailed review is provided in [8]. The definition of the traditional condition number was given in Wilkinson [16], and then used in many books and papers. The condition number is used to provide the bounds of relative errors from the perturbation of both A and all b. However, in practical applications, we only deal with a certain vectors b, and the true relative errors may be smaller, or even much smaller than the worst Cond. Such a case was studied in Chan and Faulser [2], and called the effective condition number. However, this case was first proposed in Rice [12] in 1981, but called the natural condition number. We choose the terminology, the effective condition number used in [2]. This paper is organized as follows. In Section 2, the optimal convergence rates are briefly introduced for FEM using local mesh refinements. In Sections 3 and 4, the bounds of the simplest effective condition number are derived for the homogeneous and the non-homogeneous boundary conditions, respectively. In Section 5, the relations of assumptions are discussed, and the stability analysis is improved. In Section 6, numerical experiments are carried out for Motz’s problem. Below let us describe the key results of the following sections, and their linkage. We also compare the results in this paper with those in our previous study [8,9]. To explore the local mesh refinements used in FEM, the next section is for error analysis, and Sections 3–5 are for stability analysis. Based on (1.3) and (1.5), both Cond and Cond_eff (or Cond_EE) provide the bounds of amplified factors of the solution errors, with respect to the errors in disturbance data. For the same matrix A in (1.1), the bounds of effective condition numbers may be different for different b by (1.4) and (1.6). This is distinct from Cond in (1.3), which reaches the upper bounds of amplified factors for all b. For Poisson’s equation by the finite difference method (FDM) in [8,9] and the finite elements method in this paper, the bounds of Cond_eff and Cond_EE are different for different boundary conditions. Interestingly, for the homogeneous boundary conditions, the constant bounds, Cond_eff = O(1) and Cond_EE = O(1), 1
Errors are discussed in [9] for the general perturbed system: (A + A)(x + x) = b + b.
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
1781
remain even when the minimal boundary length hmin → 0. Such an extra-ordinary feature of effective condition number is remarkably beneficial to local mesh refinements in FEM, since hmin is infinitesimal in application. Also since the case of homogeneous boundary conditions is very interesting and important, they are first studied in Section 3. In general, the homogeneous conditions cannot be satisfied for Poisson’s equation. By a transformation, however, partial homogeneous boundary conditions may be satisfied. In Section 4, we assume that the non-homogeneous boundary conditions are subjected to the domain boundary far from the corner singular point. Under some simple partition, the bounds Cond_EE = O(h−3/2 ) and Cond_EE = O(h−1/2 ) are derived in (4.19) later for the non-homogeneous Dirichlet and the nonhomogeneous Neumann boundary conditions, respectively, where h is the maximal boundary length of elements. From (1.5) and Cond_eff Cond_EE, we have
x b Cond_EE × , x b
(1.7)
to give
x b Ch−q × , x b
(1.8)
where q = 0, 32 , 12 are for the homogeneous, the non-homogeneous Dirichlet, and the non-homogeneous Neumann boundary conditions, respectively. The amplified factors O(h−q ) in (1.8) are not large, since h is not small in application. Evidently, the stability of local mesh refinements in FEM is not as severe as we thought before. Contrarily, its stability is quite well. This paper enhances the local mesh refinements which have been used in FEM for three decades. In Section 5, the analysis for assumption A2 in (4.17) is twofold. First, we may answer the question: why does the good stability happen for local mesh refinements? The reason comes from the equivalent assumption A2* in (5.6) (i.e. (5.8)) that for the PDE solutions, their portion of low frequency is not insignificant. The real amplified factors in (1.8) are much smaller 2 ), where only the high frequent portion is dominant in the PDE solutions. This fact leads than the traditional Cond = O(h− min
to the miracle of good stability analysis, if using effective condition number. Second, by assumption A2* , when the simple partition is replaced by the common quasiuniform partition in S 1 , the bounds in (4.19) are also valid. For the key results in (4.19), the analysis in Section 4 is simple and straightforward, and the analysis via Section 5 is complete and sophisticated. Sections 3–5 provide a comprehensive study for the new stability analysis of local mesh refinements used in FEM. Moreover, let us also compare this paper with our previous papers [8,9]. In papers [8,9] and this paper, the effective condition numbers are applied to the finite difference method (FDM) and the finite element method (FEM), respectively. In [8,9] the solutions of PDE are supposed to be highly smooth without singularity. In this paper, the solutions with corner singularity are discussed. Compared with [8,9], the analysis in this paper is more difficult and challenging, and its proof techniques are new and intriguing. Moreover, since the local mesh refinements are essential to FEM for singularity solutions and complicated problems, the new stability analysis in this paper is more significant and helpful to modern engineering computation. Hence, the stability analysis in this paper is the study of effective condition numbers at advanced levels. 2. Optimal convergence rates of FEM using local mesh refinements Local mesh refinements were discussed early in Gregory et al. [4] in 1978, and error analysis is provided in Wahlbin [14] in 1991. A review of local mesh refinements are provided in [14]. Hence, the local mesh refinements in finite elements have been used for three decades. This paper considers the benchmark Motz’s problem of singularity problems of Laplace’s equation, and focuses on stability analysis by the effective condition number. The stability analysis will be carried out for the local mesh refinements, to maintain the optimal convergence rates. Since the optimal convergence rates are closely related to our stability analysis, their results are introduced before an exploration of effective condition numbers is unfolded. In this section, we consider the finite element method (FEM) for solving the homogeneous boundary conditions of Poisson’s equation,
−u = f u=0
in S ,
on Γ D ,
(2.1)
∂u = 0 on Γ N , ∂n
(2.2)
where ∂∂ nu is the outward derivatives on ∂ S, and S is a polygon, having only a corner O with the interior angle Θ ∈ ( π2 , 2π ], see Fig. 1. Eqs. (2.2) are the mixed type of the Dirichlet and the Neumann boundary conditions. From Lehman [5] and Wigley [15] (also see [11], Chapter 11), the solution on the sector S 0 = {(r , θ) | 0 r r0 , 0 θ Θ} has the expansion u=
∞
di r μi cos(μi θ) (r , θ) ∈ S 0 ,
(2.3)
i =0
where (r , θ) are the polar coordinates with the origin O , di are the true coefficients, and
μi =
(i + 12 )π , Θ
i = 0, 1, . . . .
(2.4)
1782
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
(a)
(b)
Fig. 1. (a) A concave polygon, (b) the local refinements near O .
Hence when r → 0 we have
u (r , θ) = O r σ ,
(2.5)
where
σ = μ0 =
π 2Θ
1 4
(2.6)
.
Then in this paper, we consider only Θ ∈ ( π2 , 2π ] and
σ ∈ [ 14 , 1) for singularity problems.
Let S be divided into small triangles i j such that S = ments of triangulation near O are chosen as follows: r i = (ih) p 1 , where p 1 1, h =
i, j
i j . In S 0 the nodes (r i , θ j ) of triangles i j in local refine-
θ j = j θ, r0 Na
and θ =
(2.7) Θ N
. The parameter p 1 ( 1) should be chosen to reach the optimal convergence rates as
p 1 > σ1 from Theorem 2.1 given below. In the rest of the region S 1 = S \ S 0 , since the solution is smooth enough, we may still use the quasiuniform triangles i j with max(i , j )∈ S 1 h i j min(i , j )∈ S 1 h i j
C,
(2.8)
where h i j are the maximal boundary length of i j , and C is a constant independent of h. Also denote h = max{h, θ, max(i , j )∈ S 1 h i j } and assume h < 1 in application. Denote by V h0 ⊂ H 1 ( S ) the finite dimensional collection of all continuous piecewise linear interpolation functions on S h satisfying v |Γ D = 0, where H 1 ( S ) is the Sobolev space. Then the finite element method for (2.1) and (2.2) can be expressed as: To seek u h ∈ V h0 such that
∇ uh ∇ v ds = S
f v,
∀ v ∈ V h0 .
(2.9)
S
near the origin O as shown in Fig. 2, which are formulated from the four nodes A (r i , θ j ), Denote the pair of ± ij
B (r i +1 , θ j ), C (r i , θ j +1 ) and D (r i +1 , θ j +1 ). Denote u I the piecewise interpolant function on ± based on the true values ij u A , u B , u C and u D . Below, we give a lemma and an important theorem, whose simple proof is given in [10], and is omitted, since these results are well known in Wahlbin [14]. be given with (2.7) in Fig. 2. Then there exists the bound, Lemma 2.1. Let σ ∈ [ 14 , 1) and the triangles ± ij
|u − u I |21,+ ∪− C θ (r i )3 r i2σ −3 + (r i )(θ)2 r i2σ −1 , ij
(2.10)
ij
where u I is the piecewise linear interpolant of u on ± . ij Theorem 2.1. Let σ ∈ [ 14 , 1), and the local mesh refinements as (2.7) be given in S 0 and the quasiuniform triangulation in S 1 . Then the solution u h from (2.9) has the following convergence rates,
u − u h 1, S = O h τ ,
1 p1 <
1
σ
,
(2.11)
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
1783
Fig. 2. The triangles near O .
where τ = p 1 σ , and
u − uh 1, S = O h | ln h| , u − uh 1, S = O(h),
p1 >
p1 = 1
σ
1
σ
(2.12)
,
(2.13)
,
where v 1, S is the Sobolev norm. From Theorem 2.1, when p 1 σ1 and p 1 > σ1 , the reduced and the optimal convergence rates are obtained, respectively. Hence in the next sections for stability, we will study the local mesh refinements with p 1 > σ1 only. 3. Homogeneous boundary conditions Consider the corresponding eigenvalue problem with the homogeneous boundary conditions,
− w = λ w in S , w =0
∂w = 0 on Γ N . ∂n
on Γ D ,
(3.1)
From [7,11], the lowest frequent wave as the eigenfunction near O is given by
w min = w min (r , θ) = O J σ (r ) = O r σ
as r → 0,
(3.2)
where J σ (r ) is the Bessel function defined by J σ (r ) =
∞ i =0
2i +σ (−1)i r .
(i + 1) (i + σ + 1) 2
(3.3)
Hence we may have the non-homogeneous function in (2.1)
f = f (r , θ) = O r σ −2 ,
σ
1 4
as r → 0.
(3.4)
We have the following lemma. Lemma 3.1. Let (3.4) be given. Then
β¯n =
f w min ds C , S
where C is a bounded constant independent of h.
(3.5)
1784
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
Proof. Since S = S 0 ∪ S 1 , we have
β¯n =
f w min ds =
f w min ds +
S
S0
f w min ds.
(3.6)
S1
For the first term in the above equation,
Θr0 f w min ds C
S0
r0
r σ −2 r σ r dr dθ = C Θ r 2σ 0 = O(1).
(3.7)
0 0
It is easily to see
f w min ds C .
(3.8)
S1
The desired result (3.4) follows. This completes the proof of Lemma 3.1.
2
Denote the vector
un = yhmin = c 0 h . . . , w min (xi , y j ), . . .
T
(3.9)
,
where the coefficient c 0 h results from the orthogonal vector un = 1, and c 0 is a constant independent of h. From (2.9), the linear algebraic equations are obtained, Ax = b,
(3.10)
where A ∈ R is the symmetric and positive definite matrix, b results from S ˆf ds, and ˆf is the piecewise linear interT is an approximation of the eigenfunction w min corresponding to λmin of (3.1). Such an polant of f . Suppose that un = ymin n×n
assumption is supported by our numerical experiments, and more exploration is given in Remark 3.1. Hence we have
T βn = yhmin b = c 0 h
ˆf w h
min ds ≈ c 0 h
S
f w min ds.
(3.11)
S
From Lemma 3.1, we may give the following assumption A1,
β¯ =
f w min ds O(1).
(3.12)
S
The notation a b or a O(b) implies that there exists two positive constants C 1 and C 2 such that 0 < C 1 b |a| C 1 b, where b > 0. Under A1 we have
βn ≈ c 0 β¯ h O(h).
(3.13)
Besides, we have the following lemma. Lemma 3.2. Let FEM,
σ ∈ [ 14 , 1), p 1 > σ1 , and the local refinements with (2.7) and the quasiuniform partition in S 1 . Then for the linear
b = O(h).
(3.14)
Proof. Denote the linear interpolation function on
ˆf h = f 1 φ1 + f 2 φ2 + f 3 φ3 on ,
(3.15)
where 1, 2 and 3 are three vertices of , and the hat-like basic functions φi satisfy |φi | 1. Then we have
3 3 3 ˆ ( f h ) ds = f i φi ds | f i | ds = | f i | × ||.
i =1
i =1
(3.16)
i =1
Let (i , j ) = (r i , θ j ) be an interior node with a few neighboring i , j , see Fig. 3, and their control regions S i j are formulated by connecting the midpoints of edges of i , j . Hence, we obtain the discrete vector b = {b i , j } with the components,
|bi , j | c i j S i , j | f i , j |,
(3.17)
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
1785
Fig. 3. The control region S i j .
where c i j are the constants independent of h. Then we have
b 2 C
f i2j S 2i j C
r i r0
ij
f i2, j S 2i , j +
r i >r0
j
f i2, j S 2i , j .
(3.18)
j
For r ∈ (0, r0 ), there exist the bounds,
| f i , j | Cr iσ −2 , S i , j Cr i r i θ,
p r i = r i +1 − r i = (i + 1)h 1 − (ih) p 1 Ch(ih) p 1 −1 .
(3.19)
Then for the first term in (3.18) we have
r i r0
f i2, j S 2i , j C
r i2σ −4 r i2 h2 (ih)2p 1 −2
r i r0
j
= C θ
(θ)2
j
(ih)(2σ −2) p 1 h2 (ih)2p 1 −2
r i r0
= C θ
(ih)2σ p 1 r i r0
i2
(3.20)
.
Since p 1 > σ1 we obtain
(ih)2σ p 1 i2
r i r0
Ch2σ p 1
(i )2σ p 1 r i r0
i2
Ch2σ p 1 ( N a )2σ p 1 −1 = C ( N a h)2σ p 1 C where N a =
r0 h
r i r0
1 Na
1 Na
Ch,
(3.21)
= O( h1 ). Combining (3.20) and (3.21) gives f i2, j S 2i , j Ch2 .
(3.22)
j
Since the solution is smooth for r > r0 , we can show
r i >r0
f i2, j S 2i , j Ch2 .
(3.23)
j
We obtain the desired result (3.14) from (3.18), (3.22) and (3.23). This completes the proof of Lemma 3.2. From (1.6), (3.13), and Lemma 3.2, we obtain the following theorem, immediately.
2
1786
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
Theorem 3.1. Let σ ∈ [ 14 , 1), p 1 > σ1 , and the local refinements with (2.7) and the quasiuniform partition in S 1 . Under assumption A1, for the homogeneous boundary conditions in (3.1), there exists the bound, Cond_EE = O(1).
(3.24)
First, the constant in (3.24) is extra-ordinary, against our past intuition for the case of local mesh refinements before. Second, the condition p 1 > 1 in Theorem 3.1 grants the optimal convergence rate given in Theorem 2.1. From our past knowledge, the p 1 should be chosen as small as possible, but to still retain the optimal convergence rates. A natural choice is given by p 1 = σ1 + δ , where 0 < δ 1. Since the worry about instability is now removed, in computation, we may choose rather lager p 1 ,2 to provide more a flexibility in application. In theory, we may even choose p 1 → ∞, and obtain (2.13) and (3.24) simultaneously. Such a result is very interesting and important, if noting the traditional condition number
2 = O h−2p 1 → ∞ Cond = O h− min
when p 1 → ∞.
(3.25)
The renovated bounds of Cond are also provided in Remark 6.1 later. Remark 3.1. The corresponding matrix eigenvalues of (3.1) may be expressed by3 Ay = μy,
(3.26)
or Az = λBz,
(3.27) n×n
where the matrix B ∈ R is symmetric and positive. The approximation of w min in (3.1) is zmin in (3.27), but not ymin in (3.26). Strictly speaking, the equality ymin = zmin holds if and if only two matrices A and B are exchangeable. In fact, we may consider the perturbation of x and f in the following linear algebraic equations, instead of (3.10) Ax = Bf,
(3.28)
where f = {. . . , f i j , . . .} and Bf = b ∈ R . Denote the weight norm T
x B =
n
(Bx, x) = B1/2 x,
(3.29)
where matrix B is symmetric and positive definite. The eigenvectors zi of (3.27) are orthogonal in the norm · B ,
(Bzi , z j ) = δi j ,
(3.30)
where δi j = 0 if i = j and δi j = 1 if i = j. By following [8,9], we have
f B x B CondB_eff × , x B f B
(3.31)
where the effective condition number is defined by CondB_eff =
f B λhmin x B
(3.32)
,
and λhmin is the minimal eigenvalue of (3.27), to satisfy
λhmin ≈ λmin O(1),
(3.33)
and λmin is the minimal eigenvalue of (3.1). Moreover, the simplest effective condition number can also be defined by CondB_EE =
f B , βnB
where
βnB = Bf, zhmin ≈
(3.34)
f zhmin ds ≈ S
f w min ds.
(3.35)
S
By following the arguments in Lemma 3.1, we may derive In Section 6, we choose p 1 = 4 > σ1 , where σ = 12 , and good numerical results are obtained. For simplicity, we only consider the homogeneous boundary conditions in this section; the conclusions for other problems discussed in the next section are similar. 2 3
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
f B = O(1) when p 1 > Also we have βnB ≈ we obtain
1
σ
(3.36)
.
f w min ds O(1), based on A2* :
S
CondB_EE = O(1)
1787
when p 1 >
1
σ
S
u w min ds O(1), shown in (5.6) later. Hence from (3.34)–(3.36)
(3.37)
.
The same constant effective condition number as in (3.24) is achieved. Although the analysis for CondB_EE is more strict for the equality: ymin = zmin , we prefer Cond_EE to CondB_EE, because Eq. (3.10) is easier and simpler in application than (3.28). On the other hand, for the goal β¯ O(1) in assumption A1, the approximation between ymin and zmin may be greatly relaxed. Denote the real approximate function w ∗min (= w min ) for ymin , and also let
β¯ =
f w ∗min ds.
(3.38)
S
As long as S f w ∗min ds O(1), we have β¯ O(1) and βn O(h), and then Theorem 3.1 holds. Hence, the approximation: ymin ≈ zmin is not a necessary condition for β¯ O(1). In fact, we have carried some computational comparisons, to show that there does exist an approximation: ymin ≈ zmin , where their maximal relative errors are less than 4%. 4. Non-homogeneous boundary conditions Below, let us consider the non-homogeneous boundary condition far from the corner O ,
−u = f , in S , ∂ u = 0, u = u |θ=Θ = 0, ∂ n θ=0 ∂ u ∗ u = u |Γ = g , = g∗, D ∂ n Γ ∗
(4.1)
N
where Γ D∗ = Γ D ∩ ∂ S¯ 1 , Γ N∗ = Γ N ∩ ∂ S¯ 1 and ∂ S = Γ D ∪ Γ N . When ∂∂ nu |θ=Θ = g ∗ and u = u |Γ ∗ = g, we can find the harmonic D functions u¯ from [11] to satisfy
−u¯ = 0, in S 0 , ∂ u¯ = g∗, u = u¯ |θ=Θ = g . ∂ n θ=0
(4.2)
¯ we obtain the problem (4.1). Hence by a transformation u − u, The FEM solution is expressed as: To seek u h ∈ V h such that
∇ uh ∇ v ds = S
f v ds +
g ∗ v d,
∀ v ∈ V h0 ,
(4.3)
Γ∗
S
N
where V h and V h0 denote the finite collections of all continuous piecewise linear functions on S h = i j i j satisfying the non-homogeneous and the homogeneous Dirichlet boundary conditions, respectively. In S 1 , the quasiuniform i j at Γ D∗ are first chosen to be right triangles,4 as shown in Fig. 4(a). In Fig. 4(b), denote A ∈ Γ D∗ and B ∈ (Γ D∗ )h , where (Γ D∗ )h is the nearest element edge to (Γ D∗ ). For the piecewise linear function v, we have
|∇ v |2 dx dy = 1 ∪2
1 ∪2
∂v ∂n
2 +
∂v ∂s
2 dx dy ,
(4.4)
where
1 ∪2
∂v ∂n
2 dx dy =
vB − vA d
2 |1 ∪ 2 |.
In (4.5), d is the distance between Γ D∗ and (Γ D∗ )h in Fig. 4, and 4
Such a limitation is removed in Section 5.
(4.5)
1788
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
(a)
(b)
Fig. 4. Right triangles near the non-homogeneous boundary condition.
si −1 + si
|1 ∪ 2 | = d
(4.6)
,
2 with si = si +1 − si . Hence we have
1 ∪2
∂v ∂n
2
dx dy =
( v B − v A )2 d
si −1 + si
×
2
(4.7)
.
∂v 2 1 ∪2 ( ∂ s ) dx dy has no effects on b, the vector
Since the integral
b = b f + bg + bg∗ ,
(4.8)
where the bound of b f is given in (3.18), and
bg = . . . ,
si −1 + si
2d
gi, j , . . . ,
Also we have from Γ ∗ g ∗ v
bg∗ = . . . ,
si −1 + si 2
(i , j ) ∈ Γ D∗ .
g i∗, j , . . . ,
(4.9)
(i , j ) ∈ Γ N∗ .
(4.10)
Below we consider
yhmin
T
b g = c0h
gi, j
si −1 + si 2d
(i , j )∈Γ D∗
w hmin ¯i , ¯j ,
(4.11)
where (¯i , ¯j ) ∈ (Γ D∗ )h . Since w min (Γ D∗ ) = 0 we have w min ((Γ D∗ )h ) d
≈−
∂ w min , ∂n
as d → 0.
(4.12)
Then Eq. (4.11) leads to
yhmin
T
b g ≈ c0h
w min ((Γ D∗ )h )
g
d
Γ D∗
d ≈ −c 0 h Γ D∗
g
∂ w min d. ∂n
(4.13)
Also since (si −1 + si )/d2 O( h1 ), we have from (4.9)
2 si −1 + si b g = g i2, j 2d
(i , j )∈Γ D∗
−0.5 Ch
si −1 + si
(i , j )∈Γ ∗
2
g i2 Ch−0.5 g 0,Γ ∗ . D
(4.14)
D
Next, we have
yhmin
T
b g∗ = c0h
(i , j )∈Γ N∗
≈ c0h Γ N∗
si −1 + si 2
g ∗ w min d.
× g i∗, j w hmin i , j (4.15)
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
1789
Also from (4.10)
b g ∗ =
si −1 + si 2 2
(i , j )∈Γ ∗
( g ∗ )2i , j Ch0.5 g ∗ 0,ΓN∗ .
(4.16)
N
Assume the following assumption A2,
β¯ ∗ =
f w min ds −
g
w min
∂n
ΓD ∗
S
d +
g ∗ w min d O(1).
(4.17)
Γ∗ N
Then we have from (3.11), (4.13) and (4.15)
T βn = yhmin b c 0 hβ¯ ∗ O(h),
(4.18)
and from (4.8), (3.14), (4.14) and (4.16) Cond_EE =
b C h−1.5 g 0,Γ D∗ + h−0.5 g ∗ 0,ΓN∗ . |βn |
(4.19)
We have the following theorem. Theorem 4.1. Let σ ∈ [ 14 , 1), p 1 > σ1 , and the local refinements with (2.7) hold. Suppose that the triangles i j in S 1 are quasiuniform, but the i j near Γ D∗ are right triangles. Under assumption A2, for the non-homogeneous boundary conditions in (4.1), the simplest effective condition number has the bound (4.19). From Theorem 4.1 for the non-homogeneous Dirichlet boundary condition far from O , we have
Cond_EE = O h−1.5 ,
(4.20)
and for the non-homogeneous Neumann boundary condition (i.e. g = 0 and g ∗
= 0),
Cond_EE = O h−0.5 .
(4.21)
Note that the bounds of Cond_EE in (4.20) and (4.21) do not depend on p 1 (i.e. the minimal boundary length hmin ). In theory, even when p 1 → ∞, the values of Cond_EE remain to be moderate, for non-homogeneous boundary conditions. Hence, based on the simplest effective condition number Cond_EE, the instability is not severe for the FEM using the local mesh refinements near the singularity at O . 5. Intrinsic view of assumption A2 and improvements of Theorem 4.1 5.1. Intrinsic view of assumption A2 The eigenvalue problem for (4.1) is given by
− w = λ w in S , ∂ w w |Γ D = 0, ∂n
= 0.
(5.1)
ΓN
We have the following lemma. Lemma 5.1. Let u and w be the solutions of (4.1) and (5.1), respectively. Then there exists an equality
f w ds − S
g
∂w d + ∂n
ΓD
g ∗ w d = λ
u w ds.
(5.2)
S
ΓN
Proof. For the Green formulas we have
( w u − u w ) ds = S
w
(5.3)
∂S
There exist the equalities,
( w u − u w ) ds = − S
∂u ∂w −u d. ∂n ∂n w f ds + λ
S
u w ds, S
(5.4)
1790
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
and
w
∂u ∂w −u ∂n ∂n
g ∗ w d −
d =
∂S
ΓN
g
∂w d . ∂n
(5.5)
ΓD
Combining (5.3)–(5.5) gives the desired result (5.1). This completes the proof of Lemma 5.1.
2
Based on Lemma 5.1, assumption A2 in (4.17) is rewritten as A2* :
u w min O(1),
(5.6)
S
due to λmin O(1). Denote x is the solution vector,
x = . . . , u h (xi , y j ), . . .
T
(5.7)
,
assumption A2* leads to unT x O(1),
(5.8)
where (σn , un ) is the minimal eigenpair of the discrete matrix A, which is symmetric and positive definite. Note that Eq. (5.6) resembles (5.8). Eq. (5.6) implies that the solution u is rich in the eigenfunction w min , and Eq. (5.8) also implies that the solution x is rich in the eigenvector un . For the real application, the PDE solution u can only vary in a certain region. Therefore, the opposite case of (5.6) is very seldom, and the analysis of the effective condition number in this paper is useful. 5.2. Improvements of Theorem 4.1 In Section 4, the triangles i j near Γ D are to be right triangles. In the rest of this section, such a limitation can be removed. Denote the vector b g = {. . . , αi j g i j , . . .}, where
(i , j ) ∈ Γ D∗ ,
(5.9)
αi j are constants. We have the following lemma.
Lemma 5.2. For the quasiuniform i j near Γ D∗ , the constants αi j in (5.9) have the bound,
|α i j | C ,
(5.10)
where C is a constant independent of h. Proof. Define the energy I (v ) =
1
|∇ v |2 ds −
2 S
f v ds −
g ∗ v d.
Γ N∗
S
Also denote the finite dimensional collection V h∗ (⊂ H 1 ( S )) of all continuous piecewise linear functions on S h1 = satisfying v |Γ ∗ = g. Then the FEM can be rewritten as: To seek u h ∈ V h∗ such that
(5.11)
ij
i j
D
I (u h ) = min I ( v ), v ∈ V h∗
i.e.,
∂ I (v ) = 0. ∂ v i, j
(5.12)
From Zienkiewicz [17], the linear interpolant functions are given by v = ϕi v i + ϕ j v j + ϕm v m , where i, j, and m are three nodes of , and
ϕi =
ai + b i x + c i y 2| |
,
ai = x j ym − xm y j , b i = y j − ym , c i = xm − x j ,
1 1 | | = area( ) = xi 2 y
i
1 xj yj
. ym
1 xm
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
1791
Similarly, the coefficients a j , b j , c j , am , bm and cm are obtained using the counter-clock rotation, see Zienkiewicz [16]. The discrete energy of (5.11) is given by I (v ) =
1 2
∂v ∂x
2
+
∂v ∂y
2 dx dy −
f v ds −
g ∗ v d.
(5.13)
∩Γ N∗
In (5.13), only the first term on the right hand is related to the vector b g , denoted by I1(v ) =
=
1 2
1
8| |
∂v ∂x
2 +
∂v ∂y
2 dx dy
(bi v i + b j v j + bm v m )2 + (c i v i + c j v j + cm v m )2 ,
and from (5.12)
∂ I 1 ( v ) 1 2 = b i + c 2i v i + (b i b j + c i c j ) v j + (b i bm + c i cm ) v m . ∂ vi 4| |
(5.14)
Since are quasiuniform, there exist the bounds: |b i | = | y i − ym | Ch, |c i | = |xm − y i | Ch and O(h2 ). Hence the coefficients in (5.14) have the bounds
|bi b j + c i c j | C,
|bi bm + c i cm | C.
(5.15)
In (5.14), when i ∈ S 1 and j ∈ Γ D (or m ∈ Γ D ) then v i = g i (or v m = gm ) from the Dirichlet boundary conditions. Then the coefficients of g i obtained form (5.14) are O(1), based on (5.15). Since the boundary nodes (i , j ) ∈ Γ D are connected to finite quasiuniform . Hence the coefficients of g i j in (5.9) satisfy (5.10). This completes the proof of Lemma 5.2. 2 Below, let us renovate Theorem 4.1. First, assumption A2 is replaced by A2∗ in (5.6). From Lemma 5.2, we can easily show
b g Ch−0.5 g 0,Γ D∗ . To remove the limitation of right triangles i j near Γ D∗ by the common assumption: quasiuniform elements, we need only prove
yhmin
T
bg ≈ − Γ∗
∂ w min g d, ∂n
(5.16)
D
to replace (4.13), since other bounds (4.14)–(4.16) also hold. First, consider the Neumann boundary condition on Γ N∗ , which can be obtained approximately by (4.3). Denote the element node A (i , j ) ∈ Γ N∗ . The integration of elements in (4.3) related to node A are given by 3 i =1
∇ uh ∇ v ds =
3 i =1
i
g ∗ v d,
i ∩Γ N∗
which gives
Lu i j +
S i +1 + S i 2
g ∗A v A = 0,
where i +1 and i are not necessarily to be the right triangles as 1 and 2 as shown in Fig. 4(b). Since v A is arbitrary, we obtained the boundary equation u A ∈ Γ N∗ ,
S i +1 + S i ∗ S i +1 + S i ∂ v Lu i j = − gA = − . 2 2 ∂n A
Next, consider the node A ∈ Γ D∗ . We obtain from (5.17) 3 i =1
i
S i +1 + S i ∂ v ∇ uh ∇ v ds = ( Lu i j · v ) A = ( Lu i j ) · g A ≈ − · gA. 2 ∂n A
(5.17)
1792
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
Hence we obtain
yhmin
T
b g = −c 0 h
S i +1 + S i ∂ w min ∂ w min g g d. ≈ −c 0 h · 2 ∂ n ∗ ∂n (i j )∈Γ ∗
(i , j )∈Γ D
D
Γ D∗
This is the desired result (5.16). We summarize this new result as a theorem. Theorem 5.1. Let σ ∈ [ 14 , 1), p 1 > σ1 and the local mesh refinements with (2.7) hold. Suppose that the triangles i j in S 1 are quasiuniform. Under assumption A2∗ (i.e. A2), the simplified effective condition number has the bound (4.19). 6. Numerical experiments In this section, first we consider Motz’s problem (see Fig. 5)
u = 0 in S , u = 0 on O D ,
u = 500 on A B , ∂u = 0 on BC ∪ C D ∪ O A , ∂n
(6.1)
where S = {(x, y ) | −1 < x < 1, 0 < y < 1}. The true functions are found as in [6] u=
∞
1
di r i + 2 cos i +
i =0
1 2
θ
in S ,
(6.2)
with the computed expansion coefficients di . The origin O is a singular point due to the intersection between the Neumann and the Dirichlet boundary conditions. To deal with the singularity, we choose the local refinements oftriangulation near O as shown in Fig. 6. Let the solution domain S be divided into small triangles i j such that S = i , j i j . The nodes
(r i , θ j ) of the small triangles i j near O are chosen as (2.7). In the our computation, we use the larger p 1 = 4(> σ1 = 2). √ where N a and N are the division numbers along the radial direction and the semi-circle, Denote h = 2/ N a and θ = π N respectively. The numerical results are listed in Table 1. Since Cond_eff Cond, and since Eq. (1.4) is easier than (1.6) in computation, we provide only the numerical data of Cond_eff. In the table, λmax (A) and λmin (A) are the maximal and the minimal eigenvalues of A, respectively, and Ratio(Cond_eff) = Cond_eff|2N /Cond_eff|2N −1 . From Table 1 we can see
Fig. 5. Motz’s problem.
Table 1 Errors and condition numbers for the non-homogeneous Dirichlet boundary condition with p 1 = 4. N (= N a )
4
8
16
32
64
128
maxi , j |εi , j | u − uh 0, S u − uh 1, S λmax (A) λmin (A) b x Cond(A) Cond_eff Ratio(Cond_eff)
3.03(1) 2.12(1) 1.10(2) 3.97(1) 5.44(−1) 5.04(2) 4.38(2) 7.30(1) 2.12 –
8.47 6.77 5.95(1) 8.27(1) 1.50(−1) 6.88(2) 1.03(3) 5.51(2) 4.45 2.099
4.67 2.01 30.9(1) 1.68(2) 3.76(−2) 9.92(2) 2.16(3) 4.47(3) 1.22(1) 2.742
2.44 5.58(−1) 15.7(1) 3.36(2) 9.42(−3) 1.42(3) 4.35(3) 3.57(4) 3.47(1) 2.844
1.04 1.47(−1) 7.88 6.65(2) 2.39(−3) 2.03(3) 8.72(3) 2.78(5) 9.74(1) 2.807
3.87(−1) 3.72(−2) 3.95 1.32(3) 6.03(−4) 2.88(3) 1.74(4) 2.19(6) 2.74(2) 2.813
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
1793
Fig. 6. Partition for Motz’s problem with N = 32.
Table 2 Errors and condition numbers for the non-homogeneous Neumann boundary condition with p 1 = 4. N (= N a )
4
8
16
32
64
128
maxi , j |εi , j | u − uh 0, S u − uh 1, S λmax (A) λmin (A) b x Cond(A) Cond_eff Ratio(Cond_eff)
4.72(1) 4.34(1) 1.08(2) 3.97(1) 2.31(−1) 2.41(2) 7.62(2) 1.72(2) 1.37 –
1.54(1) 1.28(1) 6.00(1) 8.27(1) 7.70(−2) 2.09(2) 1.32(3) 1.07(3) 2.06 1.5036
4.52 3.62 3.12(1) 1.68(2) 1.91(−2) 1.60(2) 2.41(3) 8.79(3) 3.47 1.6845
2.26 9.66(−1) 1.58(1) 3.36(2) 5.04(−3) 1.17(2) 4.60(3) 6.67(4) 5.05 1.4553
9.93(−1) 2.49(−1) 7.90 6.65(2) 1.30(−3) 8.46(1) 8.96(3) 5.12(5) 7.26 1.4376
3.74(−1) 6.28(−2) 3.96 1.32(3) 3.30(−4) 6.05(1) 1.77(4) 4.00(6) 1.03(1) 1.4187
max |εi , j | = O h2 | ln h| , i, j
u − uh 0, S = O h2 , u − uh 1, S = O(h), −1.5 Cond_eff = O h , −1 λmax (A) = O h , λmin (A) = O h2 , Cond(A) = O h−3 ,
(6.3) (6.4) (6.5) (6.6) (6.7)
where εi , j = u (i , j ) − u h (i , j ) and u (i , j ) = u (r i , θ j ). The asymptotic rates of errors in (6.4) are the optimal, to coincide with Theorem 2.1. The effective condition number in (6.5) also verifies Theorem 5.1 for the non-homogeneous Dirichlet boundary condition. Next, we consider the non-homogeneous Neumann boundary condition (see Fig. 5),
u = 0 in S , u=0
on O D ,
∂u = g ∗ on A B , ∂n
∂u = 0 on BC ∪ C D ∪ O A , ∂n
(6.8)
where g ∗ (= 0) is the normal derivative of (6.2) on A B, to replace the non-homogeneous Dirichlet condition u = 500 on A B in (6.1). The numerical results are listed in Table 2. From Table 2 we can see clearly
max |εi , j | = O h2 | ln h| , i, j
(6.9)
1794
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
u − uh 0, S = O h2 , u − uh 1, S = O(h), −0.5 , Cond_eff = O h −1 λmax (A) = O h , λmin (A) = O h2 , Cond(A) = O h−3 .
(6.10) (6.11) (6.12) (6.13) O(h −0.5 ) again
The asymptotic rates of numerical errors are also optimal, and the effective condition number Cond_eff = verifies Theorem 5.1 for the non-homogeneous Neumann boundary condition. From the analysis and the computation in this paper, the instability is not severe for the FEM using local refinements of elements near the singularity, based on 2 ) in [13]. Moreover, we obtain the ratio from Table 2 for Cond_eff = O(1), O(h−1.5 ) or O(h−0.5 ), in contrast to Cond = O(h− min N = 128, Cond(A) Cond_eff
=
4.01(6) 1.04(1)
= 3.85(5).
(6.14)
The above facts also confirm that the effective condition number is a better criterion for numerical PDE than the traditional condition number. Remark 6.1. In Bank and Scott [1], the local mesh refinements are limited under certain restrictions, to give an improved bound of Cond,
Cond(A) Ch−2 1 + ln
hmax
hmin
,
(6.15)
where C is a constant independent of h. Evidently, Eq. (6.15) does not coincide with the numerical results in (6.7) and (6.13):
Cond(A) = O h−3 .
(6.16)
Since the mesh family based on (2.7) is degenerate, thus the assumption of non-degenerate triangulations in [1] is not complied with. Note that the local mesh refinements (2.7) with the simplest node location is much simpler than those in [1,14]. For (2.7), a further exploration displays the new bound (details appear elsewhere): Cond(A) Ch−3 | ln hmin |,
(6.17)
where C is a constant independent of h. Eq. (6.17) is consistent perfectly with (6.16). Cond in (6.17) is large, although it is acceptable in application. In summary, based on effective condition numbers, but not on the traditional condition number, the instability is not severe. This paper provides a stability justification, to support the important adaptive mesh refinements in FEM existing for three decades. Remark 6.2. Motz’s problem is a typical corner singularity with mixed boundary conditions. The other kind of corner singularity is for the Dirichlet condition, where S 0 a concave polygon with Θ ∈ (π , 2π ]. The harmonic functions in S 0 with u = 0 at θ = 0 and θ = Θ have the expansion u=
∞
∗
c i r μi sin(μ∗i θ),
(r , θ) ∈ S 0 ,
(6.18)
i =0
where c i are the true coefficients, and
μ∗i =
iπ
Θ
,
i = 1, 2, . . . .
(6.19)
Hence when r → 0 we have
∗ u (r , θ) = O r σ ,
(6.20)
where
σ ∗ = μ∗0 =
π Θ
1 2
,
(6.21)
to give σ ∗ ∈ [ 12 , 1) for π < Θ 2π of the concave polygon. Compared with (2.5) for Motz’s problem, the corner singularity of Dirichlet conditions is weaker. Moreover, the Motz’s singularity is, indeed, of the singularity of concave corners with Dirichlet conditions. Let us consider an interior crack at O D = (−1 < x < 0) ∩ ( y = 0) within S ∗ = {(x, y ) | −1 < x < 1, −1 < y < 1} with the Dirichlet condition on ∂ S ∗ where u = 0 at O D (see Fig. 7). Suppose the solution in S ∗ is symmetric with respect to the axis X , i.e., u (x, y ) =
Z.-C. Li, H.-T. Huang / Applied Numerical Mathematics 59 (2009) 1779–1795
1795
Fig. 7. The interior crack problem.
u (x, − y ). By using this symmetry, the interior crack problem may be simplified to the problem in its half domain S (see Fig. 5) with the following mixed type of Dirichlet and Neumann boundary conditions at y = 0
∂u = 0, on O A , ∂n u = 0, on O D .
(6.22) (6.23)
Eqs. (6.22) and (6.23) are exactly the same corner boundary conditions (6.1) for Motz’s problem. In this paper, we choose Motz’s problem as the corner singularity for easy comparisons with our previous study for singularity problems in [6,11]. Acknowledgements We are grateful to Professors Robert Beauwens, Peter Jimack and the referees for their valuable comments and suggestions. We are also indebted to Mr. D.T.C. Wu for the computation in Section 6. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
R.E. Bank, R. Scott, On the conditioning of finite element equations with highly refined meshes, SIAM J. Numer. Anal. 26 (1989) 1383–1394. F.C. Chan, D.E. Foulser, Effectively well-conditioned linear systems, SIAM J. Stat. Comput. 9 (1988) 963–969. G.H. Golub, C.F. van Loan, Matrix Computations, second ed., The Johns Hopkins, Baltimore and London, 1989. J.A. Gregory, D. Frishelov, B. Schiff, J.R. Whiteman, Local mesh refinement with finite elements for elliptic problems, J. Comput. Phys. 29 (1978) 133–140. R.S. Lehman, Developments near an analytic corner of solution of elliptic differential equation, J. Math. Meth. 8 (1959) 727–760. Z.C. Li, Combined Methods for Elliptic Equations with Singularities, Interfaces and Infinities, Kluwer Academic Publishers, Dordrecht, Boston, 1998. Z.C. Li, Error analysis for boundary approximation methods solving eigenvalue problems, J. Comput. Appl. Math. 200 (2007) 231–254. Z.C. Li, C.S. Chien, H.T. Huang, Effective condition number for finite difference method, J. Comput. Appl. Math. 198 (2007) 208–235. Z.C. Li, H.T. Huang, Effective condition number for numerical partial differential equations, Numer. Linear Algebra Appl. 15 (2008) 575–594. Z.C. Li, H.T. Huang, Effective condition number for the finite element method using local mesh refinements, Technical report, Department of Applied Mathematics, National Sun Yet-sen University, Kaohsiung, Taiwan, 2007. Z.C. Li, T.T. Lu, H.Y. Hu, A.H.D. Cheng, Trefftz and Collocation Methods, WIT Press, Southampton, Boston, 2008. J.R. Rice, Matrix Computations and Mathematical Software, McGraw-Hill Book Company, New York, 1981. G. Strang, G.J. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cliffs, NJ, 1973. L.B. Wahlbin, Local behavior in finite element method, in: P.G. Ciarlet, J.L. Lions (Eds.), in: Finite Element Methods, vol. I, North-Holland, 1991, pp. 352– 522. N.M. Wigley, Asymptotic expansions at a corner of solutions of mixed boundary value problems, J. Math. Meth. 4 (1964) 549–579. J.H. Wilkinson, The Algebraic Eigenvalue Problem, Clarendon Press, Oxford, 1965, 191 p. O.C. Zienkiewicz, The Finite Element Method, third ed., McGraw-Hill, London, 1977.