- Email: [email protected]
Stability analysis for the penalty plus hybrid and the direct Trefftz methods for singularity problems
ARTICLE IN PRESS
Engineering Analysis with Boundary Elements 31 (2007) 163–175 www.elsevier.com/locate/enganabound
Stability analysis for the penalty plus hybrid and the direct Trefftz methods for singularity problems Zi-Cai Lia,b,c,, Hung-Tsai Huangd, Jin Huanga,e, Leevan Linga,f a
Department of Applied Mathematics, National Sun Yat-sen University, Kaohsiung 80424, Taiwan Department of Computer Science and Engineering, National Sun Yat-sen University, Kaohsiung 80424, Taiwan c National Center for Theoretical Sciences, Taiwan d Department of Applied Mathematics, I-Shou University, Kaohsiung County 840, Taiwan e College of Applied Mathematics, University of Electronic and Science Technology of China, ChengDu, China f Department of Mathematics, Hong Kong, Baptist University, Kowloon Tong, Hong Kong
b
Received 20 October 2005; accepted 18 April 2006 Available online 11 July 2006 This paper is dedicated to S. Christiansen on the occasion of his 70th birthday
Abstract For solving the linear algebraic equations Ax ¼ b, the new stability analysis is made based on the effective condition number Cond_eff. The Cond_eff may provide a better upper bound of relative errors of x resulting from the rounding errors of b, than the traditional condition number Cond, which with too large value is, in many times, misleading. In this paper, we apply the effective condition number to the Trefftz methods (TMs) for Poisson’s equations with singularities. Two TMs, such as the penalty plus hybrid TM and the Lagrange multiple (i.e., direct) TM, are studied. We focus on the stability analysis of the solutions when the optimal superconvergence is achieved. Whenpsolving Motz’s problem, the benchmark of singularity problems, by the penalty plus hybrid TM, we have derived that Cond_eff ¼ ffiffiffi OðNð 2ÞN Þ and Cond ¼ OðN 2 2N Þ, where N is the number the singular particular functions used. For solving Motz’s problem by the pffiffiffiof N direct TM, we have derived that Cond=Cond_eff ¼ OðN 2 Þ. Numerical experiments are provided to verify the stability analysis made. In summary, the two TMs are efficient, but the penalty plus hybrid TM is more recommended, due to simplicity of algorithms without extra-variables and less limitations in applications. r 2006 Elsevier Ltd. All rights reserved. Keywords: Effective condition number; Condition number; Stability analysis; Trefftz method; Penalty plus hybrid technique; Lagrange multiplier technique; The direct Trefftz method; Singular problem; Motz’s problem
1. Introduction For solving the linear algebraic equations Ax ¼ b, the traditional condition number is defined by Cond ¼ s1 =sn in [1,2], where s1 and sn are the maximal and the minimal singular values of matrix A, respectively. The condition number is used to provide the bounds of the relative errors from the perturbation of both A and b. Such a Cond can only be reached by the worst situation of all rounding Corresponding author. Department of Applied Mathematics and Department of Computer Science and Engineering, National Sun Yat-sen University, Kaohsiung 80424, Taiwan. E-mail address: [email protected] (Z.-C. Li).
0955-7997/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2006.04.005
errors and all b. For the given b, however, the real relative errors may be smaller, and even much smaller than the Cond. Therefore, the effective condition number was introduced and studied in Chan and Foulser [3]; subsequently it was applied to boundary value problem in Christiansen and Hansen [4], and to the boundary integral equation in Christiansen and Saranen [5]. Other kinds of effective condition number are given in Banoczi et al. [6], where the effective condition number is applied to the Gaussian elimination and the QR factorization. The main references for condition number are given in [7]. Let us consider the linear algebraic equations Ax ¼ b,
(1.1)
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
164
where A 2 Rnn , x 2 Rn and b 2 Rn are the unknown and known vectors, respectively. When there occurs a perturbation of b, the errors of x satisfy (1.2)
Aðx þ DxÞ ¼ b þ Db.
The values of Cond are used to measure the relative errors of x, which for (1.2) are given by kDxk kDbk pCond , kxk kbk
(1.3)
where kxk is the Euclidean norm and the matrix norm kAk ¼ supxa0 kAxk=kxk. Let the matrix A 2 Rnn be symmetric and positive definite. The eigenvalues are arranged in a descending order, l1 Xl2 X Xln 40, and the eigenvectors ui satisfy Aui ¼ li ui , where fui g are orthogonal, with uTi uj ¼ dij , where dij ¼ 1 if i ¼ j and dij ¼ 0 if iaj. The traditional condition number is given by Cond ¼ l1 =ln . In [2,5,7], the effective condition number is defined by Cond_eff ¼
kbk ¼ ln kxk
kbk qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , Pn 2 2 ln b =l i i¼1 i
(1.4)
where bi ¼ uTi b, also to provide kDxk kDbk pCond_eff . kxk kbk
(1.5)
In this paper, we will apply the effective condition number to the Trefftz methods (TMs) for solving Poisson’s equation, in particular for Motz’s problem, the benchmark of singularity problems. When the particular solutions or the piecewise particular solutions satisfying Poisson’s equation are chosen, the approximate solutions are obtained by a linear combination of these particular solutions, and their expansion coefficients are sought by satisfying the exterior or interior boundary conditions. This is called the TM. To couple the boundary conditions, different techniques lead to different TMs. First, we use the penalty plus hybrid techniques in Li [8], to give the penalty plus hybrid TM. Next, we choose the Lagrange multiplier techniques to give the Lagrange multiplier TM (or called the direct TM). In the direct TM, the multipliers are also regarded as variables. For solving Laplace’s equation, the constants are not permitted in the particular solutions in the direct TM. Moreover, if the uniform particular solutions cannot be found on the entire domain, using the piecewise particular solutions is necessary. When the solution domain S is divided into several subdomains S i such that S ¼ [i S i and Si \ S j ¼ ; for iaj, such a limitation without constants excludes the interior subdomains Si with qSi \ qS ¼ ;. However, no such a limitation exists for the penalty plus hybrid TM. In the direct TM, there often occurs the puzzle that the Cond is extremely huge, but the approximate solutions are convincingly accurate. Based on the Cond, such numerical solutions should be distrusted, and ignored. However, one often uses them by refusing the warning of Cond. In this paper, the new effective condition number can clarify this
puzzle well. In fact, the Cond_eff from the direct TM and the penalty plus hybrid TM is reasonably large, to revalidate the numerical solutions with a doubt from huge Cond. Moreover, since the exponential convergence rates are proved in Li et al. [9] for the two TMs without integration approximation, both the direct TM and the penalty plus hybrid TM are very efficient. However, the penalty plus hybrid is more advantageous due to its simplicity and less limitation in application. For numerical methods of partial differential equations, stability analysis is as important as error analysis. The proof for the former is more challenging than the latter. This is why there have appeared numerous reports on error analysis, but not so many reports on stability analysis. The effective condition number is a new trend of stability analysis, to provide a better bound of relative errors from rounding errors. Inspired by the pioneer work of [3–5], we apply the effective condition number for different numerical methods, and find several new discoveries. The related work to this paper are [10] for Adini’s elements and [11] for penalty combinations for singularity problems. Since Adini’s elements involve the derivatives ux and uy , the estimates of the minimal eigenvalue of the stiffness matrix are new and intriguing, the new proofs for bounds of both Cond and Cond_eff are provided in [10]. Since the condition number is huge by the penalty combination of TM coupled with high FEM for singularity problems, such a penalty combination was declined. However, the effective condition number is reasonably large so that the penalty combination is reconciled. The strict analysis of stability and superconvergence is explored in [11]. The a priori estimates of upper bounds of effective condition number may enhance the stability analysis of numerical methods. This paper is organized as follows. In Sections 2 and 3, the penalty plus hybrid TM and the direct TM are described, respectively, accompanied with a brief error analysis for Motz’s problem. In Section 4, new stability analysis is made, and the bounds of Cond_eff, Cond and Cond=Cond_eff are derived in detail. In Section 5, numerical experiments are carried to verify the new stability and the error analysis made. 2. Penalty plus hybrid techniques To match the exterior boundary conditions and different particular solutions in the TM, or matching the TM with other methods, an effective coupling strategy is essential in order to yield optimal convergence rates and good stability of the numerical solutions. The coupling can be accomplished by using additional integrals along the common boundary of different particular solutions. 2.1. General algorithms Below, we provide the penalty plus hybrid couplings for combinations of the TM-FEM. Details are given in
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
[8,9,12,13]. Consider the elliptic equation q qu q qu p p þ cu ¼ f ; ðx; yÞ 2 S, qx qx qy qy
where (2.1)
ðvÞ I ða;bÞ h
with the mixed type of the Dirichlet and Robin boundary conditions u ¼ g1
on GD ,
qu þ qu ¼ g2 qn
(2.2)
on GN ,
(2.3)
where the domain S is bounded (i.e., a polygon) with the exterior boundary, G ¼ GD [ GN , and MeasðGD Þ40 for uniqueness of the solution. The functions p; c; f ; q; g1 and g2 are sufficiently smooth, and c ¼ cðx; yÞX0;
q ¼ qðx; yÞX0,
p ¼ pðx; yÞXp0 40, where p0 is a constant. The solution of problem (2.1)–(2.3) can be equivalently expressed by minimizing a quadratic functional IðvÞ:
1 ¼ 2
165
ZZ
ðprv rv þ cv2 Þ ds ZZ 1 þ ðprv rv þ cv2 Þ ds 2 S2 Z Z 1 2 2 þ qv d‘ þ w ðv g1 Þ2 d‘ 2 GN GD Z Z qv p ðv g1 Þ d‘ þ w2 ðvþ v Þ2 d‘ qn GD G0 Z qvþ qv þb p a ðvþ v Þ d‘ qn qn G0 ZZ Z fv ds g2 v d‘, ð2:6Þ S1
S
GN
where the parameters a and b satisfy a þ b ¼ 1, w2 ¼ Pc ðL þ NÞs is the penalty weight, s40, and Pc ð40Þ is a constant which is large enough but independent of L and N. Eq. (2.6) includes various combined methods, based on different a; b; Pc and s.1 When fFi g and fCi g are chosen as the particular solutions satisfying Eq. (2.1), Eq. (2.5) leads the TMs of variant combinations.
IðuÞ ¼ min IðvÞ, v2H 1 ðSÞ
2.2. Application for Motz’s problem
where the quadratic functional is ZZ 1 ½prv rv þ cv2 ds IðvÞ ¼ 2 S Z Z Z Z 1 2 qv d‘ fv ds g2 v d‘, þ 2 GN S GN
Motz’s problem is the benchmark of singularity problems. The problem was first discussed by Motz [14] in 1947 for the relaxation method. Since then many researchers have selected Motz’s problem as a prototype for verifying efficiency of numerical methods (see [8]). Motz’s problem solves the Laplace equation on the rectangle S ¼ fðx; yÞ j 1oxo1; 0oyo1g
and H 1 ðSÞ ¼ fvjv; vx ; vy 2 L2 ðSÞ; and vjGD ¼ g1 g is a subset of the Sobolev space. Let S be divided by a piecewise straight line G0 into S 1 and S2 . The admissible functions can be written as follows: 8 L P > > bi Fi in S1 ; >
v¼
i¼1
> > þ > :v ¼
N P
(2.4)
ai C i
in S 2 ;
i¼1
where fFi g and fCi g are analytic, complete and linearly independent basis functions on S1 and S2 , respectively, and ai and bi are unknown coefficients to be sought. Note that the admissible functions vþ and v in (2.4) may not satisfy the elliptic equation (2.1) exactly. Define a space H ¼ fvjv 2 L2 ðSÞ; v 2 H 1 ðS 1 Þ; v 2 1 H ðS 2 Þg. Let V h H be a finite-dimensional collection of the functions (2.4). In this section, we present the combination of methods using the penalty plus hybrid techniques: to seek an approximate solution uh 2 V h such that I ða;bÞ ðuh Þ h
¼
min I ða;bÞ ðvÞ, h v2V h
(2.5)
Du ¼
q2 u q2 u þ ¼0 qx2 qy2
in S,
(2.7)
with the mixed Neumann–Dirichlet boundary conditions (see Fig. 1): ujOD ¼ 0; ujAB ¼ 500, qu qu qu ¼ ¼ ¼ 0. qyOA qyBC qxCD
ð2:8Þ ð2:9Þ
Note that there exists a singularity at the origin (0, 0) due to the intersection point of the Neumann–Dirichlet boundary conditions. In fact, the singular solutions of (2.7)–(2.9) are found as 1 X 1 uðr; yÞ ¼ Di riþ1=2 cos i þ y, (2.10) 2 i¼0 where Di are the true expansion coefficients, and ðr; yÞ are the polar coordinates with the origin at ð0; 0Þ (see Fig. 1). Since the convergence radius, R ¼ 2, is analyzed in Rosser 1
In some case (e.g., Pc 40, and a ¼ 1 and b ¼ 0), Eq. (2.5) may not hold, but the variational equation of I hð1;0Þ ðvÞ leads to the desired equation of uh .
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
166
of S ¼ SR ¼ fðr; yÞ j 0proR; 0pyppg, it is better to choose w2 ¼ OðN þ 1Þ, i.e., w2 ¼ Pc ðN þ 1Þ, where Pc is a constant independent of N, see [12]. Define the norm
Y un = 0
C
B
kvkh ¼ fjvj21;S þ w2 kvk0;AB g1=2 .
(2.14)
We have the following theorem. un = 0
u = 500
X D
u=0
O
un = 0
A
Theorem 2.1. Suppose that there exists a constant C independent of N such that qv pCðN þ 1Þkvk0;AB 8v 2 V N . (2.15) qn 0;AB Then there exists the error bound for uPH N
Fig. 1. Motz’s problem.
ku uPH N kh pC inf ku vkh . v2V N
and Papamichael [15], the series expansions (2.10) are most suited to the entire solution domain S. Hence the admissible functions of finite terms, N X 1 uN ðr; yÞ ¼ D~ i rðiþ1=2Þ cos i þ y, (2.11) 2 i¼0 with the unknown coefficients D~ i , are most efficient for numerical Motz’s solutions, to yield the exponential convergence rates OðecN Þ, where c is a positive constant. When functions (2.11) are chosen, Eq. (2.7), ujOD ¼ 0 and qu=qyjOA ¼ 0 are satisfied automatically. Then the coefficients D~ i are sought by the collocation equations of the rest of boundary conditions in (2.8) and (2.9), which is called the TM in this paper. The very accurate solution was obtained in double precision, with the 35 leading coefficients, D~ 0 D~ 34 provided in [16], while an error at the power of D~ 31 was discovered by Lucas and Oh [17]. The very accurate solution can be regarded as a true solution for testing other numerical methods. In Georgiou et al. [18], the same singular functions as in (2.11) are chosen, but the boundary conditions are matched by Lagrange multipliers method to be discussed in Section 3. In this section, we will provide the error bounds for the penalty plus hybrid TM and the direct TM. Based on (2.5), the numerical solutions uPH by the N penalty plus hybrid TM may be sought by minimizing the following energy: TðuPH N Þ ¼ min TðvÞ, v2V N
(2.12)
where V N is the finite-dimensional collection of (2.11) and Z Z 1 qv qv v d‘ ðv 500Þ d‘ TðvÞ ¼ 2 GN qn qn AB Z 2 þw ðv 500Þ2 d‘, ð2:13Þ AB
where GN ¼ OA [ BC [ CD (see Fig. 1) and w40. Below, let us consider how to choose the weight w. First, w2 should be chosen to balance the first and the second terms on the right-hand side in (2.13). From the simple case
(2.16)
Moreover, suppose that w2 ¼ Pc ðN þ 1Þ, and that Pc is chosen large enough but independent of N. There exists the error bound for uPH N ku uPH N kh pC inf ku vkh v2V N pffiffiffiffiffi ¼ C 1 N aN ; 0oao1,
ð2:17Þ
where C 1 is also bounded and independent of N. Theorem 2.1 displays the exponential convergence rates by the penalty plus hybrid TM for Motz’s problem. Let us consider the integration approximation in (2.12). The penalty plus hybrid TM involving integral approximation is designed to seek u~ PH N 2 V N such that ~ ~ u~ PH Þ ¼ min TðvÞ, Tð N v2V N
(2.18)
where Ze Z 1 e qv qv ~ v d‘ ðv 500Þ d‘ TðvÞ ¼ 2 GN qn qn AB Ze þ w2 ðv 500Þ2 d‘,
ð2:19Þ
AB
R R where e is an approximation of by some rules. We have the following theorem. Theorem 2.2. For the solutions by the penalty plus hybrid TM, there exists the error bound, ku u~ PH N k1;S pC inf ku vk1;S v2V N ( Z Ze ! 1 qu w d‘ þ C sup qn w2V N kwk1;S GN GN ) Z ! Z e qu þ ð2:20Þ w d‘ . AB qn AB The proofs of Theorems 3.1 and 3.2 are given in [9,12]. Note that the additional errors from the integration approximation in (2.20) are analogous to those in FEM. Since the integration rules are formulated based on k-order polynomials, such as the Gaussian and Newton–Cotes
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
rules in [19], only the polynomial convergence rates can be obtained from (2.20). 3. The Lagrange multiplier techniques 3.1. General algorithms R We introduce an additional integral G0 lðvþ v Þ d‘, where l is a continuous function of Lagrange multipliers. We describe the Lagrange multiplier coupling as: to seek ðuh ; lÞ 2 V N R such that Bðuh ; v; l; mÞ ¼ f ðvÞ;
ðv; mÞ 2 V N R,
167
The variational equation, qF ðvÞ=qv ¼ 0, leads to (3.5). In Li [8], this is referred to the Lagrange multiplier method, which was just the direct TM in Jin and Cheung [21], and Kita and Kamiya [22]. The Lagrange multiplier method was first introduced by Babusˇ ka [23] to treat the constraint Dirichlet boundary condition as a natural boundary condition, and to relax the limitation on the admissible functions used. Since then the techniques of Lagrange multipliers have drawn much attention, the related references are given in [8]. 3.2. Application for Motz’s problem
(3.1)
where ðv; mÞ 2 V N R denotes that v 2 V N and m 2 R, and ZZ ðp 5 u 5v þ cuvÞ ds Bðu; v; l; mÞ ¼ S1 ZZ þ ðp 5 u 5v þ cuvÞ ds
For Motz’s problem, we may regard the Dirichlet conditionRRujAB ¼ 500 as a constraint to minimize the energy, 12 S jrvj2 ds. Then define a functional ZZ Z 1 jrvj2 ds lðv 500Þ d‘, (3.9) IðvÞ ¼ 2 S AB
S2
þ Dðu; v; l; mÞ,
ð3:2Þ
ZZ fv ds.
f ðvÞ ¼
(3.3)
S
In (3.2), the boundary integrals are Z Z quv d‘ lv d‘ Dðu; v; l; mÞ ¼ GN GD \S 2 Z Z mðu gÞ d‘ lðvþ v Þ d‘ GD \S 2 G0 Z mðuþ u Þ d‘, ð3:4Þ G0
G
G
where l ¼ qu=qn is also treated as unknown, and l; m 2 H 1=2 ðGÞ. H 1=2 ðGÞ is the negative norm in the Sobolev space defined by (see [20]) kvk1=2;G ¼ kvk20;G Z Z
2
ðvðPÞ vðQÞÞ d‘ðPÞ d‘ðQÞ kP Qk2 G G R j G uv d‘j ¼ sup . va0 kvk1=2;G þ
kuk1=2;G
S
AB
AB
8ðv; mÞ 2 H 10 ðSÞ H 1=2 ðABÞ,
ð3:10Þ
where H 10 ðSÞ ¼ fv; vx ; vy 2 L2 ðSÞ and vjOD ¼ 0g. Let lL be the L-order polynomials on AB. We choose the Chebyshev polynomials, lL ¼
L X
A~ i T i ð1 2yÞ;
0pyp1,
(3.11)
i¼0
where the Lagrange multiplier l has the true solution, l ¼ p qu=qnjG0 . In (3.2), the l is treated as an extra variable. Consider the Laplace equation with the Dirichlet condition u ¼ g on G, and no interior boundary exists. From (3.1) with p ¼ 1 and c ¼ 0 we have ZZ Z Z 5 u 5v ds mðu gÞ d‘ lv d‘ ¼ 0, (3.5) S
with the Lagrange multiplier l. From (3.9), the stationary condition for (3.9) is given by: to seek ðu; lÞ 2 H 10 ðSÞ H 1=2 ðABÞ such that ZZ Z Z ru rv ds lv d‘ mðu 500Þ d‘ ¼ 0,
In fact, we may define a functional ZZ Z 1 F ðvÞ ¼ 5 v 5v ds mðv gÞ d‘. 2 S G
1=2
where A~ i are the coefficients to be sought, and the Chebyshev polynomials are defined by T i ðxÞ ¼ cosði arccosðxÞÞ;
1pxp1.
(3.12)
Denote by V N V L the collection of finite dimensions of (2.11) and (3.11). The discrete Lagrange multiplier method (i.e., the direct TM) is given by: to seek ðu~ N ; lL Þ 2 V N V L such that Aðu~ N ; vÞ þ bðu~ N ; v; lL ; mÞ ¼ 0
8ðv; mÞ 2 V N V L ,
(3.13)
where ZZ Aðu; vÞ ¼
ru rv ds, S
,
ð3:6Þ
Z bðu; v; l; mÞ ¼
(3.7)
(3.8)
Z lv d‘
AB
mðu 500Þ d‘.
(3.14)
AB
The multiplier method for Motz’s problem was discussed in Georgiou et al. [18], where the piecewise k-order polynomials were chosen, to provide the polynomial convergence rates. Below we give a brief justification for Laplace’s equation.
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
168
From [8] we obtain the following theorem.
4. Stability analysis for Motz’s problem
Theorem 3.1. Suppose the following three assumptions hold.
4.1. The penalty plus hybrid TM
(A1) For Aðu; vÞ, there exist two positive constants C 0 and C independent of N such that,
For simplicity, we consider the TMs for Motz’s problem in Sections 2.2 and 3.2. The approaches for stability analysis can be easily extended to other Poisson’s equations by the TMs. Eq. (2.12) leads to the linear algebraic equations,
C 0 kvk21;S pAðv; vÞ 8v 2 V N , jAðu; vÞjpCkuk1;S kvk1;S
ð3:15Þ
8v 2 V N .
ð3:16Þ
R (A2) For AB mv, there exits the Ladyzhenskaya–Babuska–Brezzi (LBB) condition: 8mL 2 V L ; 9vN 2 V N ; vN a0 such that Z mL vN d‘XbkvN k1;S kmL k1=2;AB , (3.17) AB
where b40 is a constant independent of N and h. (A3) Also the following bound holds Z lv d‘pCklk1=2;AB kvk1;S 8v 2 V N . (3.18) AB
Then there exist the error bounds, ku u~ N k1;S pC inf ku vk1;S v2V N
lmax ¼ max xa0
Aðx; xÞ ; ðx; xÞ
lmin ¼ min xa0
Aðx; xÞ , ðx; xÞ
(4.2)
where Aðx; xÞ ¼ xT Ax and ðx; xÞ ¼ xT x. In this subsection, we will provide a upper bound of the effective condition number kbk , lmin kxk
(4.3)
as well as the traditional condition number
ð3:19Þ
Z2V L
Cond ¼
inf ku vk1;S
v2V N
where the matrix A 2 RðNþ1ÞðNþ1Þ is symmetric and positive definite. Denote by lmax and lmin the maximal and the minimal eigenvalues of A, respectively, defined by
Cond_eff ¼
þ inf kl Zk1=2;AB , kl lL k1=2;AB pC
(4.1)
Ax ¼ b,
lmax . lmin
(4.4)
In (4.2), the notations are given from (2.12),
þ inf kl Zk1=2;AB , Z2V L
ð3:20Þ
where C is a constant independent of N and h. Corollary 3.1. Let the conditions of Theorem 3.1 hold. Then there exist the bounds pffiffiffiffiffi ku u~ N k1;S pCf N aN þ bLþ1 g, pffiffiffiffiffi ð3:21Þ kl lL k1;S pCf N aN þ bLþ1 g, where 0oa; bo1. Corollary 3.1 implies an optimal matching between N and L obtained: pffiffiffiffiffi bLþ1 ¼ Oð N aN Þ, which leads approximately to ln a L N. ln b The direct TM involving integration approximation may follow Theorem 2.2. To close this section, let us summarize the important results: the exponential convergence rates can be obtained from Theorems 2.1 and 3.1, but only the polynomial convergence rates for the two TMs involving integration approximation.
ðx; xÞ ¼ kxk2 ¼
L X
D2i ,
(4.5)
i¼0
ZZ I ðvÞ ¼ Aðx; xÞ ¼ rv rv ds S Z Z vn v d‘ þ w2 v2 d‘ 2 AB AB Z 2 vn v d‘ þ w2 kvk20;AB , ¼ jvj1;S 2
ð4:6Þ
AB
where qv=qn ¼ qv=qn is the normal derivative, w2 ¼ Pc ðN þ 1Þ and v 2 V N . We have the following lemma. Lemma 4.1. Suppose that for v 2 V N kvk1;AB pCðN þ 1Þkvk0;AB ,
(4.7)
where C is a constant independent of N. For the penalty plus hybrid TM (2.12) for Motz’s problem, there exists the lower bound, lmin ¼ lmin ðAÞXc0 ,
(4.8)
where c0 40 is also a constant independent of N. Proof. From Li et al. [9,12], when the Pc is chosen large enough, there exists the uniformly inequality c0 kvk2h pI ðvÞ 8v 2 V N ,
(4.9)
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
where c0 ð40Þ is a constant independent of N, and the norm kvkh is defined in (2.14). Hence we have I ðvÞXc0 fjvj21;S þ w2 kvk0;AB gXc0 jvj21;S .
(4.10)
Denote the semi-disk with the radius r, S r ¼ fðr; yÞj0prpr; 0pyppg.
(4.11)
Since Sr jr¼1 S (see Fig. 1), we have jvj1;S Xjvj1;Sr jr¼1 .
(4.12)
From the Green formula, we have ZZ Z jvj21;Sr ¼ rv rv ds ¼ vr v d‘, Sr
(4.13)
‘r
where vr ¼ qv=qr, and ‘r ¼ fðr; yÞjr ¼ r; 0pyppg is the semi-circle. By calculus and the orthogonality of cosði þ 12Þy, we obtain from (2.11) ) Z Z p (X N 1 i1=2 1 vr v d‘ ¼ Di i þ r cos i þ y 2 2 0 ‘r i¼0 ( ) N X 1 Di riþ1=2 cos i þ y r dy 2 i¼0 Z N X 1 2iþ1 p 1 2 2 ¼ Di i þ r cos i þ y dy 2 2 0 i¼0 N pX 1 D2i i þ r2iþ1 , ð4:14Þ ¼ 2 i¼0 2 where Di are arbitrary coefficients. When r ¼ 1, Eq. (4.14) leads to Z N N pX 1 pX vr v d‘ ¼ D2i i þ X D2 . (4.15) 2 i¼0 2 4 i¼0 i ‘r Combining (4.10), (4.12) and (4.15) we give I ðvÞXc0
N X
D2i .
xa0
This completes the proof of Lemma 4.1.
i¼0
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2ðNþ1Þ 1 max pC 1 Nðrmax ÞN , ¼ CN 2 rmax 1
ð4:21Þ
where C 1 is a constant independent of N. This completes the proof of Lemma 4.2. & Since the vector x has the lower bound for Motz’s solution, vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N uX kxk ¼ t D2i XjD0 jX400, (4.22) i¼0
where Di are the approximate expansion coefficients, we obtain the following theorem from Lemmas 4.1 and 4.2, (4.3) and (4.22) immediately.
Cond_eff ¼
&
where w2 ¼ Pc ðN þ 1Þ. We have the following lemma. Lemma 4.2. Let w2 ¼ Pc ðN þ 1Þ. For vector b, there exists a upper bound,
where rmax of N.
i¼0
(4.17)
Then Denote Fi ¼ Fi ðr; yÞ ¼ riþ1=2 cosði þ 12Þy. PN vN ¼ i¼0 Di Fi , and the components bi of b are given from (2.12) by Z Z qFi d‘ þ 1000w2 Fi d‘, (4.18) bi ¼ 500 AB qx AB
kbkpCNðrmax ÞN ,
where we have used w2 ¼ Pc ðN þ 1Þ. Hence we obtain vffiffiffiffiffiffiffiffiffiffiffiffiffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u N uX uX 2 t kbk ¼ bi pCN t r2i max
Theorem 4.1. Let (4.7) hold. For the penalty plus hybrid TM (2.12) for Motz’s problem, there exists the bound for the effective condition number
Hence we have from (4.2) I ðvÞ Xc0 . ðx; xÞ
Proof. Since qFi =qx ¼ ði þ 12Þri1=2 cosði 12Þy, we have from (4.18) Z Z qFi d‘ þ 1000w2 jbi jp500 jFi j d‘ AB qx AB Z 1 i1=2 1 iþ r p500 cos i 2 y d‘ 2 AB Z 1 þ 1000w2 riþ1=2 j cos i þ yj d‘ 2 AB 1 i1=2 þ 1000w2 riþ1=2 p500 i þ r 2 1 pC iþ ð4:20Þ þ N rimax , 2
(4.16)
i¼0
lmin ðAÞ ¼ min
169
(4.19) pffiffiffi ¼ maxS r ¼ 2, and C is a constant independent
where rmax
kbk pCNðrmax ÞN , (4.23) lmin kxk pffiffiffi ¼ 2 and C is a constant independent of N.
Next, we will estimate a upper bound of lmax . We have the following lemma. Lemma 4.3. Let (4.7) hold. For the penalty plus hybrid TM (2.12) for Motz’s problem, when w2 ¼ Pc ðN þ 1Þ there exists a upper bound, lmax ¼ lmax ðAÞpCN 2 ðrmax Þ2N , (4.24) pffiffiffi where rmax ¼ 2 and C is a constant independent of N. Proof. From (4.6) Z 2 I ðvÞpjvj1;S þ 2 vn v d‘ þ w2 kvk20;AB , AB
(4.25)
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
170
where w2 ¼ Pc ðN þ 1Þ. We have Z vn v d‘pkvn k21=2;AB kvk21=2;AB ,
defined by (4.26)
AB
where the semi-norms and the negative norms in the Sobolev space are defined by (3.6) and (3.7). Since Dv ¼ 0 for v 2 V h , we have from Babusˇ ka and Aziz [24] kvn k1=2;AB pCkvk1;S ,
(4.27)
and from the imbedding theorem [20] kvk1=2;AB pCkvk1;S .
(4.28)
Combining (4.26)–(4.28) gives Z vn v d‘pCkvk21;S .
g¼ (4.29)
AB
I ðvÞpCð1 þ w2 Þkvk21;S .
(4.30)
Moreover, since vj1oxo0^y¼0 ¼ 0 for v 2 V N , we have from the Poincare´ inequality [20], v 2 VN.
(4.31)
Denote S rmax ¼ Sjr¼rmax defined in (4.11). We have S Srmax . Then we obtain from (4.30), (4.31) and (4.14) I ðvÞpCð1 þ w2 Þjvj21;S pCð1 þ w2 Þjvj21;Srmax pCð1 þ w
2
ÞNr2N max
N X
D2i .
Cond jlmax jkxk ¼ . Cond_eff kbk
ð4:32Þ
i¼0
The desired result (4.24) follows from (4.2) by noting w2 ¼ Pc ðN þ 1Þ, and this completes the proof of Lemma 4.3. &
First, we have the following lemma.
jlmax jpCfNðrmax Þ2N þ Lg,
Theorem 4.2. Let (4.7) hold. For the penalty plus hybrid TM (2.12) for Motz’s problem, when w2 ¼ Pc ðN þ 1Þ, there exists the bound for the traditional condition number CondpCN 2 ðrmax Þ2N , (4.33) pffiffiffi where rmax ¼ 2, and C is a constant independent of N. From Theorems 4.1 and 4.2, we can see that pffiffiffi Cond_eff ¼ OðNð 2ÞN Þ; Cond ¼ OðN 2 2N Þ, ð4:34Þ
The Cond_eff is significantly smaller than Cond. 4.2. The direct TM
Proof. From (3.9) Z ZZ 1 2 Fðx; xÞ ¼ jrvj ds lv d‘. 2 S AB
where the matrix F 2 RðNþ1ÞðNþ1Þ is symmetric but nondefinite. Denote jlmax j and jlmin j the maximal and the minimal eigenvalues in absolute values of F, respectively,
(4.40)
Since Z pklk l v d‘ 0;AB kvk0;AB AB
1 p ðklk20;AB þ kvk20;AB ÞpCðklk20;AB þ kvk21;S Þ, 2 ð4:41Þ where we have used kvk0;AB pCkvk1;S . Hence we obtain from (4.40) and (4.31) jFðx; xÞjpCðklk20;AB þ kvk21;S Þ pC 1 ðklk20;AB þ jvj21;S Þ pC 1 ðklk20;AB þ jvj21;Srmax Þ,
ð4:42Þ
where C 1 is a constant independentP of N and L. Next, from (3.11) we have l ¼ Li¼0 Ai T i ðxÞ with the coefficients Ai , where T i ðxÞ is the Chebyshev polynomials (3.12), Ai are coefficients, and x ¼ 1 2y. Since jT i ðxÞjp1, we have from the Schwarz inequality !2 Z L X 2 klk0;AB ¼ Ai T i ðxÞ d‘
Eq. (3.10) leads to the linear algebraic equations, (4.35)
(4.39)
Then we have Z ZZ 1 jrvj2 ds þ lv d‘. jFðx; xÞjp 2 S AB
AB
Fx ¼ b,
(4.38)
where C is a constant independent of N and L.
Based on Lemmas 4.1 and 4.2, we have the following theorem from (4.4).
Cond ¼ OððCond_effÞ2 Þ.
(4.37)
Lemma 4.4. For the direct TM (3.10) for Motz’s problem, there exists a upper bound,
Hence we have from (4.25) and kvk0;AB pCkvk1;S ,
kvk1;S pCjvj1;S
jFðx; xÞj jFðx; xÞj ; jlmin j ¼ min , (4.36) xa0 ðx; xÞ xa0 ðx; xÞ P PL 2 2 where ðx; xÞ ¼ N i¼0 Di þ i¼0 Ai . Hence the effective condition number (4.3) and the traditional condition number (4.4) are also valid if replacing lmax and lmin by jlmax j and jlmin j, respectively. From the definitions of Cond and Cond_eff, we have their ratio,
jlmax j ¼ max
p
Z
i¼0 L X
AB
pCL
!2 Ai jT i ðxÞj
i¼0
L X i¼0
A2i ,
d‘pjABj
L X
!2 Ai
i¼0
ð4:43Þ
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
where C is a constant independent of L. Hence from (4.42), (4.43) and (4.14) ( ) N L X X 2N 2 2 jFðx; xÞjpC Nðrmax Þ Di þ CL Ai i¼0
i¼0 N X
L X
!
171
P 2 Moreover, the solution vector kxk ¼ f N i¼0 Di þ PL 2 1=2 Oð1Þ, since the sequences fDi g and fAi g are i¼0 Ai g geometrically convergent, see Table 7 given in the next section. Then we obtain following theorem from (4.37), (4.46) and Lemma 4.4.
ð4:44Þ
Theorem 4.3. For the direct TM (3.10)for Motz’s problem, there exists the ratio
This desired result (4.38) follows from (4.36), and completes the proof of Lemma 4.4. &
Cond jlmax jkxk ¼ pCfNðrmax Þ2N þ Lg, (4.47) Cond_eff kbk pffiffiffi where rmax ¼ 2, and C is a constant independent of N and L.
2N
pCfNðrmax Þ
þ Lg
D2i
þ
i¼0
A2i
.
i¼0
The components of b from (3.10) are given by Z 1 bi ¼ 500 T i ð1 2yÞ dy.
(4.45) 5. Numerical experiments
0
Since T 0 ðxÞ ¼ 1 we have b0 ¼ 500, and vffiffiffiffiffiffiffiffiffiffiffiffiffi u L uX kbk ¼ t b2i Xjb0 jX500.
(4.46)
i¼0
From Lu et al. [25], the Gaussian rule may raise the accuracy of the leading coefficients in approximate solutions. Then we will choose the Gaussian rules of six nodes in this paper. First for Motz’s problem by the penalty plus
Table 1 The error norms, condition numbers and effective condition numbers for Motz’s problem by the penalty plus hybrid Trefftz method N M
10 240
18 240
26 240
34 240
34 1920
jj 0;AB qe qn 1;BC qe qn 1;CD jej0;S jej1;S e0 j jDD
0:143 ð1Þ 0.461
0:120 ð3Þ 0:605 ð2Þ
0:119 ð5Þ 0:835 ð4Þ
0:255 ð6Þ 0:797 ð5Þ
0:122 ð7Þ 0:121 ð5Þ
0.512
0:604 ð2Þ
0:749 ð4Þ
0:139 ð5Þ
0:948 ð6Þ
0:691 ð2Þ 0.179 0:224 ð6Þ
0:229 ð4Þ 0:510 ð2Þ 0:253 ð10Þ
0:191 ð6Þ 0:187 ð3Þ 0:481 ð12Þ
0:742 ð8Þ 0:787 ð5Þ 0:489 ð12Þ
0:149 ð8Þ 0:788 ð5Þ 0:340 ð14Þ
0:268 ð5Þ
0:418 ð9Þ
0:298 ð12Þ
0:342 ð12Þ
0:791 ð13Þ
0:219 ð4Þ 2.11 0.104 (4) 55.9
0:530 ð6Þ 2.11 0.251 (6) 739
0:133 ð9Þ 2.11 0.628 (8) 0.112 (5)
0:336 ð11Þ 2.11 0.159 (11) 0.174 (6)
0:336 ð11Þ 2.11 0.159 (11) 0.174 (6)
jD0 j e1 j jDD jD1 j lmax ðAÞ lmin ðAÞ Cond Cond_eff
Table 2 The error norms and errors of leading coefficients from the penalty plus hybrid TM for Motz’s problem by the Gaussian rule of six nodes rule as N ¼ 34, where M is the number of nodes along AB M
kek1;AB
4D0 D 0
4D1 D 1
4D2 D 2
4D3 D 3
30 60 120 240 480 960 1920 3840 7680 15360
0:155 0:361 0:102 0:255 0:604 0:117 0:122 0:134 0:138 0:137
0:305 0:764 0:191 0:489 0:101 0:419 0:340 0:159 0:262 0:639
0:258 0:634 0:158 0:342 0:178 0:396 0:791 0:182 0:110 0:249
0:157 0:391 0:987 0:246 0:574 0:847 0:593 0:617 0:206 0:542
0:123 0:314 0:786 0:198 0:501 0:135 0:700 0:795 0:692 0:296
ð4Þ ð5Þ ð5Þ ð6Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ
ð10Þ ð11Þ ð11Þ ð12Þ ð12Þ ð13Þ ð14Þ ð13Þ ð13Þ ð13Þ
ð10Þ ð11Þ ð11Þ ð12Þ ð12Þ ð13Þ ð13Þ ð12Þ ð12Þ ð12Þ
ð9Þ ð10Þ ð11Þ ð11Þ ð12Þ ð13Þ ð12Þ ð12Þ ð12Þ ð12Þ
ð8Þ ð9Þ ð10Þ ð10Þ ð11Þ ð11Þ ð13Þ ð13Þ ð12Þ ð12Þ
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
172
Table 3 The computed coefficients by penalty plus hybrid TM for Motz’s problem at N ¼ 34 by the Gaussian rule of six nodes with M ¼ 1920 along AB, where Cond ¼ 0:174ð11Þ and Cond_eff ¼ 0:350ð6Þ i
D~ i
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
401.162453745235780 87.6559201950809808 17.2379150794570393 8.07121525969756881 1.44027271701435344 0.331054885918999064 0.275437344510011672 0:869329945079974908 ð1Þ 0:336048784089699043 ð1Þ 0:153843744676367914 ð1Þ 0:730230165812159951 ð2Þ 0:318411371069231850 ð2Þ 0:122064593191541443 ð2Þ 0:530965354174686320 ð3Þ 0:271512086017575616 ð3Þ 0:120045395993149926 ð3Þ 0:505391338611362311 ð4Þ 0:231664115217506241 ð4Þ 0:115349881701277696 ð4Þ 0:529376122315752069 ð5Þ 0:229002261909030397 ð5Þ 0:106257700723052036 ð5Þ 0:530877129212484370 ð6Þ 0:245424903141403556 ð6Þ 0:108802035839763755 ð6Þ 0:511393128467733109 ð7Þ 0:254875722596865296 ð7Þ 0:111517612609778815 ð7Þ 0:497789230128346209 ð8Þ 0:235669761775463569 ð8Þ 0:117383353377804351 ð8Þ 0:357052180887959222 ð9Þ 0:157652196649668681 ð9Þ 0:752044900590329238 ð10Þ 0:374844370146114913 ð10Þ
hybrid TM, the errors, condition numbers and effective condition numbers are listed in Table 1, where ¼ u uN , and the notations are sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ZZ ZZ 2 jj0;S ¼ ds; jj1;S ¼ jrj2 ds, S
S
q q ¼ max , qn BC qy 1;BC
jj1;AB ¼ max jj; AB q q ¼ max , qn CD qx 1;CD ~ DDi ¼ Di Di ,
where Di and D~ i are the true and the approximate coefficients, respectively. We can see from Table 1, jj0;S ¼ Oðð0:67ÞN Þ;
jj1;S ¼ Oðð0:67ÞN Þ,
(5.1)
which coincide with Theorem 2.1. Next we can see from Table 1, lmax ðAÞ ¼ Oðð1:22ÞN Þ, pffiffiffi Cond ¼ Oð2N Þ; Cond_eff ¼ Oðð 2ÞN Þ. lmin ðAÞ ¼ C;
ð5:2Þ
The last asymptote of Cond_eff in the above equation can be observed, based on the data in Table 1 with the same M ¼ 240, Cond_effjN¼34 0:174ð6Þ ¼ 0:112ð5Þ Cond_effjN¼26
pffiffiffi ¼ 15:5 16 ¼ ð 2Þ8 ¼ ðrmax Þ8 ,
ð5:3Þ
where M is the total number of nodes in AB used in the Gaussian rule. Eqs. (5.2) are consistent with Lemma 4.1 and Theorems 4.1 and 4.2. In order to seek a better D0 , for N ¼ 34 we choose different M. The relative errors of Di ; i ¼ 0; 1; 2; 3 are listed
Table 4 The error norms, condition numbers and effective condition numbers for Motz’s problem by the direct Trefftz method, where the true solution jlj1;AB ¼ 357 N; L M
10; 4 240
18; 5 240
26; 7 240
jl lL j1;AB jl lL j0;AB jj AB qe qn 1;BC qe qn
0:292 0.135 0:166 ð1Þ 0.413
0:103 0:160 ð1Þ 0:163 ð3Þ 0:789 ð2Þ
0:944 0:112 0:158 0:121
0.535
0:582 ð2Þ
0:536 ð4Þ
0:211 ð5Þ
0:541 ð6Þ
jej0;S jej1;S e0 j jDD
0:706 ð2Þ 0.175 0:266 ð6Þ
0:292 ð4Þ 0:493 ð2Þ 0:356 ð10Þ
0:197 ð6Þ 0:183 ð3Þ 0:467 ð12Þ
0:463 ð8Þ 0:785 ð5Þ 0:471 ð12Þ
0:160 ð8Þ 0:777 ð5Þ 0:142 ð15Þ
0:360 ð5Þ
0:561 ð9Þ
0:362 ð12Þ
0:427 ð12Þ
0:551 ð14Þ
0.171 (4) 0:263 ð2Þ 0.650 (6) 375
0.418 0:492 0.850 0.201
0.105 (9) 0:819 ð3Þ 0.128 (13) 0.121 (5)
0.267 0:129 0.207 0.766
0.267 0:129 0.207 0.766
34; 9 240 ð2Þ ð2Þ ð5Þ ð3Þ
0:168 0:581 0:600 0:322
34; 9 7680 ð3Þ ð4Þ ð7Þ ð5Þ
0:444 0:635 0:165 0:181
ð4Þ ð5Þ ð7Þ ð5Þ
1;CD
jD0 j e1 j jDD jD1 j jlmax ðFÞj jlmin ðFÞj Cond Cond_eff
(6) ð3Þ (9) (4)
(11) ð4Þ (16) (5)
(11) ð4Þ (16) (5)
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
in Table 2. From Table 2, we have found the D0 with 15 significant digits at M ¼ 1920, and provide all leading coefficients in Table 3. Next, for Motz’s problem by the direct TM, the errors, condition numbers and effective condition numbers are listed in Tables 4–6. For N ¼ 34 and L ¼ 9 the best D0 is found in Table 6 with M ¼ 7680. For N ¼ 34, L ¼ 9 and M ¼ 7680, the leading coefficients are listed in Table 7. From Table 4, there exist the empirical asymptotes, jl lL j0;AB ¼ Oðð0:69ÞN Þ, jj0;S ¼ Oðð0:63ÞN Þ;
jj1;S ¼ Oðð0:63ÞN Þ.
ð5:4Þ
The empirical exponential convergence rates coincide with Corollary 3.1. Also from Table 4, we can also see lmin ðFÞ ¼ Oðð0:59ÞN Þ; Cond ¼ Oðð2:5ÞN Þ;
lmax ðFÞ ¼ Oðð1:5ÞN Þ, Cond_eff ¼ Oðð1:68ÞN Þ.
ð5:5Þ
173
When M ¼ 240, from Table 4 we obtain the ratios for N ¼ 34 and 26 Cond 0:207ð16Þ ¼ ¼ 0:270ð10Þ, g1 ¼ Cond_eff 0:766ð5Þ Cond 0:128ð13Þ g2 ¼ ¼ ¼ 0:106ð8Þ, ð5:6Þ Cond_eff 0:121ð5Þ respectively, to give their ratio g1 0:270ð10Þ ¼ 255:4 256 ¼ 0:106ð8Þ g2 pffiffiffi ¼ 28 ¼ ð 2Þ16 ¼ ðrmax Þ28 .
ð5:7Þ
Eq. (5.7) indicates the numerical asymptote Cond ¼ Oððrmax Þ2N Þ, Cond_eff
(5.8)
which agrees with Theorem 4.3.
Table 5 The error norms, condition numbers and effective condition numbers for Motz’s problem by the direct Trefftz method with N ¼ 34 and M ¼ 240 L
11
9
jl lL j1;AB jl lL j0;AB jj AB qe qn 1;BC qe qn 1;CD jej0;S jej1;S e0 j jDD
1.88 0.181 0:835 ð8Þ 0:288 ð4Þ
0:168 0:581 0:600 0:322
0:346 ð4Þ
0:211 ð5Þ
0:146 ð3Þ
0:399 ð2Þ
0:927 ð2Þ
0:860 ð7Þ 0:750 ð5Þ 0:471 ð12Þ
0:463 ð8Þ 0:785 ð5Þ 0:471 ð12Þ
0:760 ð6Þ 0:295 ð4Þ 0:445 ð12Þ
0:652 ð4Þ 0:166 ð2Þ 0:141 ð9Þ
0:157 ð2Þ 0:283 ð1Þ 0:403 ð7Þ
0:424 ð12Þ
0:427 ð12Þ
0:460 ð12Þ
0:235 ð8Þ
0:705 ð6Þ
0.267 0:251 0.107 0.996
0.267 0:129 0.207 0.766
0.267 0:130 0.206 0.764
0.267 (11) 0:715 ð2Þ 0.373 (13) 0.138 (3)
0.267 (11) 0:221 ð1Þ 0.121 (13) 44.6
jD0 j e1 j jDD jD1 j jlmax ðFÞj jlmin ðFÞj Cond Cond_eff
(11) ð10Þ (22) (10)
7 ð3Þ ð4Þ ð7Þ ð5Þ
(11) ð4Þ (16) (5)
5
0:587 0:674 0:104 0:231
ð2Þ ð3Þ ð4Þ ð3Þ
(11) ð2Þ (14) (3)
3
0:432 0:612 0:754 0:909
ð1Þ ð2Þ ð3Þ ð2Þ
0:538 0:705 0:171 0:692
ð1Þ ð2Þ ð1Þ ð1Þ
Table 6 The error norms and errors of leading coefficients from the direct TM for Motz’s problem by the Gaussian rule of six nodes rule as N ¼ 34 and L ¼ 9, where M is the number of nodes along AB M
kek1;AB
4D0 D
4D1 D
4D2 D
4D3 D
0:301 0:753 0:188 0:471 0:118 0:290 0:708 0:255 0:142 0:142 0:213 0:694 0:353
0:275 0:687 0:172 0:427 0:108 0:263 0:292 0:324 0:551 0:616 0:616 0:139 0:778
0:131 0:324 0:810 0:202 0:506 0:115 0:340 0:179 0:868 0:989 0:284 0:704 0:233
0:118 0:298 0:746 0:186 0:467 0:115 0:325 0:742 0:464 0:104 0:128 0:350 0:530
0
30 60 120 240 480 960 1920 3840 7680 15360 30720 61440 122880
0:314 0:973 0:242 0:600 0:138 0:158 0:163 0:164 0:165 0:164 0:164 0:164 0:163
ð5Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ
ð10Þ ð11Þ ð11Þ ð12Þ ð12Þ ð13Þ ð14Þ ð14Þ ð15Þ ð14Þ ð14Þ ð14Þ ð14Þ
1
ð10Þ ð11Þ ð11Þ ð12Þ ð12Þ ð13Þ ð14Þ ð15Þ ð14Þ ð14Þ ð14Þ ð13Þ ð14Þ
2
ð9Þ ð10Þ ð11Þ ð11Þ ð12Þ ð12Þ ð13Þ ð13Þ ð14Þ ð14Þ ð13Þ ð13Þ ð13Þ
3
ð10Þ ð10Þ ð10Þ ð10Þ ð11Þ ð11Þ ð12Þ ð13Þ ð13Þ ð12Þ ð12Þ ð13Þ ð13Þ
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
174
Table 7 The computed coefficients by the direct TM for Motz’s problem at N ¼ 34 and L ¼ 9 by the Gaussian rule of six nodes with M ¼ 7680 along AB i
D~ i
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
401.162453745234473 87.6559201950883988 17.2379150794469602 8.07121525969850850 1.44027271702260307 0.331054885920439745 0.275437344511575033 0:869329945182076569 ð1Þ 0:336048784013096499 ð1Þ 0:153843744655217714 ð1Þ 0:730230171201515840 ð2Þ 0:318411377642627540 ð2Þ 0:122064583787855107 ð2Þ 0:530965339810580820 ð3Þ 0:271512403959266711 ð3Þ 0:120045603702773163 ð3Þ 0:505388574434537121 ð4Þ 0:231663774629916405 ð4Þ 0:115357376458830408 ð4Þ 0:529403268205044528 ð5Þ 0:228968049742712711 ð5Þ 0:106253353627942843 ð5Þ 0:531674010537534927 ð6Þ 0:245588935639265154 ð6Þ 0:108599572158342144 ð6Þ 0:511072049985628876 ð7Þ 0:258961666286623336 ð7Þ 0:111973290170299084 ð7Þ 0:492162098344059451 ð8Þ 0:234519394304798108 ð8Þ 0:127222786035216330 ð8Þ 0:361778166440610682 ð9Þ 0:151771358294162104 ð9Þ 0:737166229537951134 ð10Þ 0:463940700433669773 ð10Þ
i
A~ i 0 1 2 3 4 5 6 7 8 9
340:470847534819598 18:1576584805583181 0:490581518606606148 2:19556646178029169 0:152479090907213927 0:638314780462041764 ð1Þ 0:143070661282430310 ð1Þ 0:298153855980607253 ð3Þ 0:609121543712076350 ð3Þ 0:118299765484016262 ð3Þ
pffiffiffi Since rmax ¼ 2, the term Nðrmax Þ2N is much larger than OðLÞ. Hence when the L is changed, lmax and Cond=Cond_eff retain almost the same, based on Lemma 4.4. Such a conclusion is confirmed from the data in Table 5. Interestingly, lmin and Cond_eff will depend on L significantly. When N ¼ 34 and L ¼ 9, we can see from Table 4 that the condition numbers for M ¼ 240 and 7680 are the same. Moreover, when N ¼ 34, L ¼ 9, and M ¼ 240 and 7680, the Cond ¼ 0:207ð16Þ is so huge that the solutions with
Table 8 The number of significant decimal digits of leading coefficients for the three BAMs at N ¼ 34, as well as those in Georgiou et al. [18] at N ¼ 74, and N l ¼ 33, where N l is the number of the Lagrange multiplier used i
Direct TM
Georgiou et al. [18]
P-H TM
Li et al. [9]
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
16 14 14 13 13 12 10 9 9 8 7 7 6 6 6 6 5 4 4
13 12 12 11 11 10 10 9 9 9 8 8 8 7 7 6 5 5 3
15 13 12 13 11 11 11 9 9 9 8 7 7 6 6 5 5 5 5
13 12 12 11 11 10 10 9 10 9 8 7 7 6 6 5 5 5 3
the coefficients in Table 7 are skeptical. In contrast, Cond_eff ¼ 0.766(5) is moderate so that this solution with coefficients in Table 7 is valid and trustworthy. Hence, the effective condition number in this paper is important to the stability analysis of TMs. Finally, let us compare the two TMs for Motz’s problem. The numbers of significant digits of the leading coefficients Di are listed in Table 8. Note that the accuracy of D0 in this paper is higher than that in Georgiou et al. [18] and Li et al. [9]. Also we can see from Table 8, the direct TM and the penalty plus hybrid TM provide 16 and 15 significant digits with M ¼ 7680 and 1920, respectively. Overall, both TMs provide almost the same accuracy and stability, based on Cond_eff but not on Cond, and they are both efficient. However, the penalty plus hybrid TM is more advantageous. The reasons are two-fold: (1) the simplicity in algorithms without extra-variables; (2) wide application. For the direct TM for Poisson’s equation, the constants are excluded into the particular solutions for Laplace’s equation. When the solution domain S is divided by G0 into several subdomains S i , for the direct TM, the interior subdomains are not permitted for Laplace’s equation, see [12]. However, for the penalty plus hybrid TM, such a limitation does not exist. 6. Concluding remarks To close this paper, let us give a few remarks. 1. The effective condition number Cond_eff may provide a better upper bound of relative errors for the given b than the traditional Cond does. The new stability analysis is made based on Cond_eff, but not on Cond. Note that in many cases, the huge Cond is misleading, because the numerical solutions look convincingly accurate. Such
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
conclusions in practical computation have been observed for a long time. The analysis in this paper for the first time provides the theoretical arguments. It is also well known that the stability is a serious issue for the TMs, the method of fundamental solutions and the spectral methods. The study on the effective condition number in this paper may enhance their stability analysis. 2. For penalty plus hybrid TM, frompTheorems 4.1 and ffiffiffi 4.2, the bounds of Cond_eff ¼ OðN 2ÞN and Cond ¼ OðN 2 2N Þ are proved for Motz’s problem, where N is the number of the singular functions used. Hence the Cond_eff is significantly smaller than Cond. Moreover, for the direct TM, the bounds Cond=Cond_eff ¼ OðN2N Þ explicitly display the advantages of Cond_eff over Cond. 3. Based on the error analysis, and in particular the new stability analysis, accompanied with the numerical experiments, both the penalty plus hybrid TM and the direct TM are very efficient. However, the penalty plus hybrid TM is more recommended due to simplicity and less limitation in applications. Acknowledgments We are grateful to S. Christiansen, P.C. Hansen, Alexander H-D Cheng and the referees for their valuable comments and suggestions. References [1] Golub GH, van Loan CF. Matrix computations. 2nd ed. Baltimore and London: The Johns Hopkins; 1989. [2] Higham NJ. Accuracy and stability of numerical algorithms. Philadelphia: SIAM; 1996. [3] Chan FC, Foulser DE. Effectively well-conditioned linear systems. SIAM J Stat Comput 1988;9:963–9. [4] Christiansen S, Hansen PC. The effective condition number applied to error analysis of certain boundary collocation methods. J Comp Appl Math 1994;54(1):15–36. [5] Christiansen S, Saranen J. The conditioning of some numerical methods for first kind boundary integral equations. J Comp Appl Math 1996;67(1):43–58. [6] Banoczi JM, Chiu NC, Cho GE, Ipsen ICP. The lack of influence of the right-hand side on the accuracy of linear system solution. SIAM J Stat Comput 1998;20:203–27. [7] Huang HT, Li ZC. Effective condition number and superconvergence of the Trefftz method coupled with high order FEM for singularity problems. Eng Anal Boundary Elements 2006;30:270–83.
175
[8] Li ZC. Combined methods for elliptic equations with singularities, interfaces and infinities. Boston, London: Kluwer Academic Publishers; 1998. [9] Li ZC, Chen YL, Georgiou GG, Xenohontos C. Special boundary approximation methods for Laplace equation problems with boundary singularities—applications to the Motz problem. Inter Comp Math Appl 2006;51:115–42. [10] Li ZC, Huang HT, J. Huang, New stability analysis of Adini’s elements for Poisson’s equation. Technical Report. Kaohsiung, Taiwan: Department of Applied Mathematics, National Sun Yatsen University; 2005. [11] Li ZC, Huang HT, Huang J. Stability analysis and superconvergence for the Trefftz penalty method coupled with high order FEM for singularity problems. Technical Report. Kaohsiung, Taiwan: Department of Applied Mathematics, National Sun Yat-sen University; 2005. [12] Li ZC, Lu TT, Hu HY, Cheng AHD. Trefftz and collocation methods. Southampton: WIT Press; 2006. [13] Li ZC, Lu TT, Huang HT, Cheng AH-D. Trefftz, collocation and other boundary methods—a comparison. Numer Methods PDEs, accepted for publication. [14] Motz M. The treatment of singularities of partial differential equations by relaxation methods. Quart Appl Math 1947;4:371–7. [15] Rosser JB, Papamichael N. A power series solution of a harmonic mixed boundary value problem. MRC, Technical Report. University of Wisconsin; 1975. [16] Li ZC, Mathon R, Sermer P. Boundary methods for solving elliptic problem with singularities and interfaces. SIAM J Numer Anal 1987;24:487–98. [17] Lucas TR, Oh HS. The method of auxiliary mapping for the finite element solutions of elliptic problems containing singularities. J Comp Phys 1993;1108:327–52. [18] Georgiou GC, Olson LG, Smyrlis YS. A singular function boundary integral method for the Laplace equation. Commun Numer Methods Eng 1996;12:127–34. [19] Atkinson KE. An introduction to numerical analysis. 2nd ed. New York: Wiley; 1989. [20] Ciarlet PG. Basic error estimates for elliptic problems. In: Ciarlet PG, Lions TL, editors. Finite element methods (Part 1). Amsterdam: North-Holland; 1992. p. 17–352. [21] Jin WG, Cheung YK. Trefftz direct method. Adv Eng Software 1995;24:65–9. [22] Kita E, Kamiya N. Trefftz method: an overview. Adv Eng Software 1995;24:3–12. [23] Babusˇ ka I. The finite element method with Lagrangian multipliers. Numer Math 1973;20:179–92. [24] Babusˇ ka I, Kziz AK. Survey lectures on the mathematical foundations on the finite element method. In: Aziz AK, editor. The mathematical foundations of the finite element method with applications to partial differential equations. New York and London: Academic Press; 1972. p. 3–349. [25] Lu TT, Hu HY, Li ZC. Highly accurate solutions of Motz’s and the cracked beam problems. Eng Anal Boundary Elements 2004;28: 1387–403.
Engineering Analysis with Boundary Elements 31 (2007) 163–175 www.elsevier.com/locate/enganabound
Stability analysis for the penalty plus hybrid and the direct Trefftz methods for singularity problems Zi-Cai Lia,b,c,, Hung-Tsai Huangd, Jin Huanga,e, Leevan Linga,f a
Department of Applied Mathematics, National Sun Yat-sen University, Kaohsiung 80424, Taiwan Department of Computer Science and Engineering, National Sun Yat-sen University, Kaohsiung 80424, Taiwan c National Center for Theoretical Sciences, Taiwan d Department of Applied Mathematics, I-Shou University, Kaohsiung County 840, Taiwan e College of Applied Mathematics, University of Electronic and Science Technology of China, ChengDu, China f Department of Mathematics, Hong Kong, Baptist University, Kowloon Tong, Hong Kong
b
Received 20 October 2005; accepted 18 April 2006 Available online 11 July 2006 This paper is dedicated to S. Christiansen on the occasion of his 70th birthday
Abstract For solving the linear algebraic equations Ax ¼ b, the new stability analysis is made based on the effective condition number Cond_eff. The Cond_eff may provide a better upper bound of relative errors of x resulting from the rounding errors of b, than the traditional condition number Cond, which with too large value is, in many times, misleading. In this paper, we apply the effective condition number to the Trefftz methods (TMs) for Poisson’s equations with singularities. Two TMs, such as the penalty plus hybrid TM and the Lagrange multiple (i.e., direct) TM, are studied. We focus on the stability analysis of the solutions when the optimal superconvergence is achieved. Whenpsolving Motz’s problem, the benchmark of singularity problems, by the penalty plus hybrid TM, we have derived that Cond_eff ¼ ffiffiffi OðNð 2ÞN Þ and Cond ¼ OðN 2 2N Þ, where N is the number the singular particular functions used. For solving Motz’s problem by the pffiffiffiof N direct TM, we have derived that Cond=Cond_eff ¼ OðN 2 Þ. Numerical experiments are provided to verify the stability analysis made. In summary, the two TMs are efficient, but the penalty plus hybrid TM is more recommended, due to simplicity of algorithms without extra-variables and less limitations in applications. r 2006 Elsevier Ltd. All rights reserved. Keywords: Effective condition number; Condition number; Stability analysis; Trefftz method; Penalty plus hybrid technique; Lagrange multiplier technique; The direct Trefftz method; Singular problem; Motz’s problem
1. Introduction For solving the linear algebraic equations Ax ¼ b, the traditional condition number is defined by Cond ¼ s1 =sn in [1,2], where s1 and sn are the maximal and the minimal singular values of matrix A, respectively. The condition number is used to provide the bounds of the relative errors from the perturbation of both A and b. Such a Cond can only be reached by the worst situation of all rounding Corresponding author. Department of Applied Mathematics and Department of Computer Science and Engineering, National Sun Yat-sen University, Kaohsiung 80424, Taiwan. E-mail address: [email protected] (Z.-C. Li).
0955-7997/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2006.04.005
errors and all b. For the given b, however, the real relative errors may be smaller, and even much smaller than the Cond. Therefore, the effective condition number was introduced and studied in Chan and Foulser [3]; subsequently it was applied to boundary value problem in Christiansen and Hansen [4], and to the boundary integral equation in Christiansen and Saranen [5]. Other kinds of effective condition number are given in Banoczi et al. [6], where the effective condition number is applied to the Gaussian elimination and the QR factorization. The main references for condition number are given in [7]. Let us consider the linear algebraic equations Ax ¼ b,
(1.1)
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
164
where A 2 Rnn , x 2 Rn and b 2 Rn are the unknown and known vectors, respectively. When there occurs a perturbation of b, the errors of x satisfy (1.2)
Aðx þ DxÞ ¼ b þ Db.
The values of Cond are used to measure the relative errors of x, which for (1.2) are given by kDxk kDbk pCond , kxk kbk
(1.3)
where kxk is the Euclidean norm and the matrix norm kAk ¼ supxa0 kAxk=kxk. Let the matrix A 2 Rnn be symmetric and positive definite. The eigenvalues are arranged in a descending order, l1 Xl2 X Xln 40, and the eigenvectors ui satisfy Aui ¼ li ui , where fui g are orthogonal, with uTi uj ¼ dij , where dij ¼ 1 if i ¼ j and dij ¼ 0 if iaj. The traditional condition number is given by Cond ¼ l1 =ln . In [2,5,7], the effective condition number is defined by Cond_eff ¼
kbk ¼ ln kxk
kbk qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , Pn 2 2 ln b =l i i¼1 i
(1.4)
where bi ¼ uTi b, also to provide kDxk kDbk pCond_eff . kxk kbk
(1.5)
In this paper, we will apply the effective condition number to the Trefftz methods (TMs) for solving Poisson’s equation, in particular for Motz’s problem, the benchmark of singularity problems. When the particular solutions or the piecewise particular solutions satisfying Poisson’s equation are chosen, the approximate solutions are obtained by a linear combination of these particular solutions, and their expansion coefficients are sought by satisfying the exterior or interior boundary conditions. This is called the TM. To couple the boundary conditions, different techniques lead to different TMs. First, we use the penalty plus hybrid techniques in Li [8], to give the penalty plus hybrid TM. Next, we choose the Lagrange multiplier techniques to give the Lagrange multiplier TM (or called the direct TM). In the direct TM, the multipliers are also regarded as variables. For solving Laplace’s equation, the constants are not permitted in the particular solutions in the direct TM. Moreover, if the uniform particular solutions cannot be found on the entire domain, using the piecewise particular solutions is necessary. When the solution domain S is divided into several subdomains S i such that S ¼ [i S i and Si \ S j ¼ ; for iaj, such a limitation without constants excludes the interior subdomains Si with qSi \ qS ¼ ;. However, no such a limitation exists for the penalty plus hybrid TM. In the direct TM, there often occurs the puzzle that the Cond is extremely huge, but the approximate solutions are convincingly accurate. Based on the Cond, such numerical solutions should be distrusted, and ignored. However, one often uses them by refusing the warning of Cond. In this paper, the new effective condition number can clarify this
puzzle well. In fact, the Cond_eff from the direct TM and the penalty plus hybrid TM is reasonably large, to revalidate the numerical solutions with a doubt from huge Cond. Moreover, since the exponential convergence rates are proved in Li et al. [9] for the two TMs without integration approximation, both the direct TM and the penalty plus hybrid TM are very efficient. However, the penalty plus hybrid is more advantageous due to its simplicity and less limitation in application. For numerical methods of partial differential equations, stability analysis is as important as error analysis. The proof for the former is more challenging than the latter. This is why there have appeared numerous reports on error analysis, but not so many reports on stability analysis. The effective condition number is a new trend of stability analysis, to provide a better bound of relative errors from rounding errors. Inspired by the pioneer work of [3–5], we apply the effective condition number for different numerical methods, and find several new discoveries. The related work to this paper are [10] for Adini’s elements and [11] for penalty combinations for singularity problems. Since Adini’s elements involve the derivatives ux and uy , the estimates of the minimal eigenvalue of the stiffness matrix are new and intriguing, the new proofs for bounds of both Cond and Cond_eff are provided in [10]. Since the condition number is huge by the penalty combination of TM coupled with high FEM for singularity problems, such a penalty combination was declined. However, the effective condition number is reasonably large so that the penalty combination is reconciled. The strict analysis of stability and superconvergence is explored in [11]. The a priori estimates of upper bounds of effective condition number may enhance the stability analysis of numerical methods. This paper is organized as follows. In Sections 2 and 3, the penalty plus hybrid TM and the direct TM are described, respectively, accompanied with a brief error analysis for Motz’s problem. In Section 4, new stability analysis is made, and the bounds of Cond_eff, Cond and Cond=Cond_eff are derived in detail. In Section 5, numerical experiments are carried to verify the new stability and the error analysis made. 2. Penalty plus hybrid techniques To match the exterior boundary conditions and different particular solutions in the TM, or matching the TM with other methods, an effective coupling strategy is essential in order to yield optimal convergence rates and good stability of the numerical solutions. The coupling can be accomplished by using additional integrals along the common boundary of different particular solutions. 2.1. General algorithms Below, we provide the penalty plus hybrid couplings for combinations of the TM-FEM. Details are given in
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
[8,9,12,13]. Consider the elliptic equation q qu q qu p p þ cu ¼ f ; ðx; yÞ 2 S, qx qx qy qy
where (2.1)
ðvÞ I ða;bÞ h
with the mixed type of the Dirichlet and Robin boundary conditions u ¼ g1
on GD ,
qu þ qu ¼ g2 qn
(2.2)
on GN ,
(2.3)
where the domain S is bounded (i.e., a polygon) with the exterior boundary, G ¼ GD [ GN , and MeasðGD Þ40 for uniqueness of the solution. The functions p; c; f ; q; g1 and g2 are sufficiently smooth, and c ¼ cðx; yÞX0;
q ¼ qðx; yÞX0,
p ¼ pðx; yÞXp0 40, where p0 is a constant. The solution of problem (2.1)–(2.3) can be equivalently expressed by minimizing a quadratic functional IðvÞ:
1 ¼ 2
165
ZZ
ðprv rv þ cv2 Þ ds ZZ 1 þ ðprv rv þ cv2 Þ ds 2 S2 Z Z 1 2 2 þ qv d‘ þ w ðv g1 Þ2 d‘ 2 GN GD Z Z qv p ðv g1 Þ d‘ þ w2 ðvþ v Þ2 d‘ qn GD G0 Z qvþ qv þb p a ðvþ v Þ d‘ qn qn G0 ZZ Z fv ds g2 v d‘, ð2:6Þ S1
S
GN
where the parameters a and b satisfy a þ b ¼ 1, w2 ¼ Pc ðL þ NÞs is the penalty weight, s40, and Pc ð40Þ is a constant which is large enough but independent of L and N. Eq. (2.6) includes various combined methods, based on different a; b; Pc and s.1 When fFi g and fCi g are chosen as the particular solutions satisfying Eq. (2.1), Eq. (2.5) leads the TMs of variant combinations.
IðuÞ ¼ min IðvÞ, v2H 1 ðSÞ
2.2. Application for Motz’s problem
where the quadratic functional is ZZ 1 ½prv rv þ cv2 ds IðvÞ ¼ 2 S Z Z Z Z 1 2 qv d‘ fv ds g2 v d‘, þ 2 GN S GN
Motz’s problem is the benchmark of singularity problems. The problem was first discussed by Motz [14] in 1947 for the relaxation method. Since then many researchers have selected Motz’s problem as a prototype for verifying efficiency of numerical methods (see [8]). Motz’s problem solves the Laplace equation on the rectangle S ¼ fðx; yÞ j 1oxo1; 0oyo1g
and H 1 ðSÞ ¼ fvjv; vx ; vy 2 L2 ðSÞ; and vjGD ¼ g1 g is a subset of the Sobolev space. Let S be divided by a piecewise straight line G0 into S 1 and S2 . The admissible functions can be written as follows: 8 L P > > bi Fi in S1 ; >
v¼
i¼1
> > þ > :v ¼
N P
(2.4)
ai C i
in S 2 ;
i¼1
where fFi g and fCi g are analytic, complete and linearly independent basis functions on S1 and S2 , respectively, and ai and bi are unknown coefficients to be sought. Note that the admissible functions vþ and v in (2.4) may not satisfy the elliptic equation (2.1) exactly. Define a space H ¼ fvjv 2 L2 ðSÞ; v 2 H 1 ðS 1 Þ; v 2 1 H ðS 2 Þg. Let V h H be a finite-dimensional collection of the functions (2.4). In this section, we present the combination of methods using the penalty plus hybrid techniques: to seek an approximate solution uh 2 V h such that I ða;bÞ ðuh Þ h
¼
min I ða;bÞ ðvÞ, h v2V h
(2.5)
Du ¼
q2 u q2 u þ ¼0 qx2 qy2
in S,
(2.7)
with the mixed Neumann–Dirichlet boundary conditions (see Fig. 1): ujOD ¼ 0; ujAB ¼ 500, qu qu qu ¼ ¼ ¼ 0. qyOA qyBC qxCD
ð2:8Þ ð2:9Þ
Note that there exists a singularity at the origin (0, 0) due to the intersection point of the Neumann–Dirichlet boundary conditions. In fact, the singular solutions of (2.7)–(2.9) are found as 1 X 1 uðr; yÞ ¼ Di riþ1=2 cos i þ y, (2.10) 2 i¼0 where Di are the true expansion coefficients, and ðr; yÞ are the polar coordinates with the origin at ð0; 0Þ (see Fig. 1). Since the convergence radius, R ¼ 2, is analyzed in Rosser 1
In some case (e.g., Pc 40, and a ¼ 1 and b ¼ 0), Eq. (2.5) may not hold, but the variational equation of I hð1;0Þ ðvÞ leads to the desired equation of uh .
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
166
of S ¼ SR ¼ fðr; yÞ j 0proR; 0pyppg, it is better to choose w2 ¼ OðN þ 1Þ, i.e., w2 ¼ Pc ðN þ 1Þ, where Pc is a constant independent of N, see [12]. Define the norm
Y un = 0
C
B
kvkh ¼ fjvj21;S þ w2 kvk0;AB g1=2 .
(2.14)
We have the following theorem. un = 0
u = 500
X D
u=0
O
un = 0
A
Theorem 2.1. Suppose that there exists a constant C independent of N such that qv pCðN þ 1Þkvk0;AB 8v 2 V N . (2.15) qn 0;AB Then there exists the error bound for uPH N
Fig. 1. Motz’s problem.
ku uPH N kh pC inf ku vkh . v2V N
and Papamichael [15], the series expansions (2.10) are most suited to the entire solution domain S. Hence the admissible functions of finite terms, N X 1 uN ðr; yÞ ¼ D~ i rðiþ1=2Þ cos i þ y, (2.11) 2 i¼0 with the unknown coefficients D~ i , are most efficient for numerical Motz’s solutions, to yield the exponential convergence rates OðecN Þ, where c is a positive constant. When functions (2.11) are chosen, Eq. (2.7), ujOD ¼ 0 and qu=qyjOA ¼ 0 are satisfied automatically. Then the coefficients D~ i are sought by the collocation equations of the rest of boundary conditions in (2.8) and (2.9), which is called the TM in this paper. The very accurate solution was obtained in double precision, with the 35 leading coefficients, D~ 0 D~ 34 provided in [16], while an error at the power of D~ 31 was discovered by Lucas and Oh [17]. The very accurate solution can be regarded as a true solution for testing other numerical methods. In Georgiou et al. [18], the same singular functions as in (2.11) are chosen, but the boundary conditions are matched by Lagrange multipliers method to be discussed in Section 3. In this section, we will provide the error bounds for the penalty plus hybrid TM and the direct TM. Based on (2.5), the numerical solutions uPH by the N penalty plus hybrid TM may be sought by minimizing the following energy: TðuPH N Þ ¼ min TðvÞ, v2V N
(2.12)
where V N is the finite-dimensional collection of (2.11) and Z Z 1 qv qv v d‘ ðv 500Þ d‘ TðvÞ ¼ 2 GN qn qn AB Z 2 þw ðv 500Þ2 d‘, ð2:13Þ AB
where GN ¼ OA [ BC [ CD (see Fig. 1) and w40. Below, let us consider how to choose the weight w. First, w2 should be chosen to balance the first and the second terms on the right-hand side in (2.13). From the simple case
(2.16)
Moreover, suppose that w2 ¼ Pc ðN þ 1Þ, and that Pc is chosen large enough but independent of N. There exists the error bound for uPH N ku uPH N kh pC inf ku vkh v2V N pffiffiffiffiffi ¼ C 1 N aN ; 0oao1,
ð2:17Þ
where C 1 is also bounded and independent of N. Theorem 2.1 displays the exponential convergence rates by the penalty plus hybrid TM for Motz’s problem. Let us consider the integration approximation in (2.12). The penalty plus hybrid TM involving integral approximation is designed to seek u~ PH N 2 V N such that ~ ~ u~ PH Þ ¼ min TðvÞ, Tð N v2V N
(2.18)
where Ze Z 1 e qv qv ~ v d‘ ðv 500Þ d‘ TðvÞ ¼ 2 GN qn qn AB Ze þ w2 ðv 500Þ2 d‘,
ð2:19Þ
AB
R R where e is an approximation of by some rules. We have the following theorem. Theorem 2.2. For the solutions by the penalty plus hybrid TM, there exists the error bound, ku u~ PH N k1;S pC inf ku vk1;S v2V N ( Z Ze ! 1 qu w d‘ þ C sup qn w2V N kwk1;S GN GN ) Z ! Z e qu þ ð2:20Þ w d‘ . AB qn AB The proofs of Theorems 3.1 and 3.2 are given in [9,12]. Note that the additional errors from the integration approximation in (2.20) are analogous to those in FEM. Since the integration rules are formulated based on k-order polynomials, such as the Gaussian and Newton–Cotes
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
rules in [19], only the polynomial convergence rates can be obtained from (2.20). 3. The Lagrange multiplier techniques 3.1. General algorithms R We introduce an additional integral G0 lðvþ v Þ d‘, where l is a continuous function of Lagrange multipliers. We describe the Lagrange multiplier coupling as: to seek ðuh ; lÞ 2 V N R such that Bðuh ; v; l; mÞ ¼ f ðvÞ;
ðv; mÞ 2 V N R,
167
The variational equation, qF ðvÞ=qv ¼ 0, leads to (3.5). In Li [8], this is referred to the Lagrange multiplier method, which was just the direct TM in Jin and Cheung [21], and Kita and Kamiya [22]. The Lagrange multiplier method was first introduced by Babusˇ ka [23] to treat the constraint Dirichlet boundary condition as a natural boundary condition, and to relax the limitation on the admissible functions used. Since then the techniques of Lagrange multipliers have drawn much attention, the related references are given in [8]. 3.2. Application for Motz’s problem
(3.1)
where ðv; mÞ 2 V N R denotes that v 2 V N and m 2 R, and ZZ ðp 5 u 5v þ cuvÞ ds Bðu; v; l; mÞ ¼ S1 ZZ þ ðp 5 u 5v þ cuvÞ ds
For Motz’s problem, we may regard the Dirichlet conditionRRujAB ¼ 500 as a constraint to minimize the energy, 12 S jrvj2 ds. Then define a functional ZZ Z 1 jrvj2 ds lðv 500Þ d‘, (3.9) IðvÞ ¼ 2 S AB
S2
þ Dðu; v; l; mÞ,
ð3:2Þ
ZZ fv ds.
f ðvÞ ¼
(3.3)
S
In (3.2), the boundary integrals are Z Z quv d‘ lv d‘ Dðu; v; l; mÞ ¼ GN GD \S 2 Z Z mðu gÞ d‘ lðvþ v Þ d‘ GD \S 2 G0 Z mðuþ u Þ d‘, ð3:4Þ G0
G
G
where l ¼ qu=qn is also treated as unknown, and l; m 2 H 1=2 ðGÞ. H 1=2 ðGÞ is the negative norm in the Sobolev space defined by (see [20]) kvk1=2;G ¼ kvk20;G Z Z
2
ðvðPÞ vðQÞÞ d‘ðPÞ d‘ðQÞ kP Qk2 G G R j G uv d‘j ¼ sup . va0 kvk1=2;G þ
kuk1=2;G
S
AB
AB
8ðv; mÞ 2 H 10 ðSÞ H 1=2 ðABÞ,
ð3:10Þ
where H 10 ðSÞ ¼ fv; vx ; vy 2 L2 ðSÞ and vjOD ¼ 0g. Let lL be the L-order polynomials on AB. We choose the Chebyshev polynomials, lL ¼
L X
A~ i T i ð1 2yÞ;
0pyp1,
(3.11)
i¼0
where the Lagrange multiplier l has the true solution, l ¼ p qu=qnjG0 . In (3.2), the l is treated as an extra variable. Consider the Laplace equation with the Dirichlet condition u ¼ g on G, and no interior boundary exists. From (3.1) with p ¼ 1 and c ¼ 0 we have ZZ Z Z 5 u 5v ds mðu gÞ d‘ lv d‘ ¼ 0, (3.5) S
with the Lagrange multiplier l. From (3.9), the stationary condition for (3.9) is given by: to seek ðu; lÞ 2 H 10 ðSÞ H 1=2 ðABÞ such that ZZ Z Z ru rv ds lv d‘ mðu 500Þ d‘ ¼ 0,
In fact, we may define a functional ZZ Z 1 F ðvÞ ¼ 5 v 5v ds mðv gÞ d‘. 2 S G
1=2
where A~ i are the coefficients to be sought, and the Chebyshev polynomials are defined by T i ðxÞ ¼ cosði arccosðxÞÞ;
1pxp1.
(3.12)
Denote by V N V L the collection of finite dimensions of (2.11) and (3.11). The discrete Lagrange multiplier method (i.e., the direct TM) is given by: to seek ðu~ N ; lL Þ 2 V N V L such that Aðu~ N ; vÞ þ bðu~ N ; v; lL ; mÞ ¼ 0
8ðv; mÞ 2 V N V L ,
(3.13)
where ZZ Aðu; vÞ ¼
ru rv ds, S
,
ð3:6Þ
Z bðu; v; l; mÞ ¼
(3.7)
(3.8)
Z lv d‘
AB
mðu 500Þ d‘.
(3.14)
AB
The multiplier method for Motz’s problem was discussed in Georgiou et al. [18], where the piecewise k-order polynomials were chosen, to provide the polynomial convergence rates. Below we give a brief justification for Laplace’s equation.
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
168
From [8] we obtain the following theorem.
4. Stability analysis for Motz’s problem
Theorem 3.1. Suppose the following three assumptions hold.
4.1. The penalty plus hybrid TM
(A1) For Aðu; vÞ, there exist two positive constants C 0 and C independent of N such that,
For simplicity, we consider the TMs for Motz’s problem in Sections 2.2 and 3.2. The approaches for stability analysis can be easily extended to other Poisson’s equations by the TMs. Eq. (2.12) leads to the linear algebraic equations,
C 0 kvk21;S pAðv; vÞ 8v 2 V N , jAðu; vÞjpCkuk1;S kvk1;S
ð3:15Þ
8v 2 V N .
ð3:16Þ
R (A2) For AB mv, there exits the Ladyzhenskaya–Babuska–Brezzi (LBB) condition: 8mL 2 V L ; 9vN 2 V N ; vN a0 such that Z mL vN d‘XbkvN k1;S kmL k1=2;AB , (3.17) AB
where b40 is a constant independent of N and h. (A3) Also the following bound holds Z lv d‘pCklk1=2;AB kvk1;S 8v 2 V N . (3.18) AB
Then there exist the error bounds, ku u~ N k1;S pC inf ku vk1;S v2V N
lmax ¼ max xa0
Aðx; xÞ ; ðx; xÞ
lmin ¼ min xa0
Aðx; xÞ , ðx; xÞ
(4.2)
where Aðx; xÞ ¼ xT Ax and ðx; xÞ ¼ xT x. In this subsection, we will provide a upper bound of the effective condition number kbk , lmin kxk
(4.3)
as well as the traditional condition number
ð3:19Þ
Z2V L
Cond ¼
inf ku vk1;S
v2V N
where the matrix A 2 RðNþ1ÞðNþ1Þ is symmetric and positive definite. Denote by lmax and lmin the maximal and the minimal eigenvalues of A, respectively, defined by
Cond_eff ¼
þ inf kl Zk1=2;AB , kl lL k1=2;AB pC
(4.1)
Ax ¼ b,
lmax . lmin
(4.4)
In (4.2), the notations are given from (2.12),
þ inf kl Zk1=2;AB , Z2V L
ð3:20Þ
where C is a constant independent of N and h. Corollary 3.1. Let the conditions of Theorem 3.1 hold. Then there exist the bounds pffiffiffiffiffi ku u~ N k1;S pCf N aN þ bLþ1 g, pffiffiffiffiffi ð3:21Þ kl lL k1;S pCf N aN þ bLþ1 g, where 0oa; bo1. Corollary 3.1 implies an optimal matching between N and L obtained: pffiffiffiffiffi bLþ1 ¼ Oð N aN Þ, which leads approximately to ln a L N. ln b The direct TM involving integration approximation may follow Theorem 2.2. To close this section, let us summarize the important results: the exponential convergence rates can be obtained from Theorems 2.1 and 3.1, but only the polynomial convergence rates for the two TMs involving integration approximation.
ðx; xÞ ¼ kxk2 ¼
L X
D2i ,
(4.5)
i¼0
ZZ I ðvÞ ¼ Aðx; xÞ ¼ rv rv ds S Z Z vn v d‘ þ w2 v2 d‘ 2 AB AB Z 2 vn v d‘ þ w2 kvk20;AB , ¼ jvj1;S 2
ð4:6Þ
AB
where qv=qn ¼ qv=qn is the normal derivative, w2 ¼ Pc ðN þ 1Þ and v 2 V N . We have the following lemma. Lemma 4.1. Suppose that for v 2 V N kvk1;AB pCðN þ 1Þkvk0;AB ,
(4.7)
where C is a constant independent of N. For the penalty plus hybrid TM (2.12) for Motz’s problem, there exists the lower bound, lmin ¼ lmin ðAÞXc0 ,
(4.8)
where c0 40 is also a constant independent of N. Proof. From Li et al. [9,12], when the Pc is chosen large enough, there exists the uniformly inequality c0 kvk2h pI ðvÞ 8v 2 V N ,
(4.9)
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
where c0 ð40Þ is a constant independent of N, and the norm kvkh is defined in (2.14). Hence we have I ðvÞXc0 fjvj21;S þ w2 kvk0;AB gXc0 jvj21;S .
(4.10)
Denote the semi-disk with the radius r, S r ¼ fðr; yÞj0prpr; 0pyppg.
(4.11)
Since Sr jr¼1 S (see Fig. 1), we have jvj1;S Xjvj1;Sr jr¼1 .
(4.12)
From the Green formula, we have ZZ Z jvj21;Sr ¼ rv rv ds ¼ vr v d‘, Sr
(4.13)
‘r
where vr ¼ qv=qr, and ‘r ¼ fðr; yÞjr ¼ r; 0pyppg is the semi-circle. By calculus and the orthogonality of cosði þ 12Þy, we obtain from (2.11) ) Z Z p (X N 1 i1=2 1 vr v d‘ ¼ Di i þ r cos i þ y 2 2 0 ‘r i¼0 ( ) N X 1 Di riþ1=2 cos i þ y r dy 2 i¼0 Z N X 1 2iþ1 p 1 2 2 ¼ Di i þ r cos i þ y dy 2 2 0 i¼0 N pX 1 D2i i þ r2iþ1 , ð4:14Þ ¼ 2 i¼0 2 where Di are arbitrary coefficients. When r ¼ 1, Eq. (4.14) leads to Z N N pX 1 pX vr v d‘ ¼ D2i i þ X D2 . (4.15) 2 i¼0 2 4 i¼0 i ‘r Combining (4.10), (4.12) and (4.15) we give I ðvÞXc0
N X
D2i .
xa0
This completes the proof of Lemma 4.1.
i¼0
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2ðNþ1Þ 1 max pC 1 Nðrmax ÞN , ¼ CN 2 rmax 1
ð4:21Þ
where C 1 is a constant independent of N. This completes the proof of Lemma 4.2. & Since the vector x has the lower bound for Motz’s solution, vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N uX kxk ¼ t D2i XjD0 jX400, (4.22) i¼0
where Di are the approximate expansion coefficients, we obtain the following theorem from Lemmas 4.1 and 4.2, (4.3) and (4.22) immediately.
Cond_eff ¼
&
where w2 ¼ Pc ðN þ 1Þ. We have the following lemma. Lemma 4.2. Let w2 ¼ Pc ðN þ 1Þ. For vector b, there exists a upper bound,
where rmax of N.
i¼0
(4.17)
Then Denote Fi ¼ Fi ðr; yÞ ¼ riþ1=2 cosði þ 12Þy. PN vN ¼ i¼0 Di Fi , and the components bi of b are given from (2.12) by Z Z qFi d‘ þ 1000w2 Fi d‘, (4.18) bi ¼ 500 AB qx AB
kbkpCNðrmax ÞN ,
where we have used w2 ¼ Pc ðN þ 1Þ. Hence we obtain vffiffiffiffiffiffiffiffiffiffiffiffiffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u N uX uX 2 t kbk ¼ bi pCN t r2i max
Theorem 4.1. Let (4.7) hold. For the penalty plus hybrid TM (2.12) for Motz’s problem, there exists the bound for the effective condition number
Hence we have from (4.2) I ðvÞ Xc0 . ðx; xÞ
Proof. Since qFi =qx ¼ ði þ 12Þri1=2 cosði 12Þy, we have from (4.18) Z Z qFi d‘ þ 1000w2 jbi jp500 jFi j d‘ AB qx AB Z 1 i1=2 1 iþ r p500 cos i 2 y d‘ 2 AB Z 1 þ 1000w2 riþ1=2 j cos i þ yj d‘ 2 AB 1 i1=2 þ 1000w2 riþ1=2 p500 i þ r 2 1 pC iþ ð4:20Þ þ N rimax , 2
(4.16)
i¼0
lmin ðAÞ ¼ min
169
(4.19) pffiffiffi ¼ maxS r ¼ 2, and C is a constant independent
where rmax
kbk pCNðrmax ÞN , (4.23) lmin kxk pffiffiffi ¼ 2 and C is a constant independent of N.
Next, we will estimate a upper bound of lmax . We have the following lemma. Lemma 4.3. Let (4.7) hold. For the penalty plus hybrid TM (2.12) for Motz’s problem, when w2 ¼ Pc ðN þ 1Þ there exists a upper bound, lmax ¼ lmax ðAÞpCN 2 ðrmax Þ2N , (4.24) pffiffiffi where rmax ¼ 2 and C is a constant independent of N. Proof. From (4.6) Z 2 I ðvÞpjvj1;S þ 2 vn v d‘ þ w2 kvk20;AB , AB
(4.25)
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
170
where w2 ¼ Pc ðN þ 1Þ. We have Z vn v d‘pkvn k21=2;AB kvk21=2;AB ,
defined by (4.26)
AB
where the semi-norms and the negative norms in the Sobolev space are defined by (3.6) and (3.7). Since Dv ¼ 0 for v 2 V h , we have from Babusˇ ka and Aziz [24] kvn k1=2;AB pCkvk1;S ,
(4.27)
and from the imbedding theorem [20] kvk1=2;AB pCkvk1;S .
(4.28)
Combining (4.26)–(4.28) gives Z vn v d‘pCkvk21;S .
g¼ (4.29)
AB
I ðvÞpCð1 þ w2 Þkvk21;S .
(4.30)
Moreover, since vj1oxo0^y¼0 ¼ 0 for v 2 V N , we have from the Poincare´ inequality [20], v 2 VN.
(4.31)
Denote S rmax ¼ Sjr¼rmax defined in (4.11). We have S Srmax . Then we obtain from (4.30), (4.31) and (4.14) I ðvÞpCð1 þ w2 Þjvj21;S pCð1 þ w2 Þjvj21;Srmax pCð1 þ w
2
ÞNr2N max
N X
D2i .
Cond jlmax jkxk ¼ . Cond_eff kbk
ð4:32Þ
i¼0
The desired result (4.24) follows from (4.2) by noting w2 ¼ Pc ðN þ 1Þ, and this completes the proof of Lemma 4.3. &
First, we have the following lemma.
jlmax jpCfNðrmax Þ2N þ Lg,
Theorem 4.2. Let (4.7) hold. For the penalty plus hybrid TM (2.12) for Motz’s problem, when w2 ¼ Pc ðN þ 1Þ, there exists the bound for the traditional condition number CondpCN 2 ðrmax Þ2N , (4.33) pffiffiffi where rmax ¼ 2, and C is a constant independent of N. From Theorems 4.1 and 4.2, we can see that pffiffiffi Cond_eff ¼ OðNð 2ÞN Þ; Cond ¼ OðN 2 2N Þ, ð4:34Þ
The Cond_eff is significantly smaller than Cond. 4.2. The direct TM
Proof. From (3.9) Z ZZ 1 2 Fðx; xÞ ¼ jrvj ds lv d‘. 2 S AB
where the matrix F 2 RðNþ1ÞðNþ1Þ is symmetric but nondefinite. Denote jlmax j and jlmin j the maximal and the minimal eigenvalues in absolute values of F, respectively,
(4.40)
Since Z pklk l v d‘ 0;AB kvk0;AB AB
1 p ðklk20;AB þ kvk20;AB ÞpCðklk20;AB þ kvk21;S Þ, 2 ð4:41Þ where we have used kvk0;AB pCkvk1;S . Hence we obtain from (4.40) and (4.31) jFðx; xÞjpCðklk20;AB þ kvk21;S Þ pC 1 ðklk20;AB þ jvj21;S Þ pC 1 ðklk20;AB þ jvj21;Srmax Þ,
ð4:42Þ
where C 1 is a constant independentP of N and L. Next, from (3.11) we have l ¼ Li¼0 Ai T i ðxÞ with the coefficients Ai , where T i ðxÞ is the Chebyshev polynomials (3.12), Ai are coefficients, and x ¼ 1 2y. Since jT i ðxÞjp1, we have from the Schwarz inequality !2 Z L X 2 klk0;AB ¼ Ai T i ðxÞ d‘
Eq. (3.10) leads to the linear algebraic equations, (4.35)
(4.39)
Then we have Z ZZ 1 jrvj2 ds þ lv d‘. jFðx; xÞjp 2 S AB
AB
Fx ¼ b,
(4.38)
where C is a constant independent of N and L.
Based on Lemmas 4.1 and 4.2, we have the following theorem from (4.4).
Cond ¼ OððCond_effÞ2 Þ.
(4.37)
Lemma 4.4. For the direct TM (3.10) for Motz’s problem, there exists a upper bound,
Hence we have from (4.25) and kvk0;AB pCkvk1;S ,
kvk1;S pCjvj1;S
jFðx; xÞj jFðx; xÞj ; jlmin j ¼ min , (4.36) xa0 ðx; xÞ xa0 ðx; xÞ P PL 2 2 where ðx; xÞ ¼ N i¼0 Di þ i¼0 Ai . Hence the effective condition number (4.3) and the traditional condition number (4.4) are also valid if replacing lmax and lmin by jlmax j and jlmin j, respectively. From the definitions of Cond and Cond_eff, we have their ratio,
jlmax j ¼ max
p
Z
i¼0 L X
AB
pCL
!2 Ai jT i ðxÞj
i¼0
L X i¼0
A2i ,
d‘pjABj
L X
!2 Ai
i¼0
ð4:43Þ
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
where C is a constant independent of L. Hence from (4.42), (4.43) and (4.14) ( ) N L X X 2N 2 2 jFðx; xÞjpC Nðrmax Þ Di þ CL Ai i¼0
i¼0 N X
L X
!
171
P 2 Moreover, the solution vector kxk ¼ f N i¼0 Di þ PL 2 1=2 Oð1Þ, since the sequences fDi g and fAi g are i¼0 Ai g geometrically convergent, see Table 7 given in the next section. Then we obtain following theorem from (4.37), (4.46) and Lemma 4.4.
ð4:44Þ
Theorem 4.3. For the direct TM (3.10)for Motz’s problem, there exists the ratio
This desired result (4.38) follows from (4.36), and completes the proof of Lemma 4.4. &
Cond jlmax jkxk ¼ pCfNðrmax Þ2N þ Lg, (4.47) Cond_eff kbk pffiffiffi where rmax ¼ 2, and C is a constant independent of N and L.
2N
pCfNðrmax Þ
þ Lg
D2i
þ
i¼0
A2i
.
i¼0
The components of b from (3.10) are given by Z 1 bi ¼ 500 T i ð1 2yÞ dy.
(4.45) 5. Numerical experiments
0
Since T 0 ðxÞ ¼ 1 we have b0 ¼ 500, and vffiffiffiffiffiffiffiffiffiffiffiffiffi u L uX kbk ¼ t b2i Xjb0 jX500.
(4.46)
i¼0
From Lu et al. [25], the Gaussian rule may raise the accuracy of the leading coefficients in approximate solutions. Then we will choose the Gaussian rules of six nodes in this paper. First for Motz’s problem by the penalty plus
Table 1 The error norms, condition numbers and effective condition numbers for Motz’s problem by the penalty plus hybrid Trefftz method N M
10 240
18 240
26 240
34 240
34 1920
jj 0;AB qe qn 1;BC qe qn 1;CD jej0;S jej1;S e0 j jDD
0:143 ð1Þ 0.461
0:120 ð3Þ 0:605 ð2Þ
0:119 ð5Þ 0:835 ð4Þ
0:255 ð6Þ 0:797 ð5Þ
0:122 ð7Þ 0:121 ð5Þ
0.512
0:604 ð2Þ
0:749 ð4Þ
0:139 ð5Þ
0:948 ð6Þ
0:691 ð2Þ 0.179 0:224 ð6Þ
0:229 ð4Þ 0:510 ð2Þ 0:253 ð10Þ
0:191 ð6Þ 0:187 ð3Þ 0:481 ð12Þ
0:742 ð8Þ 0:787 ð5Þ 0:489 ð12Þ
0:149 ð8Þ 0:788 ð5Þ 0:340 ð14Þ
0:268 ð5Þ
0:418 ð9Þ
0:298 ð12Þ
0:342 ð12Þ
0:791 ð13Þ
0:219 ð4Þ 2.11 0.104 (4) 55.9
0:530 ð6Þ 2.11 0.251 (6) 739
0:133 ð9Þ 2.11 0.628 (8) 0.112 (5)
0:336 ð11Þ 2.11 0.159 (11) 0.174 (6)
0:336 ð11Þ 2.11 0.159 (11) 0.174 (6)
jD0 j e1 j jDD jD1 j lmax ðAÞ lmin ðAÞ Cond Cond_eff
Table 2 The error norms and errors of leading coefficients from the penalty plus hybrid TM for Motz’s problem by the Gaussian rule of six nodes rule as N ¼ 34, where M is the number of nodes along AB M
kek1;AB
4D0 D 0
4D1 D 1
4D2 D 2
4D3 D 3
30 60 120 240 480 960 1920 3840 7680 15360
0:155 0:361 0:102 0:255 0:604 0:117 0:122 0:134 0:138 0:137
0:305 0:764 0:191 0:489 0:101 0:419 0:340 0:159 0:262 0:639
0:258 0:634 0:158 0:342 0:178 0:396 0:791 0:182 0:110 0:249
0:157 0:391 0:987 0:246 0:574 0:847 0:593 0:617 0:206 0:542
0:123 0:314 0:786 0:198 0:501 0:135 0:700 0:795 0:692 0:296
ð4Þ ð5Þ ð5Þ ð6Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ
ð10Þ ð11Þ ð11Þ ð12Þ ð12Þ ð13Þ ð14Þ ð13Þ ð13Þ ð13Þ
ð10Þ ð11Þ ð11Þ ð12Þ ð12Þ ð13Þ ð13Þ ð12Þ ð12Þ ð12Þ
ð9Þ ð10Þ ð11Þ ð11Þ ð12Þ ð13Þ ð12Þ ð12Þ ð12Þ ð12Þ
ð8Þ ð9Þ ð10Þ ð10Þ ð11Þ ð11Þ ð13Þ ð13Þ ð12Þ ð12Þ
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
172
Table 3 The computed coefficients by penalty plus hybrid TM for Motz’s problem at N ¼ 34 by the Gaussian rule of six nodes with M ¼ 1920 along AB, where Cond ¼ 0:174ð11Þ and Cond_eff ¼ 0:350ð6Þ i
D~ i
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
401.162453745235780 87.6559201950809808 17.2379150794570393 8.07121525969756881 1.44027271701435344 0.331054885918999064 0.275437344510011672 0:869329945079974908 ð1Þ 0:336048784089699043 ð1Þ 0:153843744676367914 ð1Þ 0:730230165812159951 ð2Þ 0:318411371069231850 ð2Þ 0:122064593191541443 ð2Þ 0:530965354174686320 ð3Þ 0:271512086017575616 ð3Þ 0:120045395993149926 ð3Þ 0:505391338611362311 ð4Þ 0:231664115217506241 ð4Þ 0:115349881701277696 ð4Þ 0:529376122315752069 ð5Þ 0:229002261909030397 ð5Þ 0:106257700723052036 ð5Þ 0:530877129212484370 ð6Þ 0:245424903141403556 ð6Þ 0:108802035839763755 ð6Þ 0:511393128467733109 ð7Þ 0:254875722596865296 ð7Þ 0:111517612609778815 ð7Þ 0:497789230128346209 ð8Þ 0:235669761775463569 ð8Þ 0:117383353377804351 ð8Þ 0:357052180887959222 ð9Þ 0:157652196649668681 ð9Þ 0:752044900590329238 ð10Þ 0:374844370146114913 ð10Þ
hybrid TM, the errors, condition numbers and effective condition numbers are listed in Table 1, where ¼ u uN , and the notations are sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ZZ ZZ 2 jj0;S ¼ ds; jj1;S ¼ jrj2 ds, S
S
q q ¼ max , qn BC qy 1;BC
jj1;AB ¼ max jj; AB q q ¼ max , qn CD qx 1;CD ~ DDi ¼ Di Di ,
where Di and D~ i are the true and the approximate coefficients, respectively. We can see from Table 1, jj0;S ¼ Oðð0:67ÞN Þ;
jj1;S ¼ Oðð0:67ÞN Þ,
(5.1)
which coincide with Theorem 2.1. Next we can see from Table 1, lmax ðAÞ ¼ Oðð1:22ÞN Þ, pffiffiffi Cond ¼ Oð2N Þ; Cond_eff ¼ Oðð 2ÞN Þ. lmin ðAÞ ¼ C;
ð5:2Þ
The last asymptote of Cond_eff in the above equation can be observed, based on the data in Table 1 with the same M ¼ 240, Cond_effjN¼34 0:174ð6Þ ¼ 0:112ð5Þ Cond_effjN¼26
pffiffiffi ¼ 15:5 16 ¼ ð 2Þ8 ¼ ðrmax Þ8 ,
ð5:3Þ
where M is the total number of nodes in AB used in the Gaussian rule. Eqs. (5.2) are consistent with Lemma 4.1 and Theorems 4.1 and 4.2. In order to seek a better D0 , for N ¼ 34 we choose different M. The relative errors of Di ; i ¼ 0; 1; 2; 3 are listed
Table 4 The error norms, condition numbers and effective condition numbers for Motz’s problem by the direct Trefftz method, where the true solution jlj1;AB ¼ 357 N; L M
10; 4 240
18; 5 240
26; 7 240
jl lL j1;AB jl lL j0;AB jj AB qe qn 1;BC qe qn
0:292 0.135 0:166 ð1Þ 0.413
0:103 0:160 ð1Þ 0:163 ð3Þ 0:789 ð2Þ
0:944 0:112 0:158 0:121
0.535
0:582 ð2Þ
0:536 ð4Þ
0:211 ð5Þ
0:541 ð6Þ
jej0;S jej1;S e0 j jDD
0:706 ð2Þ 0.175 0:266 ð6Þ
0:292 ð4Þ 0:493 ð2Þ 0:356 ð10Þ
0:197 ð6Þ 0:183 ð3Þ 0:467 ð12Þ
0:463 ð8Þ 0:785 ð5Þ 0:471 ð12Þ
0:160 ð8Þ 0:777 ð5Þ 0:142 ð15Þ
0:360 ð5Þ
0:561 ð9Þ
0:362 ð12Þ
0:427 ð12Þ
0:551 ð14Þ
0.171 (4) 0:263 ð2Þ 0.650 (6) 375
0.418 0:492 0.850 0.201
0.105 (9) 0:819 ð3Þ 0.128 (13) 0.121 (5)
0.267 0:129 0.207 0.766
0.267 0:129 0.207 0.766
34; 9 240 ð2Þ ð2Þ ð5Þ ð3Þ
0:168 0:581 0:600 0:322
34; 9 7680 ð3Þ ð4Þ ð7Þ ð5Þ
0:444 0:635 0:165 0:181
ð4Þ ð5Þ ð7Þ ð5Þ
1;CD
jD0 j e1 j jDD jD1 j jlmax ðFÞj jlmin ðFÞj Cond Cond_eff
(6) ð3Þ (9) (4)
(11) ð4Þ (16) (5)
(11) ð4Þ (16) (5)
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
in Table 2. From Table 2, we have found the D0 with 15 significant digits at M ¼ 1920, and provide all leading coefficients in Table 3. Next, for Motz’s problem by the direct TM, the errors, condition numbers and effective condition numbers are listed in Tables 4–6. For N ¼ 34 and L ¼ 9 the best D0 is found in Table 6 with M ¼ 7680. For N ¼ 34, L ¼ 9 and M ¼ 7680, the leading coefficients are listed in Table 7. From Table 4, there exist the empirical asymptotes, jl lL j0;AB ¼ Oðð0:69ÞN Þ, jj0;S ¼ Oðð0:63ÞN Þ;
jj1;S ¼ Oðð0:63ÞN Þ.
ð5:4Þ
The empirical exponential convergence rates coincide with Corollary 3.1. Also from Table 4, we can also see lmin ðFÞ ¼ Oðð0:59ÞN Þ; Cond ¼ Oðð2:5ÞN Þ;
lmax ðFÞ ¼ Oðð1:5ÞN Þ, Cond_eff ¼ Oðð1:68ÞN Þ.
ð5:5Þ
173
When M ¼ 240, from Table 4 we obtain the ratios for N ¼ 34 and 26 Cond 0:207ð16Þ ¼ ¼ 0:270ð10Þ, g1 ¼ Cond_eff 0:766ð5Þ Cond 0:128ð13Þ g2 ¼ ¼ ¼ 0:106ð8Þ, ð5:6Þ Cond_eff 0:121ð5Þ respectively, to give their ratio g1 0:270ð10Þ ¼ 255:4 256 ¼ 0:106ð8Þ g2 pffiffiffi ¼ 28 ¼ ð 2Þ16 ¼ ðrmax Þ28 .
ð5:7Þ
Eq. (5.7) indicates the numerical asymptote Cond ¼ Oððrmax Þ2N Þ, Cond_eff
(5.8)
which agrees with Theorem 4.3.
Table 5 The error norms, condition numbers and effective condition numbers for Motz’s problem by the direct Trefftz method with N ¼ 34 and M ¼ 240 L
11
9
jl lL j1;AB jl lL j0;AB jj AB qe qn 1;BC qe qn 1;CD jej0;S jej1;S e0 j jDD
1.88 0.181 0:835 ð8Þ 0:288 ð4Þ
0:168 0:581 0:600 0:322
0:346 ð4Þ
0:211 ð5Þ
0:146 ð3Þ
0:399 ð2Þ
0:927 ð2Þ
0:860 ð7Þ 0:750 ð5Þ 0:471 ð12Þ
0:463 ð8Þ 0:785 ð5Þ 0:471 ð12Þ
0:760 ð6Þ 0:295 ð4Þ 0:445 ð12Þ
0:652 ð4Þ 0:166 ð2Þ 0:141 ð9Þ
0:157 ð2Þ 0:283 ð1Þ 0:403 ð7Þ
0:424 ð12Þ
0:427 ð12Þ
0:460 ð12Þ
0:235 ð8Þ
0:705 ð6Þ
0.267 0:251 0.107 0.996
0.267 0:129 0.207 0.766
0.267 0:130 0.206 0.764
0.267 (11) 0:715 ð2Þ 0.373 (13) 0.138 (3)
0.267 (11) 0:221 ð1Þ 0.121 (13) 44.6
jD0 j e1 j jDD jD1 j jlmax ðFÞj jlmin ðFÞj Cond Cond_eff
(11) ð10Þ (22) (10)
7 ð3Þ ð4Þ ð7Þ ð5Þ
(11) ð4Þ (16) (5)
5
0:587 0:674 0:104 0:231
ð2Þ ð3Þ ð4Þ ð3Þ
(11) ð2Þ (14) (3)
3
0:432 0:612 0:754 0:909
ð1Þ ð2Þ ð3Þ ð2Þ
0:538 0:705 0:171 0:692
ð1Þ ð2Þ ð1Þ ð1Þ
Table 6 The error norms and errors of leading coefficients from the direct TM for Motz’s problem by the Gaussian rule of six nodes rule as N ¼ 34 and L ¼ 9, where M is the number of nodes along AB M
kek1;AB
4D0 D
4D1 D
4D2 D
4D3 D
0:301 0:753 0:188 0:471 0:118 0:290 0:708 0:255 0:142 0:142 0:213 0:694 0:353
0:275 0:687 0:172 0:427 0:108 0:263 0:292 0:324 0:551 0:616 0:616 0:139 0:778
0:131 0:324 0:810 0:202 0:506 0:115 0:340 0:179 0:868 0:989 0:284 0:704 0:233
0:118 0:298 0:746 0:186 0:467 0:115 0:325 0:742 0:464 0:104 0:128 0:350 0:530
0
30 60 120 240 480 960 1920 3840 7680 15360 30720 61440 122880
0:314 0:973 0:242 0:600 0:138 0:158 0:163 0:164 0:165 0:164 0:164 0:164 0:163
ð5Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ ð7Þ
ð10Þ ð11Þ ð11Þ ð12Þ ð12Þ ð13Þ ð14Þ ð14Þ ð15Þ ð14Þ ð14Þ ð14Þ ð14Þ
1
ð10Þ ð11Þ ð11Þ ð12Þ ð12Þ ð13Þ ð14Þ ð15Þ ð14Þ ð14Þ ð14Þ ð13Þ ð14Þ
2
ð9Þ ð10Þ ð11Þ ð11Þ ð12Þ ð12Þ ð13Þ ð13Þ ð14Þ ð14Þ ð13Þ ð13Þ ð13Þ
3
ð10Þ ð10Þ ð10Þ ð10Þ ð11Þ ð11Þ ð12Þ ð13Þ ð13Þ ð12Þ ð12Þ ð13Þ ð13Þ
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
174
Table 7 The computed coefficients by the direct TM for Motz’s problem at N ¼ 34 and L ¼ 9 by the Gaussian rule of six nodes with M ¼ 7680 along AB i
D~ i
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
401.162453745234473 87.6559201950883988 17.2379150794469602 8.07121525969850850 1.44027271702260307 0.331054885920439745 0.275437344511575033 0:869329945182076569 ð1Þ 0:336048784013096499 ð1Þ 0:153843744655217714 ð1Þ 0:730230171201515840 ð2Þ 0:318411377642627540 ð2Þ 0:122064583787855107 ð2Þ 0:530965339810580820 ð3Þ 0:271512403959266711 ð3Þ 0:120045603702773163 ð3Þ 0:505388574434537121 ð4Þ 0:231663774629916405 ð4Þ 0:115357376458830408 ð4Þ 0:529403268205044528 ð5Þ 0:228968049742712711 ð5Þ 0:106253353627942843 ð5Þ 0:531674010537534927 ð6Þ 0:245588935639265154 ð6Þ 0:108599572158342144 ð6Þ 0:511072049985628876 ð7Þ 0:258961666286623336 ð7Þ 0:111973290170299084 ð7Þ 0:492162098344059451 ð8Þ 0:234519394304798108 ð8Þ 0:127222786035216330 ð8Þ 0:361778166440610682 ð9Þ 0:151771358294162104 ð9Þ 0:737166229537951134 ð10Þ 0:463940700433669773 ð10Þ
i
A~ i 0 1 2 3 4 5 6 7 8 9
340:470847534819598 18:1576584805583181 0:490581518606606148 2:19556646178029169 0:152479090907213927 0:638314780462041764 ð1Þ 0:143070661282430310 ð1Þ 0:298153855980607253 ð3Þ 0:609121543712076350 ð3Þ 0:118299765484016262 ð3Þ
pffiffiffi Since rmax ¼ 2, the term Nðrmax Þ2N is much larger than OðLÞ. Hence when the L is changed, lmax and Cond=Cond_eff retain almost the same, based on Lemma 4.4. Such a conclusion is confirmed from the data in Table 5. Interestingly, lmin and Cond_eff will depend on L significantly. When N ¼ 34 and L ¼ 9, we can see from Table 4 that the condition numbers for M ¼ 240 and 7680 are the same. Moreover, when N ¼ 34, L ¼ 9, and M ¼ 240 and 7680, the Cond ¼ 0:207ð16Þ is so huge that the solutions with
Table 8 The number of significant decimal digits of leading coefficients for the three BAMs at N ¼ 34, as well as those in Georgiou et al. [18] at N ¼ 74, and N l ¼ 33, where N l is the number of the Lagrange multiplier used i
Direct TM
Georgiou et al. [18]
P-H TM
Li et al. [9]
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
16 14 14 13 13 12 10 9 9 8 7 7 6 6 6 6 5 4 4
13 12 12 11 11 10 10 9 9 9 8 8 8 7 7 6 5 5 3
15 13 12 13 11 11 11 9 9 9 8 7 7 6 6 5 5 5 5
13 12 12 11 11 10 10 9 10 9 8 7 7 6 6 5 5 5 3
the coefficients in Table 7 are skeptical. In contrast, Cond_eff ¼ 0.766(5) is moderate so that this solution with coefficients in Table 7 is valid and trustworthy. Hence, the effective condition number in this paper is important to the stability analysis of TMs. Finally, let us compare the two TMs for Motz’s problem. The numbers of significant digits of the leading coefficients Di are listed in Table 8. Note that the accuracy of D0 in this paper is higher than that in Georgiou et al. [18] and Li et al. [9]. Also we can see from Table 8, the direct TM and the penalty plus hybrid TM provide 16 and 15 significant digits with M ¼ 7680 and 1920, respectively. Overall, both TMs provide almost the same accuracy and stability, based on Cond_eff but not on Cond, and they are both efficient. However, the penalty plus hybrid TM is more advantageous. The reasons are two-fold: (1) the simplicity in algorithms without extra-variables; (2) wide application. For the direct TM for Poisson’s equation, the constants are excluded into the particular solutions for Laplace’s equation. When the solution domain S is divided by G0 into several subdomains S i , for the direct TM, the interior subdomains are not permitted for Laplace’s equation, see [12]. However, for the penalty plus hybrid TM, such a limitation does not exist. 6. Concluding remarks To close this paper, let us give a few remarks. 1. The effective condition number Cond_eff may provide a better upper bound of relative errors for the given b than the traditional Cond does. The new stability analysis is made based on Cond_eff, but not on Cond. Note that in many cases, the huge Cond is misleading, because the numerical solutions look convincingly accurate. Such
ARTICLE IN PRESS Z.-C. Li et al. / Engineering Analysis with Boundary Elements 31 (2007) 163–175
conclusions in practical computation have been observed for a long time. The analysis in this paper for the first time provides the theoretical arguments. It is also well known that the stability is a serious issue for the TMs, the method of fundamental solutions and the spectral methods. The study on the effective condition number in this paper may enhance their stability analysis. 2. For penalty plus hybrid TM, frompTheorems 4.1 and ffiffiffi 4.2, the bounds of Cond_eff ¼ OðN 2ÞN and Cond ¼ OðN 2 2N Þ are proved for Motz’s problem, where N is the number of the singular functions used. Hence the Cond_eff is significantly smaller than Cond. Moreover, for the direct TM, the bounds Cond=Cond_eff ¼ OðN2N Þ explicitly display the advantages of Cond_eff over Cond. 3. Based on the error analysis, and in particular the new stability analysis, accompanied with the numerical experiments, both the penalty plus hybrid TM and the direct TM are very efficient. However, the penalty plus hybrid TM is more recommended due to simplicity and less limitation in applications. Acknowledgments We are grateful to S. Christiansen, P.C. Hansen, Alexander H-D Cheng and the referees for their valuable comments and suggestions. References [1] Golub GH, van Loan CF. Matrix computations. 2nd ed. Baltimore and London: The Johns Hopkins; 1989. [2] Higham NJ. Accuracy and stability of numerical algorithms. Philadelphia: SIAM; 1996. [3] Chan FC, Foulser DE. Effectively well-conditioned linear systems. SIAM J Stat Comput 1988;9:963–9. [4] Christiansen S, Hansen PC. The effective condition number applied to error analysis of certain boundary collocation methods. J Comp Appl Math 1994;54(1):15–36. [5] Christiansen S, Saranen J. The conditioning of some numerical methods for first kind boundary integral equations. J Comp Appl Math 1996;67(1):43–58. [6] Banoczi JM, Chiu NC, Cho GE, Ipsen ICP. The lack of influence of the right-hand side on the accuracy of linear system solution. SIAM J Stat Comput 1998;20:203–27. [7] Huang HT, Li ZC. Effective condition number and superconvergence of the Trefftz method coupled with high order FEM for singularity problems. Eng Anal Boundary Elements 2006;30:270–83.
175
[8] Li ZC. Combined methods for elliptic equations with singularities, interfaces and infinities. Boston, London: Kluwer Academic Publishers; 1998. [9] Li ZC, Chen YL, Georgiou GG, Xenohontos C. Special boundary approximation methods for Laplace equation problems with boundary singularities—applications to the Motz problem. Inter Comp Math Appl 2006;51:115–42. [10] Li ZC, Huang HT, J. Huang, New stability analysis of Adini’s elements for Poisson’s equation. Technical Report. Kaohsiung, Taiwan: Department of Applied Mathematics, National Sun Yatsen University; 2005. [11] Li ZC, Huang HT, Huang J. Stability analysis and superconvergence for the Trefftz penalty method coupled with high order FEM for singularity problems. Technical Report. Kaohsiung, Taiwan: Department of Applied Mathematics, National Sun Yat-sen University; 2005. [12] Li ZC, Lu TT, Hu HY, Cheng AHD. Trefftz and collocation methods. Southampton: WIT Press; 2006. [13] Li ZC, Lu TT, Huang HT, Cheng AH-D. Trefftz, collocation and other boundary methods—a comparison. Numer Methods PDEs, accepted for publication. [14] Motz M. The treatment of singularities of partial differential equations by relaxation methods. Quart Appl Math 1947;4:371–7. [15] Rosser JB, Papamichael N. A power series solution of a harmonic mixed boundary value problem. MRC, Technical Report. University of Wisconsin; 1975. [16] Li ZC, Mathon R, Sermer P. Boundary methods for solving elliptic problem with singularities and interfaces. SIAM J Numer Anal 1987;24:487–98. [17] Lucas TR, Oh HS. The method of auxiliary mapping for the finite element solutions of elliptic problems containing singularities. J Comp Phys 1993;1108:327–52. [18] Georgiou GC, Olson LG, Smyrlis YS. A singular function boundary integral method for the Laplace equation. Commun Numer Methods Eng 1996;12:127–34. [19] Atkinson KE. An introduction to numerical analysis. 2nd ed. New York: Wiley; 1989. [20] Ciarlet PG. Basic error estimates for elliptic problems. In: Ciarlet PG, Lions TL, editors. Finite element methods (Part 1). Amsterdam: North-Holland; 1992. p. 17–352. [21] Jin WG, Cheung YK. Trefftz direct method. Adv Eng Software 1995;24:65–9. [22] Kita E, Kamiya N. Trefftz method: an overview. Adv Eng Software 1995;24:3–12. [23] Babusˇ ka I. The finite element method with Lagrangian multipliers. Numer Math 1973;20:179–92. [24] Babusˇ ka I, Kziz AK. Survey lectures on the mathematical foundations on the finite element method. In: Aziz AK, editor. The mathematical foundations of the finite element method with applications to partial differential equations. New York and London: Academic Press; 1972. p. 3–349. [25] Lu TT, Hu HY, Li ZC. Highly accurate solutions of Motz’s and the cracked beam problems. Eng Anal Boundary Elements 2004;28: 1387–403.