- Email: [email protected]
Effective condition number for weighted linear least squares problems and applications to the Trefftz method
Engineering Analysis with Boundary Elements 36 (2012) 53–62
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Effective condition number for weighted linear least squares problems and applications to the Trefftz method Yi-min Wei a,b,1, Tzon-Tzer Lu c, Hung-Tsai Huang d, Zi-Cai Li e,f, a
School of Mathematical Science, Fudan University, PR China Shanghai Key Laboratory of Contemporary Applied Mathematics, Shanghai, 200433, PR China c Department of Applied Mathematics, National Sun Yat-sen University, Kaohsiung, Taiwan d Department of Applied Mathematics, I-Shou University, Kaohsiung, Taiwan e Department of Applied Mathematics, National Sun Yat-sen University, Kaohsiung, Taiwan f Department of Computer Science and Engineering, National Sun Yat-sen University, Kaohsiung, Taiwan b
a r t i c l e i n f o
abstract
Article history: Received 14 March 2011 Accepted 14 July 2011 Available online 23 August 2011
In [27], the effective condition number Cond_eff is developed for the linear least squares problem. In this paper, we extend the effective condition number for weighted linear least squares problem with both full rank and rank-deficient cases. We apply the effective condition number to the collocation Trefftz pffiffiffi Laplace’s equation with a crack singularity, to prove that Cond_eff pffiffiffimethod (CTM) [29] for ¼ Oð LÞ and Cond ¼ OðL1=2 ð 2ÞL Þ, where L is the number of singular pffiffiffi particular solutions used. The Cond grows exponentially as L increases, but Cond_eff is only Oð LÞ. The small effective condition number explains well the high accuracy of the TM solution, but the huge Cond cannot. & 2011 Elsevier Ltd. All rights reserved.
Keywords: Condition number Effective condition number Perturbation Weighted linear least squares problem Collocation Trefftz method Singularity problem
1. Introduction Condition number measures the sensitivity of the solutions of a problem with respect to small perturbations of the given data. In [27], the effective condition number is developed for the linear least squares problem by following Rice [35], Chan and Foulser [11], Christiansen and Hansen [12], Axelsson and Kaporin [2,3], componentwise linear least squares [14,36], and structured Kenny, et al. [23], matrix [44,45]. The traditional condition number was given in Wilkinson [43] in 1965 for linear algebraic equations, and has been used in many textbooks and papers. In Rice [35], the natural condition number was first proposed for specific non-homogenous term; we choose the effective condition number used [11] for our study [26,27]. More references for condition number and effective condition number are provided in [27]. In this paper, we extend the effective condition number to Corresponding author at: Department of Applied Mathematics, National Sun Yat-sen University, Kaohsiung, Taiwan. E-mail addresses: [email protected], [email protected] (Y.-m. Wei), [email protected] (T.-T. Lu), [email protected] (H.-T. Huang), [email protected] (Z.-C. Li). 1 This author is supported by the National Natural Science Foundation of China under Grant 10871051, Doctoral Program of the Ministry of Education under Grant 20090071110003, Shanghai Science and Technology Committee under Grant 09DZ2272900 and 973 Program Project under Grant 2010CB327900.
0955-7997/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2011.07.010
weighted linear least squares problem with both full rank and rank-deficient cases. For the over-determined system Fx ¼b with F A Rmn ðm Z nÞ and rankðFÞ ¼ r r n, and the perturbed system Fðx þ DxÞ ¼ b þ Db, there exist the bounds of relative errors, JDxJ JDbJ rCond , JxJ JbJ
JDxJ JDbJ rCond_eff , JxJ JbJ
ð1:1Þ
where JxJð ¼ JxJ2 Þ is the 2-norm, the traditional condition number and the effective condition number are defined by (see [4,15,26]) Cond ¼
s1 JbJ , , Cond_eff ¼ sr sr JxJ
ð1:2Þ
and sr is the r-th singular value of matrix F. For ðF þ DFÞ ðx þ DxÞ ¼ b þ Db with the nonsingular F ¼ A A Rnn , we also derive in [27] the following computational formula: JDxJ 1 JDAJ JDbJ r þCond_eff , ð1:3Þ Cond JxJ 1d JAJ JbJ where d ¼ JAJ=sn o1. Based on (1.3), we provide in [27] a theoretical analysis that for solving the linear system by Gaussian elimination or the QR factorization, the direction of right hand (i.e., b) is insignificant for solution errors. Such a conclusion was discovered in Banoczi et al. [7] purely by the numerical experiments. We also apply the effective condition number for the finite
54
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
difference method (FDM) for Poisson’s equations, where the discretization errors are dominant, and the direction of b does influence the solution errors. In this paper, we consider the weighted linear least squares problem. Let M A Rmm be symmetric and positive definite matrix. For Fx ¼b, we may solve the weighted linear least squares problem, M1=2 Fx ¼ M1=2 b:
ð1:4Þ
For any symmetric and positive definite matrix N A R may be reduced to
nn
Gy ¼ d,
Hence, the generalized effective condition number in (1.7) can be regarded as one of the effective condition number for weighted linear least squares problem explored in this paper. This paper is organized as follows. In Section 2, the computational formulas are derived for effective condition number of weighted linear least squares problem. In Section 3, the Trefftz method for Laplace’s equation is discussed, and stability analysis is made by the effective condition number, and a linkage to Section 2 is explored. In last section, a few final remarks are made.
, Eq. (1.4) ð1:5Þ
where
2. Effective condition number of weighted linear least squares problem 2.1. Effective condition number
G ¼ M1=2 FN1=2 A Rmn , y ¼ N1=2 xA Rn ,
d ¼ M1=2 b A Rm :
ð1:6Þ
In this paper, we will extend the effective condition number to the weighted linear least squares problem (1.5). Note that Eq. (1.3) in [27] is valid only for A A Rnn . Hence, the error bounds like (1.3) of the weighted linear least squares problem are significant in application. Furthermore, more applications of effective condition number are provided for linear algebraic equations and numerical partial differential equations (PDE), such as the Trefftz method for Poisson’s equation with singularities. For solving Ax¼b with the nonsingular matrix AA Rnn , the effective condition number has been used in Axelsson and Kaporin [2,3] and Banoczi et al. [7]. Denote the norm JxJD ¼ ðxT DxÞ1=2 , where D is the symmetric and positive definite matrix. The generalized effective condition number is defined in Axelsson and Kaporin [2] as JAn xJD n kn ðA,xÞD ¼ JA JD , JxJD
ð1:7Þ
where n a0. For the preconditioned conjugate gradient methods, both the errors and the convergence rates may be estimated, by means of kn ðA,xÞD . In fact, the generalized effective condition number in (1.7) is of the effective condition number for weighted linear least squares problem explored in this paper. For Ax¼ b with the symmetric and positive definite matrix AA Rnn , we may solve the following weighted linear equations with M¼ N¼D, D1=2 Ax ¼ D1=2 b:
ð1:8Þ
The generalized effective condition number of (1.7) is given by
k1 ðA,xÞD ¼
JAxJD 1 JbJD 1 JA JD ¼ JA JD ¼ Cond_effðAÞ: JxJD JxJD
A A Cmn ,
ð2:1Þ
JBJNM ¼ max JByJN ,
B A Cnm ,
ð2:2Þ
JxJN ¼ 1
JyJM ¼ 1
as in [40]. Then we obtain the following relations: JAJMN ¼ JM1=2 AN1=2 J2 ,
A A Cmn ,
ð2:3Þ
JBJNM ¼ JN1=2 BM1=2 J2 ,
B A Cnm :
ð2:4Þ
Lemma 2.1 (van Loan [37]). Let AA Cmn ðm Z nÞ with rankðAÞ ¼ r, and M and N be Hermitian positive definite matrices of order m and n, respectively. There exist U A Cmm and V A Cnn , satisfying UH MU ¼ I and VH N1 V ¼ I, such that D 0 ð2:5Þ A¼U VH , 0 0 and the weighted Moore–Penrose inverse AyMN can be represented by ! D1 0 ð2:6Þ AyMN ¼ N1 V UH M: 0 0 In (2.5) D ¼ diagðm1 , m2 , . . . , mr Þ, m1 Z m2 Z Z mr 4 0 and m2i is the nonzero eigenvalue of N1 AH MA, where mi is called the weighted singular values, and JAJMN ¼ m1 ,
ð1:10Þ
JAyMN JNM ¼
1
mr
:
ð2:8Þ
x
n
Gn x ¼ d ,
n
mn
ð1:11Þ
where matrix A A C and m Z n. Suppose that the rank of A is r with r r n r m. Then the least-squares (M) solution in the minimum-norm (N) is given by
ð1:12Þ
x ¼ AyMN b:
where d ¼ An1 b:
ð2:7Þ
Consider the weighted linear least squares problem with rankdeficient [8,21,42], minJAxbJM ,
and then
Gn ¼ An ,
JAJMN ¼ max JAxJM ,
ð1:9Þ
Moreover, from Ax¼b we have the following weighted linear equations: An x ¼ An1 b,
Consider the weighted norm, and for generality, we use the complex matrices. Given two vectors x A Cn , y A Cm and two Hermitian positive definite matrices, MA Cmm and N A Cnn . The weighted vector norms may be defined by JyJM ¼ JM1=2 yJ2 and JxJN ¼ JN1=2 xJ2 , and the weighted matrix norms by
Hence, the generalized effective condition number (1.9) for (1.11) is given by
Assume that there exists a perturbation of A and b such that rankðAÞ ¼ rankðA þ DAÞ ¼ r,
n
JGn xJD Jd JD k1 ðG ,xÞD ¼ JðGn Þ1 JD ¼ JðGn Þ1 JD JxJD JxJD
we have the perturbed weighted linear least squares problem [6,8,9,16–20,22,31,38,39,41,42],
n
¼ Cond-effðGn Þ ¼
JAn xJD n JA JD ¼ kn ðA,xÞD : JxJD
ð1:13Þ
min JðA þ DAÞðxþ DxÞðb þ DbÞJM :
x þ Dx
ð2:9Þ
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
Then we can find the least-squares (M) solution in the minimumnorm (N) of (2.9) by
55
Remark 2.1. If vector M1=2 b is just parallel to u1 , i.e.,
b2 ¼ ¼ bn ¼ 0, we have JbJM ¼ 9b1 9 and JxJN ¼ 9b1 9=m1 . Hence,
ð2:10Þ
the effective condition number (2.15) just leads to the traditional condition number in (2.11).
Let matrix A be decomposed by the weighted singular value decomposition (2.5) D 0 A¼U VH , 0 0
Remark 2.2. When AA Rmn , and M and N are the identity matrices, Eqs. (2.14) and (2.15) lead to (1.1) and (1.2), respectively.
x þ Dx ¼ ðAþ DAÞyMN ðb þ DbÞ:
where matrices U A Cmm and V A Cnn are weighted unitary, and D A Rrr is a diagonal matrix with the positive weighted singular values mi as m1 Z m2 Z Z mr 4 0: The traditional condition number is defined by van Loan [37] CondMN ðAÞ ¼
JAJMN JAyMN JNM
¼ m1 =mr :
ð2:11Þ
First let us consider (2.9) for simple case of DA ¼ 0. If DAa 0, then the perturbation bounds from effective condition number are given in Theorem 2.1 later. Denote M1=2 U ¼ ½u1 ,u2 , . . . ,um , we have M1=2 b ¼
m X
bi ui , M1=2 Db ¼
i¼1
m X
ai ui ,
ð2:12Þ
First we cite a lemma from [31,40]. Lemma 2.2. Let matrices A, DA A Cmn with rankðAÞ ¼ rankðAþ DAÞ ¼ r r n rm and denote d ¼ JAyMN JNM JDAJMN o1. Then there exist the upper bounds pffiffiffi y 1 þ 5 JAMN JNM JDAJMN JAy JNM , JðAþ DAÞyMN AyMN JNM r 2 1JAyMN JNM JDAJMN MN JðAþ DAÞyMN JNM r
i¼1
1=2 1=2 where bi ¼ uH b and ai ¼ uH Db. It is easy to check that i M i M ! D1 0 x ¼ AyMN b ¼ N1 V UH Mb, 0 0
Dx ¼ AyMN Db ¼ N1 V
D1
0
0
0
! UH MDb:
Since M1=2 U and N1=2 V are unitary, we can obtain ! D1 0 1=2 1=2 1=2 H 1=2 VÞ JxJN ¼ JN xJ2 ¼ ðN ðM UÞ M b 0 0 2 vffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 r r X uX b bi b i ¼ u ¼t Z r, mi i m2i mr i¼1
2.2. Perturbation bounds
1JAyMN JNM JDAJMN
:
We have the following theorem. Theorem 2.1. Let matrices A, DA A Cmn with rankðAÞ ¼ rankðAþ DAÞ ¼ r rn rm and denote d ¼ JAyMN JNM JDAJMN o 1. Then we have " # pffiffiffi JDxJN Cond_eff 1þ 5 JDAJMN JDbJM CondMN ðAÞ r þ : 1d 2 JxJN JAJMN JbJM ð2:16Þ Proof. Since x ¼ AyMN b and x þ Dx ¼ ðAþ DAÞyMN ðb þ DbÞ, we have
Dx ¼ ðA þ DAÞyMN ðb þ DbÞAyMN b ¼ ½ðA þ DAÞyMN AyMN b þ ðAþ DAÞyMN Db,
i¼1
2
JAyMN JNM
and then
and ! D1 0 JDxJN ¼ JN1=2 DxJ2 ¼ ðN1=2 VÞ ðM1=2 UÞH M1=2 Db 0 0 2 vffiffiffiffiffiffiffiffiffiffiffiffiffiffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u X r r r 2 X X u u a ai t 1 i r t ¼ u ¼ a2i rJAyMN JNM JDbJM : mi i m2i mr i¼1
2
i¼1
JDxJN ¼ JN1=2 f½ðAþ DAÞyMN AyMN b þ ðA þ DAÞyMN DbgJ2 ¼ JN1=2 ½ðAþ DAÞyMN AyMN M1=2 M1=2 b þN1=2 ðA þ DAÞyMN M1=2 M1=2 DbJ2 r JðAþ DAÞyMN AyMN JNM JbJM þ JðA þ DAÞyMN JNM JDbJM :
i¼1
ð2:13Þ Hence we have JDx9N JDbJM 1 JDbJM r sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ Cond_eff , 2 JxJN mr JbJM P b r i¼1
ð2:14Þ
i 2 i
þ
m
JbJM JbJM sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2ffi ¼ m JxJN ¼ r b1 br mr þ þ
m1
JAyMN JNM 1JAyMN JNM JDAJMN
JDbJM ,
to give
where the effective condition number is defined by Cond_eff ¼
It follows from Lemma 2.2 pffiffiffi y 1þ 5 JAMN JNM JDAJMN JAy JNM JbJM JDxJN p 2 1JAyMN JNM JDAJMN MN
JAyMN JNM JbJM JxJN
:
mr
ð2:15Þ When br a 0, JDxJN JDbJM mr JbJM JDbJM r ¼ : JxJN mr br br JbJM So the simplest effective condition number can be defined by Cond_EE ¼ JbJM =br from (2.15), see [26]. Clearly, Cond_eff is a tighter bound than Cond_EE.
JDxJN r
pffiffiffi 1þ 5 d JbJM 1 JDbJM þ : 2 1d 1d mr mr
Since d ¼ CondMN ðAÞ JDAJMN =JAJMN , the desired result (2.16) is obtained and the proof of Theorem 2.1 is completed. & If rankðAÞ ¼ rankðA þ DAÞ ¼ n, then the sharp error bound is given by (see [40]) JDxJN Cond_eff pffiffiffi JDAJMN JDbJM r þ 2CondMN ðAÞ : ð2:17Þ 1d JxJN JAJMN JbJM Besides, by following [10], we may also derive similar bounds of Cond_eff, under the same assumptions in this paper.
56
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
2.3. Applications and comparisons Consider the real matrix AA Rmn ðm Z nÞ and the real vectors x A Rn and b A Rm . Let MA Rmm and NA Rnn be the identity matrices. Then we have from Theorem 2.1, " # pffiffiffi JDxJ Cond_eff 1þ 5 JDAJ JDbJ r Cond þ , ð2:18Þ JxJ 1d 2 JAJ JbJ y
where d ¼ JA JJDAJ o1, and JxJ ¼ JxJ2 . Moreover, when rankðAÞ ¼ rankðA þ DAÞ ¼ n, the following bound is given from (2.17), JDxJ Cond_eff pffiffiffi JDAJ JDbJ r þ : ð2:19Þ 2CondðAÞ JxJ 1d JAJ JbJ Let us compare our new error bounds with those in existing literature. Although the bound of (2.19) is larger than that in (1.3), which is only for the nonsingular A A Rnn , Theorem 2.1 is valid for weighted linear least squares problem. In Golub and van Loan [15, p. 228], there is the error bound for A A Rmn ðmZ nÞ JDxJ 2 r e Cond Cond tan y þ þ Oðe2 Þ, ð2:20Þ JxJ cos y where JDAJ JDbJ sn 1 , o ¼ JAJ JbJ s1 Cond
e ¼ max
ð2:21Þ
and sin y ¼
rLS JbJ
,
rLS ¼ JAxLS bJ:
ð2:22Þ
Obviously, the bound in (2.18) is smaller than (2.20), by noting 2 OðCond Þ in (2.20). Moreover, in Lawson and Hanson [24], for mn AA R ðm Z nÞ with rankðA þ DAÞ ¼ rankðAÞ r n, there exists the error bound (see (9.10) in [24, p.51]) JDxJ Cond n rLS o JDAJ JbJ JDbJ r þ 2 þ Cond JxJ 1d JAJ JAJJxJ JbJ JAJJxJ n 1 rLS o JDAJ JDbJ þ Cond_eff , ¼ Cond 2þ Cond 1d JAJ JbJ JAJJxJ ð2:23Þ where d ¼ Cond JDAJ=JAJ and rLS is the same as (2.22). The bound in (2.18) is also smaller than that in (2.23), when Cond is large. Hence (2.16), (2.18) and (2.19) are new contributions of effective condition number in this paper. Note that when AA Rnn , then rLS ¼ 0. Eq. (2.23) leads to JDxJ 1 JDAJ JDbJ r þ Cond_eff , ð2:24Þ 2 Cond JxJ 1d JAJ JbJ which is similar to (1.3). In (2.16), (2.18), and (2.24), if the term JDbJ=JbJ is dominant in the brackets, we have JDxJ Cond_eff JDbJ r , JxJ 1d JbJ
ð2:25Þ
which is very close to the right hand side of (1.2). This is the case occurring in numerical PDE. The same conclusion is drawn as in [27]. Since Eqs. (2.16), (2.18), and (2.24) may lead to (2.25), this section provides a theoretical foundation of effective condition number for weighted linear least squares problems, which is an advance level of [26,27]. The collocation Trefftz method for Laplace’s equation in Section 3 can be classified into weighted linear least squares problems, since different weights are used in the central and the Gaussian rules. The example in Section 3 displays a significance of effective condition number for numerical partial differential equations (PDE). Furthermore, a close linkage of Section 3 and this section is explored in Section 3.3. The techniques of weights are essential in science and numerical
methods, due to variety and complication of the problems studied. Therefore, the effective condition number in this paper is a development of our previous study (see [26,27]).
3. Trefftz method for Laplace’s equations In this section, the effective condition number is applied to the weighted linear least squares problem from numerical PDE. We take the Trefftz method [29] for solving Laplace’s equation with singularities for example. 3.1. Numerical algorithms and computed results In this section, we consider the crack problem (see Fig. 1)
Du ¼ 0 in S, u¼0
on OD,
u ¼ 1 on AB [ BC , un ¼ 0
u ¼ y on CD,
on OA,
ð3:1Þ
where S ¼ fðx,yÞ91 o x o1,0 o yo 1g, and un ¼ @u=@n is the outward normal derivative to @S. To solve (3.1), we may use the collocation Trefftz method (CTM) involving integration approximation (see [29,30,33]). When the Dirichlet boundary conditions involve noise, reader may refer the numerical examples in Lu et al. [32]. The exact solutions can be expanded as 1 X 1 u¼ di r i þ 1=2 cos i þ y in S, ð3:2Þ 2 i¼0 where di are the true coefficients. Choose the admissible solutions as L X 1 Di r i þ 1=2 cos i þ y in S, ð3:3Þ uL ¼ 2 i¼0 where Di are the unknown coefficients to be sought. Since the expansions (3.2) satisfy the Laplace’s equation and the boundary conditions at y ¼0 already, the coefficients Di should be chosen to satisfy the rest of the boundary conditions in (3.1), u9AB[BC ¼ 1,
u9CD ¼ y,
ð3:4Þ
as best as possible, where AB, BC and CD are shown in Fig. 1. Hence, the linear least squares method (LSM) may be used to seek the coefficients Di. Denote Z ½u,v ¼ uv d‘, ð3:5Þ Gn
where Gn ¼ AB [ BC [ CD, and VL denotes the set of finite dimensional functions in (3.3). Then, we may seek uL A VL such that Z Z ½uL ,v ¼ v d‘ þ yv d‘, 8v A VL : ð3:6Þ AB[BC
CD
Fig. 1. The crack problem.
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
Denote the energy Z Z ðv1Þ2 d‘ þ IðvÞ ¼ AB[BC
ðvyÞ2 d‘:
ð3:7Þ
CD
The solution (i.e., the coefficients Di) of (3.6) can also be expressed by: To seek uL A VL such that IðuL Þ ¼ minIðvÞ:
ð3:8Þ
v A VL
When the integrals in (3.7) involve numerical approximation, the energy is denoted by Z Z c c ^ ¼ IðvÞ ðv1Þ2 d‘ þ ðvyÞ2 d‘, ð3:9Þ AB[BC
CD
R R where b is an approximation of integral by some quadrature rule, such as the central or the Gaussian rule. Then the solution u~ L A VL is sought by ^ ^ u~ L Þ ¼ minIðvÞ: Ið
ð3:10Þ
v A VL
On the other hand, one may establish the boundary collocation equations directly from the boundary conditions (3.4), L X 1 i þ 1=2 Di r j cos i þ y ¼ 1, ðrj , yj Þ A AB [ BC , ð3:11Þ 2 j i¼0 L X i¼0
i þ 1=2 Di r j
1 cos i þ y ¼ rj sinðyj Þ, 2 j
57
nodes ðrj , yj Þ of small pffiffiffi subsections. Eqs. (3.11) and (3.12) multiplied by the weight h lead to L pffiffiffi X pffiffiffi 1 i þ 1=2 h Di r j cos iþ yj ¼ h, ðrj , yj Þ A AB [ BC , ð3:13Þ 2 i¼0 L pffiffiffi pffiffiffi X 1 i þ 1=2 Di r j cos iþ yj ¼ hrj sinðyj Þ, h 2 i¼0
ð3:12Þ
where ðr, yÞ are the polar coordinates with the origin (0, 0). Let ðrj , yj Þ denote the integration nodes. For the central rule, we choose uniform sections with the mesh-spacing h and the middle
ð3:14Þ
which are equivalent to (3.10). Such an equivalence is also valid for the Gaussian rule with suitable weights, see [33]. In computation, we choose the number m of nodes ðrj , yj Þ to be larger than (Lþ 1). Eqs. (3.13) and (3.14) lead to Fx ¼ b,
ð3:15Þ
where x A RL þ 1 is the unknown vector consisting of coefficients Di ði ¼ 0,1, . . . ,LÞ, b A Rm , and the stiffness matrix, F A RmðL þ 1Þ . The algorithms of (3.13) and (3.14) are called the collocation Trefftz methods (CTM) in [29,33], or the boundary approximation method in [25,30]. Once coefficients Di (i.e., x) are obtained, the errors on AB [ BC [ CD Z 1=2 Z JeJB ¼ JuuL JB ¼ ð1uL Þ2 d‘ þ ðyuL Þ2 d‘ ð3:16Þ AB[BC
ðrj , yj Þ A CD,
ðrj , yj Þ A CD,
CD
pffiffiffiffiffiffiffiffiffi are computable, where the boundary norm JvJB ¼ ½v,v. We will compute JeJB and the maximal errors on Gn , denoted by JeJ1, Gn ¼ maxGn 9e9. Note that Eq. (3.15) (i.e., (3.13) and (3.14)) is the weighted linear least squares problem of (3.11) and (3.12). We simply use (3.15) for simplicity in discussion.
Table 1 The error norms, condition numbers and errors of leading coefficients from the CTM for the crack problem by the Gaussian rule of six nodes rule as m ¼ 42, where 0n denotes the error less than the computer rounding errors in double precision. L
JeJB
JeJ1, Gn
DD0
D
DD1
D
DD2
D
DD3
D
Cond
Cond_eff
Cond_EE
0.176( 6) 0.908( 9) 0.795( 11) 0.667( 13) 0.566( 15) 0n 0n 0n 0n 0n 0n 0n 0n 0n
0.689( 8) 0.268( 10) 0.243( 12) 0.109( 14) 0.871( 15) 0.871( 15) 0.871( 15) 0.871( 15) 0.871( 15) 0.871( 15) 0.871( 15) 0.871( 15) 0.871( 15) 0.871( 15)
0.546( 6) 0.650( 8) 0.788( 10) 0.970( 12) 0.141( 13) 0.200( 14) 0.206( 14) 0.206( 14) 0.206( 14) 0.206( 14) 0.206( 14) 0.206( 14) 0.206( 14) 0.206( 14)
0.791( 3) 0.161( 4) 0.298( 6) 0.521( 8) 0.870( 10) 0.140( 11) 0.202( 13) 0.155( 13) 0.194( 13) 0.202( 13) 0.202( 13) 0.202( 13) 0.202( 13) 0.202( 13)
15.4 52.5 186 674 0.249(4) 0.927(4) 0.348(5) 0.132(6) 0.502(6) 0.192(7) 0.737(7) 0.286(8) 0.110(9) 0.429(9)
1.22 1.23 1.24 1.24 1.24 1.24 1.24 1.24 1.25 1.25 1.25 1.25 1.25 1.25
1.34 1.41 1.45 1.48 1.50 1.52 1.53 1.54 1.55 1.56 1.56 1.57 1.57 1.58
0
10 14 18 22 26 30 34 38 42 46 50 54 58 62
0.769( 4) 0.549( 5) 0.438( 6) 0.372( 7) 0.330( 8) 0.304( 9) 0.286( 10) 0.274( 11) 0.266( 12) 0.263( 13) 0.265( 14) 0.479( 15) 0.374( 15) 0.379( 15)
0.153( 3) 0.115( 4) 0.967( 6) 0.854( 7) 0.768( 8) 0.743( 9) 0.718( 10) 0.706( 11) 0.704( 12) 0.715( 13) 0.755( 14) 0.178( 14) 0.155( 14) 0.211( 14)
1
2
3
Table 2 The error norms, condition numbers and errors of leading coefficients from the CTM for the crack problem by the Gaussian rule of six nodes rule as L ¼54, where 0n denotes the error less than the computer rounding errors in double precision. m
JeJB
JeJ1, Gn
DD0
D
DD1
D
DD2
D
DD3
D
Cond
Cond_eff
Cond_EE
0.566( 15) 0.754( 15) 0.754( 15) 0.566( 15) 0n 0.377( 15) 0.132( 14) 0.377( 15) 0.189( 15)
0.654( 15) 0.654( 15) 0.436( 15) 0.871( 15) 0.871( 15) 0.654( 15) 0.436( 15) 0.436( 15) 0n
0.117( 14) 0.294( 14) 0.441( 14) 0.117( 14) 0.206( 14) 0.132( 14) 0.176( 14) 0.250( 14) 0.147( 14)
0.168( 13) 0.129( 13) 0.244( 13) 0.129( 14) 0.194( 14) 0.200( 14) 0.129( 15) 0.777( 15) 0.518( 15)
0.286(8) 0.286(8) 0.286(8) 0.286(8) 0.286(8) 0.286(8) 0.286(8) 0.286(8) 0.286(8)
1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25
1.57 1.57 1.57 1.57 1.57 1.57 1.57 1.57 1.57
0
18 24 30 36 42 48 54 60 66
0.676( 15) 0.105( 14) 0.925( 15) 0.750( 15) 0.479( 15) 0.640( 15) 0.152( 14) 0.651( 15) 0.529( 15)
0.102( 13) 0.311( 14) 0.244( 14) 0.244( 14) 0.178( 14) 0.200( 14) 0.333( 14) 0.200( 14) 0.178( 14)
1
2
3
58
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
In our computation, the Gaussian rule with six nodes is chosen for computing the integrals in (3.7), and m denotes the number of nodes along AB. Then m is multiples of six, and m ¼ 4m. The QR factorization is used for solving (3.15). The errors, condition numbers and effective condition numbers are listed in Tables 1 and 2. Note that all computations in this paper are carried out by Fortran programs in double precision. The highly accurate solution is found at L¼ 54 and m ¼ 42, with JeJB ¼ 0:479ð15Þ,
Cond ¼ 0:284ð8Þ,
Cond_eff ¼ 1:25,
DD0
¼ 0n ,
ð3:17Þ
Cond_EE ¼ 1:57,
ð3:18Þ
JeJ1, Gn ¼ 0:178ð14Þ,
D0
Table 3 The leading coefficients Di from the CTM for the crack problem by the Gaussian rule with six nodes as L¼ 54 and m ¼ 42 along AB. i
Di
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 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
1.17769302091071149 0.254824842629162551 0.472464287866036919( 1) 0.268074643613298927( 1) 0.207096003479169876( 2) 0.285158772208820362( 3) 0.183013497194766252( 3) 0.533051704416331353( 3) 0.157452569473614989( 4) 0.537931615260617409( 6) 0.282586859356730478( 5) 0.172189962001964248( 4) 0.991762428046167717( 7) 0.182802186655744245( 7) 0.753824634122721559( 7) 0.681235665190457081( 6) 0.148739422491196200( 8) 0.333015969702268427( 9) 0.248256678147026780( 8) 0.300174031113286268( 7) 0.781012645562570253( 11) 0.104558489837163971( 10) 0.933291671159730157( 10) 0.141328444671448315( 8) 0.311705144580195162( 12) 0.338335131879225414( 12) 0.382332681328341618( 11) 0.696105306905609050( 10) 0.248911941008673783( 14) 0.124843211054098715( 13) 0.166517138300667560( 12) 0.354249175259574404( 11) 0.179980257225143477( 15) 0.490941653852575528( 15) 0.759286687343576104( 14) 0.184733137546747423( 12) 0.851604256940787594( 17) 0.204119729589908422( 16) 0.360262461151616336( 15) 0.977970794533835831( 14) 0.181405679095255387( 17) 0.735513396111446796( 18) 0.181360402136695322( 16) 0.509563344013018680( 15) 0.436358147836334786( 18) 0.196117522501444700( 19) 0.101927109194609567( 17) 0.232534810678914514( 16) 0.623638419504983828( 19) 0.808600297512558928( 20) 0.582074741343905307( 19) 0.654386595624466984( 18) 0.381783351320010671( 20) 0.568217460684078295( 21) 0.207475181771925468( 20)
Table 4 The singular values si and the coefficients bi for the matrix resulting from the CTM solution in Table 3, where Cond¼ 0.286(8), Cond_eff ¼ 1.25 and Cond_EE ¼1.57. i
si
bi
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 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
0.199(8) 0.104(8) 0.101(8) 0.531(7) 0.351(6) 0.246(6) 0.186(6) 0.131(6) 0.178(5) 0.131(5) 0.995(4) 0.735(4) 0.166(4) 0.127(4) 0.983(3) 0.752(3) 0.252(3) 0.198(3) 0.158(3) 0.125(3) 0.567(2) 0.461(2) 0.380(2) 0.311(2) 0.176(2) 0.148(2) 0.126(2) 0.106(2) 0.707(1) 0.608(1) 0.533(1) 0.461(1) 0.348(1) 0.306(1) 0.276(1) 0.244(1) 0.202(1) 0.181(1) 0.168(1) 0.152(1) 0.135(1) 0.123(1) 0.117(1) 0.108(1) 0.101(1) 0.936 0.911 0.858 0.837 0.785 0.775 0.753 0.749 0.715 0.702
0.657( 1) 0.195( 1) 0.163 0.407( 1) 0.380( 1) 0.262( 1) 0.990( 1) 0.554( 1) 0.320( 1) 0.271( 1) 0.864( 1) 0.565( 1) 0.300( 1) 0.278( 1) 0.836( 1) 0.570( 1) 0.302( 1) 0.294( 1) 0.871( 1) 0.588( 1) 0.323( 1) 0.322( 1) 0.962( 1) 0.624( 1) 0.362( 1) 0.364( 1) 0.112 0.676( 1) 0.422( 1) 0.424( 1) 0.135 0.743( 1) 0.514( 1) 0.505( 1) 0.172 0.826( 1) 0.651( 1) 0.606( 1) 0.226 0.907( 1) 0.820( 1) 0.700( 1) 0.292 0.927( 1) 0.931( 1) 0.747( 1) 0.350 0.806( 1) 0.836( 1) 0.730( 1) 0.357 0.495( 1) 0.498( 1) 0.100 0.672
where Cond_EE ¼JbJ=bn , and 0n denotes the error less than the computer rounding errors in double precision. Since the maximal errors of harmonic functions occur only on the boundary Gn , we conclude that JeJ1,S r 0:178ð14Þ. The coefficients Di are listed in Table 3, and the singular values si and the expansion coefficients bi in Table 4. The leading coefficient D0 ¼ 1:17769302091071149 in Table 3 is exact in the sense that its error is less than the computer rounding error t ¼ 0:5ð17Þ, and D1 ,D2 and D3 have 16, 15 and 15 significant digits, respectively.2 Table 5 lists the values of smax , smin , etc. for different L used in Table 2. This very accurate
2 Such significant digits are found by comparisons with the more accurate solution at L¼100 and by Mathematics using 100 working digits.
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
Table 5 The maximal and the minimal singular values and condition numbers in Table 1 as m ¼ 42. L
smax
smin
JbJ
bn
Cond
Cond_eff
Cond_EE
10 14 18 22 26 30 34 38 42 46 50 54 58 62
11.2 37.2 131 475 0.175(4) 0.651(4) 0.245(5) 0.926(5) 0.352(6) 0.135(7) 0.517(7) 0.199(8) 0.768(8) 0.297(9)
0.716 0.710 0.704 0.713 0.704 0.703 0.688 0.702 0.702 0.702 0.702 0.702 0.701 0.701
1.05 1.05 1.05 1.05 1.05 1.05 1.05 1.05 1.05 1.05 1.05 1.05 1.05 1.05
0.782 0.749 0.728 0.713 0.702 0.894 0.688 0.683 0.679 0.676 0.674 0.672 0.670 0.669
15.4 52.5 186 674 0.249(4) 0.927(4) 0.348(5) 0.132(6) 0.502(6) 0.192(7) 0.737(7) 0.284(8) 0.110(9) 0.429(9)
1.22 1.23 1.24 1.24 1.24 1.24 1.24 1.24 1.25 1.25 1.25 1.25 1.25 1.25
1.34 1.41 1.45 1.48 1.50 1.52 1.53 1.54 1.55 1.56 1.56 1.57 1.57 1.58
59
Suppose that the integration rules are chosen so that to satisfy the following equivalence relations: JvJ 0, Gn OðJvJ0, Gn Þ,
ð3:27Þ
i.e., there exist two constants c0 ð 4 0Þ and C1 such that c0 JvJ0, Gn rJvJ 0, Gn r C1 JvJ0, Gn :
ð3:28Þ
An analysis in [33] shows that the integration rules to satisfy (3.28) are not severe, even by the simplest central rule. In the following, the constants c0 ð 4 0Þ,c1 ð 40Þ,C and C1 are independent of L, but their values may be different in different places. We have the following lemma. Lemma 3.1. Suppose that for v A VL there exists a positive constant p 4 0 such that JvJ1, Gn r CLp JvJ0, Gn :
ð3:29Þ
Under (3.27), there exists the lower bound,
smin ðFÞ Z c0 Lp=2 : solution can be used to test other methods, in particular, those for the boundary integral equation (BIE) of first kind.3 3.2. Bounds of effective condition number and condition number Below, we provide the stability analysis for the crack problem in Section 3.1, mainly based on the effective condition number. Eqs. (3.13) and (3.14) are rewritten as (3.15). We will derive an upper bound of (1.2): Cond_eff ¼
smax , smin
Ax ¼ d,
lmin ðAÞ ¼ min xa0
ðAx,xÞ : ðx,xÞ
ð3:23Þ
ð3:24Þ
Z Z c ^ ¼2 v2 d‘ 2IðvÞ ¼ 2 v2 d‘ ¼ 2JvJ20, Gn : ðAx,xÞ ¼ 2IðvÞ Gn
ð3:25Þ
Gn
Denote 2
Z c
ð3:33Þ
Combining (3.32) and (3.33) yields IðvÞ Zc0 Lp JvJ21,S :
ð3:34Þ
Denote the semi-disk with the radius r, Sr ¼ fðr, yÞ90 rr r r,0 r y r pg: Since Sr 9r ¼ 1 S, we have JvJ1,S ZJvJ1,Sr 9
r¼1
Z 9v91,Sr 9
r¼1
:
From the Green formula, ZZ Z 2 9v91,Sr ¼ rvrv ds ¼ vn v d‘, Sr
ð3:35Þ
ð3:36Þ
‘r
‘r ¼ fðr, yÞ9r ¼ r,0 r y r pg: D2i ,
i¼0
JvJ 0, Gn ¼
JvJ21,S rCJvJ21=2, Gn :
ð3:32Þ
where ‘r is the semi-circle
In (3.23), the notations are L X
Hence we obtain from (3.27) and (3.31),
On the other hand, since Dv ¼ 0 for v A VL , we have from Oden and Reddy [34, p.189],
where lmax and lmin are the maximal and the minimal eigenvalues of A, respectively, defined by
ðx,xÞ ¼ JxJ2 ¼
In (3.31), the semi-norms in the Sobolev space is defined by ( )1=2 Z Z ðvðPÞvðQ ÞÞ2 2 d‘ðPÞ d‘ðQ Þ : JvJ1=2, G ¼ JvJ0, G þ JPQ J2 G G
ð3:20Þ
where the matrix A ¼ FT F A RðL þ 1ÞðL þ 1Þ is symmetric and positive definite, and d ¼ FT b. Then there exist the relations pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi smax ðFÞ ¼ lmax ðAÞ, smin ðFÞ ¼ lmin ðAÞ, ð3:22Þ
ðAx,xÞ , ðx,xÞ
ð3:31Þ
ðAx,xÞ ¼ 2JvJ 0, Gn Zc0 JvJ20, Gn Z c1 Lp JvJ21=2, Gn :
ð3:21Þ
xa0
JvJ1=2, Gn rCðJvJ1, Gn JvJ0, Gn Þ1=2 r C1 Lp=2 JvJ0, Gn :
ð3:19Þ
where smax and smin are the maximal and the minimal singular values of F, respectively. In this application, rankðFÞ ¼ n ¼ Lþ 1, smax ¼ s1 and smin ¼ sn 4 0. The normal equation of (3.15) is given by
lmax ðAÞ ¼ max
Proof. We have from (3.29) and Babuska and Aziz [5, p. 21],
2
JbJ , smin JxJ
and that of Cond ¼
ð3:30Þ
v2 d‘:
ð3:26Þ
Gn
3 Motz’s problem in [33] with mixed type of Dirichlet and Neumann boundary conditions cannot be used as the test model for the BIE of first kind, since the latter results from only the Dirichlet boundary condition.
By calculus, from the orthogonality of cosðiþ 12Þy we obtain from (3.3) ) Z p (X Z L 1 i1=2 1 vn v d‘ ¼ Di iþ r cos i þ y 2 2 ‘r 0 i¼0 ( ) L X 1 Di ri þ 1=2 cos iþ y r dy 2 i¼0 Z L X 1 2i þ 1 p 1 ¼ D2i i þ r cos2 iþ y dy 2 2 0 i¼0 L pX 1 2i þ 1 ¼ D2i i þ r : ð3:37Þ 2 2i¼0
60
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
Then when r ¼ 1, Z L L pX 1 pX Z vn v d‘ ¼ D2i i þ D2 : 2 2i¼0 4i¼0 i ‘r
ð3:38Þ
Moreover, since S Sr 9r ¼ pffiffi2 , we have from (3.37)
Z L pX 1 2i þ 1
2 2 vn v d‘ ¼ D2i i þ r 9v91,S r 9v91,Spffi ¼
2
2 2i¼0 ‘pffi2
Combining (3.34)–(3.36) and (3.38) yields IðvÞ Z c0 Lp
L X
r D2i :
ð3:39Þ
i¼0
^ IðvÞ IðvÞ Z c0 Z c1 Lp , ðx,xÞ x a 0 ðx,xÞ
lmin ðAÞ ¼ min
L 1 pffiffiffi 2L þ 1 X ð 2Þ Lþ D2i : 2 2 i¼0
p
L pffiffiffi X D2i : IðvÞ rCLð 2Þ2L
ð3:53Þ
smin ¼ smin ðFÞ Zc1 Lp=2 :
Hence from (3.28) and (3.53), the maximal eigenvalue lmax ðAÞ has the following bound: ð3:41Þ
This is the desired result (3.30). This completes the proof of Lemma 3.1. & Theorem 3.1. Let (3.27) and (3.29) hold. Then for (3.15) there exists the following bound of the effective condition number :
ð3:54Þ
i¼0
ð3:40Þ
and from (3.22)
Cond_eff r CL
ð3:52Þ
Combining (3.51) and (3.52) yields
Hence we have from (3.27) and (3.39)
p=2
pffiffi
r¼ 2
ð3:42Þ
lmax ðAÞ ¼ max xa0
^ IðvÞ JxJ2
r C max xa0
IðvÞ JxJ2
pffiffiffi r C1 Lð 2Þ2L :
Then we have from (3.22), pffiffiffi pffiffiffi
smax ¼ smax ðAÞ r C Lð 2ÞL :
ð3:55Þ
ð3:56Þ
This is the desired result (3.47), and this completes the proof of Lemma 3.2. & Based on Lemmas 3.1 and 3.2, we have the following theorem.
Proof. The vector b in (3.15) has the components pffiffiffi pffiffiffi T b ¼ f. . . , ai h, . . . , ai yi h, . . .g,
ð3:43Þ
where 0 o yi o1 and ai is the integration weights with ai ¼ Oð1Þ. Since h r C1=m, where m is the dimension of b, we have sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X X pffiffiffi ffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffi JbJ r ðai hÞ2 r a2i h r max9ai 9 mh r C: ð3:44Þ i
i
i
Also since the lower bounds of the true coefficients are known in Table 3, and the vector x satisfies vffiffiffiffiffiffiffiffiffiffiffiffiffiffi u L uX D2i ZD0 Z 1: ð3:45Þ JxJ ¼ t i¼0
Theorem 3.2. Let (3.27) and (3.29) hold. Then for (3.15) there exists the bound for the traditional condition number pffiffiffi smax Cond ¼ r CLðp þ 1Þ=2 ð 2ÞL : ð3:57Þ
smin
For the sectorial S, p ¼1 is given in [25, p. 112]. Hence we have the following corollary from Lemmas 3.1 and 3.2. Corollary 3.1. Let (3.27) and (3.29) hold. Also assume p ¼1. Then for (3.15), there exist the bounds
smin Z c0 L1=2 ,
ð3:58Þ
pffiffiffi pffiffiffi smax rC Lð 2ÞL :
ð3:59Þ
Hence we have from Lemma 3.1, (3.44) and (3.45), Cond_eff ¼
JbJ rCLp=2 : smin JxJ
From Theorems 3.1 and 3.2, we have the following corollary. ð3:46Þ
This is the desired result (3.42), and completes the proof of Theorem 3.1. & Next, we will estimate the upper bound of smax , to have the following lemma. Lemma 3.2. Let (3.27) hold. Then for (3.15) there exists the upper bound, pffiffiffi pffiffiffi smax ¼ smax ðFÞ rC Lð 2ÞL : ð3:47Þ Proof. From the embedding theorem in Ciarlet [13], IðvÞ ¼ JvJ20, Gn
r CJvJ21,S :
ð3:48Þ
From (3.27) and (3.48) we have 2
^ ¼ 2JvJ n r2CIðvÞ ¼ CJvJ2 n rC1 JvJ2 : ðAx,xÞ ¼ 2IðvÞ 0, G 1,S 0, G
ð3:49Þ
Since v9y ¼ 041 o x o 0 ¼ 0 for v A VL , there exists the bound from the Poincare inequality [13], JvJ1,S rC9v91,S ,
ð3:50Þ
where 9v91,S is the semi-norm of v on S. Hence we have from (3.48) and (3.50), 2
IðvÞ r C9v91,S :
ð3:51Þ
Corollary 3.2. Let (3.27) and (3.29) hold. Also assume p ¼1. Then for (3.15) there exist the bounds Cond_eff r CL1=2 ,
ð3:60Þ
pffiffiffi Condr CLð 2ÞL :
ð3:61Þ
Eqs. (3.59) and display that when L-1, Cond_eff and pffiffi(3.60) ffi Cond grow as Oð LÞ and exponentially, respectively. Now, let us examine again the data in Table 1. We can see JeJB ¼ Oðð0:56ÞL Þ,
JeJ1, Gn ¼ Oðð0:56ÞL Þ,
ð3:62Þ
Cond_eff ¼ Oð1Þ,
pffiffiffi Cond ¼ Oðð1:4ÞL Þ ¼ Oðð 2ÞL Þ:
ð3:63Þ
Eqs. (3.63) are consistent with Corollary 3.2. From Table 5, we can see pffiffiffi smin 0:7, smax ¼ ð1:4ÞL Oðð 2ÞL Þ, ð3:64Þ JbJ ¼ 1:0,
bn 0:67, Cond_EE 1:57:
ð3:65Þ
Eqs. (3.64) also agree with Corollary 3.1. Besides, the error analysis in [33] proves that the CTM has the exponentially convergence rates, which is also validated by (3.62). From this section, the CTM has an excellent stability, based on Cond_eff. Hence, the CTM becomes a very efficient boundary method, to
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
compete with other boundary methods, such as the boundary element method (BEM) and the BIE method, see [28,29]. 3.3. Linkage to Section 2 Finally, let us link the numerical computation in this section to the theoretical estimation (2.25) in Section 2. Take the leading coefficient D0 in Table 3 for example, which only has an error less than the unit of rounding error. Such an extreme accuracy can be explained by the effective condition number, but not by the traditional condition number. Denote by g the rounding error of the computer used. Since the double precision is used (say, g ¼ 1017 ), coefficient D0 has the relative error 9DD0 9 r g ¼ 1017 : 9D0 9
ð3:66Þ
For the coefficients Di in Table 3, we find vffiffiffiffiffiffiffiffiffiffiffiffiffiffi u L uX D2i ¼ mJxJ, 9D0 9 ¼ mt where L¼54 and the factor m ¼ 1:2. Then we have qP ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi L 2 9DD0 9 9DD0 9 JDxJ i ¼ 0 DDi ¼ r ¼ : mJxJ mJxJ mJxJ 9D0 9
ð3:67Þ
ð3:68Þ
ð3:70Þ
Since the errors JDAJ=JAJ from the solution method will satisfy JDAJ 1 1 5 ¼ , JAJ Cond 0:286ð8Þ
ð3:71Þ
in double precision, we have d 51. Therefore, we have from (3.69) JDxJ JDbJ r 1:25 : JxJ JbJ
ð3:72Þ
Combining (3.67) and (3.71) gives 9DD0 9 1:25 JDbJ 1:25 JDbJ JDbJ ¼ ¼ 1:02 : r JbJ 1:2 JbJ JbJ m 9D0 9
ð3:73Þ
Since the relative error of b results from rounding error, we have 9DD0 9 ð r 1017 Þ r 1:02 g ¼ 1:02ð17Þ: 9D0 9
ð3:74Þ
How perfect coincidence between the numerical computation and the theoretical estimation! Eq. (3.74) strongly verifies the estimation (2.25), which indicates a key role of effective condition number in (2.16), (2.18) and (2.19). This explores a linkage of this section to Section 2. For ðF þ DFÞðxþ DxÞ ¼ b þ Db, we may reform it as n
Fðx þ DxÞ ¼ b þ Db ,
ð3:75Þ
where
Dbn ¼ DbDFðx þ DxÞ:
ð3:78Þ
we also obtain (3.74). This again exhibits the significance of effective condition number for stability analysis. On the other hand, for the solution x (i.e., coefficients Di) in Table 3, from (1.1) we have the error bounds JDxJ JDbJ JDbJ rCond ¼ 0:286ð8Þ JxJ JbJ JbJ r0:286ð8Þ g ¼ 0:286ð9Þ:
ð3:79Þ
Combining (3.67) and (3.78) yields 9DD0 9 JDbJ ð r 1017 Þ rCond mJbJ 9D0 9 ¼ Cond
g 0:286ð8Þ 1017 ¼ m m
0:286ð9Þ ¼ 2:38ð10Þ: 1:2
ð3:80Þ
In (3.80) there exists a big gap (about a factor 107) between the actual error and the estimation of D0.4 This fact illustrates inefficacy of Cond. Similarly for ðF þ DFÞðxþ DxÞ ¼ b þ Db, we have JDxJ JDb J rCond , JxJ JbJ
ð3:81Þ
and also reach (3.80), which is misleading, indeed.
where JDAJ : JAJ
n
JDb J r g ¼ 1017 , JbJ
n
From Table 1 at L¼54, Cond_eff¼1.25 and Cond¼0.286 (8). The relative errors of x are given from (2.25), JDxJ 1 JDbJ 1 JDbJ r ¼ , ð3:69Þ Cond_eff 1:25 JxJ 1d JbJ 1d JbJ
d ¼ Cond
Since the minimal bound is also given by the rounding error
¼
i¼0
61
ð3:76Þ
4. Concluding remarks To close this paper, let us address the new contributions in this paper. (1) For weighted linear least squares problem with both full rank and rank-deficient cases, we develop the effective condition number defined in (2.15), and derive in Theorem 2.1 the new error bounds from the effective condition number. These are the new computational formulas, compared to the existing literature including our previous study in [26,27]. (2) In Section 3.3, the extremely accurate leading coefficient D0 is explained by the very small effective condition number. In contrast, the huge condition number is misleading for stability analysis. Moreover, the close linkage between Sections 2 and 3 is explored, and Eq. (2.25) is verified perfectly by numerical computation of Motz’s problem by the CTM. This indicates a key role of effective condition number in (2.16), (2.18) and (2.19). (3) In [29,28], the error analysis is explored for Trefftz and collocation Trefftz methods. Although the numerical results of Cond_eff are provided, but no strict stability analysis was given wherein. This paper is devoted to the new stability analysis based on Cond_eff, to show an advantage of the collocation Trefftz method (CTM) over other boundary methods. Moreover, the analysis and numerical examples in this paper display that the effective condition number is particularly significant to the CTM and the spectral method, where the huge condition number is misleading. (4) By Axelsson’s work [1–3], the effective condition number is beneficial to studying the iteration solution methods. From [26,27] and this paper, the effective condition number is important to numerical PDE, and it may become a new trend of stability analysis.
From (1.1) we have n
JDxJ JDb J r Cond_eff : JxJ JbJ
ð3:77Þ
4 Eq. (3.80) indicates that D0 has at most 10 significant digits. However, the numerical D0 in Table 3 (also see many numerical values of D0 in Table 1) has 17 significant digits!
62
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
References [1] Axelsson O. Iterative solution methods. Cambridge: Cambridge University Press; 1994. [2] Axelsson O, Kaporin I. On the sublinear and superlinear rate of convergence of conjugate gradient methods. Numer Algorithms 2000;25:1–22. [3] Axelsson O, Kaporin I. Error norm estimation and stopping criteria in preconditioned conjugate gradient iterations. Numer Linear Algebra Appl 2001;8:265–86. [4] Atkinson KE. An introduction to numerical analysis. New York: John Wiley & Sons; 1989. [5] Babuska I, Aziz AK. Survey lectures in the mathematical foundations of the finite elements. In: Aziz AK, editor. The mathematical foundation of the finite element methods with applications to partial differential equations. New York, London: Academic Press; 1973. p. 3–369. [6] Barlow JL, Handy SL. The direct solution of weighted and equality constrained least squares problems. SIAM J Sci Stat Comput 1988;9:704–16. [7] Banoczi JM, Chiu M, Cho GE, Ipsen CF. The lack of influence of right-hand side on accuracy of linear system solution. SIAM J Sci Comput 1998;20: 203–27. ˚ Numerical methods for least squares problems. Philadelphia: SIAM; ¨ [8] Bjorck A. 1996. [9] Bobrovnikova EY, Vavasis SA. Accurate solution of weighted least squares by iterative methods. SIAM J Matrix Anal Appl 2001;22:1153–74. [10] Campbell S, Meyer C. Generalized inverses of linear transformations. New York: Pitman, Dover; 1979; Campbell S, Meyer C. Generalized Inverses of Linear Transformations. New York: Dover; 1991. [11] Chan TF, Foulser DE. Effectively well-conditioned linear systems. SIAM J Sci Stat Comput 1988;9:963–9. [12] Christiansen S, Hansen PC. The effective condition number applied to error analysis of certain boundary collocation methods. J Comput Appl Math 1994;54:15–36. [13] Ciarlet PG. Basic error estimates for elliptic problems. In: Ciarlet PG, Lions JL, editors. Finite element methods (Part I). North-Holland; 1991. p. 17–352. [14] Cucker F, Diao H, Wei Y. On mixed and componentwise condition numbers for Moore–Penrose inverse and linear least squares problems. Math Comput 2007;78(258):947–63. [15] Golub GH, van Loan CF. Matrix computation.3rd ed. Baltimore, London: The John Hopkins University Press; 1996. [16] Gulliksson M. Iterative refinement for constrained and weighted linear least squares. BIT 1994;34:239–53. [17] Gulliksson M. On the modified Gram–Schmidt algorithm for weighted and constrained linear least squares problems. BIT 1995;35:453–68. [18] Gulliksson M. Backward error analysis for the constrained and weighted linear least squares problem when using the weighted QR decomposition. SIAM J Matrix Anal Appl 1995;16:675–87. [19] Gulliksson ME, Wedin PA, Wei Y. Perturbation identities for regularized Tikhonov inverses and weighted pseudoinverses. BIT 2000;40:513–23. [20] Gulliksson ME, Jin X, Wei Y. Perturbation bound for constrained and weighted least squares problem. Linear Algebra Appl 2000;349:221–32.
[21] Higham NJ. Accuracy and stability of numerical algorithm.2nd ed. Philadelphia: SIAM; 2002. [22] Hough PD, Vavasis SA. Complete orthogonal decomposition for weighted least squares. SIAM J Matrix Anal Appl 1997;18:369–72. [23] Kenney CS, Laub AJ, Reese MS. Statistical condition estimation for linear least squares. SIAM J Matrix Anal Appl 1998;19:906–23. [24] Lawson CL, Hanson RJ. Solving least squares problems. Englewood Cliffs, NJ: Prentice-Hall, Inc.; 1974. Reprint, SIAM, Philadelphia, 1995. [25] Li ZC. Combined methods for elliptic equations with singularities, interfaces and infinities. Dordrecht, Boston: Kluwer Academic Publishers; 1998. [26] Li ZC, Chien CS, Huang HT. Effective condition number for finite difference method. J Comput Appl Math 2007;198:208–35. [27] Li ZC, Huang HT. Effective condition number for numerical partial differential equations. Numer Linear Algebra Appl 2008;15:575–94. [28] Li ZC, Lu TT, Huang HT, Cheng AH-D. Trefftz, collocation and other boundary methods, a comparison. Numer Meth PDE 2007;23:93–144. [29] Li ZC, Lu TT, Hu HY, Cheng AH-D. Trefftz and collocation methods. Southampton, Boston: WIT Publishers; 2008. [30] Li ZC, Mathon R, Sermer P. Boundary methods for solving elliptic problem with singularities and interfaces. SIAM J Numer Anal 1987;24:487–98. [31] Lin L, Lu TT, Wei Y. On level-2 condition number for the weighted Moore– Penrose inverse. Comput Math Appl 2008;55:788–800. [32] Lu TT, Chang CM, Huang HT, Li ZC. Stability analysis of Trefftz methods for the stick-slip problem. Eng Anal Boundary Elem 2009;33:474–84. [33] Lu TT, Hu HY, Li ZC. Highly accurate solutions of Motz’s and the cracked beam problems. Eng Anal Boundary Elem 2004;28:1387–403. [34] Oden JT, Reddy JN. An introduction to the mathematical theory of finite elements. New York: John Wiley & Sons; 1976. [35] Rice JR. Matrix computations and mathematical software. New York: McGraw-Hill Book Company; 1981. [36] Rump S. Ill-conditionedness need not be the componentwise near to illposedness for least squares problems. BIT 1999;39:143–51. [37] van Loan CF. Generalizing the singular value decomposition. SIAM J Numer Anal 1976;13:76–83. [38] van Loan CF. On the method of weighting for equality-constrained least squares problems. SIAM J Numer Anal 1985;22:851–64. [39] Vavasis SA. Stable numerical algorithms for equilibrium systems. SIAM J Matrix Anal Appl 1994;15:1108–31. [40] Wang G, Wei Y, Qiao S. Generalized inverses: theory and computations. Beijing: Science Press; 2004. [41] Wei M, De Pierro AR. Upper perturbation bounds of weighted projections, weighted and constrained least squares problem. SIAM J Matrix Anal Appl 2000;21:931–51. [42] Wei M. Supremum and stability of weighted pseudoinverses and weighted least squares problems: analysis and computations. New York: Nova Science Publishers; 2001. [43] Wilkinson JH. The algebraic eigenvalue problems. Oxford University Press; 1965. [44] Xiang H, Wei Y. Structured mixed and componentwise condition numbers of some structured matrices. J Comput Appl Math 2007;202:217–29. [45] Xu W, Wei Y, Qiao S. Condition numbers for structured least squares problems. BIT 2006;46:203–25.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Effective condition number for weighted linear least squares problems and applications to the Trefftz method Yi-min Wei a,b,1, Tzon-Tzer Lu c, Hung-Tsai Huang d, Zi-Cai Li e,f, a
School of Mathematical Science, Fudan University, PR China Shanghai Key Laboratory of Contemporary Applied Mathematics, Shanghai, 200433, PR China c Department of Applied Mathematics, National Sun Yat-sen University, Kaohsiung, Taiwan d Department of Applied Mathematics, I-Shou University, Kaohsiung, Taiwan e Department of Applied Mathematics, National Sun Yat-sen University, Kaohsiung, Taiwan f Department of Computer Science and Engineering, National Sun Yat-sen University, Kaohsiung, Taiwan b
a r t i c l e i n f o
abstract
Article history: Received 14 March 2011 Accepted 14 July 2011 Available online 23 August 2011
In [27], the effective condition number Cond_eff is developed for the linear least squares problem. In this paper, we extend the effective condition number for weighted linear least squares problem with both full rank and rank-deficient cases. We apply the effective condition number to the collocation Trefftz pffiffiffi Laplace’s equation with a crack singularity, to prove that Cond_eff pffiffiffimethod (CTM) [29] for ¼ Oð LÞ and Cond ¼ OðL1=2 ð 2ÞL Þ, where L is the number of singular pffiffiffi particular solutions used. The Cond grows exponentially as L increases, but Cond_eff is only Oð LÞ. The small effective condition number explains well the high accuracy of the TM solution, but the huge Cond cannot. & 2011 Elsevier Ltd. All rights reserved.
Keywords: Condition number Effective condition number Perturbation Weighted linear least squares problem Collocation Trefftz method Singularity problem
1. Introduction Condition number measures the sensitivity of the solutions of a problem with respect to small perturbations of the given data. In [27], the effective condition number is developed for the linear least squares problem by following Rice [35], Chan and Foulser [11], Christiansen and Hansen [12], Axelsson and Kaporin [2,3], componentwise linear least squares [14,36], and structured Kenny, et al. [23], matrix [44,45]. The traditional condition number was given in Wilkinson [43] in 1965 for linear algebraic equations, and has been used in many textbooks and papers. In Rice [35], the natural condition number was first proposed for specific non-homogenous term; we choose the effective condition number used [11] for our study [26,27]. More references for condition number and effective condition number are provided in [27]. In this paper, we extend the effective condition number to Corresponding author at: Department of Applied Mathematics, National Sun Yat-sen University, Kaohsiung, Taiwan. E-mail addresses: [email protected], [email protected] (Y.-m. Wei), [email protected] (T.-T. Lu), [email protected] (H.-T. Huang), [email protected] (Z.-C. Li). 1 This author is supported by the National Natural Science Foundation of China under Grant 10871051, Doctoral Program of the Ministry of Education under Grant 20090071110003, Shanghai Science and Technology Committee under Grant 09DZ2272900 and 973 Program Project under Grant 2010CB327900.
0955-7997/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2011.07.010
weighted linear least squares problem with both full rank and rank-deficient cases. For the over-determined system Fx ¼b with F A Rmn ðm Z nÞ and rankðFÞ ¼ r r n, and the perturbed system Fðx þ DxÞ ¼ b þ Db, there exist the bounds of relative errors, JDxJ JDbJ rCond , JxJ JbJ
JDxJ JDbJ rCond_eff , JxJ JbJ
ð1:1Þ
where JxJð ¼ JxJ2 Þ is the 2-norm, the traditional condition number and the effective condition number are defined by (see [4,15,26]) Cond ¼
s1 JbJ , , Cond_eff ¼ sr sr JxJ
ð1:2Þ
and sr is the r-th singular value of matrix F. For ðF þ DFÞ ðx þ DxÞ ¼ b þ Db with the nonsingular F ¼ A A Rnn , we also derive in [27] the following computational formula: JDxJ 1 JDAJ JDbJ r þCond_eff , ð1:3Þ Cond JxJ 1d JAJ JbJ where d ¼ JAJ=sn o1. Based on (1.3), we provide in [27] a theoretical analysis that for solving the linear system by Gaussian elimination or the QR factorization, the direction of right hand (i.e., b) is insignificant for solution errors. Such a conclusion was discovered in Banoczi et al. [7] purely by the numerical experiments. We also apply the effective condition number for the finite
54
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
difference method (FDM) for Poisson’s equations, where the discretization errors are dominant, and the direction of b does influence the solution errors. In this paper, we consider the weighted linear least squares problem. Let M A Rmm be symmetric and positive definite matrix. For Fx ¼b, we may solve the weighted linear least squares problem, M1=2 Fx ¼ M1=2 b:
ð1:4Þ
For any symmetric and positive definite matrix N A R may be reduced to
nn
Gy ¼ d,
Hence, the generalized effective condition number in (1.7) can be regarded as one of the effective condition number for weighted linear least squares problem explored in this paper. This paper is organized as follows. In Section 2, the computational formulas are derived for effective condition number of weighted linear least squares problem. In Section 3, the Trefftz method for Laplace’s equation is discussed, and stability analysis is made by the effective condition number, and a linkage to Section 2 is explored. In last section, a few final remarks are made.
, Eq. (1.4) ð1:5Þ
where
2. Effective condition number of weighted linear least squares problem 2.1. Effective condition number
G ¼ M1=2 FN1=2 A Rmn , y ¼ N1=2 xA Rn ,
d ¼ M1=2 b A Rm :
ð1:6Þ
In this paper, we will extend the effective condition number to the weighted linear least squares problem (1.5). Note that Eq. (1.3) in [27] is valid only for A A Rnn . Hence, the error bounds like (1.3) of the weighted linear least squares problem are significant in application. Furthermore, more applications of effective condition number are provided for linear algebraic equations and numerical partial differential equations (PDE), such as the Trefftz method for Poisson’s equation with singularities. For solving Ax¼b with the nonsingular matrix AA Rnn , the effective condition number has been used in Axelsson and Kaporin [2,3] and Banoczi et al. [7]. Denote the norm JxJD ¼ ðxT DxÞ1=2 , where D is the symmetric and positive definite matrix. The generalized effective condition number is defined in Axelsson and Kaporin [2] as JAn xJD n kn ðA,xÞD ¼ JA JD , JxJD
ð1:7Þ
where n a0. For the preconditioned conjugate gradient methods, both the errors and the convergence rates may be estimated, by means of kn ðA,xÞD . In fact, the generalized effective condition number in (1.7) is of the effective condition number for weighted linear least squares problem explored in this paper. For Ax¼ b with the symmetric and positive definite matrix AA Rnn , we may solve the following weighted linear equations with M¼ N¼D, D1=2 Ax ¼ D1=2 b:
ð1:8Þ
The generalized effective condition number of (1.7) is given by
k1 ðA,xÞD ¼
JAxJD 1 JbJD 1 JA JD ¼ JA JD ¼ Cond_effðAÞ: JxJD JxJD
A A Cmn ,
ð2:1Þ
JBJNM ¼ max JByJN ,
B A Cnm ,
ð2:2Þ
JxJN ¼ 1
JyJM ¼ 1
as in [40]. Then we obtain the following relations: JAJMN ¼ JM1=2 AN1=2 J2 ,
A A Cmn ,
ð2:3Þ
JBJNM ¼ JN1=2 BM1=2 J2 ,
B A Cnm :
ð2:4Þ
Lemma 2.1 (van Loan [37]). Let AA Cmn ðm Z nÞ with rankðAÞ ¼ r, and M and N be Hermitian positive definite matrices of order m and n, respectively. There exist U A Cmm and V A Cnn , satisfying UH MU ¼ I and VH N1 V ¼ I, such that D 0 ð2:5Þ A¼U VH , 0 0 and the weighted Moore–Penrose inverse AyMN can be represented by ! D1 0 ð2:6Þ AyMN ¼ N1 V UH M: 0 0 In (2.5) D ¼ diagðm1 , m2 , . . . , mr Þ, m1 Z m2 Z Z mr 4 0 and m2i is the nonzero eigenvalue of N1 AH MA, where mi is called the weighted singular values, and JAJMN ¼ m1 ,
ð1:10Þ
JAyMN JNM ¼
1
mr
:
ð2:8Þ
x
n
Gn x ¼ d ,
n
mn
ð1:11Þ
where matrix A A C and m Z n. Suppose that the rank of A is r with r r n r m. Then the least-squares (M) solution in the minimum-norm (N) is given by
ð1:12Þ
x ¼ AyMN b:
where d ¼ An1 b:
ð2:7Þ
Consider the weighted linear least squares problem with rankdeficient [8,21,42], minJAxbJM ,
and then
Gn ¼ An ,
JAJMN ¼ max JAxJM ,
ð1:9Þ
Moreover, from Ax¼b we have the following weighted linear equations: An x ¼ An1 b,
Consider the weighted norm, and for generality, we use the complex matrices. Given two vectors x A Cn , y A Cm and two Hermitian positive definite matrices, MA Cmm and N A Cnn . The weighted vector norms may be defined by JyJM ¼ JM1=2 yJ2 and JxJN ¼ JN1=2 xJ2 , and the weighted matrix norms by
Hence, the generalized effective condition number (1.9) for (1.11) is given by
Assume that there exists a perturbation of A and b such that rankðAÞ ¼ rankðA þ DAÞ ¼ r,
n
JGn xJD Jd JD k1 ðG ,xÞD ¼ JðGn Þ1 JD ¼ JðGn Þ1 JD JxJD JxJD
we have the perturbed weighted linear least squares problem [6,8,9,16–20,22,31,38,39,41,42],
n
¼ Cond-effðGn Þ ¼
JAn xJD n JA JD ¼ kn ðA,xÞD : JxJD
ð1:13Þ
min JðA þ DAÞðxþ DxÞðb þ DbÞJM :
x þ Dx
ð2:9Þ
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
Then we can find the least-squares (M) solution in the minimumnorm (N) of (2.9) by
55
Remark 2.1. If vector M1=2 b is just parallel to u1 , i.e.,
b2 ¼ ¼ bn ¼ 0, we have JbJM ¼ 9b1 9 and JxJN ¼ 9b1 9=m1 . Hence,
ð2:10Þ
the effective condition number (2.15) just leads to the traditional condition number in (2.11).
Let matrix A be decomposed by the weighted singular value decomposition (2.5) D 0 A¼U VH , 0 0
Remark 2.2. When AA Rmn , and M and N are the identity matrices, Eqs. (2.14) and (2.15) lead to (1.1) and (1.2), respectively.
x þ Dx ¼ ðAþ DAÞyMN ðb þ DbÞ:
where matrices U A Cmm and V A Cnn are weighted unitary, and D A Rrr is a diagonal matrix with the positive weighted singular values mi as m1 Z m2 Z Z mr 4 0: The traditional condition number is defined by van Loan [37] CondMN ðAÞ ¼
JAJMN JAyMN JNM
¼ m1 =mr :
ð2:11Þ
First let us consider (2.9) for simple case of DA ¼ 0. If DAa 0, then the perturbation bounds from effective condition number are given in Theorem 2.1 later. Denote M1=2 U ¼ ½u1 ,u2 , . . . ,um , we have M1=2 b ¼
m X
bi ui , M1=2 Db ¼
i¼1
m X
ai ui ,
ð2:12Þ
First we cite a lemma from [31,40]. Lemma 2.2. Let matrices A, DA A Cmn with rankðAÞ ¼ rankðAþ DAÞ ¼ r r n rm and denote d ¼ JAyMN JNM JDAJMN o1. Then there exist the upper bounds pffiffiffi y 1 þ 5 JAMN JNM JDAJMN JAy JNM , JðAþ DAÞyMN AyMN JNM r 2 1JAyMN JNM JDAJMN MN JðAþ DAÞyMN JNM r
i¼1
1=2 1=2 where bi ¼ uH b and ai ¼ uH Db. It is easy to check that i M i M ! D1 0 x ¼ AyMN b ¼ N1 V UH Mb, 0 0
Dx ¼ AyMN Db ¼ N1 V
D1
0
0
0
! UH MDb:
Since M1=2 U and N1=2 V are unitary, we can obtain ! D1 0 1=2 1=2 1=2 H 1=2 VÞ JxJN ¼ JN xJ2 ¼ ðN ðM UÞ M b 0 0 2 vffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 r r X uX b bi b i ¼ u ¼t Z r, mi i m2i mr i¼1
2.2. Perturbation bounds
1JAyMN JNM JDAJMN
:
We have the following theorem. Theorem 2.1. Let matrices A, DA A Cmn with rankðAÞ ¼ rankðAþ DAÞ ¼ r rn rm and denote d ¼ JAyMN JNM JDAJMN o 1. Then we have " # pffiffiffi JDxJN Cond_eff 1þ 5 JDAJMN JDbJM CondMN ðAÞ r þ : 1d 2 JxJN JAJMN JbJM ð2:16Þ Proof. Since x ¼ AyMN b and x þ Dx ¼ ðAþ DAÞyMN ðb þ DbÞ, we have
Dx ¼ ðA þ DAÞyMN ðb þ DbÞAyMN b ¼ ½ðA þ DAÞyMN AyMN b þ ðAþ DAÞyMN Db,
i¼1
2
JAyMN JNM
and then
and ! D1 0 JDxJN ¼ JN1=2 DxJ2 ¼ ðN1=2 VÞ ðM1=2 UÞH M1=2 Db 0 0 2 vffiffiffiffiffiffiffiffiffiffiffiffiffiffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u X r r r 2 X X u u a ai t 1 i r t ¼ u ¼ a2i rJAyMN JNM JDbJM : mi i m2i mr i¼1
2
i¼1
JDxJN ¼ JN1=2 f½ðAþ DAÞyMN AyMN b þ ðA þ DAÞyMN DbgJ2 ¼ JN1=2 ½ðAþ DAÞyMN AyMN M1=2 M1=2 b þN1=2 ðA þ DAÞyMN M1=2 M1=2 DbJ2 r JðAþ DAÞyMN AyMN JNM JbJM þ JðA þ DAÞyMN JNM JDbJM :
i¼1
ð2:13Þ Hence we have JDx9N JDbJM 1 JDbJM r sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ Cond_eff , 2 JxJN mr JbJM P b r i¼1
ð2:14Þ
i 2 i
þ
m
JbJM JbJM sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2ffi ¼ m JxJN ¼ r b1 br mr þ þ
m1
JAyMN JNM 1JAyMN JNM JDAJMN
JDbJM ,
to give
where the effective condition number is defined by Cond_eff ¼
It follows from Lemma 2.2 pffiffiffi y 1þ 5 JAMN JNM JDAJMN JAy JNM JbJM JDxJN p 2 1JAyMN JNM JDAJMN MN
JAyMN JNM JbJM JxJN
:
mr
ð2:15Þ When br a 0, JDxJN JDbJM mr JbJM JDbJM r ¼ : JxJN mr br br JbJM So the simplest effective condition number can be defined by Cond_EE ¼ JbJM =br from (2.15), see [26]. Clearly, Cond_eff is a tighter bound than Cond_EE.
JDxJN r
pffiffiffi 1þ 5 d JbJM 1 JDbJM þ : 2 1d 1d mr mr
Since d ¼ CondMN ðAÞ JDAJMN =JAJMN , the desired result (2.16) is obtained and the proof of Theorem 2.1 is completed. & If rankðAÞ ¼ rankðA þ DAÞ ¼ n, then the sharp error bound is given by (see [40]) JDxJN Cond_eff pffiffiffi JDAJMN JDbJM r þ 2CondMN ðAÞ : ð2:17Þ 1d JxJN JAJMN JbJM Besides, by following [10], we may also derive similar bounds of Cond_eff, under the same assumptions in this paper.
56
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
2.3. Applications and comparisons Consider the real matrix AA Rmn ðm Z nÞ and the real vectors x A Rn and b A Rm . Let MA Rmm and NA Rnn be the identity matrices. Then we have from Theorem 2.1, " # pffiffiffi JDxJ Cond_eff 1þ 5 JDAJ JDbJ r Cond þ , ð2:18Þ JxJ 1d 2 JAJ JbJ y
where d ¼ JA JJDAJ o1, and JxJ ¼ JxJ2 . Moreover, when rankðAÞ ¼ rankðA þ DAÞ ¼ n, the following bound is given from (2.17), JDxJ Cond_eff pffiffiffi JDAJ JDbJ r þ : ð2:19Þ 2CondðAÞ JxJ 1d JAJ JbJ Let us compare our new error bounds with those in existing literature. Although the bound of (2.19) is larger than that in (1.3), which is only for the nonsingular A A Rnn , Theorem 2.1 is valid for weighted linear least squares problem. In Golub and van Loan [15, p. 228], there is the error bound for A A Rmn ðmZ nÞ JDxJ 2 r e Cond Cond tan y þ þ Oðe2 Þ, ð2:20Þ JxJ cos y where JDAJ JDbJ sn 1 , o ¼ JAJ JbJ s1 Cond
e ¼ max
ð2:21Þ
and sin y ¼
rLS JbJ
,
rLS ¼ JAxLS bJ:
ð2:22Þ
Obviously, the bound in (2.18) is smaller than (2.20), by noting 2 OðCond Þ in (2.20). Moreover, in Lawson and Hanson [24], for mn AA R ðm Z nÞ with rankðA þ DAÞ ¼ rankðAÞ r n, there exists the error bound (see (9.10) in [24, p.51]) JDxJ Cond n rLS o JDAJ JbJ JDbJ r þ 2 þ Cond JxJ 1d JAJ JAJJxJ JbJ JAJJxJ n 1 rLS o JDAJ JDbJ þ Cond_eff , ¼ Cond 2þ Cond 1d JAJ JbJ JAJJxJ ð2:23Þ where d ¼ Cond JDAJ=JAJ and rLS is the same as (2.22). The bound in (2.18) is also smaller than that in (2.23), when Cond is large. Hence (2.16), (2.18) and (2.19) are new contributions of effective condition number in this paper. Note that when AA Rnn , then rLS ¼ 0. Eq. (2.23) leads to JDxJ 1 JDAJ JDbJ r þ Cond_eff , ð2:24Þ 2 Cond JxJ 1d JAJ JbJ which is similar to (1.3). In (2.16), (2.18), and (2.24), if the term JDbJ=JbJ is dominant in the brackets, we have JDxJ Cond_eff JDbJ r , JxJ 1d JbJ
ð2:25Þ
which is very close to the right hand side of (1.2). This is the case occurring in numerical PDE. The same conclusion is drawn as in [27]. Since Eqs. (2.16), (2.18), and (2.24) may lead to (2.25), this section provides a theoretical foundation of effective condition number for weighted linear least squares problems, which is an advance level of [26,27]. The collocation Trefftz method for Laplace’s equation in Section 3 can be classified into weighted linear least squares problems, since different weights are used in the central and the Gaussian rules. The example in Section 3 displays a significance of effective condition number for numerical partial differential equations (PDE). Furthermore, a close linkage of Section 3 and this section is explored in Section 3.3. The techniques of weights are essential in science and numerical
methods, due to variety and complication of the problems studied. Therefore, the effective condition number in this paper is a development of our previous study (see [26,27]).
3. Trefftz method for Laplace’s equations In this section, the effective condition number is applied to the weighted linear least squares problem from numerical PDE. We take the Trefftz method [29] for solving Laplace’s equation with singularities for example. 3.1. Numerical algorithms and computed results In this section, we consider the crack problem (see Fig. 1)
Du ¼ 0 in S, u¼0
on OD,
u ¼ 1 on AB [ BC , un ¼ 0
u ¼ y on CD,
on OA,
ð3:1Þ
where S ¼ fðx,yÞ91 o x o1,0 o yo 1g, and un ¼ @u=@n is the outward normal derivative to @S. To solve (3.1), we may use the collocation Trefftz method (CTM) involving integration approximation (see [29,30,33]). When the Dirichlet boundary conditions involve noise, reader may refer the numerical examples in Lu et al. [32]. The exact solutions can be expanded as 1 X 1 u¼ di r i þ 1=2 cos i þ y in S, ð3:2Þ 2 i¼0 where di are the true coefficients. Choose the admissible solutions as L X 1 Di r i þ 1=2 cos i þ y in S, ð3:3Þ uL ¼ 2 i¼0 where Di are the unknown coefficients to be sought. Since the expansions (3.2) satisfy the Laplace’s equation and the boundary conditions at y ¼0 already, the coefficients Di should be chosen to satisfy the rest of the boundary conditions in (3.1), u9AB[BC ¼ 1,
u9CD ¼ y,
ð3:4Þ
as best as possible, where AB, BC and CD are shown in Fig. 1. Hence, the linear least squares method (LSM) may be used to seek the coefficients Di. Denote Z ½u,v ¼ uv d‘, ð3:5Þ Gn
where Gn ¼ AB [ BC [ CD, and VL denotes the set of finite dimensional functions in (3.3). Then, we may seek uL A VL such that Z Z ½uL ,v ¼ v d‘ þ yv d‘, 8v A VL : ð3:6Þ AB[BC
CD
Fig. 1. The crack problem.
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
Denote the energy Z Z ðv1Þ2 d‘ þ IðvÞ ¼ AB[BC
ðvyÞ2 d‘:
ð3:7Þ
CD
The solution (i.e., the coefficients Di) of (3.6) can also be expressed by: To seek uL A VL such that IðuL Þ ¼ minIðvÞ:
ð3:8Þ
v A VL
When the integrals in (3.7) involve numerical approximation, the energy is denoted by Z Z c c ^ ¼ IðvÞ ðv1Þ2 d‘ þ ðvyÞ2 d‘, ð3:9Þ AB[BC
CD
R R where b is an approximation of integral by some quadrature rule, such as the central or the Gaussian rule. Then the solution u~ L A VL is sought by ^ ^ u~ L Þ ¼ minIðvÞ: Ið
ð3:10Þ
v A VL
On the other hand, one may establish the boundary collocation equations directly from the boundary conditions (3.4), L X 1 i þ 1=2 Di r j cos i þ y ¼ 1, ðrj , yj Þ A AB [ BC , ð3:11Þ 2 j i¼0 L X i¼0
i þ 1=2 Di r j
1 cos i þ y ¼ rj sinðyj Þ, 2 j
57
nodes ðrj , yj Þ of small pffiffiffi subsections. Eqs. (3.11) and (3.12) multiplied by the weight h lead to L pffiffiffi X pffiffiffi 1 i þ 1=2 h Di r j cos iþ yj ¼ h, ðrj , yj Þ A AB [ BC , ð3:13Þ 2 i¼0 L pffiffiffi pffiffiffi X 1 i þ 1=2 Di r j cos iþ yj ¼ hrj sinðyj Þ, h 2 i¼0
ð3:12Þ
where ðr, yÞ are the polar coordinates with the origin (0, 0). Let ðrj , yj Þ denote the integration nodes. For the central rule, we choose uniform sections with the mesh-spacing h and the middle
ð3:14Þ
which are equivalent to (3.10). Such an equivalence is also valid for the Gaussian rule with suitable weights, see [33]. In computation, we choose the number m of nodes ðrj , yj Þ to be larger than (Lþ 1). Eqs. (3.13) and (3.14) lead to Fx ¼ b,
ð3:15Þ
where x A RL þ 1 is the unknown vector consisting of coefficients Di ði ¼ 0,1, . . . ,LÞ, b A Rm , and the stiffness matrix, F A RmðL þ 1Þ . The algorithms of (3.13) and (3.14) are called the collocation Trefftz methods (CTM) in [29,33], or the boundary approximation method in [25,30]. Once coefficients Di (i.e., x) are obtained, the errors on AB [ BC [ CD Z 1=2 Z JeJB ¼ JuuL JB ¼ ð1uL Þ2 d‘ þ ðyuL Þ2 d‘ ð3:16Þ AB[BC
ðrj , yj Þ A CD,
ðrj , yj Þ A CD,
CD
pffiffiffiffiffiffiffiffiffi are computable, where the boundary norm JvJB ¼ ½v,v. We will compute JeJB and the maximal errors on Gn , denoted by JeJ1, Gn ¼ maxGn 9e9. Note that Eq. (3.15) (i.e., (3.13) and (3.14)) is the weighted linear least squares problem of (3.11) and (3.12). We simply use (3.15) for simplicity in discussion.
Table 1 The error norms, condition numbers and errors of leading coefficients from the CTM for the crack problem by the Gaussian rule of six nodes rule as m ¼ 42, where 0n denotes the error less than the computer rounding errors in double precision. L
JeJB
JeJ1, Gn
DD0
D
DD1
D
DD2
D
DD3
D
Cond
Cond_eff
Cond_EE
0.176( 6) 0.908( 9) 0.795( 11) 0.667( 13) 0.566( 15) 0n 0n 0n 0n 0n 0n 0n 0n 0n
0.689( 8) 0.268( 10) 0.243( 12) 0.109( 14) 0.871( 15) 0.871( 15) 0.871( 15) 0.871( 15) 0.871( 15) 0.871( 15) 0.871( 15) 0.871( 15) 0.871( 15) 0.871( 15)
0.546( 6) 0.650( 8) 0.788( 10) 0.970( 12) 0.141( 13) 0.200( 14) 0.206( 14) 0.206( 14) 0.206( 14) 0.206( 14) 0.206( 14) 0.206( 14) 0.206( 14) 0.206( 14)
0.791( 3) 0.161( 4) 0.298( 6) 0.521( 8) 0.870( 10) 0.140( 11) 0.202( 13) 0.155( 13) 0.194( 13) 0.202( 13) 0.202( 13) 0.202( 13) 0.202( 13) 0.202( 13)
15.4 52.5 186 674 0.249(4) 0.927(4) 0.348(5) 0.132(6) 0.502(6) 0.192(7) 0.737(7) 0.286(8) 0.110(9) 0.429(9)
1.22 1.23 1.24 1.24 1.24 1.24 1.24 1.24 1.25 1.25 1.25 1.25 1.25 1.25
1.34 1.41 1.45 1.48 1.50 1.52 1.53 1.54 1.55 1.56 1.56 1.57 1.57 1.58
0
10 14 18 22 26 30 34 38 42 46 50 54 58 62
0.769( 4) 0.549( 5) 0.438( 6) 0.372( 7) 0.330( 8) 0.304( 9) 0.286( 10) 0.274( 11) 0.266( 12) 0.263( 13) 0.265( 14) 0.479( 15) 0.374( 15) 0.379( 15)
0.153( 3) 0.115( 4) 0.967( 6) 0.854( 7) 0.768( 8) 0.743( 9) 0.718( 10) 0.706( 11) 0.704( 12) 0.715( 13) 0.755( 14) 0.178( 14) 0.155( 14) 0.211( 14)
1
2
3
Table 2 The error norms, condition numbers and errors of leading coefficients from the CTM for the crack problem by the Gaussian rule of six nodes rule as L ¼54, where 0n denotes the error less than the computer rounding errors in double precision. m
JeJB
JeJ1, Gn
DD0
D
DD1
D
DD2
D
DD3
D
Cond
Cond_eff
Cond_EE
0.566( 15) 0.754( 15) 0.754( 15) 0.566( 15) 0n 0.377( 15) 0.132( 14) 0.377( 15) 0.189( 15)
0.654( 15) 0.654( 15) 0.436( 15) 0.871( 15) 0.871( 15) 0.654( 15) 0.436( 15) 0.436( 15) 0n
0.117( 14) 0.294( 14) 0.441( 14) 0.117( 14) 0.206( 14) 0.132( 14) 0.176( 14) 0.250( 14) 0.147( 14)
0.168( 13) 0.129( 13) 0.244( 13) 0.129( 14) 0.194( 14) 0.200( 14) 0.129( 15) 0.777( 15) 0.518( 15)
0.286(8) 0.286(8) 0.286(8) 0.286(8) 0.286(8) 0.286(8) 0.286(8) 0.286(8) 0.286(8)
1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25
1.57 1.57 1.57 1.57 1.57 1.57 1.57 1.57 1.57
0
18 24 30 36 42 48 54 60 66
0.676( 15) 0.105( 14) 0.925( 15) 0.750( 15) 0.479( 15) 0.640( 15) 0.152( 14) 0.651( 15) 0.529( 15)
0.102( 13) 0.311( 14) 0.244( 14) 0.244( 14) 0.178( 14) 0.200( 14) 0.333( 14) 0.200( 14) 0.178( 14)
1
2
3
58
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
In our computation, the Gaussian rule with six nodes is chosen for computing the integrals in (3.7), and m denotes the number of nodes along AB. Then m is multiples of six, and m ¼ 4m. The QR factorization is used for solving (3.15). The errors, condition numbers and effective condition numbers are listed in Tables 1 and 2. Note that all computations in this paper are carried out by Fortran programs in double precision. The highly accurate solution is found at L¼ 54 and m ¼ 42, with JeJB ¼ 0:479ð15Þ,
Cond ¼ 0:284ð8Þ,
Cond_eff ¼ 1:25,
DD0
¼ 0n ,
ð3:17Þ
Cond_EE ¼ 1:57,
ð3:18Þ
JeJ1, Gn ¼ 0:178ð14Þ,
D0
Table 3 The leading coefficients Di from the CTM for the crack problem by the Gaussian rule with six nodes as L¼ 54 and m ¼ 42 along AB. i
Di
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 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
1.17769302091071149 0.254824842629162551 0.472464287866036919( 1) 0.268074643613298927( 1) 0.207096003479169876( 2) 0.285158772208820362( 3) 0.183013497194766252( 3) 0.533051704416331353( 3) 0.157452569473614989( 4) 0.537931615260617409( 6) 0.282586859356730478( 5) 0.172189962001964248( 4) 0.991762428046167717( 7) 0.182802186655744245( 7) 0.753824634122721559( 7) 0.681235665190457081( 6) 0.148739422491196200( 8) 0.333015969702268427( 9) 0.248256678147026780( 8) 0.300174031113286268( 7) 0.781012645562570253( 11) 0.104558489837163971( 10) 0.933291671159730157( 10) 0.141328444671448315( 8) 0.311705144580195162( 12) 0.338335131879225414( 12) 0.382332681328341618( 11) 0.696105306905609050( 10) 0.248911941008673783( 14) 0.124843211054098715( 13) 0.166517138300667560( 12) 0.354249175259574404( 11) 0.179980257225143477( 15) 0.490941653852575528( 15) 0.759286687343576104( 14) 0.184733137546747423( 12) 0.851604256940787594( 17) 0.204119729589908422( 16) 0.360262461151616336( 15) 0.977970794533835831( 14) 0.181405679095255387( 17) 0.735513396111446796( 18) 0.181360402136695322( 16) 0.509563344013018680( 15) 0.436358147836334786( 18) 0.196117522501444700( 19) 0.101927109194609567( 17) 0.232534810678914514( 16) 0.623638419504983828( 19) 0.808600297512558928( 20) 0.582074741343905307( 19) 0.654386595624466984( 18) 0.381783351320010671( 20) 0.568217460684078295( 21) 0.207475181771925468( 20)
Table 4 The singular values si and the coefficients bi for the matrix resulting from the CTM solution in Table 3, where Cond¼ 0.286(8), Cond_eff ¼ 1.25 and Cond_EE ¼1.57. i
si
bi
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 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
0.199(8) 0.104(8) 0.101(8) 0.531(7) 0.351(6) 0.246(6) 0.186(6) 0.131(6) 0.178(5) 0.131(5) 0.995(4) 0.735(4) 0.166(4) 0.127(4) 0.983(3) 0.752(3) 0.252(3) 0.198(3) 0.158(3) 0.125(3) 0.567(2) 0.461(2) 0.380(2) 0.311(2) 0.176(2) 0.148(2) 0.126(2) 0.106(2) 0.707(1) 0.608(1) 0.533(1) 0.461(1) 0.348(1) 0.306(1) 0.276(1) 0.244(1) 0.202(1) 0.181(1) 0.168(1) 0.152(1) 0.135(1) 0.123(1) 0.117(1) 0.108(1) 0.101(1) 0.936 0.911 0.858 0.837 0.785 0.775 0.753 0.749 0.715 0.702
0.657( 1) 0.195( 1) 0.163 0.407( 1) 0.380( 1) 0.262( 1) 0.990( 1) 0.554( 1) 0.320( 1) 0.271( 1) 0.864( 1) 0.565( 1) 0.300( 1) 0.278( 1) 0.836( 1) 0.570( 1) 0.302( 1) 0.294( 1) 0.871( 1) 0.588( 1) 0.323( 1) 0.322( 1) 0.962( 1) 0.624( 1) 0.362( 1) 0.364( 1) 0.112 0.676( 1) 0.422( 1) 0.424( 1) 0.135 0.743( 1) 0.514( 1) 0.505( 1) 0.172 0.826( 1) 0.651( 1) 0.606( 1) 0.226 0.907( 1) 0.820( 1) 0.700( 1) 0.292 0.927( 1) 0.931( 1) 0.747( 1) 0.350 0.806( 1) 0.836( 1) 0.730( 1) 0.357 0.495( 1) 0.498( 1) 0.100 0.672
where Cond_EE ¼JbJ=bn , and 0n denotes the error less than the computer rounding errors in double precision. Since the maximal errors of harmonic functions occur only on the boundary Gn , we conclude that JeJ1,S r 0:178ð14Þ. The coefficients Di are listed in Table 3, and the singular values si and the expansion coefficients bi in Table 4. The leading coefficient D0 ¼ 1:17769302091071149 in Table 3 is exact in the sense that its error is less than the computer rounding error t ¼ 0:5ð17Þ, and D1 ,D2 and D3 have 16, 15 and 15 significant digits, respectively.2 Table 5 lists the values of smax , smin , etc. for different L used in Table 2. This very accurate
2 Such significant digits are found by comparisons with the more accurate solution at L¼100 and by Mathematics using 100 working digits.
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
Table 5 The maximal and the minimal singular values and condition numbers in Table 1 as m ¼ 42. L
smax
smin
JbJ
bn
Cond
Cond_eff
Cond_EE
10 14 18 22 26 30 34 38 42 46 50 54 58 62
11.2 37.2 131 475 0.175(4) 0.651(4) 0.245(5) 0.926(5) 0.352(6) 0.135(7) 0.517(7) 0.199(8) 0.768(8) 0.297(9)
0.716 0.710 0.704 0.713 0.704 0.703 0.688 0.702 0.702 0.702 0.702 0.702 0.701 0.701
1.05 1.05 1.05 1.05 1.05 1.05 1.05 1.05 1.05 1.05 1.05 1.05 1.05 1.05
0.782 0.749 0.728 0.713 0.702 0.894 0.688 0.683 0.679 0.676 0.674 0.672 0.670 0.669
15.4 52.5 186 674 0.249(4) 0.927(4) 0.348(5) 0.132(6) 0.502(6) 0.192(7) 0.737(7) 0.284(8) 0.110(9) 0.429(9)
1.22 1.23 1.24 1.24 1.24 1.24 1.24 1.24 1.25 1.25 1.25 1.25 1.25 1.25
1.34 1.41 1.45 1.48 1.50 1.52 1.53 1.54 1.55 1.56 1.56 1.57 1.57 1.58
59
Suppose that the integration rules are chosen so that to satisfy the following equivalence relations: JvJ 0, Gn OðJvJ0, Gn Þ,
ð3:27Þ
i.e., there exist two constants c0 ð 4 0Þ and C1 such that c0 JvJ0, Gn rJvJ 0, Gn r C1 JvJ0, Gn :
ð3:28Þ
An analysis in [33] shows that the integration rules to satisfy (3.28) are not severe, even by the simplest central rule. In the following, the constants c0 ð 4 0Þ,c1 ð 40Þ,C and C1 are independent of L, but their values may be different in different places. We have the following lemma. Lemma 3.1. Suppose that for v A VL there exists a positive constant p 4 0 such that JvJ1, Gn r CLp JvJ0, Gn :
ð3:29Þ
Under (3.27), there exists the lower bound,
smin ðFÞ Z c0 Lp=2 : solution can be used to test other methods, in particular, those for the boundary integral equation (BIE) of first kind.3 3.2. Bounds of effective condition number and condition number Below, we provide the stability analysis for the crack problem in Section 3.1, mainly based on the effective condition number. Eqs. (3.13) and (3.14) are rewritten as (3.15). We will derive an upper bound of (1.2): Cond_eff ¼
smax , smin
Ax ¼ d,
lmin ðAÞ ¼ min xa0
ðAx,xÞ : ðx,xÞ
ð3:23Þ
ð3:24Þ
Z Z c ^ ¼2 v2 d‘ 2IðvÞ ¼ 2 v2 d‘ ¼ 2JvJ20, Gn : ðAx,xÞ ¼ 2IðvÞ Gn
ð3:25Þ
Gn
Denote 2
Z c
ð3:33Þ
Combining (3.32) and (3.33) yields IðvÞ Zc0 Lp JvJ21,S :
ð3:34Þ
Denote the semi-disk with the radius r, Sr ¼ fðr, yÞ90 rr r r,0 r y r pg: Since Sr 9r ¼ 1 S, we have JvJ1,S ZJvJ1,Sr 9
r¼1
Z 9v91,Sr 9
r¼1
:
From the Green formula, ZZ Z 2 9v91,Sr ¼ rvrv ds ¼ vn v d‘, Sr
ð3:35Þ
ð3:36Þ
‘r
‘r ¼ fðr, yÞ9r ¼ r,0 r y r pg: D2i ,
i¼0
JvJ 0, Gn ¼
JvJ21,S rCJvJ21=2, Gn :
ð3:32Þ
where ‘r is the semi-circle
In (3.23), the notations are L X
Hence we obtain from (3.27) and (3.31),
On the other hand, since Dv ¼ 0 for v A VL , we have from Oden and Reddy [34, p.189],
where lmax and lmin are the maximal and the minimal eigenvalues of A, respectively, defined by
ðx,xÞ ¼ JxJ2 ¼
In (3.31), the semi-norms in the Sobolev space is defined by ( )1=2 Z Z ðvðPÞvðQ ÞÞ2 2 d‘ðPÞ d‘ðQ Þ : JvJ1=2, G ¼ JvJ0, G þ JPQ J2 G G
ð3:20Þ
where the matrix A ¼ FT F A RðL þ 1ÞðL þ 1Þ is symmetric and positive definite, and d ¼ FT b. Then there exist the relations pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi smax ðFÞ ¼ lmax ðAÞ, smin ðFÞ ¼ lmin ðAÞ, ð3:22Þ
ðAx,xÞ , ðx,xÞ
ð3:31Þ
ðAx,xÞ ¼ 2JvJ 0, Gn Zc0 JvJ20, Gn Z c1 Lp JvJ21=2, Gn :
ð3:21Þ
xa0
JvJ1=2, Gn rCðJvJ1, Gn JvJ0, Gn Þ1=2 r C1 Lp=2 JvJ0, Gn :
ð3:19Þ
where smax and smin are the maximal and the minimal singular values of F, respectively. In this application, rankðFÞ ¼ n ¼ Lþ 1, smax ¼ s1 and smin ¼ sn 4 0. The normal equation of (3.15) is given by
lmax ðAÞ ¼ max
Proof. We have from (3.29) and Babuska and Aziz [5, p. 21],
2
JbJ , smin JxJ
and that of Cond ¼
ð3:30Þ
v2 d‘:
ð3:26Þ
Gn
3 Motz’s problem in [33] with mixed type of Dirichlet and Neumann boundary conditions cannot be used as the test model for the BIE of first kind, since the latter results from only the Dirichlet boundary condition.
By calculus, from the orthogonality of cosðiþ 12Þy we obtain from (3.3) ) Z p (X Z L 1 i1=2 1 vn v d‘ ¼ Di iþ r cos i þ y 2 2 ‘r 0 i¼0 ( ) L X 1 Di ri þ 1=2 cos iþ y r dy 2 i¼0 Z L X 1 2i þ 1 p 1 ¼ D2i i þ r cos2 iþ y dy 2 2 0 i¼0 L pX 1 2i þ 1 ¼ D2i i þ r : ð3:37Þ 2 2i¼0
60
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
Then when r ¼ 1, Z L L pX 1 pX Z vn v d‘ ¼ D2i i þ D2 : 2 2i¼0 4i¼0 i ‘r
ð3:38Þ
Moreover, since S Sr 9r ¼ pffiffi2 , we have from (3.37)
Z L pX 1 2i þ 1
2 2 vn v d‘ ¼ D2i i þ r 9v91,S r 9v91,Spffi ¼
2
2 2i¼0 ‘pffi2
Combining (3.34)–(3.36) and (3.38) yields IðvÞ Z c0 Lp
L X
r D2i :
ð3:39Þ
i¼0
^ IðvÞ IðvÞ Z c0 Z c1 Lp , ðx,xÞ x a 0 ðx,xÞ
lmin ðAÞ ¼ min
L 1 pffiffiffi 2L þ 1 X ð 2Þ Lþ D2i : 2 2 i¼0
p
L pffiffiffi X D2i : IðvÞ rCLð 2Þ2L
ð3:53Þ
smin ¼ smin ðFÞ Zc1 Lp=2 :
Hence from (3.28) and (3.53), the maximal eigenvalue lmax ðAÞ has the following bound: ð3:41Þ
This is the desired result (3.30). This completes the proof of Lemma 3.1. & Theorem 3.1. Let (3.27) and (3.29) hold. Then for (3.15) there exists the following bound of the effective condition number :
ð3:54Þ
i¼0
ð3:40Þ
and from (3.22)
Cond_eff r CL
ð3:52Þ
Combining (3.51) and (3.52) yields
Hence we have from (3.27) and (3.39)
p=2
pffiffi
r¼ 2
ð3:42Þ
lmax ðAÞ ¼ max xa0
^ IðvÞ JxJ2
r C max xa0
IðvÞ JxJ2
pffiffiffi r C1 Lð 2Þ2L :
Then we have from (3.22), pffiffiffi pffiffiffi
smax ¼ smax ðAÞ r C Lð 2ÞL :
ð3:55Þ
ð3:56Þ
This is the desired result (3.47), and this completes the proof of Lemma 3.2. & Based on Lemmas 3.1 and 3.2, we have the following theorem.
Proof. The vector b in (3.15) has the components pffiffiffi pffiffiffi T b ¼ f. . . , ai h, . . . , ai yi h, . . .g,
ð3:43Þ
where 0 o yi o1 and ai is the integration weights with ai ¼ Oð1Þ. Since h r C1=m, where m is the dimension of b, we have sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X X pffiffiffi ffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffi JbJ r ðai hÞ2 r a2i h r max9ai 9 mh r C: ð3:44Þ i
i
i
Also since the lower bounds of the true coefficients are known in Table 3, and the vector x satisfies vffiffiffiffiffiffiffiffiffiffiffiffiffiffi u L uX D2i ZD0 Z 1: ð3:45Þ JxJ ¼ t i¼0
Theorem 3.2. Let (3.27) and (3.29) hold. Then for (3.15) there exists the bound for the traditional condition number pffiffiffi smax Cond ¼ r CLðp þ 1Þ=2 ð 2ÞL : ð3:57Þ
smin
For the sectorial S, p ¼1 is given in [25, p. 112]. Hence we have the following corollary from Lemmas 3.1 and 3.2. Corollary 3.1. Let (3.27) and (3.29) hold. Also assume p ¼1. Then for (3.15), there exist the bounds
smin Z c0 L1=2 ,
ð3:58Þ
pffiffiffi pffiffiffi smax rC Lð 2ÞL :
ð3:59Þ
Hence we have from Lemma 3.1, (3.44) and (3.45), Cond_eff ¼
JbJ rCLp=2 : smin JxJ
From Theorems 3.1 and 3.2, we have the following corollary. ð3:46Þ
This is the desired result (3.42), and completes the proof of Theorem 3.1. & Next, we will estimate the upper bound of smax , to have the following lemma. Lemma 3.2. Let (3.27) hold. Then for (3.15) there exists the upper bound, pffiffiffi pffiffiffi smax ¼ smax ðFÞ rC Lð 2ÞL : ð3:47Þ Proof. From the embedding theorem in Ciarlet [13], IðvÞ ¼ JvJ20, Gn
r CJvJ21,S :
ð3:48Þ
From (3.27) and (3.48) we have 2
^ ¼ 2JvJ n r2CIðvÞ ¼ CJvJ2 n rC1 JvJ2 : ðAx,xÞ ¼ 2IðvÞ 0, G 1,S 0, G
ð3:49Þ
Since v9y ¼ 041 o x o 0 ¼ 0 for v A VL , there exists the bound from the Poincare inequality [13], JvJ1,S rC9v91,S ,
ð3:50Þ
where 9v91,S is the semi-norm of v on S. Hence we have from (3.48) and (3.50), 2
IðvÞ r C9v91,S :
ð3:51Þ
Corollary 3.2. Let (3.27) and (3.29) hold. Also assume p ¼1. Then for (3.15) there exist the bounds Cond_eff r CL1=2 ,
ð3:60Þ
pffiffiffi Condr CLð 2ÞL :
ð3:61Þ
Eqs. (3.59) and display that when L-1, Cond_eff and pffiffi(3.60) ffi Cond grow as Oð LÞ and exponentially, respectively. Now, let us examine again the data in Table 1. We can see JeJB ¼ Oðð0:56ÞL Þ,
JeJ1, Gn ¼ Oðð0:56ÞL Þ,
ð3:62Þ
Cond_eff ¼ Oð1Þ,
pffiffiffi Cond ¼ Oðð1:4ÞL Þ ¼ Oðð 2ÞL Þ:
ð3:63Þ
Eqs. (3.63) are consistent with Corollary 3.2. From Table 5, we can see pffiffiffi smin 0:7, smax ¼ ð1:4ÞL Oðð 2ÞL Þ, ð3:64Þ JbJ ¼ 1:0,
bn 0:67, Cond_EE 1:57:
ð3:65Þ
Eqs. (3.64) also agree with Corollary 3.1. Besides, the error analysis in [33] proves that the CTM has the exponentially convergence rates, which is also validated by (3.62). From this section, the CTM has an excellent stability, based on Cond_eff. Hence, the CTM becomes a very efficient boundary method, to
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
compete with other boundary methods, such as the boundary element method (BEM) and the BIE method, see [28,29]. 3.3. Linkage to Section 2 Finally, let us link the numerical computation in this section to the theoretical estimation (2.25) in Section 2. Take the leading coefficient D0 in Table 3 for example, which only has an error less than the unit of rounding error. Such an extreme accuracy can be explained by the effective condition number, but not by the traditional condition number. Denote by g the rounding error of the computer used. Since the double precision is used (say, g ¼ 1017 ), coefficient D0 has the relative error 9DD0 9 r g ¼ 1017 : 9D0 9
ð3:66Þ
For the coefficients Di in Table 3, we find vffiffiffiffiffiffiffiffiffiffiffiffiffiffi u L uX D2i ¼ mJxJ, 9D0 9 ¼ mt where L¼54 and the factor m ¼ 1:2. Then we have qP ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi L 2 9DD0 9 9DD0 9 JDxJ i ¼ 0 DDi ¼ r ¼ : mJxJ mJxJ mJxJ 9D0 9
ð3:67Þ
ð3:68Þ
ð3:70Þ
Since the errors JDAJ=JAJ from the solution method will satisfy JDAJ 1 1 5 ¼ , JAJ Cond 0:286ð8Þ
ð3:71Þ
in double precision, we have d 51. Therefore, we have from (3.69) JDxJ JDbJ r 1:25 : JxJ JbJ
ð3:72Þ
Combining (3.67) and (3.71) gives 9DD0 9 1:25 JDbJ 1:25 JDbJ JDbJ ¼ ¼ 1:02 : r JbJ 1:2 JbJ JbJ m 9D0 9
ð3:73Þ
Since the relative error of b results from rounding error, we have 9DD0 9 ð r 1017 Þ r 1:02 g ¼ 1:02ð17Þ: 9D0 9
ð3:74Þ
How perfect coincidence between the numerical computation and the theoretical estimation! Eq. (3.74) strongly verifies the estimation (2.25), which indicates a key role of effective condition number in (2.16), (2.18) and (2.19). This explores a linkage of this section to Section 2. For ðF þ DFÞðxþ DxÞ ¼ b þ Db, we may reform it as n
Fðx þ DxÞ ¼ b þ Db ,
ð3:75Þ
where
Dbn ¼ DbDFðx þ DxÞ:
ð3:78Þ
we also obtain (3.74). This again exhibits the significance of effective condition number for stability analysis. On the other hand, for the solution x (i.e., coefficients Di) in Table 3, from (1.1) we have the error bounds JDxJ JDbJ JDbJ rCond ¼ 0:286ð8Þ JxJ JbJ JbJ r0:286ð8Þ g ¼ 0:286ð9Þ:
ð3:79Þ
Combining (3.67) and (3.78) yields 9DD0 9 JDbJ ð r 1017 Þ rCond mJbJ 9D0 9 ¼ Cond
g 0:286ð8Þ 1017 ¼ m m
0:286ð9Þ ¼ 2:38ð10Þ: 1:2
ð3:80Þ
In (3.80) there exists a big gap (about a factor 107) between the actual error and the estimation of D0.4 This fact illustrates inefficacy of Cond. Similarly for ðF þ DFÞðxþ DxÞ ¼ b þ Db, we have JDxJ JDb J rCond , JxJ JbJ
ð3:81Þ
and also reach (3.80), which is misleading, indeed.
where JDAJ : JAJ
n
JDb J r g ¼ 1017 , JbJ
n
From Table 1 at L¼54, Cond_eff¼1.25 and Cond¼0.286 (8). The relative errors of x are given from (2.25), JDxJ 1 JDbJ 1 JDbJ r ¼ , ð3:69Þ Cond_eff 1:25 JxJ 1d JbJ 1d JbJ
d ¼ Cond
Since the minimal bound is also given by the rounding error
¼
i¼0
61
ð3:76Þ
4. Concluding remarks To close this paper, let us address the new contributions in this paper. (1) For weighted linear least squares problem with both full rank and rank-deficient cases, we develop the effective condition number defined in (2.15), and derive in Theorem 2.1 the new error bounds from the effective condition number. These are the new computational formulas, compared to the existing literature including our previous study in [26,27]. (2) In Section 3.3, the extremely accurate leading coefficient D0 is explained by the very small effective condition number. In contrast, the huge condition number is misleading for stability analysis. Moreover, the close linkage between Sections 2 and 3 is explored, and Eq. (2.25) is verified perfectly by numerical computation of Motz’s problem by the CTM. This indicates a key role of effective condition number in (2.16), (2.18) and (2.19). (3) In [29,28], the error analysis is explored for Trefftz and collocation Trefftz methods. Although the numerical results of Cond_eff are provided, but no strict stability analysis was given wherein. This paper is devoted to the new stability analysis based on Cond_eff, to show an advantage of the collocation Trefftz method (CTM) over other boundary methods. Moreover, the analysis and numerical examples in this paper display that the effective condition number is particularly significant to the CTM and the spectral method, where the huge condition number is misleading. (4) By Axelsson’s work [1–3], the effective condition number is beneficial to studying the iteration solution methods. From [26,27] and this paper, the effective condition number is important to numerical PDE, and it may become a new trend of stability analysis.
From (1.1) we have n
JDxJ JDb J r Cond_eff : JxJ JbJ
ð3:77Þ
4 Eq. (3.80) indicates that D0 has at most 10 significant digits. However, the numerical D0 in Table 3 (also see many numerical values of D0 in Table 1) has 17 significant digits!
62
Y.-m. Wei et al. / Engineering Analysis with Boundary Elements 36 (2012) 53–62
References [1] Axelsson O. Iterative solution methods. Cambridge: Cambridge University Press; 1994. [2] Axelsson O, Kaporin I. On the sublinear and superlinear rate of convergence of conjugate gradient methods. Numer Algorithms 2000;25:1–22. [3] Axelsson O, Kaporin I. Error norm estimation and stopping criteria in preconditioned conjugate gradient iterations. Numer Linear Algebra Appl 2001;8:265–86. [4] Atkinson KE. An introduction to numerical analysis. New York: John Wiley & Sons; 1989. [5] Babuska I, Aziz AK. Survey lectures in the mathematical foundations of the finite elements. In: Aziz AK, editor. The mathematical foundation of the finite element methods with applications to partial differential equations. New York, London: Academic Press; 1973. p. 3–369. [6] Barlow JL, Handy SL. The direct solution of weighted and equality constrained least squares problems. SIAM J Sci Stat Comput 1988;9:704–16. [7] Banoczi JM, Chiu M, Cho GE, Ipsen CF. The lack of influence of right-hand side on accuracy of linear system solution. SIAM J Sci Comput 1998;20: 203–27. ˚ Numerical methods for least squares problems. Philadelphia: SIAM; ¨ [8] Bjorck A. 1996. [9] Bobrovnikova EY, Vavasis SA. Accurate solution of weighted least squares by iterative methods. SIAM J Matrix Anal Appl 2001;22:1153–74. [10] Campbell S, Meyer C. Generalized inverses of linear transformations. New York: Pitman, Dover; 1979; Campbell S, Meyer C. Generalized Inverses of Linear Transformations. New York: Dover; 1991. [11] Chan TF, Foulser DE. Effectively well-conditioned linear systems. SIAM J Sci Stat Comput 1988;9:963–9. [12] Christiansen S, Hansen PC. The effective condition number applied to error analysis of certain boundary collocation methods. J Comput Appl Math 1994;54:15–36. [13] Ciarlet PG. Basic error estimates for elliptic problems. In: Ciarlet PG, Lions JL, editors. Finite element methods (Part I). North-Holland; 1991. p. 17–352. [14] Cucker F, Diao H, Wei Y. On mixed and componentwise condition numbers for Moore–Penrose inverse and linear least squares problems. Math Comput 2007;78(258):947–63. [15] Golub GH, van Loan CF. Matrix computation.3rd ed. Baltimore, London: The John Hopkins University Press; 1996. [16] Gulliksson M. Iterative refinement for constrained and weighted linear least squares. BIT 1994;34:239–53. [17] Gulliksson M. On the modified Gram–Schmidt algorithm for weighted and constrained linear least squares problems. BIT 1995;35:453–68. [18] Gulliksson M. Backward error analysis for the constrained and weighted linear least squares problem when using the weighted QR decomposition. SIAM J Matrix Anal Appl 1995;16:675–87. [19] Gulliksson ME, Wedin PA, Wei Y. Perturbation identities for regularized Tikhonov inverses and weighted pseudoinverses. BIT 2000;40:513–23. [20] Gulliksson ME, Jin X, Wei Y. Perturbation bound for constrained and weighted least squares problem. Linear Algebra Appl 2000;349:221–32.
[21] Higham NJ. Accuracy and stability of numerical algorithm.2nd ed. Philadelphia: SIAM; 2002. [22] Hough PD, Vavasis SA. Complete orthogonal decomposition for weighted least squares. SIAM J Matrix Anal Appl 1997;18:369–72. [23] Kenney CS, Laub AJ, Reese MS. Statistical condition estimation for linear least squares. SIAM J Matrix Anal Appl 1998;19:906–23. [24] Lawson CL, Hanson RJ. Solving least squares problems. Englewood Cliffs, NJ: Prentice-Hall, Inc.; 1974. Reprint, SIAM, Philadelphia, 1995. [25] Li ZC. Combined methods for elliptic equations with singularities, interfaces and infinities. Dordrecht, Boston: Kluwer Academic Publishers; 1998. [26] Li ZC, Chien CS, Huang HT. Effective condition number for finite difference method. J Comput Appl Math 2007;198:208–35. [27] Li ZC, Huang HT. Effective condition number for numerical partial differential equations. Numer Linear Algebra Appl 2008;15:575–94. [28] Li ZC, Lu TT, Huang HT, Cheng AH-D. Trefftz, collocation and other boundary methods, a comparison. Numer Meth PDE 2007;23:93–144. [29] Li ZC, Lu TT, Hu HY, Cheng AH-D. Trefftz and collocation methods. Southampton, Boston: WIT Publishers; 2008. [30] Li ZC, Mathon R, Sermer P. Boundary methods for solving elliptic problem with singularities and interfaces. SIAM J Numer Anal 1987;24:487–98. [31] Lin L, Lu TT, Wei Y. On level-2 condition number for the weighted Moore– Penrose inverse. Comput Math Appl 2008;55:788–800. [32] Lu TT, Chang CM, Huang HT, Li ZC. Stability analysis of Trefftz methods for the stick-slip problem. Eng Anal Boundary Elem 2009;33:474–84. [33] Lu TT, Hu HY, Li ZC. Highly accurate solutions of Motz’s and the cracked beam problems. Eng Anal Boundary Elem 2004;28:1387–403. [34] Oden JT, Reddy JN. An introduction to the mathematical theory of finite elements. New York: John Wiley & Sons; 1976. [35] Rice JR. Matrix computations and mathematical software. New York: McGraw-Hill Book Company; 1981. [36] Rump S. Ill-conditionedness need not be the componentwise near to illposedness for least squares problems. BIT 1999;39:143–51. [37] van Loan CF. Generalizing the singular value decomposition. SIAM J Numer Anal 1976;13:76–83. [38] van Loan CF. On the method of weighting for equality-constrained least squares problems. SIAM J Numer Anal 1985;22:851–64. [39] Vavasis SA. Stable numerical algorithms for equilibrium systems. SIAM J Matrix Anal Appl 1994;15:1108–31. [40] Wang G, Wei Y, Qiao S. Generalized inverses: theory and computations. Beijing: Science Press; 2004. [41] Wei M, De Pierro AR. Upper perturbation bounds of weighted projections, weighted and constrained least squares problem. SIAM J Matrix Anal Appl 2000;21:931–51. [42] Wei M. Supremum and stability of weighted pseudoinverses and weighted least squares problems: analysis and computations. New York: Nova Science Publishers; 2001. [43] Wilkinson JH. The algebraic eigenvalue problems. Oxford University Press; 1965. [44] Xiang H, Wei Y. Structured mixed and componentwise condition numbers of some structured matrices. J Comput Appl Math 2007;202:217–29. [45] Xu W, Wei Y, Qiao S. Condition numbers for structured least squares problems. BIT 2006;46:203–25.