Adaptive least squares finite integration method for higher-dimensional singular perturbation problems with multiple boundary layers

Applied Mathematics and Computation 271 (2015) 232–250

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Adaptive least squares finite integration method for higher-dimensional singular perturbation problems with multiple boundary layers D.F. Yun a, Z.H. Wen b, Y.C. Hon a,b,∗ a b

Department of Mathematics, City University of Hong Kong, Kowloon Tong, Hong Kong Special Administrative Region College of Mathematics, Taiyuan University of Technology, Taiyuan, China

a r t i c l e

i n f o

Keywords: Singular perturbation Boundary layer Least squares Finite integration method

a b s t r a c t Based on the recently developed finite integration method for solving one-dimensional partial differential equation, we extend in this paper the method by using the technique of least squares to tackle higher-dimensional singular perturbation problems with multiple boundary layers. Theoretical convergence and numerical stability tests indicate that, even with the most simple numerical trapezoidal integration rule, the proposed method provides a stable, efficient, and highly accurate approximate solutions to the singular perturbation problems. An adaptive scheme on the refinement of integration points is also devised to better capture the stiff boundary layers. Illustrative examples are given in both 1D and 2D with comparison among some existing numerical methods. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Singular perturbation problems arise in many physical modeling processes of fluid flows. Due to the stiffness and complexity of the boundary layers exerted by the singularly perturbed equations, it is very difficult to obtain exact solutions for these equations. The rapid advancement of computing technology allows the simulation of the boundary layers through seeking numerical approximate solutions to the mathematical models. During the last decades, these singular perturbation problems have been extensively studied by many researchers using various kinds of numerical methods [1–5] among which the finite difference method (FDM) and finite element method (FEM) are mostly well established. Despite their many attractive features, the FDM and FEM are mesh-dependent which require an appropriate generation of good mesh for accurate and stable approximation to the solution. The generation of good mesh, however, is not a trivial task for problems defined under irregular domain. For stiff problems with boundary layers to be investigated in this paper, some kinds of adaptive schemes on mesh refinement have to be established. In fact, when the singularly perturbed equation is perturbed at a very small parameter  attached to a higher derivative term in the equation, the solution normally forms a layer where the solution derivatives change dramatically within a very narrow interval. This causes most standard computational methods fail to produce satisfactory approximation. It had been proven by Miller et al. [6] that the numerical methods composed of upwind- and central-difference operators on uniform mesh were defective for solving the singularly perturbed differential equations. In fact, Miller et al. had proven that the numerical

∗ Corresponding author at.: Department of Mathematics, City University of Hong Kong, Kowloon Tong, Hong Kong Special Administrative Region. Tel.: +(852) 3442 8675. E-mail addresses: [email protected] (D.F. Yun), [email protected] (Z.H. Wen), [email protected], [email protected] (Y.C. Hon).

http://dx.doi.org/10.1016/j.amc.2015.08.116 0096-3003/© 2015 Elsevier Inc. All rights reserved.

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250

233

oscillation errors are unbounded as  → 0 on uniform meshes. A good review of the numerical methods for these singularly perturbed problems can be found in [7,8]. Recent advances in numerical approximation for the singularly perturbed problems are concentrated on performing the computations on specially designed meshes. Miller et al. [6] reported that the Shishkin meshes [9] provide a better numerical results by using upwind- and central-difference schemes. Similar result was also obtained by Sun and Stynes [10] in applying the FEM using piecewise polynomial test/trial functions on the Shishkin meshes. The use of multiquadric collocation method (MQCM), which is a kind of kernel-based approximation methods, was first proposed by Hon [11] to solve a one-dimensional singular perturbation problem with one boundary layer. The meshless advantage of MQCM was shown to provide an efficient and accurate approximation of the boundary layer. Ling and Trummer [12] later combined the kernel-based approximation method with an adaptive scheme to obtain an improved highly accurate solution for solving the one-dimensional boundary layer problem. For multiple boundary layers problems, Zhang [13] gave a theoretical approach from which numerical approximation can be constructed. In this paper, we extend the recently developed finite integration method (FIM) [14,15] for solving higher-dimensional singularly perturbed equations with multiple boundary layers. Unlike FDM which uses finite quotient formula, the FIM uses numerical quadratic integration rule and hence avoids the wellknown roundoff-discretization error problem in using FDM. Both theoretical and numerical results given in this paper show that, even with the most simple numerical trapezoidal integration rule, the FIM provides a very stable, efficient, and highly accurate approximation to the solutions of these higher-dimensional singular perturbation problems. Numerical comparisons with the MQCM and FDM with Shishkin grids show that the FIM achieves better accuracy at lower computational cost. This paper is organized as follow. The FIM is briefly introduced in Section 2. The formulation of the methodology to solve a singularly perturbed equation with multiple boundary layers is given in Section 3. Theoretical convergence analysis is discussed in Section 4. Section 5 demonstrates the improved accuracy of the method in solving the singularly perturbed equation with one single layer by comparing with the exact solution and some recent works. Numerical stability test is performed to verify the convergence order in one-dimensional case. An adaptive scheme is proposed in Section 6 with illustrative example given in Section 7 on its superior stability and higher convergence to solve the perturbed equation with double layers. In Section 8, we combine the least squares technique with FIM to give a new least-squares finite integration method (LSFIM) to solve a twodimensional singular perturbed problem with stiff boundary layer. Stability and convergence order are derived with numerical examples verified for two-dimensional case. Some concluding remarks with suggested future works are provided in the final conclusion section. 2. Finite integration method Based on our recently developed finite integration method for solving one-dimensional partial differential equations [15] and the technique of least squares, we will extend in this section the method to solve higher-dimensional singular perturbation problems with multiple boundary layers. For simplicity, we first consider 1D case and define F(1) (x) to be the integration of an integrable function f(x)

F (1) (x) =



x

f (ξ ) dξ ,

a

(1)

where x ∈ [a, b] ⊂ R. Assume a partition a = x1 < x1 < . . . < xN = b of the interval [a, b] into subintervals, [a, b] = ∪N−1 [xi , xi+1 ], the above integrai=1 tion when x = xk can be written as

F (1) (xk ) =



xk

a

k−1  

f (ξ ) dξ =

i=1

xi+1 xi

f (ξ ) dξ .

(2)

In each subinterval, the use of numerical trapezoidal rule to the integration gives:



xi+1 xi

f (ξ ) dξ = 

f (xi ) + f (xi+1 ) f  (ξi ) (xi+1 − xi ) − (xi+1 − xi )3 2 6

i 2

fi +

i 2

fi+1 −

f  (ξi ) 3 i , 6

(3)

where fi = f (xi ), fi+1 = f (xi+1 ), i = xi+1 − xi , i = 1, 2, . . . , k − 1 and ξi ∈ (xi−1 , xi ). Combining Eqs. (2) and (3), we have

F (1) (x

k

)=

k−1  i=1

= 

1 2 k  i=1



i 2

f1 +

fi + k−1 

i 2

fi+1

−

2 k−1  i=1

k−1   f (ξi ) 3 i 6 i=1

i−1 + i

i=2

a(ki1) fi −



fi +

f  (ξi ) 3 i , 6

k−1 2

fk −

k−1   f (ξi ) 3 i 6 i=1

(4)

234

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250

(1)

where a11 = 0 if k = 1, and when k ≥ 2:

⎧ 1 ⎪ , ⎪ ⎪ ⎪ ⎨ 2 i−1 + i a(ki1) = ⎪ 2 ⎪ ⎪ ⎪ ⎩ i−1 ,

i = 1, i = 2, 3, . . . , k − 1, i = k.

2

Since







k−1  k−1  maxx | f  (x)| f (ξi ) 3 maxx | f  (x)|

2 i ≤ ( max i ) i ≤ (b − a)( max i 2 ),

6 6 i i

i=1 6

i=1

(5)

f  (ξi ) 3 + +  the error term k−1 i=1 6 i −→ 0 as max1≤i≤N−1 i → 0 providing f is bounded. If we ignore this small magnitude error term, Eq. (4) can be expressed in matrix form as

A(1) f ≈ F(1) , (1)

(1) T

(1)

(1)

where F(1) = [F1 , F2 , . . . , FN ] , Fk

T

integration is given by

⎛

0

0

···

0

1

0

0

···

0

1 +2

2

0

···

0

1 +2

2 +3

3

···

0

.. .

.. .

.. .

..

.. .

1 +2

2 +3

3 +4

···

0

0

⎜ 1 ⎜2 ⎜ 1 ⎜2 ⎜ (1) A = ⎜ 1 ⎜2 ⎜ . ⎜ .. ⎝

2 2

2

2

1 2

2

2

(1)

= F (1) (xk ), f = [ f1 , f2 , . . . , fN ] . The matrix A(1) = (ai j )

2

2

2

.

N−1

N×N

of the first order numerical

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎟ ⎠

(6)

2

For double-layer integration, we let

F (2) (x) =



x

 η a

a

f (ξ ) dξ dη.

Computing the above integration at xk with composite trapezoidal rule, we have

F (2) (xk ) =



xk a

 η a

f (ξ ) dξ dη



f (ξ ) dξ

a x1 i=1 i=1 x=ηi 

k−1 3  x 

3 k i i−1



 (1)  (1)  j  i 

f (x)

= aki ai j f j − − f (ξ ) dξ

6 6 x1  k  = a(ki1)

i=1

xi

k−1  i 3 f (ξ ) dξ − 6

j=1



x

x=ξi j

j=1

i=1 j=1

i=1

j=1

k k i−1 3

   j  (2) (1)  aki fi − aki f (x)

i=1

i=1

j=1

6

− x=ξi j

i=1



k  k i i−1 3

   j  f (x)

= a(ki1) a(i j1) f j − a(ki1) 6

x=ξi j

k−1  i=1

k−1  i 3 − 6 i=1

i 6

 3 

x x1



x=ηi

x

x1



f (ξ ) dξ



f (ξ ) dξ

x=ηi

,

(7)

x=ηi

where ηi ∈ (xi−1 , xi ), ξi j ∈ (x j−1 , x j ). Since



 k i−1 3



j  (1) f (x)

aki 6

i=1 j=1 x=ξ

ij

  N−1

k   maxx | f  (x)|

(1) aki  j ( max i 2 )

≤ 6 i

i=1 j=1 ≤

maxx | f  (x)| (b − a)2 ( max i 2 ), 6 i

(8)

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250

235



  x  k−1

 k−1  (x)|



 maxx | f  (x)| max | f

i 3 x f (ξ ) dξ

i ( max i 2 ) ≤ (b − a)( max i 2 ),

≤ 6 6 i i

i=1 6

x1 i=1 x=η

(9)

i

k

3 (1) i−1  j  f (x)|

k−1 i 3  x

i=1

aki

j=1 6

x=ξi j

−→ 0 and

i=1

6

(



f (ξ ) dξ ) |x=ηi −→ 0 as max1≤i≤N−1 i → 0+ providing that   f , f are both bounded. Ignoring these two small magnitude terms, we can now express the double-layer integration in the matrix form as the error term

x1

A(2) f ≈ F(2) ,

(10)

(2) (2) (2) T (2) [F1 , F2 , . . . , FN ] , Fk

(1) 2

where F(2) = = F (2) (xk ), f = [ f1 , f2 , . . . , fN ] , and A(2) ≈ A(1) A(1) = (A Similarly, for multi-layer integration, we let

F (m) (x) =



x

···

a

 ξ3  ξ2 a

a

T

) .

f (ξ1 )dξ1 dξ2 · · · dξm .

Applying the composite trapezoidal rule to the above integration at xk gives:

F (m) (xk ) =



xk

···

a

 ξ3  ξ2 a

a

f (ξ1 )dξ1 dξ2 · · · dξm ≈

k  im =1

···

i1 i2  k   a(ki1) · · · a(i 1i) a(i 1j) f j  a(kim) fi , i1 =1 j=1

m

2 1

1

i=1

whose matrix form is given by

A(m) f ≈ F(m) , (m)

where F(m) = [F1

(m)

, F2

(m) T

, . . . , FN

(m)

] , Fk

m

(1) (1) (1) (1) = F (m) (xk ), f = [ f1 , f2 , . . . , fN ] , and A(m) ≈ A  A · · · A  = (A ) . T

m

From Eqs. (5) (8) and (9), and by the method of mathematical induction, we then have the following error estimate: Lemma 2.1. Using the numerical trapezoidal rule, we have the error estimate

A(m) f − F(m)  = O(h2 ), m = 1, 2, · · · , where h = maxi |xi+1 − xi |, providing f is smooth enough. We remark here that if the interval [a, b] is divided uniformly by the partition {xi }i=1 , the integration matrix A(1) is lower triangular with all elements predetermined based on the quadrature rule [15]. In other words, the cost of inversion of the resultant matrix from using this FIM is extremely low in comparing with the traditional methods such as FDM, FEM, and MQCM. Since numerical integration is much more stable than numerical differentiation, the FIM turns out to be a very stable and highly accurate numerical method. These advantages persist in higher dimensional problems because the resultant multiple integration matrix A(m) is also lower triangular. For non-uniformly distributed partition, the use of radial basis functions in the integration [15] will keep the FIM highly accurate and mesh-independent but at the loss of the lower triangular matrix advantage. The extension of FIM to higher dimension is routine. Consider a two-dimensional rectangular domain [a, b] × [c, d] with N1 × N2 nodes. Define F(1) (x, y) to be the integration of an integrable function f(x, y) along the x direction: N

F (1) (xk , yk ) =



xk

a

f (ξ , yk ) dξ ,

(11)

where (xk , yk ) denotes the node at the ith column and jth row as indexed in Fig. 1(I). Here, k = ( j − 1)N1 + i and the nodal values of integration (11) can be represented in matrix form as

F(I 1) ≈ A(x,I1) fI ,

(12)

(1)

where FI = [F (1) (x1 , y1 ), F (1) (x2 , y2 ), . . . , F (1) (xM , yM )], fI = [ f (x1 , y1 ), f (x2 , y2 ), . . . , f (xM , yM )] with nodes arranged as shown in Fig. 1(I). The first order integration matrix in x direction is given by

⎛

A(1) ⎜ O ⎜ O A(x,I1) = ⎜ ⎜ . ⎝ .. O

O A(1) O .. . O

O O A(1) .. . O

··· ··· ··· .. . ···

⎞

O O ⎟ ⎟ O ⎟, .. ⎟ ⎠ . ( 1) A

where each sub-matrix A(1) has the same form of (6) when N = N1 and the elements 1 , 2 , · · · , N1 −1 are the nodal distances i = xi+1 − xi in x direction. Similarly, define G(1) (x, y) to be the integration of f(x, y) along the y direction:

G(1) (xk , yk ) =



yk c

f (xk , η) dη,

(13)

236

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250

(I) 17

(II)

18

19

20

14

15

16

9

10

11

12

5

6

7

8

13

1

2

5 P

P

1

15

20

4

9

14

19

3

8

13

18

7

12

17

2

1

4

3

10

6

11

16

Fig. 1. Two different index arrangement of nodes.

where (xk , yk ) denotes the node at the ith column and jth row as indexed in Fig. 1(II). Here, k = (i − 1)N2 + j and the nodal values of integration (13) can be represented in matrix form as 1) G(II1) ≈ A(y,II fII ,

(14)

(1)

where GII = [G(1) (x1 , y1 ), G(1) (x2 , y2 ), . . . , G(1) (xM , yM )], fII = [ f (x1 , y1 ), f (x2 , y2 ), . . . , f (xM , yM )] with nodes arranged as shown in Fig. 1(II). The first order integration matrix in y direction is given by

⎛

A(1) ⎜ O ⎜ O 1) A(y,II =⎜ ⎜ . ⎝ .. O

O A(1) O .. . O

O O A(1) .. . O

··· ··· ··· .. . ···

⎞

O O ⎟ ⎟ O ⎟, .. ⎟ ⎠ . ( 1) A

where each sub-matrix A(1) the same form of (6) when N = N2 and the elements 1 , 2 , · · · , N2 −1 are the nodal distances i = yi+1 − yi in y direction. We note here that the two different index arrangement of nodes can be transformed into each other through the following non-singular matrix P = ( pi, j )M×M :

fI = P fII .

(15)

For instance, if k = ( j − 1)N1 + i, then the kth element of fI is the same as the ((i − 1)N2 + j)th element of fII . This implies that pk,(i−1)N2 + j = 1 and all remaining elements of the kth row are zeroes. As an example, if N1 = 3, N2 = 4, we have M = 12, and

⎛1

⎜0 ⎜0 ⎜ ⎜0 ⎜0 ⎜ ⎜0 P=⎜ ⎜0 ⎜ ⎜0 ⎜0 ⎜ ⎜0 ⎝ 0 0

0 0 0 1 0 0 0 0 0 0 0 0

0 0 0 0 0 0 1 0 0 0 0 0

0 0 0 0 0 0 0 0 0 1 0 0

0 1 0 0 0 0 0 0 0 0 0 0

0 0 0 0 1 0 0 0 0 0 0 0

0 0 0 0 0 0 0 1 0 0 0 0

0 0 0 0 0 0 0 0 0 0 1 0

0 0 1 0 0 0 0 0 0 0 0 0

0 0 0 0 0 1 0 0 0 0 0 0

0 0 0 0 0 0 0 0 1 0 0 0

⎞

0 0⎟ 0⎟ ⎟ 0⎟ ⎟ 0⎟ 0⎟ ⎟, 0⎟ ⎟ 0⎟ 0⎟ ⎟ 0⎟ ⎠ 0 1

Evaluating Eq. (12) in respectively the x and the y directions and using (15), we obtain the integration matrix

F(II1)  P −1 F(I 1) = P −1 A(x,I1) P fII , In the following, we set (II) as the default index of nodes, and denote

A(x1)  P −1 A(x,I1) P,

1) A(y1)  A(y,II ,

to be the first order integration matrices along respectively the x and the y axes.

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250

For multi-layer integration, we define F (2) (x, y) 

a

a

f (ξ1 , y) dξ1 dξ2 and denote

T , FM(2) ] ,

F(2) = [F1(2) , F2(2) , · · · where

Fk(2) = F (2) (xk , yk ), (1)

 x  ξ2

237

 ξ2

F (1) (ξ2 , y) 

a

T

F(1) = [F1(1) , F2(1) , · · · , FM(1) ] .

f (ξ1 , y) dξ1 ,

= F (1) (xk , yk ). We then obtain the following formula for the integration matrices:

Here, Fk

F(2) ≈ A(x1) F(1) ≈ A(x1) A(x1) f, xξ is the integration matrix corresponding to the double-layer integration a a 2 f (ξ1 , y)dξ1 dξ2 . The yx yη (1) (1) (1) (1) integration matrices corresponding to c c 2 f (x, η1 )dη1 dη2 and c a f (ξ , η)dξ dη are respectively Ay Ay and Ay Ax . This implies that all higher order integration matrices can be obtained from the products of lower order integration matrices. (1) (1)

where the product Ax Ax

3. FIM for one-dimensional boundary layer problem To demonstrate the advantage of the FIM for solving singular perturbation problems, we consider the following two-point boundary value problem (BVP):

− u (x) + p(x)u (x) + q(x)u(x) = f (x), u(a) = α ,

u(b) =

x ∈ (a, b),

(16)

β,

(17)

where a, b, α and β are constants with a < b. The coefficient  > 0 is called singularly perturbed parameter whose value is usually very small ( 1). The functions p(x), q(x), and f(x) are assumed to be sufficiently smooth on [a, b]. It is well known that if  1, the solution of the problem (16)–(17) normally forms a boundary layer near the end points a or b or an interior point at which the solution’s derivatives change dramatically within a very narrow interval containing the layer. Based on the last section on the introduction of the FIM, we assume a partition a = x1 < x2 < . . . < xN = b of the interval [a, b]. Applying the integration operation twice on both sides of Eq. (16), we obtain the following integral equation

−



x

 η a

a

u (ξ )dξ dη +





xk

a

 η a

 η

u (ξ ) dξ dη +

p(ξ )u (ξ ) dξ dη +

a

a

At each point x = xk , we have

−

x



 η

xk a

a



 η

x

a

a

p(ξ )u (ξ ) dξ dη +



xk

a

q(ξ )u(ξ ) dξ dη =

 η a



x

 η a

a

q(ξ )u(ξ ) dξ dη =



xk a

f (ξ ) dξ dη.  η a

f (ξ ) dξ dη.

(18)

(19)

The following numerical approximations are obtained from simple integration by part:



xk x1

 η x1

u (ξ ) dξ dη =



xk x1





u (η) − u (x1 ) dη

= u(xk ) − u(x1 ) − u (x1 )(xk − x1 ), 

xk x1

 η x1

p(ξ )u (ξ ) dξ dη =



xk x1

 = ≈

xk x1





η  η

p(ξ )u(ξ )

− p (ξ )u(ξ )dξ dη x1 x1

p(η)u(η) dη −



xk x1

xk x1



xk x1

 η x1

 η x1

q(ξ )u(ξ ) dξ dη ≈

p(x1 )u(x1 ) dη −



xk

 η

x1

x1

p (ξ )u(ξ ) dξ dη

k k   a(ki1) p(xi )u(xi ) − p(x1 )u(x1 )(xk − x1 ) − a(ki2) p (xi )u(xi ), i=1



(20)

(21)

i=1

k  a(ki2) q(xi )u(xi ),

(22)

i=1

f (ξ ) dξ dη ≈

k  a(ki2) f (xi ).

(23)

i=1

Combining Eqs. (19)–(23), we obtain the following matrix system for the approximate solution of Eq. (18):

− u + AP0 u − A2 P1 u + A2 Qu = c0 x + c1 i + A2 f,

(24)

238

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250

where c0 , c1 are integral constants, u = [u1 , u2 , . . . , uN ] , f = [ f1 , f2 , . . . , fN ] , i = [1, 1, . . . , 1] , A = A(1) , and T

T

⎛

p(x1 ) ⎜ 0 ⎜ 0 P0 = ⎜ ⎜ . ⎝ .. 0

⎛

0 p(x2 ) 0 .. . 0

p (x1 ) ⎜ 0 ⎜ 0 P1 = ⎜ ⎜ . ⎝ .. 0

⎛

q(x1 ) ⎜ 0 ⎜ 0 Q=⎜ ⎜ . ⎝ .. 0

0

p (x

2

)

0 q(x2 ) 0 .. . 0

··· ··· ··· .. . ···

0 0 q(x3 ) .. . 0

0 0 ⎟ ⎟ 0 ⎟, .. ⎟ . ⎠ p(xN )

⎞

··· ··· ··· .. . ···

0 0 p (x3 ) .. . 0

0 .. . 0

⎞

··· ··· ··· .. . ···

0 0 p(x3 ) .. . 0

T

0 0 ⎟ ⎟ 0 ⎟, .. ⎟ . ⎠  p (xN )

⎞

0 0 ⎟ ⎟ 0 ⎟. .. ⎟ . ⎠ q(xN )

The boundary conditions (17) now determine the values of u1 and uN as:

u1 = α ,

uN = β .

(25)

As a result, the approximation of the solution u can be obtained from solving Eqs. (24) and (25), which is a linear system of and c0 , c1 : N + 2 equations with N + 2 unknowns {ui }N i=1

G u = F,



(26)



− I + AP0 − A2 P1 + A2 Q E1 E2

−x −i 0 0 0 0

u c0 c1

 =

A2 f

α β

,

(27)

where I is an identity matrix of size N, E1 = [1, 0, . . . , 0]1×N , E2 = [0, 0, . . . , 1]1×N . As mentioned, except the two columns induced from the boundary conditions, the resultant system (27) is a lower triangular matrix whose numerical solution can be efficiently obtained at comparatively very low cost of computation. Fig. 2 shows the sparsity pattern of the coefficient matrix of the linear system by FIM. In our numerical examples, we use the MATLAB built-in command backslash to solve the linear system (27). The computational complexity of Gaussian elimination method is O(N3 ) for general matrix. By using the index arrangement of nodes as shown in Fig. 1, the resultant matrix by FIM is nearly lower triangular as shown in Fig. 2. A total of O(N2 ) operations are needed to transform the resultant matrix into a lower triangular matrix, whose computational complexity is O(N2 ). Therefore, the overall computational complexity of solving (27) is O(N2 ). 4. Numerical error estimate From the derivation of the resultant integration matrix system for the solution, it can be observed that the source of error comes mainly from the numerical quadrature formula, whose convergence order estimation has been well established. In this section, we will derive the convergence error estimate in using the FIM for solving the singular perturbation problems (16) (17) as follow: in each subinterval (xi , xi+1 ), we have



xi+1

xi

p(ξ )u(ξ ) dξ =

p(xi )u(xi ) + p(xi+1 )u(xi+1 ) 3 i − i ( p(x)u(x)) |x=ξi , 2 6

where ξi ∈ (xi , xi+1 ), i = xk+1 − xk , i = 1, 2, . . . , N − 1. From (21) we have



xk

x1

p(ξ )u(ξ ) dξ =

k−1   i=1

=

xi+1 xi

p(ξ )u(ξ ) dξ

k k−1   i 3  a(ki1) p(xi )u(xi ) − ( p(x)u(x)) |x=ξi , 6 i=1

i=1

where ξi ∈ (xi , xi+1 ), i = xi+1 − xi , and



xk

x1

 η x1

p (ξ )u(ξ ) dξ dη =

k−1   i=1

xi+1 xi

 η x1

p (ξ )u(ξ ) dξ dη

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250

239

0 5 10 15 20 25 30 35 40 45 50 0

10

20

30 nz = 1376

40

50

Fig. 2. Visualize sparsity pattern.

=

 k  a(ki1) i=1

xi x1

p (ξ )u(ξ ) dξ −

k−1  i 3 6



x1

i=1



x

p (ξ )u(ξ ) dξ



k i i−1     j3  = a(ki1) a(i j1) p (x j )u(x j ) − ( p (x)u(x)) |x=ξi j 6 i=1

−

j=1

k−1  i=1

=

i

 3 

6



j=1

x x1

p (ξ )u(ξ ) dξ



|x=ηi

i i−1 k k      j3  a(ki1) a(i j1) p (x j )u(x j ) − a(ki1) ( p (x)u(x)) |x=ξi j 6 i=1

j=1

k−1  i 3 − 6



i=1

x x1

i=1

=

|x=ηi

p (ξ )u(ξ ) dξ



j=1

|x=ηi

k k i−1     j3  a(ki2) p (xi )u(xi ) − a(ki1) ( p (x)u(x)) |x=ξi j 6 i=1

−

k−1  i=1

i

3



6

i=1

x x1

j=1

p (ξ )u(ξ ) dξ



|x=ηi ,

where  j = x j+1 − x j , i = xi+1 − xi , and ξi j ∈ (x j , x j+1 ), ηi ∈ (xi , xi+1 ) are dummy variables. Similarly, from (22) we have



xk x1

 η x1

q(ξ )u(ξ ) dξ dη =

k−1   i=1

=

 η

xi

x1

 k  a(ki1) i=1

=

xi+1

xi x1

q(ξ )u(ξ ) dξ dη

q(ξ )u(ξ ) dξ −

k−1  i 3 6 i=1



x x1

q(ξ )u(ξ ) dξ



|x=ηi

k k i i−1      j3  a(ki1) a(i j1) q(x j )u(x j ) − a(ki1) (q(x)u(x)) |x=ξi j 6 i=1

j=1

i=1

j=1

240

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250

−

k−1  i 3 6



x1

i=1

=

x

q(ξ )u(ξ ) dξ



|x=ηi

k k k−1 i−1      j3 i 3  a(ki2) q(xi )u(xi ) − a(ki1) (q(x)u(x)) |x=ξi j − 6 6 i=1

i=1

j=1

xk

 η

x1

f (ξ ) dξ dη =

x1



k k k−1 i−1      j 3  i 3 f (x)|x=ξi j − a(ki2) f (xi ) − a(ki1) 6 6 i=1

i=1

j=1

x x1

i=1

f (ξ ) dξ



x

x1

i=1

and from (23) we have





q(ξ )u(ξ ) dξ



|x=ηi ,

|x=ξi .

From the error estimate on the following small magnitude terms emerged during the integration process:



 k−1 



k−1  i 2

i 3

  ( p(x)u(x)) |x=ξi ≤ ( max |( p(x)u(x)) |) max i

x 6 i

i=1 6

i=1   k i−1   j3  i 2   ≤ (b − a)( max |( p(x)u(x)) |) max , a(ki1) ( p (x)u(x)) |x=ξi j x

≤ ( max |( p (x)u(x)) x



x

p (ξ )u(ξ ) dξ

x1

i=1



|x=ξi

k i−1    j3 a(ki1) (q(x)u(x)) |x=ξi j 6 i=1

j=1

k−1  i 3 6



i=1

x x1

q(ξ )u(ξ ) dξ



|x=ξi

k i−1    j3 a(ki1) ( f (x)) |x=ξi 6 i=1

k−1  i 3 6

j=1



i=1

x x1

f (ξ ) dξ



|x=ξi



  |) max i i

 2

6

i=1

N  a(ki1)



j=1 N−1 

i=1

6



i

j=1

  2 ( max |( p (x)u(x)) |) max i , x 6 i

   k−1 



x  i 2



 ≤ max

p (ξ )u(ξ ) dξ

max i x

6 i

x1 i=1

  

 

x i 2



 ≤ (b − a) max

p (ξ )u(ξ ) dξ

max , x

6 i

x1   i 2 2  ≤ 2(b − a) ( max |(q(x)u(x)) |) max , ≤ 2(b − a)

k−1  i 3 6

6

i

2

x

i

6

 

 

x i 2



≤ (b − a) max

q(ξ )u(ξ ) dξ

max , x

6 i

x1   i 2 2  ≤ 2(b − a) ( max |( f (x)) |) max , 

x

i

6

 

 

x i 2



≤ (b − a) max

f (ξ ) dξ

max , x

6 i

x1 

we obtain that the discretized scheme (24) admits a small magnitude term with the convergence order of O(h2 ), where h = maxi |xi+1 − xi |. Let u denote the exact solution of Eq. (19), i.e., G u = F + δ where δ represents the small magnitude terms during the procedure of integration by parts. If u denote the numerical solution of Eq. (26), the difference u − u can be calculated by

G (u − u) = δ. From Eq. (27) the coefficient matrix, − I + AP0 − A2 P1 + A2 Q, is strictly lower triangular, in which the two rightmost columns 0 0 ) in (27) becomes a triangular matrix as both [−x − i] can be eliminated by elementary column operations. The sub-matrix ( 0 0 E1 , E2 have only one non-zero element. Since the elements on the main diagonal are non-zero, G is non-singular [20]. Finally, we have u − u ≤ G−1  δ = O(h2 ). Remark 4.1. Different from one-dimensional case, where the invertibility of the matrix G can be achieved by elementary row and column operations, the non-singularity of the resultant matrix for two- or higher-dimensional problems is not obvious due to the inclusion of general boundary conditions and free undetermined function terms. Moreover, numerical tests indicate that the resultant matrix is not fully ranked. To avoid singularities under these higher-dimensional cases, we use least square technique in this paper to get the error estimate

 −1 u − u ≤  GT G  GT  δ = O(h2 ),

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250

241

Table 1 Numerical convergence order of FIM for Example 1.

 = 10(−0.5)

=1

N

10 20 40 80 160

 = 10(−1)

 · ∞

Order

 · ∞

Order

 · ∞

Order

1.01E − 04 2.51E − 05 6.29E − 06 1.57E − 06 3.93E − 07

– 2.00 2.00 2.00 2.00

2.48E − 03 6.16E − 04 1.54E − 04 3.85E − 05 9.63E − 06

– 2.01 2.00 2.00 2.00

3.45E − 02 7.87E − 03 1.93E − 03 4.79E − 04 1.20E − 04

– 2.13 2.03 2.01 2.00

where G is an over-determined matrix. 5. Numerical 1D example with single boundary layer Example 1. We first consider the following singular perturbation problem:

 u + u = 0, x ∈ (0, 1),

(28)

u(0) = 0,

(29)

u(1) = 1.

which will form a single boundary layer near the left endpoint when  is small. The exact solution is given by x

1 − e− 

u(x) =

1

1 − e− 

.

We assume a uniform partition 0 = x0 < x1 < . . . < xN = 1, where xi = twice on both sides of Eq. (28), we have

i N,i

= 1, 2, . . . , N. Applying the integration operation

 u + A u = c0 x + c1 i,

(30)

T

T

where u = [u0 , u1 , . . . , uN ] , c0 and c1 are two integral constants and i = [1, 1, . . . , 1] . From the boundary conditions (29), we obtain:

u0 = 0,

uN = 1.

(31)

From Eqs. (30) and (31), we obtain a linear system of N + 3 algebraic equations with N + 3 unknowns u0 , u1 , . . . , uN , c0 and c1 , whose coefficient matrix has the pattern as displayed in Fig. 2 when N = 50. The total number of non-zero elements is 1376 whose values are determined by using the built-in function SPY in MATLAB. It is noted here that the nearly lower-triangular structure of the coefficient matrix ensures the simplicity and efficiency of the FIM. In this 1D numerical computation, we use MATLAB (Gaussian elimination with partial pivoting) for the solution of the linear system of Eqs. (30) and (31). To verify the convergence order of O(h2 ) in our numerical Example 1, we consider the case when  = 1, 10(−0.5) and 10(−1) . Table 1 shows the numerical convergence order when the total number N equals to 10, 20, 40, 80 and 160, respectively. Here, the E log E1

convergence order is computed by log 22 , where E1 and E2 denote the maximum pointwise errors u − u when the number of nodes are equal to N and 2N, respectively. As shown in Table 1, the numerical results agree well with the analytical analysis given in the last section. When the total number of points is N = 100 and the perturbed parameter is  = 10−4 . The numerical oscillation error plot displayed in Fig. 3 indicates the well known difficulty in solving these kinds of singular perturbation problems by numerical methods such as FDM and FEM using uniform grids/meshes. To overcome this large oscillation near boundary layer, some kinds of optimal point placement strategies were developed [16,17] among which the one introduced by Shishkin [9] is to subdivide the interval [0, 1] into two regions [0, x N ] and [x N , 1], 2

where x N is the transition point determined by 2

x N = min

1

2

2

,

2

  ln N .

The total number of N points are uniformly located in these two subintervals as:

⎧ ⎪ ⎨ (i − 1) h1 , xi =  ⎪ ⎩ x N + i − N h2 , 2

i = 1, 2, . . . , i=

2

where h1 , h2 are defined by

h1 =

xN 2

N 2

−1

,

h2 =

1 − xN 2

N 2

.

N , 2

N + 1, . . . , N, 2

(32)

242

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250

2 FIM Exact solution

1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

0

0.2

0.4

0.6

0.8

1

Fig. 3. FIM on uniform points.  = 10−4 , N = 100.

1.4 FIM Exact solutions 1.2

1

0.8

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

Fig. 4. Numerical solution by FIM on Shishkin points when  = 10−4 , N = 100.

We then recompute the solution of the boundary layer problem (28) and (29) by using FIM on these Shishkin points. As shown in Figs. 4 and 5, the numerical oscillation error has been greatly reduced. The result is compared with the solutions given by upwind central difference FDM and the multiquadric collocation method (MQCM) in the recent paper of Hon [11] where the total number of points are chosen to be 40 and 160 respectively. Numerical comparison results displayed in Table 2 indicate that the MQCM gives much better approximation than the FDM for all different values of  and slightly better than the FIM when  is large. This is due to the spectral convergence of the MQCM. However, when  is small, the FIM gives superior results than the

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250

243

1 FIM Exact solutions

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

1

2

3

4

5

6 −4

x 10 Fig. 5. Numerical solution near the boundary layer.

Table 2 Comparison errors at point x N on Shishkin points. 2



N = 40

0.1 0.05 0.025 0.0125 0.00625 0.003125 0.0015625 0.00078125 0.000390625

N = 160

FDM

MQCM

FIM

FDM

MQCM

FIM

1.89E − 02 3.48E − 02 6.48E − 02 1.19E − 01 2.11E − 01 3.45E − 01 5.10E − 01 6.69E − 01 7.94E − 01

8.16E − 05 1.43E − 04 4.08E − 03 1.29E − 02 2.90E − 02 5.91E − 02 8.13E − 02 1.38E − 01 3.12E − 01

2.87E − 04 2.90E − 04 2.90E − 04 2.90E − 04 2.90E − 04 4.39E − 04 2.24E − 03 7.22E − 03 1.33E − 02

1.47E − 03 2.39E − 03 4.20E − 03 7.78E − 03 1.49E − 02 2.88E − 02 5.54E − 02 1.04E − 01 1.89E − 01

1.86E − 08 7.08E − 06 7.31E − 06 3.41E − 05 1.07E − 04 4.27E − 04 8.81E − 04 2.62E − 03 4.58E − 03

1.11E − 05 1.09E − 05 1.09E − 05 1.09E − 05 1.09E − 05 1.09E − 05 1.09E − 05 1.09E − 05 1.12E − 05

Table 3 Comparison between results by MQCM and FIM when  = 10−4 . MQCM

FIM

N

50

100

150

100

200

300

Error at x N 2 Time (s)

1.32E − 02 4.76E − 03

3.40E − 03 1.69E − 02

6.66E − 04 3.74E − 02

3.68E − 03 3.93E − 03

9.71E − 03 8.37E − 03

2.73E − 06 2.18E − 02

MQCM because the MQCM suffers from the well known problem of ill-conditioning due to its full coefficient matrix and optimal choice of shape parameter. The superior advantage of FIM over MQCM is further illustrated by making the comparison between MQCM and FIM under different total number of nodes N and error differences at point xN/2 . The computations are implemented on Shishkin points when  = 10−4 . It can be observed from Table 3 that FIM requires less CPU time than MQCM for the same accuracy. This verifies our claim that FIM gives a stable, efficient and accurate approximation to the solutions of these kinds of singularly perturbed problems. Although the numerical results have shown to be much improved by using Shishkin points, the construction of these points requires a good priori knowledge on the location and width of the boundary layer, which in general are not known. In the

244

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250 Table 4 The number of points N and error difference  · ∞ under varied δ .



δ = 10−1

−4

10 10−6

δ = 10−2

δ = 10−3

δ = 10−4

N

 · ∞

N

 · ∞

N

 · ∞

N

 · ∞

37 73

5.91E − 03 5.42E − 03

57 100

7.50E − 04 6.31E − 04

166 177

6.96E − 05 6.16E − 05

457 487

6.13E − 06 5.10E − 06

Table 5 Iterative results by AFIM when  = 10−8 . Iterations N

20 58

200 103

210 112

230 162

390 488

700 1110

 · ∞

1.00+00

1.60E − 03

4.34E − 04

6.05E − 05

8.77E − 06

9.45E − 07

following section, we devise an adaptive scheme in locating more points near the boundary layer for solving the singular perturbation problems (28) and (29) when the location and width of the layer are not known. 6. Adaptive FIM Based on the work presented in [8], we devise in this section an adaptive finite integration method (AFIM) for solving the singularly perturbed problems when the location and width of the boundary layer are both not known. The main idea is to iteratively place more points near the layer where the error differences between the numerical approximation and the exact solution are comparatively large; and remove points where the errors are small. The procedure is given as follows. Step 0. Set values of total number of points N, coefficient of proportionality τ , error tolerance δ and initial partition P0 : 0 = N N x1 < x2 < . . . < xN = 1. Let {uk }k=1 denote approximations of {u(xk )}k=1 , where u(xk ) is the exact solution’s value at point xk . Define

Dk,1 =

uk − uk−1

k−1

,

Dk,2 =

Dk+1,1 − Dk,1 , k−1 + k

to be the approximation’s first and second divided differences, respectively, where k = xk+1 − xk . The maximum error maxk |uk − u(xk )| can be estimated by the indicator

M[xk−1 , xk , xk+1 ] = (2k−1 + 2k )(|Dk,2 | + |Dk,1 | + |uk |). Step 1. For all k = 2, 3, . . . , N − 1, compute {M[xk−1 , xk , xk+1 ]}k=2 and mean Mμ of these indicators. N

Step 2. If maxk M[xk−1 , xk , xk+1 ] is attained at k = i, then add two points xk−1 +xk xk +xk+1 , 2 2

xi−1 +xi xi +xi+1 , 2 2

on both sides of xi .

on both sides of xk . Step 3. If M[xk−1 , xk , xk+1 ] > τ Mμ , then add two points Step 4. If M[xk−1 , xk , xk+1 ] < Mμ /τ , M[xk−2 , xk−1 , xk ] < Mμ /τ (3 ≤ k ≤ N − 1) and M[xk , xk+1 , xk+2 ] < Mμ /τ (2 ≤ k ≤ N − 2), then remove the point xk . x −xk−1 x +x x −xk−1 x +x < 13 , then add a point k 2k+1 between points xk and xk+1 . If xk −x > 3, then add a point k−12 k between Step 5. If xk −x k+1

k

k+1

k

points xk−1 and xk . Stop adding points until





1

xk − xk−1

≤ ≤3 3 xk+1 − xk

is satisfied at all points xk . Step 6. Update the total number of points N to include all the added and removed points. Repeat Step 1 to 5 using the modified partition until

max M[xk−1 , xk , xk+1 ] < δ k

is achieved. Denote  · ∞ to be the maximum error difference between the numerical approximation and the exact solution over all N points {xk }k=1 . The recomputed results of the problem (28) and (29) by AFIM with initial N = 20, τ = 20 are displayed in Table 4. For comparison with the results given in [11], we display in Table 5 the  · ∞ error, total number of iterations and total number of points N when  = 10−8 . Here, we use initial N = 20, τ = 20. It can be observed from Table 5 that an order of accuracy O(10−4 ) is achieved by AFIM when N = 112, whereas a total of 300 points were needed by adaptive MQCM given in [11].

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250

245

Table 6 Comparisons between FDM and FIM when  = 10−6 . N1

100 200 400 500 1000

FDM

FIM

L −  · ∞

R −  · ∞

 · ∞

– 3.44E − 04 8.61E − 05 – 1.38E − 05

– 6.88E − 04 1.72E − 04 – 2.76E − 05

5.98E − 04 1.48E − 04 3.68E − 05 2.35E − 05 5.86E − 06

Table 7 Error in  · ∞ for different N1 and  by FIM. N1

 = 10−6

 = 10−8

 = 10−10

 = 10−12

250 500 1000

9.44E − 05 2.35E − 05 5.86E − 06

1.68E − 04 4.18E − 05 1.04E − 05

2.62E − 04 6.53E − 05 1.63E − 05

3.78E − 04 9.40E − 05 2.35E − 05

7. Numerical 1D example with double boundary layers Example 2. To further illustrate the efficiency and accuracy of the proposed FIM and AFIM, we consider the following double boundary layers problem given in [13]:

 u − u = 0, x ∈ (0, 1),

(33)

u(0) = 1,

(34)

u(1) = 2.

whose exact solution is −1

u(x) =

−1

x

−x

(e √ − 2)e √ + (2 − e √ )e √ e

−1 √



−e

√1



.

For small value of the singular perturbation parameter  , the exact solution can be approximated by

u(x) ≈ 2e as e

−1 √



≈ 0 and

−(√ 1−x)

+e



√1

2−e 

−1 √

√1

e  −e 

−x √



,

≈ 1.

We remark here that Zhang in his recent paper [13] divided this double boundary layers problem into three boundary value problems (BVPs) with approximations obtained respectively by second order FDM on piecewise uniformly distributed points. A good approximation of this approach requires a good priori estimate on the location and width of boundary layers. Here, we apply both FIM and AFIM on the whole interval to solve the double layers problem (33) and (34) without the need to solve three separate BVPs. In our computation, we divide the whole interval [0, 1] into three subintervals

[0, t1 ], √

[t1 , t2 ],

[t2 , 1],

√ where t1 =   , t2 = 1 − t1 = 1 −   . In each subinterval, the points are chosen to be uniformly distributed by placing a total of N1 points in [0, t1 ]; N2 points in [t1 , t2 ]; and N3 points in [t2 , 1]. Due to the boundary condition u(1) = 2 u(0), we let N3 = 2 N1 to allow more points placed in [t2 , 1] for the stiffer right-end boundary layer. The FIM is then applied on this piecewise uniform partition with varied N1 , fixed N2 = 10 and  = 10−6 . The comparison result displayed in Table 6 shows that the FIM is slightly better than the FDM given in [13]. Here L·∞ , R·∞ denote respectively the maximum error differences at the left- and the right-end boundary layers. For smaller singular perturbation parameter  and number of points N1 with fixed N2 = 10, it can be seen from Table 7 that the FIM still works very well for extremely small  = 10−12 whereas the FDM fails to produce acceptable approximation. Furthermore, we apply the AFIM with initial total number of points N = 20, τ = 20, and display the numerical results in Table 8 under different error tolerance δ against different perturbed parameters  . The numerical results shown in Table 8 further verifies that the proposed AFIM works very well even for double boundary layers problem with very small  . Compared with the numerical accuracy of O(10−4 ) when  = 10−12 obtained in [18] using adaptive radial based functions collocation method and sine-transformation, the AFIM gives better result at lower cost of computation.

246

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250 Table 8 Error in  · ∞ under different parameters’ values.



δ = 10−1

−6

10 10−8 10−10 10−12

δ = 10−2

δ = 10−3

δ = 10−4

N

 · ∞

N

 · ∞

N

 · ∞

N

 · ∞

43 59 62 78

5.46E − 03 5.88E − 03 5.49E − 03 5.37E − 03

94 100 119 126

5.73E − 04 7.65E − 04 7.32E − 04 6.31E − 04

244 253 269 276

6.70E − 05 6.94E − 05 6.47E − 05 5.74E − 05

729 727 741 766

7.86E − 06 6.68E − 06 5.37E − 06 7.61E − 06

Fig. 6. L-shaped domain.

8. Least Squares FIM for 2D boundary layer problem In this section, we combine the least squares technique with the FIM for solving higher-dimensional problems. Although the example is given in domain under rectangular axes, the method is readily applicable to problems defined on other type of domains, such as circle, sector, and annular domain by simple variable substitution or coordinate transformation techniques. We first consider the following BVP for illustration: Example 3. Consider the following two-dimensional Poisson problem defined in a L-shaped domain ⊂ [0, 1]2 :

uxx + uyy = f,

in ;

u|∂ = 0 on ∂ ,

(35)

where f is given so that the problem (35) admits an exact solution u(x, y) = x(x −

1 1 2 )(1 − x)y( 2

− y)(1 − y).

To apply the FIM for the solution of the above two-dimensional problem, we choose a uniform nodes distribution {(xi , yi )} j=1 as shown in Fig. 6. Applying double integral operators for x and y variables respectively on both sides of Eq. (35), we obtain the following equivalent integral form:



y 0

 η2 0

=

u(x, η1 ) dη1 dη2 +



y  η2  0

0

x  ξ2 0

0



x

 ξ2

0

0

M

u(ξ1 , y) dξ1 dξ2

f (ξ1 , η1 ) dξ1 dξ2 dη1 dη2 + xζ0 (y) + ζ1 (y) + yθ0 (x) + θ1 (x)

(36)

where ζ 0 (y), ζ 1 (y), θ 0 (x), θ 1 (x) are four one-dimensional undetermined functions arise from integration. To fix the values of these four functions, we use the standard cubic spline interpolation technique in solving the two-dimensional singularly perturbed problems (36) as follows: let {x∗1 , x∗2 , · · · , x∗r } and {y∗1 , y∗2 , · · · , y∗s } be the interpolation centers for interpolating θ 0 (x), θ 1 (x) and ζ 0 (y), ζ 1 (y):

θ0 (x) =

s  λ(j0) φ j (x), j=1

θ1 (x) =

s  λ(j1) φ j (x), j=1

ζ0 (y) =

s  ρ (j 0) φ j (y), j=1

ζ1 (y) =

s  ρ (j 1) φ j (y), j=1

(37)

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250

247

0

10

20

30

40

50

0

10

20 nz = 460

30

Fig. 7. Visualize sparsity pattern.

where φ j = r3j , r j = x − x∗j  or y − y∗j , and  ·  is the usual Euclidean norm. Eq. (37) can be written in the linear system form as:

θ0 = x (0) , θ1 = x (1) , ζ0 = y ρ(0) , ζ1 = y ρ(1) , in

which (1)

x = (φ j (xi ))M×s , y = (φ j (yi ))M×r

(1)

(1) T

(0)

are

coefficient

(0) T

(0)

(1)

(1)

matrices

and

T

(0) = [λ(10) , λ(20) , · · · , λ(s0) ] , (1) =

(1) T

[λ1 , λ2 , · · · , λs ] , ρ(0) = [ρ1 , ρ2 , · · · , ρr ] , ρ(1) = [ρ1 , ρ2 , · · · , ρr ] are vectors of unknowns to be determined. Replacing the integral operators by corresponding integration matrices, Eq. (36) becomes

(A(x1) A(x1) + A(y1) A(y1) )u + Xy ρ(0) + y ρ(1) + Yx (0) + x (1) = A(y1) A(y1) A(x1) A(x1) f,

(38)

where u = [u1 , u2 , · · · , uM ] is the unknown solution in vector form; X = diag(x1 , x2 , · · · , xM ), Y = diag(y1 , y2 , · · · , yM ) are diagT onal matrices; and f = [ f (x1 , y1 ), f (x2 , y2 ), · · · , f (xM , yM )] is the vector of given function values. Assume that there are a total of 2(N1 + N2 ) − 4 boundary nodes. There exists a matrix B of size [2(N1 + N2 ) − 4] × M with elements being either 0 or 1 such that T

B u = 0,

(39)

where 0 is a column vector with all 2(N1 + N2 ) − 4 elements equal to 0. From (38) and (39), we obtain the following linear system of equations for the numerical solution u:

⎛



A(x1) A(x1) + A(y1) A(y1) B

Xy

y O

Yx

u

⎞

 ⎜ ρ(0) ⎟   x ⎜ ρ(1) ⎟ = A(y1) A(y1) A(x1) A(x1) f . ⎟ ⎜ 0 ⎝(0) ⎠ ( 1) 

(40)

The size of the coefficient matrix of the above linear system is (M + 2(N1 + N2 ) − 4) × (M + 2(r + s)). If we choose r = N1 − 1, s = N2 − 1, the coefficient matrix will be a square matrix. As shown in Fig. 7 when N1 = N2 = 6, the sub-matrices (1) (1) (1) (1) (Ax Ax + Ay Ay ) displayed by using the SPY function has a total number of 460 non-zero elements and is clearly overB determined with column fully ranked. Without careful manipulation, the square coefficient matrix of this linear system (40) will be ill-conditioned. In our computation, we choose the values of r and s so that the resultant linear system is

248

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250 Table 9 Numerical convergence order of FIM for Example 3. r =s=N−1

N

11 21 41 81

r =s=N−2

r =s=N−3

r =s=N−4

 · ∞

Order

 · ∞

Order

 · ∞

Order

 · ∞

Order

2.04E − 04 4.09E − 05 9.95E − 06 2.52E − 06

– 2.32 2.04 1.98

3.26E − 04 5.69E − 05 1.03E − 05 2.48E − 06

– 2.52 2.46 2.06

4.68E − 04 8.46E − 05 1.31E − 05 2.51E − 06

– 2.47 2.69 2.38

7.11E − 04 1.19E − 04 1.79E − 05 2.68E − 06

– 2.58 2.73 2.74

1 9

3

6

2

5

7

1

4

7

1 d

d

0

1 d

d

1

Fig. 8. Domain partition.

over-determined and combine the least squares technique with FIM to a new least squares finite integration method (LSFIM) to obtain approximation for the solution of this over-determined system. The advantage of using this LSFIM is to allow higher flexibility to determine optimally the four free terms of integration functions that satisfy the given boundary condition. As an illustration, we let N1 = N2 = N, r = s ≤ N − 1. We simply use the backslash command in MATLAB to solve the overdetermined system (40) and obtain the least squares approximation of the solution. The comparison of the numerical solution with the analytical solution on maximum pointwise errors for different number of nodes is displayed in Table 9 from which it can be verified that the LSFIM provides a very accurate approximation to the 2D Poisson problem on L-shape domain. Finally, we apply this LSFIM to solve the following two-dimensional singular perturbation problem. Example 4. Consider



− 2

 ∂ 2u ∂ 2u + + 2u = f (x, y, ), in = (0, 1) × (0, 1), ∂ x2 ∂ y2

u = 0,

on

(41)

∂ ,

(42)

where the function f(x, y,  ) is chosen so that its exact solution is given by



u(x, y) =

1−

x

e−  + e−

1−x

1

1 + e− 





1−

y

e−  + e−

1−y

1

1 + e− 





+ x(1 − x)y(1 − y).

This is a typical 2D singularly perturbed problem with exponential boundary layers at all sides of the domain . Suppose that there are a total of N1 × N2 nodes in [0, 1] × [0, 1]. The original domain is partitioned into nine nonoverlapping sub-domains i , 1 ≤ i ≤ 9, i.e. = ∪9i=1 i as shown in Fig. 8, where d = 2| ln |. Each sub-domain i contains N

N

Nix × Niy uniformly distributed nodes. For simplicity, we take Nix = 31 , Niy = 32 and let r = N1 − 3, s = N2 − 3 so that 3r , 3s nodes lie in [0, d], [d, 1 − d] and [1 − d, 1] along x and y axes respectively. For comparison with the results obtained by FEM and RBF given in [19], we display the maximum pointwise errors in Table 10 and profiles of numerical pointwise error differences in Fig. 9 for different N1 = N2 = 24, 48 and 96, r = s = N1 − 3 = N2 − 3, and  = 10−2 .

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250

Fig. 9. Profile of pointwise errors under different nodes arrangements when  = 10−2 .

249

250

D.F. Yun et al. / Applied Mathematics and Computation 271 (2015) 232–250 Table 10 Comparison of maximum pointwise error under different number of nodes when  = 10−2 . M

FEM

RBF

LSFIM

24 × 24 48 × 48 96 × 96

2.56E − 01 1.33E − 01 6.60E − 02

9.00E − 02 6.00E − 02 2.50E − 02

6.92E − 02 9.97E − 03 4.34E − 03

9. Conclusion In this paper we combine the recently developed finite integration method (FIM) with the technique of least squares to solve higher-dimensional singular perturbation problems with multiple boundary layers. An adaptive FIM (AFIM) is devised to better approximate the solution without a priori knowledge on the location and width of the boundary layer. Theoretical analysis is investigated and numerical comparisons with the exact solution and several existing numerical methods (FDM, FEM, MQCM) are made to verify that both the LSFIM and AFIM schemes offer a very stable, highly accurate, and efficient approach for solving these kinds of singularly perturbed problems with boundary layers. The FIM is similar to the mesh-dependent FDM with the intrinsic difference in approximating the solution and its derivatives by integration instead of finite quotient formula. This completely eliminates the well-known optimal roundoff and discretization errors problem in using FDM. The use of higher quadrature rule will give a better convergent FIM scheme although the use of the simplest composite trapezoidal rule in this paper already demonstrates its superior advantage over the traditional FDM and FEM. Furthermore, the use of radial basis functions (RBFs) in the FIM will result into a meshless computational method but then the resultant matrix will be full instead of nearly lower triangular. The open problem of choosing optimal shape parameter in using the radial basis functions will still hinder the development of FIM(RBF) method. Acknowledgments The work described in this paper was partially supported by a grant from the Research Grant Council of the Hong Kong Special Administrative Region, China (Project no. CityU 101211); the Shanxi Hundred Talent Plan; and the 2013 Youth Team Startup Project of Taiyuan University of Technology (Project no. 2013T056). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20]

J.D. Lambert, Numerical Methods for Ordinary Differential Systems: The Initial Value Problem, John Wiley & Sons, New York, 1991. E. Hairer, Solving Ordinary Differential Equations, vols. 1 and 2, Springer-Verlag, Berlin, New York, 1993. W.E. Boyce, R.C. DiPrima, Numerical Methods for Ordinary Differential Systems: The Initial Value Problem, John Wiley & Sons, New York, 2001. S. Roberts, J. Shipman, Two-point Boundary Value Problems: Shooting Methods, Elsevier, New York, 1972. M. Scott, H. Watts, Computational solution of linear two point boundary value problems via orthonormalization, SIAM J. Numer. Anal. 14 (1977) 40–70. J.J.H. Miller, E. O’Rioda, G.I. Shishkin, On piecewise-uniform meshes for upwind- and central-difference operators for solving singularly perturbed problems, IMA J. Numer. Anal. 15 (1995) 89–99. M.G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations Convection–Diffusion and Flow Problems, Springer Series in Computational Mathematics, Springer-Verlag, USA, 1996, p. 24. B. Kreiss, H.O. Kreiss, Numerical methods for singular perturbation problems, SIAM J. Numer. Anal. 18 (1981) 262–276. G. Shishkin, Grid approximation of singularly perturbed boundary value problems with a regular boundary layer, Sov. J. Numer. Anal. Math. Model. 4 (1989) 397–417. G. Sun, M. Stynes, Finite-element methods for singularly perturbed high-order elliptic two-point boundary value problems, II: convection–diffusion-type problems, IMA J. Numer. Anal. 15 (1995) 197–219. Y.C. Hon, Multiquadric collocation method with adaptive technique for problems with boundary layer, Int. J. Appl. Sci. Comput. 6 (3) (1999) 173–184. L. Ling, M.R. Trummer, Multiquadric collocation method with integral formulation for boundary layer problems, Comput. Math. Appl. 48 (56) (2004) 927– 941. W.Q. Zhang, Numerical solutions of singular perturbation problems with multiple boundary layers and interior layers, Int. J. Pure Appl. Math. 83 (4) (2013) 549–558. M. Li, Y.C. Hon, T. Korakianitis, P.H. Wen, Finite integration method for nonlocal elastic bar under static and dynamic loads, Eng. Anal. Boundary Elem. 37 (5) (2013) 842–849. P.H. Wen, Y.C. Hon, M. Li, T. Korakianitis, Finite integration method for partial differential equations, Appl. Math. Model. 37 (24) (2013) 10092–10106. E.C. Gartland, Graded-mesh difference scheme for singularly perturbed two-point boundary value problem, Math. Comput. 51 (1995) 631–657. R. Vulanovic, On a numerical solution of a type of singularly perturbed boundary value problem by using a special discretization mesh, Univ. u Novorn Sadu Zb. Rad. Prirod. Mat. Fak. Ser. Mat. 13 (1983) 187–201. L. Ling, M.R. Trummer, Adaptive multiquadric collocation for boundary layer problems, J. Comp. Appl. Math. 188 (2006) 265–282. X. Zhou, Y.C. Hon, J.C. Li, Overlapping domain decomposition method by radial basis functions, Appl. Numer. Math. 44 (2003) 241–255. Y. Li, M. Li, Y.C. Hon, Improved finite integration method for multi-dimensional nonlinear burgers equation with shock wave, Neural Parallel Sci. Comput. 23 (2015) 63–86.