- Email: [email protected]
A spectral collocation technique based on integrated Chebyshev polynomials for biharmonic problems in irregular domains
Available online at www.sciencedirect.com
Applied Mathematical Modelling 33 (2009) 284–299 www.elsevier.com/locate/apm
A spectral collocation technique based on integrated Chebyshev polynomials for biharmonic problems in irregular domains N. Mai-Duy a,*, H. See b, T. Tran-Cong a a
Computational Engineering and Science Research Centre (CESRC), The University of Southern Queensland, Toowoomba, QLD 4350, Australia b School of Chemical and Biomolecular Engineering, The University of Sydney, NSW 2006, Australia Received 27 February 2007; received in revised form 2 November 2007; accepted 5 November 2007 Available online 17 November 2007
Abstract In this paper, an integral collocation approach based on Chebyshev polynomials for numerically solving biharmonic equations [N. Mai-Duy, R.I. Tanner, A spectral collocation method based on integrated Chebyshev polynomials for biharmonic boundary-value problems, J. Comput. Appl. Math. 201 (1) (2007) 30–47] is further developed for the case of irregularly shaped domains. The problem domain is embedded in a domain of regular shape, which facilitates the use of tensor product grids. Two relevant important issues, namely the description of the boundary of the domain on a tensor product grid and the imposition of double boundary conditions, are handled effectively by means of integration constants. Several schemes of the integral collocation formulation are proposed, and their performances are numerically investigated through the interpolation of a function and the solution of 1D and 2D biharmonic problems. Results obtained show that they yield spectral accuracy. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Integral collocation formulation; Biharmonic problems; Complex geometries; Fictitious domains; Chebyshev polynomials
1. Introduction Many engineering problems, such as the deformation of a thin plate and the motion of a fluid, are governed by the biharmonic equations – fourth-order partial differential equations (PDEs). Generally, problems involving high-order PDEs and complex geometries are more difficult to solve than those with second-order PDEs and regular geometries, respectively. Spectral collocation/pseudo-spectral methods (cf. [1–3]) are known to have the capability to provide an exponential rate of convergence as the grid is refined or the degree of the interpolation polynomial is increased. The drawback of these techniques is that they require a computational domain be square 1 6 x; y 6 1. For solving problems with complex geometries, domain decompositions and coordinate transformations have *
Corresponding author. Tel.: +61 7 4631 1324; fax: +61 7 4631 2526. E-mail address: [email protected] (N. Mai-Duy).
0307-904X/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2007.11.002
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
285
usually been employed to convert irregular domains into regular ones (e.g. [4]). Another way, which is a subject of the present study, is based on the use of fictitious domain. It is noted that fictitious-domain/domainembedding techniques can be traced back to the early 1950s. The basic idea behind these techniques is to extend domains of complicated shapes to those of simpler shapes from which the generation of meshes is simple and well-established efficient numerical solvers can be applied. Fictitious domains have been widely used in the context of finite elements, where the boundary conditions are implemented by means of Lagrange multipliers (e.g. [5] and references therein). In the present work, the problem domain of irregular shape is embedded in the reference square. This new domain is then discretized using a tensor product grid. Clearly, the grid points do not generally lie on the boundary of the actual domain. It is thus difficult to impose boundary conditions here with conventional differential approaches. In our earlier work [6] which deals with biharmonic problems defined in rectangular domains, it has been shown that the use of integration to construct the Chebyshev approximations provides an effective implementation of double boundary conditions. Unlike conventional differential formulations, the integral collocation formulation has the capability to generate extra expansion coefficients (integration constants). These additional unknown values can be utilized to incorporate normal derivative boundary conditions into the Chebyshev approximations. This paper is concerned with the development of the integral collocation formulation for the case of irregularly shaped domains. Three schemes of the integral collocation formulation are presented, and their performances are numerically investigated by considering several 1D and 2D test problems. The remainder of the paper is organized as follows. A brief review of differential and integral collocation formulations is given in Section 2. In this section, the integral collocation formulation is analyzed and several practical schemes of the formulation are presented. The proposed numerical procedure based on integral collocation schemes and fictitious domains is then described and verified through the interpolation of a function and the solution of 1D and 2D biharmonic equations in Sections 3–5, respectively. Section 6 gives some concluding remarks. 2. Collocation formulations 2.1. Differential formulation Consider a univariate function f ðxÞ defined in [1, 1]. This function can be represented by the Chebyshev interpolant of degree N as follows: N X f ðxÞ ¼ ak T k ðxÞ; ð1Þ k¼0 N fak gk¼0
N
where are the coefficients of expansion and fT k ðxÞgk¼0 are the Chebyshev polynomials of the first kind defined as T k ðxÞ ¼ cosðk arccosðxÞÞ. Expressions of derivatives of f are then obtained through differentiation. At the Gauss–Lobatto (G–L) points, N pi fxi gNi¼0 ¼ cos ; ð2Þ N i¼0 the values of derivatives of f are simply computed by c df b ð1Þ f^ ¼ D b f^ ; ¼D dx d 2 d f b ð2Þ f^ ¼ D b 2 f^ ; ¼D dx2 p d d f b ðpÞ f^ ¼ D b p f^ ; ¼D dxp
ð3Þ ð4Þ
ð5Þ
286
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
T where the symbol ^: is used to denote a vector/matrix that is associated with a grid line, f^ ¼ ðf0 ; f1 ; . . . ; fN Þ , k k k k c d f b ð:Þ are the differentiation matrices. The entries of ¼ ðddxfk0 ; ddxfk1 ; . . . ; ddxfkN ÞT with k ¼ f1; 2; . . . ; pg, and D dxk
b D b ð1Þ Þ are given by Dð iþj
b ij ¼ ci ð1Þ ; 0 6 i; j 6 N ; i 6¼ j; D cj xi xj xi b ii ¼ D ; 1 6 i 6 N 1; 2ð1 x2i Þ 2 b 00 ¼ D b NN ¼ 2N þ 1 ; D 6
ð6Þ ð7Þ ð8Þ
b can also be where c0 ¼ cN ¼ 2 and ci ¼ 1 for i ¼ f1; 2; . . . ; N 1g. It is noted that the diagonal entries of D obtained in the way that represents exactly the derivative of a constant b ii ¼ D
N X
b ij : D
ð9Þ
j¼0;j6¼i
For the case of smooth functions, the Chebyshev approximation scheme (1)–(5) is known to be very accurate (exponential accuracy) with the error being OðN a Þ in which a depends on the regularity of a function. It should be emphasized that there is a reduction in accuracy for the approximation of derivatives and this reduction is an increasing function of derivative order (cf. [3]). 2.2. Integral formulation The integral collocation formulation uses a truncated Chebyshev series of degree N to represent a derivative of an unknown function f, e.g. a derivative of order p, N N X dp f ðxÞ X ðpÞ ¼ a T ðxÞ ¼ ak I k ðxÞ: k k dxp k¼0 k¼0
ð10Þ
Expressions for lower-order derivatives and the function itself are then obtained through integration as N dp1 f ðxÞ X ðp1Þ ¼ ak I k ðxÞ þ c1 ; dxp1 k¼0
ð11Þ
N dp2 f ðxÞ X ðp2Þ ¼ ak I k ðxÞ þ c1 x þ c2 ; dxp2 k¼0
ð12Þ
N df ðxÞ X xp2 xp3 ð1Þ ¼ þ c2 þ þ cp2 x þ cp1 ; ak I k ðxÞ þ c1 dx ðp 2Þ! ðp 3Þ! k¼0 N X
ð13Þ
xp1 xp2 þ c2 þ þ cp1 x þ cp ; ð14Þ ðp 1Þ! ðp 2Þ! k¼0 R ðpÞ R ðp1Þ R ð1Þ ðp1Þ ðp2Þ ð0Þ p where I k ðxÞ ¼ I k ðxÞdx; I k ðxÞ ¼ I k ðxÞdx; . . . ; I k ðxÞ ¼ I k ðxÞdx, and fci gi¼1 are the constants of integration. f ðxÞ ¼
ð0Þ
ak I k ðxÞ þ c1
2.3. An analysis of the integral formulation Since a truncated Chebyshev series expansion representing dp f =dxp is the interpolant of degree N, the integration process defined by (10)–(14) leads to an approximate expression for f that is the interpolation
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
287
polynomial of degree ðN þ pÞ with ðN þ pÞ coefficients. Based on this observation, in addition to (14), we also consider the following expansion: f ðxÞ ¼
N X
ak T k ðxÞ þ c1 g1 ðxÞ þ c2 g2 ðxÞ þ þ cp gp ðxÞ;
ð15Þ
k¼0
p
where fgi ðxÞgi¼1 are chosen in such a way that they are polynomials, all basis functions in (15) are linearly independent, and the resultant expansion (15) has the same degree of the polynomial and the same number of expansion coefficients as (14). Possible choices for such basis functions include
T N þ1 ðxÞ; T N þ2 ðxÞ; . . . ; T N þp ðxÞ N þ1 N þ2 x ; x ; . . . ; xN þp :
ð16Þ
and
ð17Þ
From here on, ICSI ; ICSII and ICSIII are used to represent three schemes of the integral collocation formulation, (14), (15) and (16), and (15) and (17), respectively. The value of p in (10) is regarded as the order of the integral collocation scheme, denoted by ICSp . A differential collocation scheme can be considered as a special case of ICS by letting p be zero (ICS0 ). ðiÞ k¼N ;i¼p N To make notations simple, we also use fI k ðxÞgk¼0;i¼0 to denote the basis functions associated with fak gk¼0 p p in (15) and its derivatives, and fgi ðxÞgi¼1 to represent the basis functions associated with fci gi¼1 in (14). Solution procedures for the three schemes are exactly the same. The evaluation of f and its derivatives at the G–L points leads to p d d f b ðpÞ^s; ¼I ½p dxp p1 dd f b ðp1Þ^s; ¼I ½p dxp1
ð18Þ ð19Þ
c df b ð1Þ^s; ¼I ½p dx ð0Þ b ^s; f^ ¼ I
ð20Þ ð21Þ
½p
where subscript [] and superscript () are used to indicate the orders of ICS and derivative function, respectively; ^s ¼ ð^ a; ^cÞT in which ^ a ¼ ða0 ; a1 ; . . . ; aN ÞT and ^c ¼ ðc1 ; c2 ; . . . ; cp ÞT ; 2
ðpÞ
I 0 ðx0 Þ;
6 ðpÞ 6 b ðpÞ ¼ 6 I 0 ðx1 Þ; I ½p 6 4
ðpÞ I 0 ðxN Þ;
ðpÞ
dp g 1 dxp
ðx0 Þ;
dp g2 dxp
ðx0 Þ;
;
I N ðx1 Þ;
ðpÞ
dp g 1 dxp
ðx1 Þ;
dp g2 dxp
ðx1 Þ;
;
dp gp dxp dp gp dxp
ðpÞ ; I N ðxN Þ;
dp g1 dxp
ðxN Þ;
dp g2 dxp
ðxN Þ;
;
dp gp dxp
ðpÞ
;
I N ðx0 Þ;
ðpÞ
;
I 1 ðx0 Þ; I 1 ðx1 Þ; ðpÞ I 1 ðxN Þ;
ðx0 Þ
3
7 ðx1 Þ 7 7; 7 5 ðxN Þ
288
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
2 b ðp1Þ I ½p
ðp1Þ
I0
ðx0 Þ;
6 6 6 ðp1Þ ðx1 Þ; 6I ¼6 0 6 6 4 ðp1Þ
I0
;
IN
ðx0 Þ;
dp1 g1 dxp1
;
ðp1Þ I N ðx1 Þ;
dp1 g1 dxp1
ðp1Þ
ðp1Þ
ðxN Þ; ; I N
; and 2 ð0Þ ð0Þ I ðx0 Þ; I 1 ðx0 Þ; 6 0 6 ð0Þ 6 I ðx1 Þ; I ð0Þ ðx1 Þ; 0 1 ð0Þ b ¼6 I 6 ½p 6 6 4 ð0Þ ð0Þ I 0 ðxN Þ; I 1 ðxN Þ;
dp1 g2 dxp1
ðx1 Þ;
dp1 g2 dxp1
ðxN Þ;
p1
d g1 dxp1
ðx0 Þ; ðx1 Þ;
ðxN Þ;
p1
d g2 dxp1
;
dp1 gp dxp1
;
dp1 gp dxp1
ðxN Þ; ;
g1 ðx0 Þ;
g2 ðx0 Þ;
;
I N ðx1 Þ;
ð0Þ
g1 ðx1 Þ;
g2 ðx1 Þ;
;
I N ðx0 Þ;
;
ð0Þ
I N ðxN Þ; g1 ðxN Þ; g2 ðxN Þ;
ðx0 Þ
3
7 7 7 ðx1 Þ 7 7; 7 7 5
ð0Þ
;
;
ðx0 Þ;
dp1 gp dxp1
gp ðx0 Þ
ðxN Þ 3
7 7 gp ðx1 Þ 7 7 7: 7 7 5 ; gp ðxN Þ
3. Function interpolation Consider a function f ðxÞ defined in [1, 1]. The domain of interest is discretized using the G–L points. Here, we concern the case, where given information consists of the values of the function at the grid points and some ‘‘extra” values. The latter can be the values of f and its derivatives at some points that do not coincide with the grid nodes. Let xbi and fbi ðdk fbi =dxk Þ with i ¼ f1; 2; . . .g denote the extra points and the extra information values, respectively. Unlike conventional differential formulations, the integral collocation formulation can easily incorporate extra information into the Chebyshev approximations. Two approaches are proposed below. 3.1. Approach 1 For the sake of simplicity, assume that there are p=2 extra points (p-an even number) and each point is associated with two given values, f and df =dx. One thus has p extra values T f^ extra ¼ fb1 ; dfb1 =dx; . . . ; fbp2 ; dfbp2 =dx : ð22Þ The expansion coefficients can be determined using the ICS scheme of order p ðICSp Þ ! " ð0Þ # b I f^ ½p b s; ^s ¼ C^ ¼ b f^ extra B
ð23Þ
where 2
ð0Þ
I 0 ðxb1 Þ;
6 ð1Þ 6 I 0 ðxb1 Þ; 6 b ¼6 B 6 6 ð0Þ 6 I ðxbp Þ; 4 0 2 ð1Þ
I 0 ðxbp2 Þ;
ð0Þ
I N ðxb1 Þ;
ð0Þ
g1 ðxb1 Þ;
g2 ðxb1 Þ;
I 1 ðxb1 Þ; ;
ð1Þ
I N ðxb1 Þ;
ð1Þ
dg1 dx
dg2 dx
ð0Þ I 1 ðxbp2 Þ;
;
ð0Þ I N ðxbp2 Þ;
g1 ðxbp2 Þ;
g2 ðxbp2 Þ;
;
;
I N ðxbp2 Þ;
dg1 dx
dg2 dx
;
I 1 ðxb1 Þ; ;
ð1Þ
I 1 ðxbp2 Þ;
ð1Þ
ðxb1 Þ;
ðxbp2 Þ;
;
ðxb1 Þ; ;
ðxbp2 Þ;
gp ðxb1 Þ
3
7 ðxb1 Þ 7 7 7 7; 7 gp ðxbp2 Þ 7 5
dgp dx
dgp dx
ðxbp2 Þ
b is the system matrix of dimension ðN þ 1 þ pÞ ðN þ 1 þ pÞ and other notations are defined as before. C The above expression indicates that the integral formulation takes into account the extra information values. After solving (23) for ^s, one can easily calculate the values of derivatives of f at the grid points using (18)–(20).
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
289
3.2. Approach 2 As mentioned earlier, the constants of integration are generated for the purpose of dealing with extra information. It can be seen that every grid point is associated with the same set of integration constants. The relationship between ^c and f^ extra is as follows ^ a ^ b b 2^c; b 1^ f extra ¼ B aþB ð24Þ ¼B ^c or b 1 B b 1^ b 1 f^ extra ; ^c ¼ B aþB 2 2
ð25Þ
b 1 and B b 2 are the first (N+1) and the last p columns of B, b respectively. It is noted that one can solve where B (24) for ^c in an analytical manner. Substitution of (25) into (21) yields b a þ ^k; f^ ¼ C^
ð26Þ
where
b ð0Þ B b 1 B b1 b ¼ b ð0Þ I and C I 2 ½p ½p 2 1 ^k ¼ I b ð0Þ B b 1 f^ extra ; 2 ½p 2 ð0Þ b ð0Þ are the first (N+1) and the last p columns of I b ð0Þ , respectively. The expansion b and I in which I ½p ½p ½p 1
2
^ and (25) for ^c. coefficients are then obtained through (26) for a Approach 1 and Approach 2 are equivalent from the mathematical point of view. However, Approach 1 involves solving one set of equations, while Approach 2 involves solving two smaller sets of equations. Consider a function f ¼ sinðpxÞ with 1 6 x 6 1. In addition to the grid function values, there are four extra values given (f and df =dx at x ¼ 1=3 and x ¼ 1=3). The three integral collocation schemes are employed to evaluate the values of derivatives of f at the grid points. Since there are four extra information values, one can employ ICSs of order 4. Each scheme is implemented in conjunction with Approach 1 and Approach 2. Results obtained are given in Tables 1–3. They indicate that the three schemes of the Table 1 f ¼ sinðpxÞ; 1 6 x 6 1: relative L2 errors, denoted by Ne, of derivatives of f at the grid points and the condition numbers of the system matrix A, denoted by cond(A), obtained by the ICSI4 scheme Grid (N + 1)
cond(A)
N e ðdf =dxÞ
N e ðd2 f =dx2 Þ
N e ðd3 f =dx3 Þ
N e ðd4 f =dx4 Þ
Approach 1 5 7 9 11 13 15 17
3.4e+5 2.0e+6 1.0e+8 1.6e+9 4.8e+8 5.4e+8 3.9e+9
1.1e02 3.8e04 7.6e06 1.0e07 1.1e09 9.1e12 4.7e14
1.2e01 5.5e03 1.6e04 3.1e06 4.4e08 4.7e10 2.9e12
2.7e01 2.3e02 1.0e03 2.9e05 5.7e07 7.9e09 6.3e11
1.2e+00 1.4e01 9.0e03 3.4e04 8.8e06 1.5e07 1.5e09
Approach 2 5 7 9 11 13 15 17
7.1e+3 4.8e+4 1.8e+6 2.7e+7 8.9e+6 1.2e+7 7.2e+7
1.1e02 3.8e04 7.6e06 1.0e07 1.1e09 9.1e12 6.1e14
1.2e01 5.5e03 1.6e04 3.1e06 4.4e08 4.7e10 3.9e12
2.7e01 2.3e02 1.0e03 2.9e05 5.7e07 7.9e09 8.5e11
1.2e+00 1.4e01 9.0e03 3.4e04 8.8e06 1.5e07 2.1e09
In addition to the grid function values, there are four extra values imposed (f and df =dx at x ¼ 1=3 and x ¼ 1=3).
290
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
Table 2 f ¼ sinðpxÞ; 1 6 x 6 1: relative L2 errors of derivatives of f at the grid points and the condition numbers of the system matrix A by the ICSII 4 scheme Grid (N+1)
cond(A)
N e ðdf =dxÞ
N e ðd2 f =dx2 Þ
N e ðd3 f =dx3 Þ
N e ðd4 f =dx4 Þ
Approach 1 5 7 9 11 13 15 17
3.7e+1 5.6e+1 7.1e+2 3.6e+3 3.9e+2 1.7e+2 6.9e+2
1.1e02 3.8e04 7.6e06 1.0e07 1.1e09 9.1e12 8.6e14
1.2e01 5.5e03 1.6e04 3.1e06 4.4e08 4.7e10 5.3e12
2.7e01 2.3e02 1.0e03 2.9e05 5.7e07 7.9e09 1.1e10
1.2e+00 1.4e01 9.0e03 3.4e04 8.8e06 1.5e07 2.7e09
Approach 2 5 7 9 11 13 15 17
1.8e+1 1.6e+1 1.6e+2 6.8e+2 8.6e+1 4.1e+1 1.3e+2
1.1e02 3.8e04 7.6e06 1.0e07 1.1e09 9.2e12 3.6e14
1.2e01 5.5e03 1.6e04 3.1e06 4.4e08 4.7e10 1.2e12
2.7e01 2.3e02 1.0e03 2.9e05 5.7e07 7.9e09 3.0e11
1.2e+00 1.4e01 9.0e03 3.4e04 8.8e06 1.5e07 5.5e10
In addition to the grid function values, there are 4 extra values imposed (f and df =dx at x ¼ 1=3 and x ¼ 1=3).
Table 3 f ¼ sinðpxÞ; 1 6 x 6 1: relative L2 errors of derivatives of f at the grid points and the condition numbers of the system matrix A by the ICSIII 4 scheme Grid (N+1)
cond(A)
N e ðdf =dxÞ
N e ðd2 f =dx2 Þ
N e ðd3 f =dx3 Þ
N e ðd4 f =dx4 Þ
Approach 1 5 7 9 11 13 15 17
2.5e+3 1.5e+4 8.3e+5 2.2e+7 1.3e+7 3.4e+7 5.0e+8
1.1e02 3.8e04 7.6e06 1.0e07 1.1e09 9.1e12 8.6e14
1.2e01 5.5e03 1.6e04 3.1e06 4.4e08 4.7e10 5.3e12
2.7e01 2.3e02 1.0e03 2.9e05 5.7e07 7.9e09 1.1e10
1.2e+00 1.4e01 9.0e03 3.4e04 8.8e06 1.5e07 2.7e09
Approach 2 5 7 9 11 13 15 17
2.4e+06 1.5e+08 8.7e+10 2.4e+13 1.5e+14 3.8e+15 4.1e+17
1.1e02 3.8e04 7.6e06 1.1e07 2.0e07 1.0e06 2.9e05
1.2e01 5.5e03 1.6e04 3.3e06 3.5e06 1.7e05 7.9e04
2.7e01 2.3e02 1.0e03 3.0e05 2.4e05 1.7e04 8.8e03
1.2e+00 1.4e01 9.0e03 3.5e04 2.3e04 1.9e03 1.3e01
In addition to the grid function values, there are four extra values imposed (f and df =dx at x ¼ 1=3 and x ¼ 1=3).
integral formulation yield similar degrees of accuracy on grids where their system matrices are wellconditioned. It should be pointed out that the system matrix of each integral collocation scheme has an entirely different range of the condition number. In each scheme, Approach 1 and Approach 2 also strongly affect the matrix condition number. Approach 2 is seen to be much better than Approach 1, except for ICSIII . The ICSII scheme appears to be the best one as its condition number is very low, ranged from 101 to 102 (Table 2). The reason for that is probably due to the fact that the system matrix of ICSII is composed largely of Chebyshev polynomials. The approximation scheme based on ICSII and Approach 2 is recommended for use in the interpolation of a function and its derivatives.
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
291
4. One-dimensional biharmonic problems Consider the following 1D biharmonic equation: d4 u d2 u þ ¼ bðxÞ; dx4 dx2
xb1 6 x 6 xb2 ; jxbi j 6 1;
ð27Þ
where bðxÞ is a known driving function, subject to Dirichlet boundary conditions at both ends u1 ; uðxb1 Þ ¼ uðxb2 Þ ¼ u2 ;
du d u1 ðxb1 Þ ¼ ; dx dx du d u2 ðxb2 Þ ¼ : dx dx
The problem domain is embedded in [1, 1] and the extended domain is discretized using the G–L points. Making use of (14)/(15) and its relevant derivatives with p ¼ 4, the governing Eq. (27) and the boundary conditions can be transformed into N X
ð4Þ
ak I k ðxÞ þ
k¼0 N X
N X
ð0Þ
ak I k ðxb1 Þ þ c1 ð1Þ
ak I k ðxb1 Þ þ c1
k¼0 N X
ð0Þ
ak I k ðxb2 Þ þ c1
k¼0 N X
ð28Þ
x3b1 x2 þ c2 b1 þ c3 xb1 þ c4 ¼ u1 ; 6 2
ð29Þ
x2b1 d u1 þ c2 xb1 þ c3 ¼ ; 2 dx
ð30Þ
x3b2 x2 þ c2 b2 þ c3 xb2 þ c4 ¼ u2 ; 6 2
ð31Þ
x2b2 d u2 þ c2 xb2 þ c3 ¼ : 2 dx
ð32Þ
k¼0
k¼0 N X
ð2Þ
ak I k ðxÞ þ c1 x þ c2 ¼ bðxÞ;
ð1Þ
ak I k ðxb2 Þ þ c1
k¼0
N
The evaluation of (28) at the whole set of the G–L points fxi gi¼0 plus the boundary conditions (29)–(32) leads to a determinate system of equations " ð4Þ # b þI b ð2Þ I ½4 ½4 ^s ¼ ^t; ð33Þ b B where ^s ¼ ða0 ; a1 ; . . . ; aN ; c1 ; c2 ; c3 ; c4 ÞT ;
^t ¼ ðb0 ; b1 ; . . . ; bN ; u1 ; du1 =dx; u2 ; du2 =dxÞT
and 2
ð0Þ
ð0Þ
ð0Þ
I 0 ðxb1 Þ; I 1 ðxb1 Þ; ; I N ðxb1 Þ;
6 ð1Þ 6 I ðx Þ; I ð1Þ ðx Þ; ; I ð1Þ ðx Þ; b1 b1 b1 6 N 1 b B ¼ 6 0ð0Þ ð0Þ ð0Þ 6 I ðxb2 Þ; I ðxb2 Þ; ; I ðxb2 Þ; N 1 4 0 ð1Þ ð1Þ ð1Þ I 0 ðxb2 Þ; I 1 ðxb2 Þ; ; I N ðxb2 Þ;
x3b1 =6; x2b1 =2; xb1 ; x2b1 =2; xb1 ;
1;
x3b2 =6; x2b2 =2; xb2 ; x2b2 =2; xb2 ;
1;
1
3
7 07 7 7: 17 5 0
The resultant system (33) can be solved in a direct manner (like Approach 1 in the case of function interpolation) or by splitting it into two smaller sets of equations (like Approach 2).
292
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
The three schemes are numerically verified using the following data xb1 ¼ 2=3;
xb2 ¼ þ2=3;
2
b ¼ p sinðpxÞ þ p4 sinðpxÞ; pffiffiffi u1 ¼ 3=2; d u1 =dx ¼ p=2; pffiffiffi u2 ¼ þ 3=2; d u2 =dx ¼ p=2: The exact solution of this problem is given by ue ¼ sinðpxÞ: Results obtained are presented in Tables 4–6. Relative L2 errors of u ðN e ðuÞÞ are computed at the grid points. Unlike the case of function interpolation, the construction of the system matrix here is mainly based on the approximation of derivative functions (the differential equation) rather than based on the original function. It can be seen that the first scheme involves more Chebyshev polynomials T k ðxÞ than the others. Numerical results show that ICSI yields a system matrix with the condition number much lower than those associated with ICSII and ICSIII . Its values are considerably small, especially for Approach 2 (Table 4). It is recommended that the numerical scheme based on ICSI and Approach 2 be considered for solving 1D biharmonic equations. 5. Two-dimensional biharmonic problems Consider a 2D Dirichlet biharmonic problem. The governing equation takes the form o4 u o4 u o4 u þ 2 þ ¼ bðx; yÞ; ox4 ox2 oy 2 oy 4
ð34Þ
Table 4 1D biharmonic problem, Dirichlet boundary conditions, 2=3 6 x 6 2=3: relative L2 errors of the solution u and the condition numbers of the system matrix A by the ICSI4 scheme Grid (N+1)
Approach 1
Approach 2
cond(A)
N e ðuÞ
cond(A)
N e ðuÞ
5 7 9 11 13 15 17 19 21
4.6e+1 5.3e+1 6.0e+1 6.6e+1 7.2e+1 7.7e+1 8.2e+1 8.6e+1 9.1e+1
1.0e01 2.3e03 1.2e05 2.4e07 4.9e10 1.4e11 1.5e14 2.2e15 6.8e16
2.2e+0 2.3e+0 2.2e+0 2.2e+0 2.1e+0 2.1e+0 2.1e+0 2.1e+0 2.1e+0
1.0e01 2.3e03 1.2e05 2.4e07 4.9e10 1.4e11 1.5e14 1.0e15 4.4e16
Table 5 1D biharmonic problem, Dirichlet boundary conditions, 2=3 6 x 6 2=3: relative L2 errors of the solution u and the condition numbers of the system matrix A by the ICSII 4 scheme Grid (N+1)
Approach 1
Approach 2
cond(A)
N e ðuÞ
cond(A)
N e ðuÞ
5 7 9 11 13 15 17 19 21
1.4e+5 9.4e+5 4.3e+6 1.5e+7 4.8e+7 1.2e+8 3.1e+8 6.8e+8 1.4e+9
1.0e01 2.3e03 1.2e05 2.4e07 4.9e10 1.4e11 2.3e14 4.9e15 2.0e15
1.5e+3 4.6e+3 4.1e+4 8.9e+4 4.2e+5 8.2e+5 2.5e+6 4.9e+6 1.0e+7
1.0e01 2.3e03 1.2e05 2.4e07 4.9e10 1.4e11 6.2e14 3.9e14 2.2e13
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
293
Table 6 1D biharmonic problem, Dirichlet boundary conditions, 2=3 6 x 6 2=3: relative L2 errors of the solution u and the condition numbers of the system matrix A by the ICSIII 4 scheme Grid (N+1)
Approach 1
Approach 2
cond(A)
N e ðuÞ
cond(A)
N e ðuÞ
5 7 9 11 13 15 17 19 21
1.9e+3 1.6e+4 1.6e+5 1.0e+6 4.9e+6 1.7e+7 5.3e+7 2.3e+8 1.5e+9
1.0e01 2.3e03 1.2e05 2.4e07 4.9e10 1.4e11 1.5e14 1.5e14 2.7e15
2.3e+03 2.2e+04 1.4e+05 1.5e+06 1.5e+07 1.8e+08 1.7e+09 2.0e+10 2.0e+11
1.0e01 2.3e03 1.2e05 2.4e07 4.9e10 1.5e11 1.7e12 5.0e12 4.5e12
where bðx; yÞ is a known driving function, subject to double boundary conditions (u and ou=on, n – the normal direction) along the boundary. The proposed numerical procedure is presented in detail through a domain of irregular shape X depicted in Fig. 1. This irregular domain is embedded in the reference square, which allows the use of tensor product grids ðN x þ 1Þ ðN y þ 1Þ. The present method divides the prescribed boundary conditions into two groups. The first group is made of the given values of the solution at the regular boundary points (grid points which lie on the actual boundary), while the second group is formed from the remaining boundary conditions. The latter consists of normal derivative boundary conditions at the regular boundary
Fig. 1. 2D Biharmonic equation: irregular domain, extended domain and discretization. The mark + is used to denote the interior points of the actual domain X.
294
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
points, and the boundary conditions at the irregular boundary points (the intersections of grid lines and irregular boundaries). The construction of approximate expressions for o2 u=ox2 and o4 u=ox4 is similar to that for o2 u=oy 2 and o4 u=oy 4 . Only derivatives of u with respect to y are considered here. Unlike the case of 1D biharmonic problems, the Chebyshev approximations will be expressed in terms of nodal variable values (physical space) to avoid the problem of increasing the system matrix size. Some typical cases are as follows. 5.1. Case 1 – line aa0 Along this line, one needs to impose the values of u only. The ICS0 scheme can be employed here. The values of o2 u=oy 2 and o4 u=oy 4 at the grid points are computed using (3)–(5). 5.2. Case 2 – line bb0 This line intersects the actual boundary at two points y b1 and y b2 ðy b1 < y b2 Þ. The first boundary point y b1 is a grid node (regular boundary point). There is one extra value associated with this node, namely ou1 =on. If the second boundary point y b2 is also a grid node, the treatment for y b2 will be the same as that for y b1 . There are thus two extra values in total along this line. To impose them, one can use the ICS2 scheme. The transformation of the spectral space into the physical space is based on the following system: 0 1 0 1 0 1 ^ u " ð0Þ # ^ ^ a a b I B ou1 C C bB C ½2 B ¼ ¼ C c ð35Þ c @ oy A @ 1A @ 1 A; b B o u2 c2 c2 oy b is the conversion matrix of dimension ðN y þ 3Þ ðN y þ 3Þ; ^a ¼ ða0 ; a1 ; . . . ; aN y ÞT ; ^u ¼ ðu0 ; u1 ; where C . . . ; uN y ÞT , and 2 3 ð1Þ ð1Þ ð1Þ dg2 1 I 0 ðy b1 Þ; I 1 ðy b1 Þ; ; I N y ðy b1 Þ; dg ðy Þ; ðy Þ b1 b1 dy dy b ¼4 5 : B ð1Þ ð1Þ ð1Þ dg1 dg2 I 0 ðy b2 Þ; I 1 ðy b2 Þ; ; I N y ðy b2 Þ; dy ðy b2 Þ; dy ðy b2 Þ ½2
Solving (35), in a direct manner (Approach 1), yields 0 1 0 1 ^ u ^ a B C b 1 B ou1 C @ c1 A ¼ C @ oy A: o u2 c2 oy
ð36Þ
The values of o2 u=oy 2 and o4 u=oy 4 at the grid points are then computed by 0 1 ^ u d 2 o ou B ð2Þ 1 b C b @ oyu1 C ¼I A; ½2 2 oy o u
ð37Þ
2
0
oy
1 ^ u d 4 ou o u1 C b ð4Þ C b 1 B ¼I @ oy A; ½2 oy 4 o u
ð38Þ
2
oy
where 2 b ð4Þ I ½2
d2 T 0 dy 2
ðy 0 Þ;
6 6 d2 T 0 6 ¼ 6 dy 2 ðy 1 Þ; 6 4 d2 T 0 ðy N y Þ; dy 2
d2 T 1 dy 2 d2 T 1 dy 2
ðy 0 Þ; ðy 1 Þ;
d2 T 1 dy 2
ðy N y Þ;
;
d2 T N y dy 2
ðy 0 Þ;
0;
0
3
7 7 ; ðy 1 Þ; 0; 0 7 7: 7 5 2 d T ; dy N2 y ðy N y Þ; 0; 0 d2 T N y dy 2
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
295
If the second boundary point y b2 does not coincide with any grid points, there are two extra values, namely u2 and o u2 =on, at y b2 , leading to a total of three extra values along line bb0 . They can be imposed through the ICS3 scheme. The transformation system is given by 0 1 0 1 0 1 ^ u ^ ^ a a " # ð0Þ B B ou1 C C B b I ½3 B c1 C c1 C B oy C C bB ð39Þ B C¼ B C¼C B C; u2 A @ @ c2 A b @ c2 A B o u2 c3 c3 oy b is the conversion matrix of dimension ðN y þ 4Þ ðN y þ 4Þ and where C 2 3 ð1Þ ð1Þ ð1Þ dg3 dg2 1 ðy Þ; ðy Þ; ðy Þ I 0 ðy b1 Þ; I 1 ðy b1 Þ; ; I N y ðy b1 Þ; dg b1 b1 b1 dy dy dy 6 7 7 ð0Þ ð0Þ ð0Þ b ¼6 B I ðy Þ; I ðy Þ; ; I ðy Þ; g ðy Þ; g ðy Þ; g 6 0 b2 1 b2 2 b2 3 ðy b2 Þ 7 : b2 b2 Ny 1 4 5 ð1Þ ð1Þ ð1Þ dg3 dg2 1 I 0 ðy b2 Þ; I 1 ðy b2 Þ; ; I N y ðy b2 Þ; dg ðy Þ; ðy Þ; ðy Þ b2 b2 b2 dy dy dy ½3
The remaining steps for obtaining expressions of o2 u=oy 2 and o4 u=oy 4 are similar to the previous case and therefore omitted here for brevity. 5.3. Case 3 – line cc0 This case involves four intersection points: y b1 ; y b2 ; y b3 and y b4 . The first and last points are regular boundary points. Assume that y b2 and y b3 are not grid points. There are six extra values to be imposed. The process of transforming the expansion coefficients into the nodal variable values is based on the following system: 0 1 0 1 0 1 ^ u ^ ^ a a B ou1 C Bc C Bc C B oy C B 1C B 1C B C B C B C C " B u C B B c2 C # 2 B C B C b ð0Þ B c2 C B ou2 C I C C B ½6 bB B oy C ¼ ð40Þ B c3 C ¼ C B c 3 C; B C C C B B b C B B C C B B u B 3C B c4 C B c4 C B ou3 C B C B C B oy C @ c5 A @ c5 A @ A o u4 c6 c6 oy b is the conversion matrix of dimension ðN y þ 7Þ ðN y þ 7Þ and where C 2
ð1Þ
I 0 ðy b1 Þ;
6 6 ð0Þ 6 I 0 ðy b2 Þ; 6 6 ð1Þ 6 I ðy Þ; 6 0 b2 b B¼6 6 ð0Þ 6 I 0 ðy b3 Þ; 6 6 ð1Þ 6 I ðy Þ; 6 0 b3 4 ð1Þ I 0 ðy b4 Þ;
ð1Þ
; I N y ðy b1 Þ;
dg1 dy
ðy b1 Þ;
ð0Þ
; I N y ðy b2 Þ; g1 ðy b2 Þ; ð1Þ
; I N y ðy b2 Þ;
dg1 dy
ðy b2 Þ;
ð0Þ
; I N y ðy b3 Þ; g1 ðy b3 Þ; ð1Þ
dg1 dy
ðy b3 Þ;
ð1Þ
dg1 dy
ðy b4 Þ;
; I N y ðy b3 Þ; ; I N y ðy b4 Þ;
;
dg6 dy
ðy b1 Þ
3
7 7 ; g6 ðy b2 Þ 7 7 7 dg6 ; dy ðy b2 Þ 7 7 7 : 7 ; g6 ðy b3 Þ 7 7 7 dg6 ; dy ðy b3 Þ 7 7 5 dg6 ; dy ðy b4 Þ ½6
It can be seen that the Chebyshev approximations of derivatives at a grid point are now expressed in terms of the nodal values of u along the grid line that goes through that point. As with finite-difference and finite-element techniques, one will gather these approximations together to form global matrices for the discretization of the PDE. Their final forms can be written as
296
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299 i of u ¼ oxi i of u ¼ oy i
e ix ~ D u þ ~k ix ;
ð41Þ
e iy ~ D u þ ~k iy ;
ð42Þ
e iy are e ix and D where ~: is used to denote a vector/matrix that is associated with a 2D tensor product grid, D known matrices of dimension ðN x þ 1ÞðN y þ 1Þ ðN x þ 1ÞðN y þ 1Þ, and ~k ix and ~k iy are known vectors of length ðN x þ 1ÞðN y þ 1Þ. The mixed partial fourth-order derivatives can be computed using the following relation o4 u 1 o2 o 2 u o2 o2 u ¼ þ : ð43Þ ox2 oy 2 2 ox2 oy 2 oy 2 ox2
Fig. 2. Biharmonic equation, single domain: geometry and discretization. The mark + is used to denote the interior points of the actual domain X.
Table 7 Biharmonic equation, single domain, Approach 2: relative L2 norms of the solution u at the interior points of the actual domain by the ICSI and ICSII schemes Grid
44 66 88 10 10 12 12 14 14
ICSI
ICSII
cond(A)
N e ðuÞ
cond(A)
N e ðuÞ
2.1e+2 7.7e+3 8.6e+4 6.5e+5 4.4e+6 3.4e+7
3.0e01 1.4e02 1.5e04 6.0e07 8.7e09 7.5e11
2.1e+2 7.7e+3 8.6e+4 6.5e+5 4.4e+6 3.4e+7
3.0e01 1.4e02 1.5e04 6.0e07 8.7e09 7.5e11
Results concerning the matrix condition number are also included.
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
297
E D
F
G
A
B
C
Subdomain 2
Subdomain 1
Subdomain 3
Fig. 3. Biharmonic equation, three subdomains: geometry, extended subdomains and discretization. The mark + is used to denote the interior points of the actual domain X.
298
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
In the calculation of the RHS of (43), approximate expressions (41) and (42) are used to evaluate the values of o2 u=ox2 and o2 u=oy 2 at the grid points, while second-order differential operators are simply replaced by o2 b 2 Iy ðÞ; ðÞ ¼ D ð44Þ x ox o2 b 2 ðÞ; ðÞ ¼ Ix D ð45Þ y oy b x and D b y are the differentiation matrices in the x- and ywhere denotes the Kronecker tensor product; D b directions, respectively ( D ij are defined by (6)–(8); and Iy and Ix are the identity matrices of dimension ðN y þ 1Þ ðN y þ 1Þ and ðN x þ 1Þ ðN x þ 1Þ, respectively. In (44) and (45), the grid points are numbered from bottom to top and from left to right. It is worth mentioning that approximate expressions for derivatives of u already contain information about the boundary of X (location and value). By collocating the governing equation at the grid points and then deleting rows corresponding to points that lie on the boundary, a determinate system of algebraic equations is obtained, which is solved for the approximate solution. The proposed procedure is tested using the following functions: b ¼ 4 cosðpxÞ cosðpyÞ þ cosðpxÞ þ cosðpyÞ; 1 ue ¼ 4 ½1 þ cosðpxÞ½1 þ cosðpyÞ: p
ð46Þ ð47Þ
Two cases, namely single domain and multi-domains, are studied. 5.4. Single domain A unit circular domain is considered (Fig. 2). This domain is embedded in the reference square. Results obtained by the ICSI and ICSII schemes are presented in Table 7. Unlike the cases of function approximations and 1D biharmonic equations, a numerical solution here is solved directly in the physical space. It can be seen that, in the physical space, the ICSI and ICSII schemes essentially yield the same results with respect to the condition number and the relative L2 error. An exponential rate of convergence with grid refinement is achieved. 5.5. Multi-domains An irregular domain, which is displayed in Fig. 3, is divided into three subdomains. Sudomain 1 is a simply connected domain, while subdomains 2 and 3 are multiply connected domains. Points A, B, C, D, E, F and G are located at (0, 0), (0, 1), (1, 1), (1, 1), (7/12, 1), (1, 7/12) and (1, 0), respectively. The circular hole is of radius 1/3 and centered at (1/2, 1/2), while the square hole is taken as ½1=6; 1=6 ½5=6; 5=6. Along the subdomain interfaces, the unknowns are chosen to be u and ou=on. These unknown values are determined using the continuity of o2 u=on2 and o3 u=on3 across the interfaces. Table 8 presents relative L2 errors of the solution u at the interior points of the three subdomains and of the whole domain. It can be seen that the proposed technique yields spectral accuracy. Table 8 Biharmonic equation, three subdomains: relative L2 norms of the solution u at the interior points of the three subdomains and of the whole domain by the ICS scheme Grid 55 77 99 11 11 13 13
N 1e ðuÞ
N 2e ðuÞ
N 3e ðuÞ
N e ðuÞ
1.4e03 6.3e05 2.1e07 9.6e10 4.5e12
4.5e03 4.9e05 3.6e07 2.0e09 1.1e11
1.4e03 1.5e05 1.3e07 2.5e09 9.4e12
2.8e03 4.8e05 2.5e07 1.9e09 8.7e12
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
299
6. Concluding remarks This paper reports a Chebyshev integral collocation approach for solving biharmonic equations in irregular domains. The problem domain is embedded in the reference square, and this extended domain is handled using integral collocation schemes. Boundary conditions are simply divided into two groups. The first group is made of the given values of the solution at the regular boundary points, while the second group is formed from the remaining boundary conditions. The latter consists of normal derivative boundary conditions at the regular boundary points, and the boundary conditions at the irregular boundary points. All boundary conditions in the second group can be implemented in a similar fashion, making the present numerical procedure very attractive in terms of simplicity. Very accurate results are achieved using coarse grids. Acknowledgements This research is supported by the Australian Research Council. The authors are grateful to Prof. R.I. Tanner and Prof. X.-J. Fan for fruitful discussions throughout this study. We would like to thank the referees for their helpful comments. References [1] [2] [3] [4] [5]
D. Gottlieb, S.A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, SIAM, Philadelphia, 1977. B. Fornberg, A Practical Guide to Pseudospectral Methods, Cambridge University Press, Cambridge, 1998. L.N. Trefethen, Spectral Methods in MATLAB, SIAM, Philadelphia, 2000. S.A. Orszag, Spectral methods for problems in complex geometries, J. Comput. Phys. 37 (1) (1980) 70–92. R. Glowinski, Yu. Kuznetsov, Distributed Lagrange multipliers based on fictitious domain method for second order elliptic problems, Comput. Methods Appl. Mech. Eng. 196 (8) (2007) 1498–1506. [6] N. Mai-Duy, R.I. Tanner, A spectral collocation method based on integrated Chebyshev polynomials for biharmonic boundary-value problems, J. Comput. Appl. Math. 201 (1) (2007) 30–47.
Applied Mathematical Modelling 33 (2009) 284–299 www.elsevier.com/locate/apm
A spectral collocation technique based on integrated Chebyshev polynomials for biharmonic problems in irregular domains N. Mai-Duy a,*, H. See b, T. Tran-Cong a a
Computational Engineering and Science Research Centre (CESRC), The University of Southern Queensland, Toowoomba, QLD 4350, Australia b School of Chemical and Biomolecular Engineering, The University of Sydney, NSW 2006, Australia Received 27 February 2007; received in revised form 2 November 2007; accepted 5 November 2007 Available online 17 November 2007
Abstract In this paper, an integral collocation approach based on Chebyshev polynomials for numerically solving biharmonic equations [N. Mai-Duy, R.I. Tanner, A spectral collocation method based on integrated Chebyshev polynomials for biharmonic boundary-value problems, J. Comput. Appl. Math. 201 (1) (2007) 30–47] is further developed for the case of irregularly shaped domains. The problem domain is embedded in a domain of regular shape, which facilitates the use of tensor product grids. Two relevant important issues, namely the description of the boundary of the domain on a tensor product grid and the imposition of double boundary conditions, are handled effectively by means of integration constants. Several schemes of the integral collocation formulation are proposed, and their performances are numerically investigated through the interpolation of a function and the solution of 1D and 2D biharmonic problems. Results obtained show that they yield spectral accuracy. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Integral collocation formulation; Biharmonic problems; Complex geometries; Fictitious domains; Chebyshev polynomials
1. Introduction Many engineering problems, such as the deformation of a thin plate and the motion of a fluid, are governed by the biharmonic equations – fourth-order partial differential equations (PDEs). Generally, problems involving high-order PDEs and complex geometries are more difficult to solve than those with second-order PDEs and regular geometries, respectively. Spectral collocation/pseudo-spectral methods (cf. [1–3]) are known to have the capability to provide an exponential rate of convergence as the grid is refined or the degree of the interpolation polynomial is increased. The drawback of these techniques is that they require a computational domain be square 1 6 x; y 6 1. For solving problems with complex geometries, domain decompositions and coordinate transformations have *
Corresponding author. Tel.: +61 7 4631 1324; fax: +61 7 4631 2526. E-mail address: [email protected] (N. Mai-Duy).
0307-904X/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2007.11.002
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
285
usually been employed to convert irregular domains into regular ones (e.g. [4]). Another way, which is a subject of the present study, is based on the use of fictitious domain. It is noted that fictitious-domain/domainembedding techniques can be traced back to the early 1950s. The basic idea behind these techniques is to extend domains of complicated shapes to those of simpler shapes from which the generation of meshes is simple and well-established efficient numerical solvers can be applied. Fictitious domains have been widely used in the context of finite elements, where the boundary conditions are implemented by means of Lagrange multipliers (e.g. [5] and references therein). In the present work, the problem domain of irregular shape is embedded in the reference square. This new domain is then discretized using a tensor product grid. Clearly, the grid points do not generally lie on the boundary of the actual domain. It is thus difficult to impose boundary conditions here with conventional differential approaches. In our earlier work [6] which deals with biharmonic problems defined in rectangular domains, it has been shown that the use of integration to construct the Chebyshev approximations provides an effective implementation of double boundary conditions. Unlike conventional differential formulations, the integral collocation formulation has the capability to generate extra expansion coefficients (integration constants). These additional unknown values can be utilized to incorporate normal derivative boundary conditions into the Chebyshev approximations. This paper is concerned with the development of the integral collocation formulation for the case of irregularly shaped domains. Three schemes of the integral collocation formulation are presented, and their performances are numerically investigated by considering several 1D and 2D test problems. The remainder of the paper is organized as follows. A brief review of differential and integral collocation formulations is given in Section 2. In this section, the integral collocation formulation is analyzed and several practical schemes of the formulation are presented. The proposed numerical procedure based on integral collocation schemes and fictitious domains is then described and verified through the interpolation of a function and the solution of 1D and 2D biharmonic equations in Sections 3–5, respectively. Section 6 gives some concluding remarks. 2. Collocation formulations 2.1. Differential formulation Consider a univariate function f ðxÞ defined in [1, 1]. This function can be represented by the Chebyshev interpolant of degree N as follows: N X f ðxÞ ¼ ak T k ðxÞ; ð1Þ k¼0 N fak gk¼0
N
where are the coefficients of expansion and fT k ðxÞgk¼0 are the Chebyshev polynomials of the first kind defined as T k ðxÞ ¼ cosðk arccosðxÞÞ. Expressions of derivatives of f are then obtained through differentiation. At the Gauss–Lobatto (G–L) points, N pi fxi gNi¼0 ¼ cos ; ð2Þ N i¼0 the values of derivatives of f are simply computed by c df b ð1Þ f^ ¼ D b f^ ; ¼D dx d 2 d f b ð2Þ f^ ¼ D b 2 f^ ; ¼D dx2 p d d f b ðpÞ f^ ¼ D b p f^ ; ¼D dxp
ð3Þ ð4Þ
ð5Þ
286
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
T where the symbol ^: is used to denote a vector/matrix that is associated with a grid line, f^ ¼ ðf0 ; f1 ; . . . ; fN Þ , k k k k c d f b ð:Þ are the differentiation matrices. The entries of ¼ ðddxfk0 ; ddxfk1 ; . . . ; ddxfkN ÞT with k ¼ f1; 2; . . . ; pg, and D dxk
b D b ð1Þ Þ are given by Dð iþj
b ij ¼ ci ð1Þ ; 0 6 i; j 6 N ; i 6¼ j; D cj xi xj xi b ii ¼ D ; 1 6 i 6 N 1; 2ð1 x2i Þ 2 b 00 ¼ D b NN ¼ 2N þ 1 ; D 6
ð6Þ ð7Þ ð8Þ
b can also be where c0 ¼ cN ¼ 2 and ci ¼ 1 for i ¼ f1; 2; . . . ; N 1g. It is noted that the diagonal entries of D obtained in the way that represents exactly the derivative of a constant b ii ¼ D
N X
b ij : D
ð9Þ
j¼0;j6¼i
For the case of smooth functions, the Chebyshev approximation scheme (1)–(5) is known to be very accurate (exponential accuracy) with the error being OðN a Þ in which a depends on the regularity of a function. It should be emphasized that there is a reduction in accuracy for the approximation of derivatives and this reduction is an increasing function of derivative order (cf. [3]). 2.2. Integral formulation The integral collocation formulation uses a truncated Chebyshev series of degree N to represent a derivative of an unknown function f, e.g. a derivative of order p, N N X dp f ðxÞ X ðpÞ ¼ a T ðxÞ ¼ ak I k ðxÞ: k k dxp k¼0 k¼0
ð10Þ
Expressions for lower-order derivatives and the function itself are then obtained through integration as N dp1 f ðxÞ X ðp1Þ ¼ ak I k ðxÞ þ c1 ; dxp1 k¼0
ð11Þ
N dp2 f ðxÞ X ðp2Þ ¼ ak I k ðxÞ þ c1 x þ c2 ; dxp2 k¼0
ð12Þ
N df ðxÞ X xp2 xp3 ð1Þ ¼ þ c2 þ þ cp2 x þ cp1 ; ak I k ðxÞ þ c1 dx ðp 2Þ! ðp 3Þ! k¼0 N X
ð13Þ
xp1 xp2 þ c2 þ þ cp1 x þ cp ; ð14Þ ðp 1Þ! ðp 2Þ! k¼0 R ðpÞ R ðp1Þ R ð1Þ ðp1Þ ðp2Þ ð0Þ p where I k ðxÞ ¼ I k ðxÞdx; I k ðxÞ ¼ I k ðxÞdx; . . . ; I k ðxÞ ¼ I k ðxÞdx, and fci gi¼1 are the constants of integration. f ðxÞ ¼
ð0Þ
ak I k ðxÞ þ c1
2.3. An analysis of the integral formulation Since a truncated Chebyshev series expansion representing dp f =dxp is the interpolant of degree N, the integration process defined by (10)–(14) leads to an approximate expression for f that is the interpolation
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
287
polynomial of degree ðN þ pÞ with ðN þ pÞ coefficients. Based on this observation, in addition to (14), we also consider the following expansion: f ðxÞ ¼
N X
ak T k ðxÞ þ c1 g1 ðxÞ þ c2 g2 ðxÞ þ þ cp gp ðxÞ;
ð15Þ
k¼0
p
where fgi ðxÞgi¼1 are chosen in such a way that they are polynomials, all basis functions in (15) are linearly independent, and the resultant expansion (15) has the same degree of the polynomial and the same number of expansion coefficients as (14). Possible choices for such basis functions include
T N þ1 ðxÞ; T N þ2 ðxÞ; . . . ; T N þp ðxÞ N þ1 N þ2 x ; x ; . . . ; xN þp :
ð16Þ
and
ð17Þ
From here on, ICSI ; ICSII and ICSIII are used to represent three schemes of the integral collocation formulation, (14), (15) and (16), and (15) and (17), respectively. The value of p in (10) is regarded as the order of the integral collocation scheme, denoted by ICSp . A differential collocation scheme can be considered as a special case of ICS by letting p be zero (ICS0 ). ðiÞ k¼N ;i¼p N To make notations simple, we also use fI k ðxÞgk¼0;i¼0 to denote the basis functions associated with fak gk¼0 p p in (15) and its derivatives, and fgi ðxÞgi¼1 to represent the basis functions associated with fci gi¼1 in (14). Solution procedures for the three schemes are exactly the same. The evaluation of f and its derivatives at the G–L points leads to p d d f b ðpÞ^s; ¼I ½p dxp p1 dd f b ðp1Þ^s; ¼I ½p dxp1
ð18Þ ð19Þ
c df b ð1Þ^s; ¼I ½p dx ð0Þ b ^s; f^ ¼ I
ð20Þ ð21Þ
½p
where subscript [] and superscript () are used to indicate the orders of ICS and derivative function, respectively; ^s ¼ ð^ a; ^cÞT in which ^ a ¼ ða0 ; a1 ; . . . ; aN ÞT and ^c ¼ ðc1 ; c2 ; . . . ; cp ÞT ; 2
ðpÞ
I 0 ðx0 Þ;
6 ðpÞ 6 b ðpÞ ¼ 6 I 0 ðx1 Þ; I ½p 6 4
ðpÞ I 0 ðxN Þ;
ðpÞ
dp g 1 dxp
ðx0 Þ;
dp g2 dxp
ðx0 Þ;
;
I N ðx1 Þ;
ðpÞ
dp g 1 dxp
ðx1 Þ;
dp g2 dxp
ðx1 Þ;
;
dp gp dxp dp gp dxp
ðpÞ ; I N ðxN Þ;
dp g1 dxp
ðxN Þ;
dp g2 dxp
ðxN Þ;
;
dp gp dxp
ðpÞ
;
I N ðx0 Þ;
ðpÞ
;
I 1 ðx0 Þ; I 1 ðx1 Þ; ðpÞ I 1 ðxN Þ;
ðx0 Þ
3
7 ðx1 Þ 7 7; 7 5 ðxN Þ
288
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
2 b ðp1Þ I ½p
ðp1Þ
I0
ðx0 Þ;
6 6 6 ðp1Þ ðx1 Þ; 6I ¼6 0 6 6 4 ðp1Þ
I0
;
IN
ðx0 Þ;
dp1 g1 dxp1
;
ðp1Þ I N ðx1 Þ;
dp1 g1 dxp1
ðp1Þ
ðp1Þ
ðxN Þ; ; I N
; and 2 ð0Þ ð0Þ I ðx0 Þ; I 1 ðx0 Þ; 6 0 6 ð0Þ 6 I ðx1 Þ; I ð0Þ ðx1 Þ; 0 1 ð0Þ b ¼6 I 6 ½p 6 6 4 ð0Þ ð0Þ I 0 ðxN Þ; I 1 ðxN Þ;
dp1 g2 dxp1
ðx1 Þ;
dp1 g2 dxp1
ðxN Þ;
p1
d g1 dxp1
ðx0 Þ; ðx1 Þ;
ðxN Þ;
p1
d g2 dxp1
;
dp1 gp dxp1
;
dp1 gp dxp1
ðxN Þ; ;
g1 ðx0 Þ;
g2 ðx0 Þ;
;
I N ðx1 Þ;
ð0Þ
g1 ðx1 Þ;
g2 ðx1 Þ;
;
I N ðx0 Þ;
;
ð0Þ
I N ðxN Þ; g1 ðxN Þ; g2 ðxN Þ;
ðx0 Þ
3
7 7 7 ðx1 Þ 7 7; 7 7 5
ð0Þ
;
;
ðx0 Þ;
dp1 gp dxp1
gp ðx0 Þ
ðxN Þ 3
7 7 gp ðx1 Þ 7 7 7: 7 7 5 ; gp ðxN Þ
3. Function interpolation Consider a function f ðxÞ defined in [1, 1]. The domain of interest is discretized using the G–L points. Here, we concern the case, where given information consists of the values of the function at the grid points and some ‘‘extra” values. The latter can be the values of f and its derivatives at some points that do not coincide with the grid nodes. Let xbi and fbi ðdk fbi =dxk Þ with i ¼ f1; 2; . . .g denote the extra points and the extra information values, respectively. Unlike conventional differential formulations, the integral collocation formulation can easily incorporate extra information into the Chebyshev approximations. Two approaches are proposed below. 3.1. Approach 1 For the sake of simplicity, assume that there are p=2 extra points (p-an even number) and each point is associated with two given values, f and df =dx. One thus has p extra values T f^ extra ¼ fb1 ; dfb1 =dx; . . . ; fbp2 ; dfbp2 =dx : ð22Þ The expansion coefficients can be determined using the ICS scheme of order p ðICSp Þ ! " ð0Þ # b I f^ ½p b s; ^s ¼ C^ ¼ b f^ extra B
ð23Þ
where 2
ð0Þ
I 0 ðxb1 Þ;
6 ð1Þ 6 I 0 ðxb1 Þ; 6 b ¼6 B 6 6 ð0Þ 6 I ðxbp Þ; 4 0 2 ð1Þ
I 0 ðxbp2 Þ;
ð0Þ
I N ðxb1 Þ;
ð0Þ
g1 ðxb1 Þ;
g2 ðxb1 Þ;
I 1 ðxb1 Þ; ;
ð1Þ
I N ðxb1 Þ;
ð1Þ
dg1 dx
dg2 dx
ð0Þ I 1 ðxbp2 Þ;
;
ð0Þ I N ðxbp2 Þ;
g1 ðxbp2 Þ;
g2 ðxbp2 Þ;
;
;
I N ðxbp2 Þ;
dg1 dx
dg2 dx
;
I 1 ðxb1 Þ; ;
ð1Þ
I 1 ðxbp2 Þ;
ð1Þ
ðxb1 Þ;
ðxbp2 Þ;
;
ðxb1 Þ; ;
ðxbp2 Þ;
gp ðxb1 Þ
3
7 ðxb1 Þ 7 7 7 7; 7 gp ðxbp2 Þ 7 5
dgp dx
dgp dx
ðxbp2 Þ
b is the system matrix of dimension ðN þ 1 þ pÞ ðN þ 1 þ pÞ and other notations are defined as before. C The above expression indicates that the integral formulation takes into account the extra information values. After solving (23) for ^s, one can easily calculate the values of derivatives of f at the grid points using (18)–(20).
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
289
3.2. Approach 2 As mentioned earlier, the constants of integration are generated for the purpose of dealing with extra information. It can be seen that every grid point is associated with the same set of integration constants. The relationship between ^c and f^ extra is as follows ^ a ^ b b 2^c; b 1^ f extra ¼ B aþB ð24Þ ¼B ^c or b 1 B b 1^ b 1 f^ extra ; ^c ¼ B aþB 2 2
ð25Þ
b 1 and B b 2 are the first (N+1) and the last p columns of B, b respectively. It is noted that one can solve where B (24) for ^c in an analytical manner. Substitution of (25) into (21) yields b a þ ^k; f^ ¼ C^
ð26Þ
where
b ð0Þ B b 1 B b1 b ¼ b ð0Þ I and C I 2 ½p ½p 2 1 ^k ¼ I b ð0Þ B b 1 f^ extra ; 2 ½p 2 ð0Þ b ð0Þ are the first (N+1) and the last p columns of I b ð0Þ , respectively. The expansion b and I in which I ½p ½p ½p 1
2
^ and (25) for ^c. coefficients are then obtained through (26) for a Approach 1 and Approach 2 are equivalent from the mathematical point of view. However, Approach 1 involves solving one set of equations, while Approach 2 involves solving two smaller sets of equations. Consider a function f ¼ sinðpxÞ with 1 6 x 6 1. In addition to the grid function values, there are four extra values given (f and df =dx at x ¼ 1=3 and x ¼ 1=3). The three integral collocation schemes are employed to evaluate the values of derivatives of f at the grid points. Since there are four extra information values, one can employ ICSs of order 4. Each scheme is implemented in conjunction with Approach 1 and Approach 2. Results obtained are given in Tables 1–3. They indicate that the three schemes of the Table 1 f ¼ sinðpxÞ; 1 6 x 6 1: relative L2 errors, denoted by Ne, of derivatives of f at the grid points and the condition numbers of the system matrix A, denoted by cond(A), obtained by the ICSI4 scheme Grid (N + 1)
cond(A)
N e ðdf =dxÞ
N e ðd2 f =dx2 Þ
N e ðd3 f =dx3 Þ
N e ðd4 f =dx4 Þ
Approach 1 5 7 9 11 13 15 17
3.4e+5 2.0e+6 1.0e+8 1.6e+9 4.8e+8 5.4e+8 3.9e+9
1.1e02 3.8e04 7.6e06 1.0e07 1.1e09 9.1e12 4.7e14
1.2e01 5.5e03 1.6e04 3.1e06 4.4e08 4.7e10 2.9e12
2.7e01 2.3e02 1.0e03 2.9e05 5.7e07 7.9e09 6.3e11
1.2e+00 1.4e01 9.0e03 3.4e04 8.8e06 1.5e07 1.5e09
Approach 2 5 7 9 11 13 15 17
7.1e+3 4.8e+4 1.8e+6 2.7e+7 8.9e+6 1.2e+7 7.2e+7
1.1e02 3.8e04 7.6e06 1.0e07 1.1e09 9.1e12 6.1e14
1.2e01 5.5e03 1.6e04 3.1e06 4.4e08 4.7e10 3.9e12
2.7e01 2.3e02 1.0e03 2.9e05 5.7e07 7.9e09 8.5e11
1.2e+00 1.4e01 9.0e03 3.4e04 8.8e06 1.5e07 2.1e09
In addition to the grid function values, there are four extra values imposed (f and df =dx at x ¼ 1=3 and x ¼ 1=3).
290
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
Table 2 f ¼ sinðpxÞ; 1 6 x 6 1: relative L2 errors of derivatives of f at the grid points and the condition numbers of the system matrix A by the ICSII 4 scheme Grid (N+1)
cond(A)
N e ðdf =dxÞ
N e ðd2 f =dx2 Þ
N e ðd3 f =dx3 Þ
N e ðd4 f =dx4 Þ
Approach 1 5 7 9 11 13 15 17
3.7e+1 5.6e+1 7.1e+2 3.6e+3 3.9e+2 1.7e+2 6.9e+2
1.1e02 3.8e04 7.6e06 1.0e07 1.1e09 9.1e12 8.6e14
1.2e01 5.5e03 1.6e04 3.1e06 4.4e08 4.7e10 5.3e12
2.7e01 2.3e02 1.0e03 2.9e05 5.7e07 7.9e09 1.1e10
1.2e+00 1.4e01 9.0e03 3.4e04 8.8e06 1.5e07 2.7e09
Approach 2 5 7 9 11 13 15 17
1.8e+1 1.6e+1 1.6e+2 6.8e+2 8.6e+1 4.1e+1 1.3e+2
1.1e02 3.8e04 7.6e06 1.0e07 1.1e09 9.2e12 3.6e14
1.2e01 5.5e03 1.6e04 3.1e06 4.4e08 4.7e10 1.2e12
2.7e01 2.3e02 1.0e03 2.9e05 5.7e07 7.9e09 3.0e11
1.2e+00 1.4e01 9.0e03 3.4e04 8.8e06 1.5e07 5.5e10
In addition to the grid function values, there are 4 extra values imposed (f and df =dx at x ¼ 1=3 and x ¼ 1=3).
Table 3 f ¼ sinðpxÞ; 1 6 x 6 1: relative L2 errors of derivatives of f at the grid points and the condition numbers of the system matrix A by the ICSIII 4 scheme Grid (N+1)
cond(A)
N e ðdf =dxÞ
N e ðd2 f =dx2 Þ
N e ðd3 f =dx3 Þ
N e ðd4 f =dx4 Þ
Approach 1 5 7 9 11 13 15 17
2.5e+3 1.5e+4 8.3e+5 2.2e+7 1.3e+7 3.4e+7 5.0e+8
1.1e02 3.8e04 7.6e06 1.0e07 1.1e09 9.1e12 8.6e14
1.2e01 5.5e03 1.6e04 3.1e06 4.4e08 4.7e10 5.3e12
2.7e01 2.3e02 1.0e03 2.9e05 5.7e07 7.9e09 1.1e10
1.2e+00 1.4e01 9.0e03 3.4e04 8.8e06 1.5e07 2.7e09
Approach 2 5 7 9 11 13 15 17
2.4e+06 1.5e+08 8.7e+10 2.4e+13 1.5e+14 3.8e+15 4.1e+17
1.1e02 3.8e04 7.6e06 1.1e07 2.0e07 1.0e06 2.9e05
1.2e01 5.5e03 1.6e04 3.3e06 3.5e06 1.7e05 7.9e04
2.7e01 2.3e02 1.0e03 3.0e05 2.4e05 1.7e04 8.8e03
1.2e+00 1.4e01 9.0e03 3.5e04 2.3e04 1.9e03 1.3e01
In addition to the grid function values, there are four extra values imposed (f and df =dx at x ¼ 1=3 and x ¼ 1=3).
integral formulation yield similar degrees of accuracy on grids where their system matrices are wellconditioned. It should be pointed out that the system matrix of each integral collocation scheme has an entirely different range of the condition number. In each scheme, Approach 1 and Approach 2 also strongly affect the matrix condition number. Approach 2 is seen to be much better than Approach 1, except for ICSIII . The ICSII scheme appears to be the best one as its condition number is very low, ranged from 101 to 102 (Table 2). The reason for that is probably due to the fact that the system matrix of ICSII is composed largely of Chebyshev polynomials. The approximation scheme based on ICSII and Approach 2 is recommended for use in the interpolation of a function and its derivatives.
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
291
4. One-dimensional biharmonic problems Consider the following 1D biharmonic equation: d4 u d2 u þ ¼ bðxÞ; dx4 dx2
xb1 6 x 6 xb2 ; jxbi j 6 1;
ð27Þ
where bðxÞ is a known driving function, subject to Dirichlet boundary conditions at both ends u1 ; uðxb1 Þ ¼ uðxb2 Þ ¼ u2 ;
du d u1 ðxb1 Þ ¼ ; dx dx du d u2 ðxb2 Þ ¼ : dx dx
The problem domain is embedded in [1, 1] and the extended domain is discretized using the G–L points. Making use of (14)/(15) and its relevant derivatives with p ¼ 4, the governing Eq. (27) and the boundary conditions can be transformed into N X
ð4Þ
ak I k ðxÞ þ
k¼0 N X
N X
ð0Þ
ak I k ðxb1 Þ þ c1 ð1Þ
ak I k ðxb1 Þ þ c1
k¼0 N X
ð0Þ
ak I k ðxb2 Þ þ c1
k¼0 N X
ð28Þ
x3b1 x2 þ c2 b1 þ c3 xb1 þ c4 ¼ u1 ; 6 2
ð29Þ
x2b1 d u1 þ c2 xb1 þ c3 ¼ ; 2 dx
ð30Þ
x3b2 x2 þ c2 b2 þ c3 xb2 þ c4 ¼ u2 ; 6 2
ð31Þ
x2b2 d u2 þ c2 xb2 þ c3 ¼ : 2 dx
ð32Þ
k¼0
k¼0 N X
ð2Þ
ak I k ðxÞ þ c1 x þ c2 ¼ bðxÞ;
ð1Þ
ak I k ðxb2 Þ þ c1
k¼0
N
The evaluation of (28) at the whole set of the G–L points fxi gi¼0 plus the boundary conditions (29)–(32) leads to a determinate system of equations " ð4Þ # b þI b ð2Þ I ½4 ½4 ^s ¼ ^t; ð33Þ b B where ^s ¼ ða0 ; a1 ; . . . ; aN ; c1 ; c2 ; c3 ; c4 ÞT ;
^t ¼ ðb0 ; b1 ; . . . ; bN ; u1 ; du1 =dx; u2 ; du2 =dxÞT
and 2
ð0Þ
ð0Þ
ð0Þ
I 0 ðxb1 Þ; I 1 ðxb1 Þ; ; I N ðxb1 Þ;
6 ð1Þ 6 I ðx Þ; I ð1Þ ðx Þ; ; I ð1Þ ðx Þ; b1 b1 b1 6 N 1 b B ¼ 6 0ð0Þ ð0Þ ð0Þ 6 I ðxb2 Þ; I ðxb2 Þ; ; I ðxb2 Þ; N 1 4 0 ð1Þ ð1Þ ð1Þ I 0 ðxb2 Þ; I 1 ðxb2 Þ; ; I N ðxb2 Þ;
x3b1 =6; x2b1 =2; xb1 ; x2b1 =2; xb1 ;
1;
x3b2 =6; x2b2 =2; xb2 ; x2b2 =2; xb2 ;
1;
1
3
7 07 7 7: 17 5 0
The resultant system (33) can be solved in a direct manner (like Approach 1 in the case of function interpolation) or by splitting it into two smaller sets of equations (like Approach 2).
292
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
The three schemes are numerically verified using the following data xb1 ¼ 2=3;
xb2 ¼ þ2=3;
2
b ¼ p sinðpxÞ þ p4 sinðpxÞ; pffiffiffi u1 ¼ 3=2; d u1 =dx ¼ p=2; pffiffiffi u2 ¼ þ 3=2; d u2 =dx ¼ p=2: The exact solution of this problem is given by ue ¼ sinðpxÞ: Results obtained are presented in Tables 4–6. Relative L2 errors of u ðN e ðuÞÞ are computed at the grid points. Unlike the case of function interpolation, the construction of the system matrix here is mainly based on the approximation of derivative functions (the differential equation) rather than based on the original function. It can be seen that the first scheme involves more Chebyshev polynomials T k ðxÞ than the others. Numerical results show that ICSI yields a system matrix with the condition number much lower than those associated with ICSII and ICSIII . Its values are considerably small, especially for Approach 2 (Table 4). It is recommended that the numerical scheme based on ICSI and Approach 2 be considered for solving 1D biharmonic equations. 5. Two-dimensional biharmonic problems Consider a 2D Dirichlet biharmonic problem. The governing equation takes the form o4 u o4 u o4 u þ 2 þ ¼ bðx; yÞ; ox4 ox2 oy 2 oy 4
ð34Þ
Table 4 1D biharmonic problem, Dirichlet boundary conditions, 2=3 6 x 6 2=3: relative L2 errors of the solution u and the condition numbers of the system matrix A by the ICSI4 scheme Grid (N+1)
Approach 1
Approach 2
cond(A)
N e ðuÞ
cond(A)
N e ðuÞ
5 7 9 11 13 15 17 19 21
4.6e+1 5.3e+1 6.0e+1 6.6e+1 7.2e+1 7.7e+1 8.2e+1 8.6e+1 9.1e+1
1.0e01 2.3e03 1.2e05 2.4e07 4.9e10 1.4e11 1.5e14 2.2e15 6.8e16
2.2e+0 2.3e+0 2.2e+0 2.2e+0 2.1e+0 2.1e+0 2.1e+0 2.1e+0 2.1e+0
1.0e01 2.3e03 1.2e05 2.4e07 4.9e10 1.4e11 1.5e14 1.0e15 4.4e16
Table 5 1D biharmonic problem, Dirichlet boundary conditions, 2=3 6 x 6 2=3: relative L2 errors of the solution u and the condition numbers of the system matrix A by the ICSII 4 scheme Grid (N+1)
Approach 1
Approach 2
cond(A)
N e ðuÞ
cond(A)
N e ðuÞ
5 7 9 11 13 15 17 19 21
1.4e+5 9.4e+5 4.3e+6 1.5e+7 4.8e+7 1.2e+8 3.1e+8 6.8e+8 1.4e+9
1.0e01 2.3e03 1.2e05 2.4e07 4.9e10 1.4e11 2.3e14 4.9e15 2.0e15
1.5e+3 4.6e+3 4.1e+4 8.9e+4 4.2e+5 8.2e+5 2.5e+6 4.9e+6 1.0e+7
1.0e01 2.3e03 1.2e05 2.4e07 4.9e10 1.4e11 6.2e14 3.9e14 2.2e13
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
293
Table 6 1D biharmonic problem, Dirichlet boundary conditions, 2=3 6 x 6 2=3: relative L2 errors of the solution u and the condition numbers of the system matrix A by the ICSIII 4 scheme Grid (N+1)
Approach 1
Approach 2
cond(A)
N e ðuÞ
cond(A)
N e ðuÞ
5 7 9 11 13 15 17 19 21
1.9e+3 1.6e+4 1.6e+5 1.0e+6 4.9e+6 1.7e+7 5.3e+7 2.3e+8 1.5e+9
1.0e01 2.3e03 1.2e05 2.4e07 4.9e10 1.4e11 1.5e14 1.5e14 2.7e15
2.3e+03 2.2e+04 1.4e+05 1.5e+06 1.5e+07 1.8e+08 1.7e+09 2.0e+10 2.0e+11
1.0e01 2.3e03 1.2e05 2.4e07 4.9e10 1.5e11 1.7e12 5.0e12 4.5e12
where bðx; yÞ is a known driving function, subject to double boundary conditions (u and ou=on, n – the normal direction) along the boundary. The proposed numerical procedure is presented in detail through a domain of irregular shape X depicted in Fig. 1. This irregular domain is embedded in the reference square, which allows the use of tensor product grids ðN x þ 1Þ ðN y þ 1Þ. The present method divides the prescribed boundary conditions into two groups. The first group is made of the given values of the solution at the regular boundary points (grid points which lie on the actual boundary), while the second group is formed from the remaining boundary conditions. The latter consists of normal derivative boundary conditions at the regular boundary
Fig. 1. 2D Biharmonic equation: irregular domain, extended domain and discretization. The mark + is used to denote the interior points of the actual domain X.
294
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
points, and the boundary conditions at the irregular boundary points (the intersections of grid lines and irregular boundaries). The construction of approximate expressions for o2 u=ox2 and o4 u=ox4 is similar to that for o2 u=oy 2 and o4 u=oy 4 . Only derivatives of u with respect to y are considered here. Unlike the case of 1D biharmonic problems, the Chebyshev approximations will be expressed in terms of nodal variable values (physical space) to avoid the problem of increasing the system matrix size. Some typical cases are as follows. 5.1. Case 1 – line aa0 Along this line, one needs to impose the values of u only. The ICS0 scheme can be employed here. The values of o2 u=oy 2 and o4 u=oy 4 at the grid points are computed using (3)–(5). 5.2. Case 2 – line bb0 This line intersects the actual boundary at two points y b1 and y b2 ðy b1 < y b2 Þ. The first boundary point y b1 is a grid node (regular boundary point). There is one extra value associated with this node, namely ou1 =on. If the second boundary point y b2 is also a grid node, the treatment for y b2 will be the same as that for y b1 . There are thus two extra values in total along this line. To impose them, one can use the ICS2 scheme. The transformation of the spectral space into the physical space is based on the following system: 0 1 0 1 0 1 ^ u " ð0Þ # ^ ^ a a b I B ou1 C C bB C ½2 B ¼ ¼ C c ð35Þ c @ oy A @ 1A @ 1 A; b B o u2 c2 c2 oy b is the conversion matrix of dimension ðN y þ 3Þ ðN y þ 3Þ; ^a ¼ ða0 ; a1 ; . . . ; aN y ÞT ; ^u ¼ ðu0 ; u1 ; where C . . . ; uN y ÞT , and 2 3 ð1Þ ð1Þ ð1Þ dg2 1 I 0 ðy b1 Þ; I 1 ðy b1 Þ; ; I N y ðy b1 Þ; dg ðy Þ; ðy Þ b1 b1 dy dy b ¼4 5 : B ð1Þ ð1Þ ð1Þ dg1 dg2 I 0 ðy b2 Þ; I 1 ðy b2 Þ; ; I N y ðy b2 Þ; dy ðy b2 Þ; dy ðy b2 Þ ½2
Solving (35), in a direct manner (Approach 1), yields 0 1 0 1 ^ u ^ a B C b 1 B ou1 C @ c1 A ¼ C @ oy A: o u2 c2 oy
ð36Þ
The values of o2 u=oy 2 and o4 u=oy 4 at the grid points are then computed by 0 1 ^ u d 2 o ou B ð2Þ 1 b C b @ oyu1 C ¼I A; ½2 2 oy o u
ð37Þ
2
0
oy
1 ^ u d 4 ou o u1 C b ð4Þ C b 1 B ¼I @ oy A; ½2 oy 4 o u
ð38Þ
2
oy
where 2 b ð4Þ I ½2
d2 T 0 dy 2
ðy 0 Þ;
6 6 d2 T 0 6 ¼ 6 dy 2 ðy 1 Þ; 6 4 d2 T 0 ðy N y Þ; dy 2
d2 T 1 dy 2 d2 T 1 dy 2
ðy 0 Þ; ðy 1 Þ;
d2 T 1 dy 2
ðy N y Þ;
;
d2 T N y dy 2
ðy 0 Þ;
0;
0
3
7 7 ; ðy 1 Þ; 0; 0 7 7: 7 5 2 d T ; dy N2 y ðy N y Þ; 0; 0 d2 T N y dy 2
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
295
If the second boundary point y b2 does not coincide with any grid points, there are two extra values, namely u2 and o u2 =on, at y b2 , leading to a total of three extra values along line bb0 . They can be imposed through the ICS3 scheme. The transformation system is given by 0 1 0 1 0 1 ^ u ^ ^ a a " # ð0Þ B B ou1 C C B b I ½3 B c1 C c1 C B oy C C bB ð39Þ B C¼ B C¼C B C; u2 A @ @ c2 A b @ c2 A B o u2 c3 c3 oy b is the conversion matrix of dimension ðN y þ 4Þ ðN y þ 4Þ and where C 2 3 ð1Þ ð1Þ ð1Þ dg3 dg2 1 ðy Þ; ðy Þ; ðy Þ I 0 ðy b1 Þ; I 1 ðy b1 Þ; ; I N y ðy b1 Þ; dg b1 b1 b1 dy dy dy 6 7 7 ð0Þ ð0Þ ð0Þ b ¼6 B I ðy Þ; I ðy Þ; ; I ðy Þ; g ðy Þ; g ðy Þ; g 6 0 b2 1 b2 2 b2 3 ðy b2 Þ 7 : b2 b2 Ny 1 4 5 ð1Þ ð1Þ ð1Þ dg3 dg2 1 I 0 ðy b2 Þ; I 1 ðy b2 Þ; ; I N y ðy b2 Þ; dg ðy Þ; ðy Þ; ðy Þ b2 b2 b2 dy dy dy ½3
The remaining steps for obtaining expressions of o2 u=oy 2 and o4 u=oy 4 are similar to the previous case and therefore omitted here for brevity. 5.3. Case 3 – line cc0 This case involves four intersection points: y b1 ; y b2 ; y b3 and y b4 . The first and last points are regular boundary points. Assume that y b2 and y b3 are not grid points. There are six extra values to be imposed. The process of transforming the expansion coefficients into the nodal variable values is based on the following system: 0 1 0 1 0 1 ^ u ^ ^ a a B ou1 C Bc C Bc C B oy C B 1C B 1C B C B C B C C " B u C B B c2 C # 2 B C B C b ð0Þ B c2 C B ou2 C I C C B ½6 bB B oy C ¼ ð40Þ B c3 C ¼ C B c 3 C; B C C C B B b C B B C C B B u B 3C B c4 C B c4 C B ou3 C B C B C B oy C @ c5 A @ c5 A @ A o u4 c6 c6 oy b is the conversion matrix of dimension ðN y þ 7Þ ðN y þ 7Þ and where C 2
ð1Þ
I 0 ðy b1 Þ;
6 6 ð0Þ 6 I 0 ðy b2 Þ; 6 6 ð1Þ 6 I ðy Þ; 6 0 b2 b B¼6 6 ð0Þ 6 I 0 ðy b3 Þ; 6 6 ð1Þ 6 I ðy Þ; 6 0 b3 4 ð1Þ I 0 ðy b4 Þ;
ð1Þ
; I N y ðy b1 Þ;
dg1 dy
ðy b1 Þ;
ð0Þ
; I N y ðy b2 Þ; g1 ðy b2 Þ; ð1Þ
; I N y ðy b2 Þ;
dg1 dy
ðy b2 Þ;
ð0Þ
; I N y ðy b3 Þ; g1 ðy b3 Þ; ð1Þ
dg1 dy
ðy b3 Þ;
ð1Þ
dg1 dy
ðy b4 Þ;
; I N y ðy b3 Þ; ; I N y ðy b4 Þ;
;
dg6 dy
ðy b1 Þ
3
7 7 ; g6 ðy b2 Þ 7 7 7 dg6 ; dy ðy b2 Þ 7 7 7 : 7 ; g6 ðy b3 Þ 7 7 7 dg6 ; dy ðy b3 Þ 7 7 5 dg6 ; dy ðy b4 Þ ½6
It can be seen that the Chebyshev approximations of derivatives at a grid point are now expressed in terms of the nodal values of u along the grid line that goes through that point. As with finite-difference and finite-element techniques, one will gather these approximations together to form global matrices for the discretization of the PDE. Their final forms can be written as
296
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299 i of u ¼ oxi i of u ¼ oy i
e ix ~ D u þ ~k ix ;
ð41Þ
e iy ~ D u þ ~k iy ;
ð42Þ
e iy are e ix and D where ~: is used to denote a vector/matrix that is associated with a 2D tensor product grid, D known matrices of dimension ðN x þ 1ÞðN y þ 1Þ ðN x þ 1ÞðN y þ 1Þ, and ~k ix and ~k iy are known vectors of length ðN x þ 1ÞðN y þ 1Þ. The mixed partial fourth-order derivatives can be computed using the following relation o4 u 1 o2 o 2 u o2 o2 u ¼ þ : ð43Þ ox2 oy 2 2 ox2 oy 2 oy 2 ox2
Fig. 2. Biharmonic equation, single domain: geometry and discretization. The mark + is used to denote the interior points of the actual domain X.
Table 7 Biharmonic equation, single domain, Approach 2: relative L2 norms of the solution u at the interior points of the actual domain by the ICSI and ICSII schemes Grid
44 66 88 10 10 12 12 14 14
ICSI
ICSII
cond(A)
N e ðuÞ
cond(A)
N e ðuÞ
2.1e+2 7.7e+3 8.6e+4 6.5e+5 4.4e+6 3.4e+7
3.0e01 1.4e02 1.5e04 6.0e07 8.7e09 7.5e11
2.1e+2 7.7e+3 8.6e+4 6.5e+5 4.4e+6 3.4e+7
3.0e01 1.4e02 1.5e04 6.0e07 8.7e09 7.5e11
Results concerning the matrix condition number are also included.
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
297
E D
F
G
A
B
C
Subdomain 2
Subdomain 1
Subdomain 3
Fig. 3. Biharmonic equation, three subdomains: geometry, extended subdomains and discretization. The mark + is used to denote the interior points of the actual domain X.
298
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
In the calculation of the RHS of (43), approximate expressions (41) and (42) are used to evaluate the values of o2 u=ox2 and o2 u=oy 2 at the grid points, while second-order differential operators are simply replaced by o2 b 2 Iy ðÞ; ðÞ ¼ D ð44Þ x ox o2 b 2 ðÞ; ðÞ ¼ Ix D ð45Þ y oy b x and D b y are the differentiation matrices in the x- and ywhere denotes the Kronecker tensor product; D b directions, respectively ( D ij are defined by (6)–(8); and Iy and Ix are the identity matrices of dimension ðN y þ 1Þ ðN y þ 1Þ and ðN x þ 1Þ ðN x þ 1Þ, respectively. In (44) and (45), the grid points are numbered from bottom to top and from left to right. It is worth mentioning that approximate expressions for derivatives of u already contain information about the boundary of X (location and value). By collocating the governing equation at the grid points and then deleting rows corresponding to points that lie on the boundary, a determinate system of algebraic equations is obtained, which is solved for the approximate solution. The proposed procedure is tested using the following functions: b ¼ 4 cosðpxÞ cosðpyÞ þ cosðpxÞ þ cosðpyÞ; 1 ue ¼ 4 ½1 þ cosðpxÞ½1 þ cosðpyÞ: p
ð46Þ ð47Þ
Two cases, namely single domain and multi-domains, are studied. 5.4. Single domain A unit circular domain is considered (Fig. 2). This domain is embedded in the reference square. Results obtained by the ICSI and ICSII schemes are presented in Table 7. Unlike the cases of function approximations and 1D biharmonic equations, a numerical solution here is solved directly in the physical space. It can be seen that, in the physical space, the ICSI and ICSII schemes essentially yield the same results with respect to the condition number and the relative L2 error. An exponential rate of convergence with grid refinement is achieved. 5.5. Multi-domains An irregular domain, which is displayed in Fig. 3, is divided into three subdomains. Sudomain 1 is a simply connected domain, while subdomains 2 and 3 are multiply connected domains. Points A, B, C, D, E, F and G are located at (0, 0), (0, 1), (1, 1), (1, 1), (7/12, 1), (1, 7/12) and (1, 0), respectively. The circular hole is of radius 1/3 and centered at (1/2, 1/2), while the square hole is taken as ½1=6; 1=6 ½5=6; 5=6. Along the subdomain interfaces, the unknowns are chosen to be u and ou=on. These unknown values are determined using the continuity of o2 u=on2 and o3 u=on3 across the interfaces. Table 8 presents relative L2 errors of the solution u at the interior points of the three subdomains and of the whole domain. It can be seen that the proposed technique yields spectral accuracy. Table 8 Biharmonic equation, three subdomains: relative L2 norms of the solution u at the interior points of the three subdomains and of the whole domain by the ICS scheme Grid 55 77 99 11 11 13 13
N 1e ðuÞ
N 2e ðuÞ
N 3e ðuÞ
N e ðuÞ
1.4e03 6.3e05 2.1e07 9.6e10 4.5e12
4.5e03 4.9e05 3.6e07 2.0e09 1.1e11
1.4e03 1.5e05 1.3e07 2.5e09 9.4e12
2.8e03 4.8e05 2.5e07 1.9e09 8.7e12
N. Mai-Duy et al. / Applied Mathematical Modelling 33 (2009) 284–299
299
6. Concluding remarks This paper reports a Chebyshev integral collocation approach for solving biharmonic equations in irregular domains. The problem domain is embedded in the reference square, and this extended domain is handled using integral collocation schemes. Boundary conditions are simply divided into two groups. The first group is made of the given values of the solution at the regular boundary points, while the second group is formed from the remaining boundary conditions. The latter consists of normal derivative boundary conditions at the regular boundary points, and the boundary conditions at the irregular boundary points. All boundary conditions in the second group can be implemented in a similar fashion, making the present numerical procedure very attractive in terms of simplicity. Very accurate results are achieved using coarse grids. Acknowledgements This research is supported by the Australian Research Council. The authors are grateful to Prof. R.I. Tanner and Prof. X.-J. Fan for fruitful discussions throughout this study. We would like to thank the referees for their helpful comments. References [1] [2] [3] [4] [5]
D. Gottlieb, S.A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, SIAM, Philadelphia, 1977. B. Fornberg, A Practical Guide to Pseudospectral Methods, Cambridge University Press, Cambridge, 1998. L.N. Trefethen, Spectral Methods in MATLAB, SIAM, Philadelphia, 2000. S.A. Orszag, Spectral methods for problems in complex geometries, J. Comput. Phys. 37 (1) (1980) 70–92. R. Glowinski, Yu. Kuznetsov, Distributed Lagrange multipliers based on fictitious domain method for second order elliptic problems, Comput. Methods Appl. Mech. Eng. 196 (8) (2007) 1498–1506. [6] N. Mai-Duy, R.I. Tanner, A spectral collocation method based on integrated Chebyshev polynomials for biharmonic boundary-value problems, J. Comput. Appl. Math. 201 (1) (2007) 30–47.