High-order implicit finite difference schemes for the two-dimensional Poisson equation

Applied Mathematics and Computation 309 (2017) 222–244

Contents lists available at ScienceDirect

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

High-order implicit finite difference schemes for the two-dimensional Poisson equation Miguel Uh Zapata a,∗, Reymundo Itzá Balam b a b

CONACYT - Centro de Investigación en Matemáticas A. C., CIMAT Unidad Mérida, Merida, Yucatan, Mexico Facultad de Ciencias, Universidad Autónoma de Mexico (UNAM), Ciudad Universitaria-México D.F., Mexico City, Mexico

a r t i c l e

i n f o

Keywords: 2D Poisson equation High-order accuracy Implicit finite difference Compact difference scheme Wave plane theory

a b s t r a c t In this paper, a new family of high-order finite difference schemes is proposed to solve the two-dimensional Poisson equation by implicit finite difference formulas of (2M + 1 ) operator points. The implicit formulation is obtained from Taylor series expansion and wave plane theory analysis, and it is constructed from a few modifications to the standard finite difference schemes. The approximations achieve (2M + 4 )-order accuracy for the inner grid points and up to eighth-order accuracy for the boundary grid points. Using a Successive Over-Relaxation method, the high-order implicit schemes have faster convergence as M is increased, compensating the additional computation of more operator points. Thus, the proposed solver results in an attractive method, easy to implement, with higher order accuracy but nearly the same computation cost as those of explicit or compact formulation. In addition, particular case M = 1 yields a new compact finite difference schemes of sixth-order of accuracy. Numerical experiments are presented to verify the feasibility of the proposed method and the high accuracy of these difference schemes. © 2017 Elsevier Inc. All rights reserved.

1. Introduction In this paper, we seek high-order accuracy numerical solutions of the Poisson equation

 p( x ) = f ( x ) ,

x∈

(1)

where  is a two-dimensional rectangular domain with Dirichlet boundary conditions defined on ∂ , and  is the Laplacian operator. The solution p and the forcing function f are assumed to be sufficiently smooth and to have the required continuous partial derivatives. When f (x, y ) = 0, Eq. (1) becomes the Laplace equation. The Poisson and Laplace equations are partial differential equations with broad applications in different fields, such as computational fluid dynamics, Structural Mechanics, theoretical physics, etc. Eq. (1) is often used to describe equilibrium phenomena for many variables such as pressure, water surface elevation, temperature and concentration [1] and its numerical solution is fundamental for the computational simulation of the corresponding applied problems. A large amount of work has been done in the past decades to develop numerical methods for the solution of the Poisson equation producing accurate results. It is known that the standard second-order finite difference method has to take small mesh sizes for obtaining desirable accuracy. Therefore, in order to obtain satisfactory numerical results with reasonable computational ∗

Corresponding author at: CIMAT, PCTY, Sierra Papacal, CP 97302, Mérida, Yucatán, México. E-mail addresses: [email protected] (M. Uh Zapata), [email protected] (R. Itzá Balam).

http://dx.doi.org/10.1016/j.amc.2017.04.006 0 096-30 03/© 2017 Elsevier Inc. All rights reserved.

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

223

cost, a reasonable approach is to develop a higher-order finite difference method. As one of the most effective numerical implementations, the fourth-order compact nine-point-finite difference schemes have received much attention due to their advantages in solution accuracy and the time efficiency. We refer to [2–7], among others. On the contrary, few sixth-order schemes have been developed to compute the solution. Previously, Chu and Fan [8] proposed a three point combined compact difference scheme with a Hermitian polynomial approximation for solving a special two dimensional convection diffusion equation. Sun and Zhang [9] proposed a sixth-order explicit finite difference discretization strategy for solving the 2D convection diffusion equation based on Alternating Direction Implicit Method and the Richardson Extrapolation Technique. Nabavi et al. [10] proposed a new sixth-order accurate compact finite difference method for solving the Helmholtz equation. The work of Wang and Zhang [11] derived an implicit sixth-order compact finite difference scheme for 2D Poisson equation using Taylor series expansions. One of the most recent sixth-order discretization method was proposed by Zhai et al. [12]. They choose a special dual partition and employ Lagrange interpolation and Simpson integral formula to derive difference schemes. In this work, we developed high-order finite difference schemes (fourth-, sixth-, eighth-order, etc.) by using the idea of implicit schemes and any number of stencil operator points for resolving the Poisson Eq. (1). The scheme is referred to as an implicit high order compact scheme because the first derivative and the second derivative of the dependent variables are computed at the same time. The algorithms that we propose have several advantages. First, the implicit finite difference scheme is accurate up to the eight-order and thus provides precise numerical results. Moreover, the scheme is more accurate if the solution of the problem is not significantly affected by the boundary conditions. The general algorithm is not a compact scheme, however the proposed high-order implicit schemes have a better rate of convergence as more operator points are considered. As a consequence, the smaller number of iteration steps compensates for the required computational time generated by the arithmetic operations associated with the additional nodes. Thereby making our algorithm competitive in comparison with the compact schemes for solving the Poisson equation. Furthermore, compact schemes of fourth- and sixth-order can be derived as special cases. Lastly, the proposal scheme is written such that an exact knowledge of the forcing function and its corresponding partial derivatives is not required, as other authors’ formulas [12]. Instead, the scheme requires the evaluation of the forcing function at specific mesh points, as it happens in most practical cases. This paper is organized as follows. The implicit finite difference formulas will be derived first. Subsequently, the accuracy of the scheme will be analyzed. The high-order implicit finite difference scheme for the Poisson equation will be described in Section 3. In Section 4, two numerical examples will be presented to examine the accuracy and effectiveness of the proposed scheme. The results will be compared with those from the explicit standard difference and compact nine-pointdifference schemes. Finally, conclusions will be provided in Section 5. 2. Implicit finite difference formulas This section shows the approximations of the second-order derivative of univariate functions p by implicit finite differences of high orders based on classical plane wave theory. These approximations will be used to obtain finite difference schemes of high-order for the one- and two-dimensional Poisson equation. 2.1. Explicit finite difference formula Let us consider the classic one-dimensional finite difference approximation for a second derivative of a real-valued function p. Given a small value h > 0, the second-order derivative at x can be approximated by the following centered finite difference operator

δx2 p =

1 [ p(x + h ) − 2 p(x ) + p(x − h )], h2

(2)

where the approximation error is to the second-order of accuracy. In order to improve the accuracy of (2), we define an approximation to the second derivative by taking more points with appropriate weighting. Thus for a fixed positive integer M, the second-order derivative of a function p can be approximated by a formula of 2M + 1 nodes as

∂2 p δ2 p ( x ) = 2 ( x ) + O ( h2M ) , ∂ x2 δx

(3)

where the explicit finite difference (EFD) of (2M)th-order of accuracy is given by M δ2 p 1  (x ) = 2 cm p(x + mh ), 2 δx h m=−M

(4)

where values cm are finite difference coefficients which can be easily determined using wave plane theory [13]. For each number m ∈ {1, 2, . . . , M}, coefficients cm and c−m correspond to the same values, resulting in a symmetric finite difference formula (4).

224

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

2.2. Implicit finite difference formula Following the work of Claerbout [14], the accuracy of centered finite difference approximation (2) can be improved not only by considering more points but also by adding a higher-order term to the operator as follows



1 + bh2

 ∂ 2 ∂ 2 p δ2 p ≈ , ∂ x2 ∂ x2 δ x2

(5)

where b is an adjustable constant and operator δδx2 is given by (4). This formulation is called implicit because approximations to the second-order derivative appear in both sides of Eq. (5). Thus, using centered operators (2) at the left-hand side of (5) yields an implicit finite difference (IFD) formula of (2M + 2 )th-order of accuracy given by equation 2



1 + bh2 δx2

M  ∂2 p 1  (x ) = 2 cm p(x + mh ) + O(h2M+2 ). 2 ∂x h m=−M

(6)

The value b and coefficients cm are also determinate by wave plane theory [13]. Notice that the accuracy of IFD formula (6) is increased by 2 compared with EFD formula (4) using the same number of operator nodes. 2.3. High-order implicit finite difference formula Motivated by the idea of increasing the operator order, this paper also proposes an implicit formulation given by the equation



1 + bh2

 2 4 ∂2 ∂ p δ2 p 4 ∂ + dh ≈ , 2 4 ∂x ∂ x ∂ x2 δ x2

(7)

where the new term is added to obtain even higher order of accuracy. Thus, a new high-order implicit finite difference (HIFD) formula of β th-order of accuracy is given by



1 + bh2 δx2 + dh4 δx4

M ∂2 p 1  ( x ) = cm p(x + mh ) + O(hβ ), ∂ x2 h2 m=−M

(8)

where operator δx2 is defined by (2) and δx4 = δx2 δx2 . It can be seen that δx4 is equivalent to the centered finite difference approximation

δx4 p =

1 [ p(x + 2h ) − 4 p(x + h ) + 6 p(x ) − 4 p(x − h ) + p(x − 2h )]. h4

(9)

Not unlike the IFD formula, the finite difference coefficients and order of accuracy of the method will be determined according to the plane wave theory as described below. It will also be shown that β = 2M + 4. It improves the order of accuracy of the EFD formula by 4. However the new HIFD formula requires two more nodes on the left-hand side in contrast to the IFD formula for which those nodes are not required. 2.3.1. HIFD coefficients The first step in deriving the coefficients of HIFD method (8)–(9) is to consider the function p as a Fourier mode

p(x ) = p0 eikx ,

√ where p0 is a constant value, i = −1 and k represents the wavenumber. Substituting (10) into Eq. (8), we obtain



−k



1 − (4b − 16d ) sin

2

2

kh 2







2 − 4d sin (kh ) p(x ) = 2 h 2

(10)



 1 c0 + cm cos (mkh ) p(x ), 2 M

(11)

m=1

Comparing the coefficients of function p at (11) and using Taylor series expansions of the trigonometric functions yields

−2α 2 +

∞ 

M ∞ M       1 ψn−1 −4b − 4d 22n−2 − 4 α 2n = c0 + cm + ψn m2n cm α 2n ,

2

n=2

m=1

n=1

(12)

m=1

where α = kh/2 and

ψn =

(−1 )n 22n . ( 2 n )!

In order for Eq. (12) to be satisfied, it is necessary that M  m=1

2n

m cm =

⎧ 1 ⎪ ⎨− 2 c0 ,

n = 0,

1,

n = 1,

−λn b − μn d,

n ≥ 2,

⎪ ⎩

(13)

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

225

Table 1 Coefficients and order of accuracy of the EFD, IFD and HIFD method.

EFD

IFD

HIFD

where

M

Order

c0

c1

1 2 3 4 1 2 3 4 1 2 3 4

2 4 6 8 4 6 8 10 6 8 10 12

−2 −5/2 −49/18 −205/72 −2 −17/10 −225/151 −319/237 −2 −53/42 −185/187 −281/340

1 4/3 3/2 8/5 1 4/5 21/32 124/223 1 32/63 85/257 67/294

c2

c3

c4

−1/12 −3/20 −1/5

1/90 8/315

−1/560

1/20 51/560 121/992

−1/438 −1/192

31/252 49/304 101/564

1/5274

1/367 1/157

−1/11554

b

d

1/12 2/15 9/56 8/45 1/12 13/63 6/23 75/251

−1/240 1/164 9/674 5/257

  λn = −2n(2n − 1 ), and μn = −2n(2n − 1 ) 22n−2 − 4 .

Finally, using Eq. (13) for n = {1, 2, . . . , M + 2} the problem results in the following linear system

⎡

( 12 )1 ⎢ ( 12 )2 ⎢ ( 12 )3 ⎢ ⎢ . ⎢ .. ⎢ ⎢ ( 12 )n ⎢ ⎢ .. ⎢ . ⎢ 2 M ⎢ (1 ) ⎢ ⎢(12 )M+1 ⎣ (12 )M+2

( 22 )1 ( 22 )2 ( 22 )3

··· ··· ···

( 22 )n

···

( 22 )M

.. .

( M 2 )1 ( M 2 )2 ( M 2 )3

···

.. . aM ( M 2 )n .. . ( M 2 )M

2 M+1

(2 )

(22 )M+2

.. .

0

0

λ1 λ2

μ1 μ2

λn

μn

λM

μM

.. . .. .

.. . .. .

⎤⎡ ⎤

⎡ ⎤

c1 1 ⎥⎢ c2 ⎥ ⎢0⎥ ⎥⎢ c ⎥ ⎢0⎥ ⎥⎢ 3 ⎥ ⎢ ⎥ ⎥⎢ . ⎥ ⎢ . ⎥ ⎥⎢ .. ⎥ ⎢ .. ⎥ ⎥⎢ ⎥ ⎢ ⎥ ⎥⎢ cn ⎥ ⎢0⎥ ⎥⎢ ⎥ = ⎢ ⎥. ⎥⎢ .. ⎥ ⎢ .. ⎥ ⎥⎢ . ⎥ ⎢ . ⎥ ⎥⎢ ⎥ ⎢ ⎥ ⎥⎢cM ⎥ ⎢0⎥

···

2 M+1

(M )

λM+1

⎥⎢ ⎥ ⎢ ⎥ μM+1 ⎥ ⎦⎣ b ⎦

···

(M2 )M+2

λM+2

μM+2

d

(14)

⎢ ⎥ ⎢0 ⎥ ⎣ ⎦ 0

Once coefficients cm (m = 1, 2, . . . , M) and values b, d are obtained by solving system (14), c0 is obtained using (13). The corresponding coefficients for the first four values M of the EFD, IFD and HIFD formulas are given in Table 1. Note that aside from the case M = 1, the coefficients cm and values b, d are different for the EFD, IFD and HIFD formulas.

2.3.2. Accuracy To analyze the accuracy of HIFD formula (8), we consider the absolute error obtained from Eq. (11) as follows



     M  kh  2 1  2 2 2 εM = − 2 c0 + cm cos (mkh ) + k 1 − (4b − 16d ) sin − 4d sin (kh ) . 2  h 2  m=1

(15)

This error can be calculated substituting trigonometric Taylor series into (15) and using (13) for n = {1, 2, . . . , M + 2}. Thus, the truncation error of the HIFD method is

 

   ∞ M       2 n εM =  ψn m2n+1 cm − ψn−1 −4b − 4d 22n−2 − 4 (k/2 ) h2n−2 . n=M+3  m=1

(16)

The minimum power of h in the error function is given by the term h2(M+3 )−2 ; therefore, high-order implicit finite difference formula (8) has (2M + 4 )th-order accuracy. The corresponding order of accuracy for the first four values M of all formulas are also given in Table 1.

3. Implicit finite difference schemes for the Poisson equation In this section, we show how to apply the EFD, IFD and HIFD formulas to obtain schemes with high-order numerical solutions of the Poisson equation. First, the one-dimensional method is presented to describe the characteristics of the implicit scheme and later the two-dimensional case is studied.

226

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

3.1. One-dimensional schemes Let us consider the one-dimensional Poisson equation with Dirichlet boundary conditions

d2 p ( x ) = f ( x ), x ∈ ( x0 , xL ), dx2 p( x 0 ) = g 0 , p( x L ) = g L .

(17)

where g0 and gL are fixed and known real values. The numerical domain is divided into Nx number of subintervals, h = (xL − x0 )/Nx is the mesh size, and xi = x0 + (i − 1 )h with i ∈ {1, 2, . . . , Nx + 1}. We will use the notation pi = p(xi ) and fi = f (xi ) to denote an unknown value and a forcing function value at the ith point of the grid, respectively. In accordance with formulas (8)–(9), an approximation to one-dimensional Poisson Eq. (17) at xi is given by the onedimensional (2M + 4 )th-order implicit finite difference scheme HIFD(M): M   1  cm pi+m = 1 + bh2 δx2 + dh4 δx4 fi 2 h

(18)

m=−M

where pi+m = p(xi + mh ) and the forcing function approximation is given by





1 + bh2 δx2 + dh4 δx4 fi = [d, b − 4d, 1 − 2b + 6d, b − 4d, d] · [ fi−2 , fi−1 , fi , fi+1 , fi+2 ].

(19)

When only b is taken into account in (19), Eq. (18) represents an implicit finite difference scheme IFD(M) of (2M + 2 )th-order of accuracy corresponding to formula (6). When b = d = 0, Eq. (18) represents an explicit finite difference scheme EFD(M) of (2M)th-order of accuracy corresponding to formula (3). Notice that all schemes differ only in the approximation of the forcing function. Thus, we obtain high-order schemes with the same coefficient matrix. For M = 1, the EFD(1) scheme is reduced to the classical second-order centered finite difference approximation; taking b = 1/12, d = 0, a fourth-order IFD(1) scheme is obtained; and considering b = 1/12 and d = −1/240, a sixth-order HIFD(1) scheme is obtained. The resulting matrix equation of the EFD(1) scheme is tridiagonal which is easy to solve using the Thomas algorithm. Otherwise, it is necessary to solve a more general matrix system. In reason of its efficiency and simplicity in implementation, we consider the Successive Over-Relaxation (SOR) method to obtain the solution of the system. This iterative algorithm is widely used and has shown good results for large linear system. Another important remark is that values cm , b and d do not necessarily correspond to the same number M for each point xi ; in other words we should denote HIFD(Mi ) for each xi . Notice that the HIFD(M) scheme (18)–(19) cannot be applied over grid points close to the boundary for an arbitrary number M. This is because those nodes would depend on some evaluations outside the domain. Hence, the total amount of operator points used in (18) depends on the maximum number of available nodes on the grid to apply a centered formula. In general, the HIFD(Mi ) scheme adopted in this paper follows



i−1 Mi = M Lx − i

i = 2, . . . M i = M + 1, . . . Lx − M . i = Lx − M + 1, . . . Lx − 1

(20)

For simplicity, we keep notation HIFD(M) for any point xi . Finally, the HIFD(M) scheme, applied over the whole domain, always requires right-hand side values f (x0 − h ) and f (xL + h ) at x0 and xL , respectively. If function f is not known outside the numerical domain [x0 , xL ], then these values can be obtained by extrapolation. 3.2. Two-dimensional schemes Let us consider the two-dimensional Poisson equation with Dirichlet boundary conditions

∂2 p ∂2 p (x, y ) + 2 (x, y ) = f (x, y ), (x, y ) ∈  2 ∂x ∂y p(x, y ) = g(x, y ), (x, y ) ∈ ∂ ,

(21)

where  = (x0 , xL ) × (y0 , yL ). The numerical domain is divided into Nx and Ny number of subintervals in x and y directions with x and y as mesh size, respectively. We will use the notation pi j = p(xi , y j ) and fi j = f (xi , y j ) where (xi , y j ) = (x0 + (i − 1 )x, y0 + ( j − 1 )y ) is the ijth point of the grid. Applying formula (8) to the second derivatives in x and y at Eq. (21) yields

∂y4

δ2 p δ2 p + ∂x4 2 = ∂x4 ∂y4 f + O(x2Mx +4 + y2My +4 ) 2 δx δy

(22)

where Mx and My are the fixed numbers related with the number of nodes used in the x- and y-direction respectively; and



 4 ∂2 4 ∂ ∂ = 1 + bx 2 + dx 4 ≈ 1 + bx δx2 + dx δx4 , ∂x ∂x   2 ∂ ∂4 ∂y4 = 1 + by2 2 + dy4 4 ≈ 1 + by δy2 + dy δy4 . ∂y ∂y 4 x

2

(23)

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

227

Fig. 1. EFD, IFD and HIFD stencils for the two-dimensional Poisson equation for different values of Mx = My .

The subscripts at b and d are not derivatives, they only indicate the dependance of the x and y direction respectively. Thus the two-dimensional high-order implicit finite difference scheme HIFD(Mx , My ) is given by My Mx  1 1  c D p + cyn Dx pi, j+n = Dy Dx fi j xm y i + m, j x 2 y 2 m=−Mx

(24)

n=−My

where the application of the operators are reduced to

Dx pis = [ pi−2,s , pi−1,s , pis , pi+1,s , pi+2,s ] ux T ,

(25)

Dy ps j = [ ps, j−2 , ps, j−1 , ps j , ps, j+1 , ps, j+2 ] uy T , for any index s = i + m or s = j + n with m ∈ {−My , −Mx + 1, . . . , 0, . . . , Mx − 1, Mx } 1, . . . , 0, . . . , My − 1, My }. The vectors related to the implicit coefficients are given by

ux = [dx , bx − 4dx , 1 − 2bx + 6dx , bx − 4dx , dx ],

and

with

n ∈ {−My , −My +

(26)

uy = [dy , by − 4dy , 1 − 2by + 6dy , by − 4dy , dy ]. The forcing function is approximated by

⎡f

i+2, j+2

⎢ fi+2, j+1 Dy Dx fi j = uy ⎢ fi+2, j ⎣ fi+2, j−1 fi+2, j−2

fi+1, j+2 fi+1, j+1 fi+1, j fi+1, j−1 fi+1, j−2

fi, j+2 fi, j+1 fi, j fi, j−1 fi, j−2

fi−1, j+2 fi−1, j+1 fi−1, j fi−1, j−1 fi−1, j−2

⎤

fi+2, j+2 fi+2, j+1 ⎥ fi+2, j ⎥ux T . ⎦ fi+2, j−1 fi+2, j−2

(27)

An EFD(Mx , Ny ) scheme is obtained in the absence of terms with b and d. In the case of only taking implicit term b, scheme (24)-(27) is reduced to an implicit difference scheme IFD(Mx , My ). The order of accuracy of these schemes is given by the minimum value between the order of accuracy of direction x and y. When Mx = My , scheme HIFD(Mx , My ) will be denoted as HIFD(M), and the same notation also applies for all schemes throughout the remainder of this paper. Scheme (24) is not a nine point compact scheme for Mx = My = 1 because its stencil requires unknown variables ps, j+2 , ps, j−2 , pi+2,s and pi−2,s . Thus, the discretization results in more diagonals to the final matrix system of our problem. However, IFD(Mx , My ) and HIFD(Mx , My ) are already compact for Mx ≥ 2, My ≥ 2. In other words, for each M they both have the same stencil and no additional points are required. In the following subsections we will describe the resulting coefficients for the case Mx = My = 1 and propose a compact scheme, named CHIFD(1) in order to avoid the problem previously described. The stencil of explicit and implicit schemes for Mx = My = 1 and Mx = My = 2 are illustrated in Fig. 1.

228

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

Fig. 2. Labeling of the 2D grid points.

3.2.1. Finite difference scheme for M=1 For convenience of exposition, we introduce some notations as follows

φV = φi, j+1 + φi, j−1 , φH = φi+1, j + φi−1, j , φC = φi+1, j+1 + φi+1, j−1 + φi−1, j+1 + φi−1, j−1 , φOV = φi, j+2 + φi, j−2 , φOH = φi+2, j + φi−2, j , φOCV = φi−1, j−2 + φi+1, j−2 + φi−1, j+2 + φi+1, j+2 , φOCH = φi−2, j−1 + φi−2, j+1 + φi+2, j−1 + φi+2, j+1 , where φ represents the unknowns or forcing function evaluated at different points of the stencil, as shown in Fig. 2. For Mx = My = 1, second-order explicit finite difference scheme EFD(1) is the classical five point centered finite difference approximation of the Poisson equation given by

4 pi, j + pi, j+1 + pi, j−1 + pi+1, j + pi−1, j = fi, j .

(28)

Fourth-order finite difference scheme IFD(1) is given by the following nine point compact scheme

  ω1 pi, j + ω2 pV + ω3 pH + ω4 pC = x2 ξ1 fi, j + ξ2 ( fV + fH ) + ξ3 fC ,

(29)

where

ω1 = −10(1 + γ 2 ), ω2 = 5γ 2 − 1, ω3 = 5 − γ 2 , ω4 = 12 (1 + γ 2 ), ξ1 =

25 , 6

ξ2 =

5 , 12

ξ3 =

1 , 24

and γ = x/y. In this case, the proposed implicit family recovers the same scheme as the fourth-order compact scheme obtained by several authors [4,12] using different techniques such as direct Taylor expansions [4], dual partitions [12], etc. In the case of the HIFD(1), the scheme is given by the eighteen point stencil

  1 pi, j + ω 2 pV + ω  3 pH + ω 4 pC + d pOV − 2 pOCV + γ 2 ( pOH − 2 pOCH ) ω   = x2 ξ1 fi, j + ξ2 ( fV + fH ) + ξ3 fC + ξ4 ( fOV + fOH ) + ξ5 ( fOCV + fOCH ) ,

where

1 2

1 = −2a2 (1 + γ 2 ), ω  2 = a2 − a1 γ 2 , ω 3 = −a1 + a2 γ 2 , ω  4 = a1 ( 1 + γ 2 ), ω 1 6

ξ1 = a22 , ξ2 =

1 a1 a2 , 12

a1 = 1 − 48d,

ξ3 =

1 2 a , 24 1

a2 = 5 + 36d,

1 2

ξ4 = da2 , ξ5 = da1 , d=−

1 . 240

Notice that the compact IFD(1) scheme is recovered if d = 0; otherwise scheme (30) is not a compact scheme.

(30)

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

229

3.2.2. Compact high-order implicit finite difference scheme As already indicated above, scheme (30) is not compact because it does not use the minimum nine grid points in the two-dimensional discretization formula. However, if we consider the same grid size in both x and y directions, the result is a compact scheme of the 2D Poisson equation. Thanks to this supposition, the fourth-order derivatives terms can be moved from the left-hand side to the right-hand side as described below. Let us consider Mx = My = 1, bx = by = 1/12, d = dx = dy = −1/240 and h = x = y. Then Eq. (22) yields

∂y2 (∂x4 pxx ) + ∂x2 (∂y4 pyy ) = −dh4 ( pyyyyxx + pxxxxyy ) + ∂y4 ∂x4 f + O(h6 ), where

(31)

    ∂2 ∂2 ∂x2 = 1 + bh2 2 , ∂y2 = 1 + bh2 2 , ∂x ∂y

and ∂x4 , ∂y4 are given by (23). Now using relation

−dh4 ( pyyyyxx + pxxxxyy ) = −dh4 fxxyy =

1 4 h fxxyy , 240

and

∂y4 ∂x4 f = f +

1 4 1 2 1 4 h ( fxx + fyy ) − h ( fxxxx + fyyyy ) + h fxxyy + O(h6 ), 12 240 144

Eq. (31) is simplified as



6 1 + bδy2

 δ2 p   δ2 p 1 4 1 1 4 + 6 1 + bδx2 = 6 f + h2 ( fxx + fyy ) − h ( fxxxx + fyyyy ) + h fxxyy + O(h6 ) 2 2 40 15 δx δ y2

(32)

Other authors present compact finite difference with similar (but not the same) right-hand side formulas [2,3,11,12]. However, various second- and fourth-order partial derivatives of the forcing function must be known exactly, in order to compute the solution. Instead, ∂y4 ∂x4 fi j is approximated as Dy Dx fij . Finally, the compact sixth-order implicit finite difference scheme CHIFD(1) for a 2D Poisson equation is given by













−20 pi, j + 4( pV + pH ) + pC = h2 [ ξ1 − 24d fi, j + ξ2 + 12d ( fV + fH ) + ξ3 − 6d fC +ξ4 ( fOV + fOH ) + ξ5 ( fOCV + fOCH )],

(33)

where ξ1 , ξ2 and ξ3 are defined the same as Eq. (30). Now, scheme (33) only requires the evaluation of the forcing function at several mesh points, as happens in most practical cases. 4. Numerical examples In this section, the proposed high-order finite difference methods for the Poisson equation are tested. In the first example, the accuracy of the method for a one-dimensional case is investigated; later we study the order of the algorithms using a two-dimensional rectangular domain and finally we study the convergence using the SOR iterative method. In those simulations, we test the Poisson equation with a given right-hand side forcing function and compare it with its analytic solution. In the following examples, the simulations have been performed using the SOR iterative linear solver with a relaxation value of ω = 1.3 and a tolerance value of ε = 10−12 . The order of accuracy is calculated as

log(·N1 /·N2 )

Order =

log (N2 /N1 )

,

where ·N and ·N are the L2 norms of fields with resolutions associated with grid size N1 and N2 respectively. The L2 1 2 norm of the error is calculated as



e  =

Ny Nx  

1/2

( pi j − p¯ i j ) xy 2

,

i=1 j=1

where Nx = Ny is the number of sub-intervals, pij is the result from the numerical solution, and p¯ i j is the result from the analytic solution respectively. 4.1. Example 1: one-dimensional problems Let us consider the boundary value problem for the 1D Poisson Eq. (17) where the right-hand forcing function f is assumed given. Provided that the solution p of the problem is sufficiently smooth, the approximation (18) is (2M )th-order accurate for the EFD(M), (2M + 2 )th-order accurate for the IFD(M) and (2M + 4 )th-order accurate for the HIFD(M). More precisely, to obtain the desired order of accuracy, the EFD(M), IFD(M) and HIFD(M) schemes theoretically need that the

230

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244 Table 2 L2 -norm and order of accuracy of the EFD, IFD and HIFD schemes for Example 1.1. EFD(1) Nx 16 32 64 128

EFD(2)

2

L -norm

Order −2

5.058×10 1.035×10−2 2.510×10−3 6.228×10−4

— 2.39 2.09 2.03

IFD(1) Nx 16 32 64 128

16 32 64 128

L -norm

Order −2

1.815×10 1.095×10−3 7.340×10−5 4.688×10−6

— 4.23 3.99 4.01

IFD(2)

2

L -norm

Order −2

1.179×10 4.802×10−4 2.874×10−5 1.778×10−6

— 4.83 4.15 4.06

HIFD(1) Nx

EFD(3)

2

L -norm

Order −3

7.275×10 8.398×10−5 1.330×10−6 2.087×10−8

— 6.73 6.11 6.06

L -norm −2

1.088×10 2.243×10−4 4.316×10−6 7.140×10−8

Order

L2 -norm

Order

— 5.85 5.83 5.98

9.570×10−3 6.638×10−5 3.853×10−7 1.679×10−9

— 7.49 7.60 7.93

IFD(3)

2

L -norm

Order −3

6.417×10 5.689×10−5 8.180×10−7 1.250×10−8

— 7.12 6.26 6.10

HIFD(2)

2

EFD(4)

2

IFD(4)

2

L -norm −3

4.734×10 1.216×10−5 4.402×10−8 1.684×10−10

Order

L2 -norm

Order

— 8.99 8.29 8.12

3.934×10−3 3.756×10−6 3.605×10−9 3.540×10−12

— 10.48 10.25 10.10

HIFD(3)

2

L -norm

Order −3

4.237×10 8.043×10−6 2.520×10−8 9.283×10−11

— 9.45 8.51 8.18

HIFD(4)

2

L -norm −3

3.313×10 1.830×10−6 1.269×10−9 1.133×10−12

Order

L2 -norm

Order

— 11.31 10.73 10.24

2.845×10−3 6.034×10−7 9.811×10−11 2.949×10−14

— 12.75 12.87 11.83

(2M + 2 )th, (2M + 4 )th and (2M + 6 )th derivatives of the solution, respectively, are bounded over the domain . This is because at least the first term of the truncation error must exist, as shown in Eq. (16) and the error formulas described by Liu and Sen [13]. In this case, if these schemes are stable, then they will converge with the corresponding rate. However, it is expected that the convergence rate of the proposed schemes are reduced for functions that are less smooth than those minimally required. Thus, the regularity of the solution p to the Poisson problem is immediately determined by the forcing function f. It will therefore be convenient to use different forcing functions with different degrees of regularity and to investigate directly the convergence properties of scheme (18). 4.1.1. Example 1.1 Let us first consider the following infinitely smooth and differentiable forcing function and exact solution of the onedimensional Poisson problem:

f (x ) =



2

−4π 2 (x + 1/4 )



2 2 − 8π 2 e−4π (x+1/4 ) ,

2 2 p(x ) = e−4π (x+1/4 ) ,

(34)

and the Dirichlet boundary conditions taken from the exact solution. Given the smoothness of the function p and f, we expect the Taylor series analysis of the truncation error to hold. In order to verify this claim in a numerical simulation, we lay down a computational grid on the interval −1 ≤ x ≤ 1 of constant grid spacing h = 2/Nx for Nx =16, 32, 64, 128 and M up to 4. Table 2 shows the comparison of the errors between the numerical solutions and their corresponding analytic solutions. The numerical approximation is computed using the explicit, implicit, and high-order implicit finite difference approximations. The errors are measured in the L2 -norm at all grid points. The results show a good approximation of the numerical solution, and in fact the solution of HIFD(4) with Nx = 128 already reaches values close to the machine precision. The convergence analysis indicates that the EFD(M) method has (2M)th-order accuracy, the IFD(M) method has (2M+2)th-order accuracy, and the HIFD(M) method has (2M+4)th-order accuracy, as we expected. We remark that this problem was chosen so that there is no effect of the boundary approximations in the results of the accuracy of our proposed methods. 4.2. Example 1.2 Now we will analyze the case when the regularity is formally insufficient to guarantee the corresponding order of convergence. We consider the numerical approximation to a function with finite differentiability. First, we analyze the case with a discontinuous forcing function and later we consider two examples of functions f ∈ C k [−1, 1], a continuous kth derivative on the interval −1 ≤ x ≤ 1, but not f ∈ C k+1 [−1, 1]. For simplicity, we will denote C k = C k [−1, 1]. Finally, an example with bounded kth derivative on the interval −1 ≤ x ≤ 1 is presented. We remark that the solution p always has two additional derivatives with respect to f, i.e. pk+2 ≡ f k . Note that the number k now has a different meaning from that used in Section 2.3. In this paper, we present three different cases where the function f (k+1 ) is not continuous at x = 0: the case in which one-sided limits are not equal, that in which the function does not approach a finite value, and that in which the function does not approach a particular value (oscillation). Thus, the forcing functions and exact solutions for this example

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

1

Example 1.2 (b)

Example 1.2 (c)

1

0.8

0.8

0.6

0.6

0.4

0.4

6

231

Example 1.2 (d)

4

f

2 0 -2 0.2

0.2

-4

f(x) p(x)

0 -1

-0.5

0

0.5

f(x) p(x)

f(x) p(x)

1

x

0 -1

-0.5

0

0.5

-6 -1

1

-0.5

0

0.5

1

x

x Fig. 3. Forcing functions and solutions presented in Example 1.2 and k = 0.

are selected as follows:

 Example 1.2.(a ) : f (x ) =

−1 1



x≤0 , x>0

p( x ) =

−x2 /2 x2 /2

x≤0 , x>0

f k (x ) f (x ) , k ≥ 0, p(x ) = k+2 , k ≥ 0, f k (1 ) f k (1 )  x f 0 ( x ) = |x|, f k ( x ) = fk−1 (t )dt , k ≥ 1,

Example 1.2.(b) : f (x ) =

−1

Example 1.2.(c ) : f (x ) =

f k (x ) , k ≥ 0, f k (1 )

f 0 ( x ) = |x|1/2 , f k ( x ) =

Example 1.2.(d ) : f (x ) =

⎧ ⎨0

1

p( x ) = 

x −1

fk+2 (x ) , k ≥ 0, f k (1 )

(35)

fk−1 (t )dt , k ≥ 1,

1

⎩rk (x ) sin + sk (x ) cos x x ⎧ 0 x = 0 , ⎨   p( x ) = , ⎩xk+4 sin 1 otherwise

x=0 otherwise

,

x

where rk (x ) = (k + 2 )(k + 3 )xk+2 − xk ,

sk (x ) = −2(k + 3 )xk+1 ,

and the Dirichlet boundary conditions are taken from the exact solution. Fig. 3 shows the plots of the last three examples for k = 0. For Example 1.2 (a) the forcing function is discontinuous at x = 0; thus the second derivative of p is only bounded for [−1, 1]. For Example 1.2 (b), f0 (x ) = |x| is a continuous function with discontinuous first derivative such that the right- and left-hand limits both exist, but are not equal. For Example 1.2 (c), f 0 (x ) = |x|1/2 has discontinuous first derivative for which the one-sided limits exist and both of them are −∞. Notice in examples (b) and (c) that the function f always takes values between 0 and 1. By the Fundamental Theorem of the Calculus, it is easy to prove that functions (b) and (c) satisfy f ∈ Ck and f ∈ / C k+1 . Let us now consider Example 1.2 (d). For even numbers k = 2m (m = 0, 1, 2, . . . ), it can be shown that f(m) is bounded, denoted here by f ∈ C∗m and p ∈ C ∗m+2 ; for odd numbers k = 2m + 1 (m = 0, 1, 2, . . . ), it can be shown that f(m) is continuous, f ∈ Cm and p ∈ C m+2 . In all these cases we defined p(m ) (0 ) = 0. For both odd and even values of m, the next derivative f (m+1 ) is not continuous at x = 0; even more this is an essential discontinuity, i.e. the one-sided limits does not exist (not even as ± ∞). For example, let us consider the case k = 0. The solution p(x ) = x4 sin(1/x ) (x = 0) has a bounded second derivative p (x ) = f (x ) for x ∈ [−1, 1], but this function is not continuous at x = 0. Thus we have that f ∈ C∗0 and p ∈ C∗2 . Let us now

232

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244 Table 3 Convergence rate of the proposed schemes using different test functions for Example 1.2. Example 1.2 (a)

EFD(M)

IFD(M)

HIFD(M)

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

1.03

1.03

1.03

1.03

1.03

1.03

1.03

1.03

1.03

1.03

1.03

1.03

Example 1.2 (b)

EFD(M)

k=

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

2.05 2.05 2.05 2.05 2.05 2.05 2.05 2.05 2.05 2.05 2.04

2.05 4.05 4.02 4.06 4.02 4.00 3.97 3.94 3.91 3.88 3.85

2.05 4.05 4.03 4.06 4.03 4.00 3.98 3.95 3.92 3.88 3.85

2.05 4.05 4.03 4.06 4.02 4.00 3.98 3.95 3.91 3.88 3.85

2.05 exact 4.10 4.10 4.10 4.10 4.10 4.10 4.10 4.10 4.10

2.05 4.07 4.11 6.11 6.08 6.12 6.09 6.06 6.04 6.01 5.98

2.05 4.09 4.11 6.11 6.08 6.12 6.08 6.06 6.04 6.01 5.98

2.05 4.09 4.11 6.13 6.10 6.12 6.09 6.06 6.04 6.01 5.98

2.05 4.06 4.12 6.10 6.16 6.16 6.16 6.15 6.15 6.15 6.15

2.05 4.09 4.11 5.75 5.60 5.74 6.60 7.13 7.70 7.98 8.06

2.05 4.09 4.11 6.18 6.44 6.39 7.35 7.98 8.30 8.15 8.10

2.05 4.09 4.11 5.98 6.14 5.92 6.84 7.43 8.61 8.23 8.12

0 1 2 3 4 5 6 7 8 9 10

f∈ 0

C C1 C2 C3 C4 C5 C6 C7 C8 C9 C10

IFD(M)

HIFD(M)

Example 1.2 (c)

EFD(M)

k=

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

1.57 2.04 2.05 2.05 2.05 2.05 2.05 2.05 2.05 2.05 2.04

1.54 4.15 4.08 4.06 4.02 4.00 3.97 3.94 3.91 3.88 3.85

1.54 4.01 4.09 4.06 4.02 4.00 3.97 3.94 3.91 3.88 3.85

1.54 4.01 4.09 4.06 4.02 4.00 3.97 3.94 3.91 3.88 3.85

1.54 3.44 3.70 4.10 4.10 4.10 4.10 4.10 4.10 4.10 4.10

1.54 3.56 3.59 6.16 6.13 6.12 6.08 6.06 6.03 6.00 5.97

1.54 3.57 3.59 6.07 6.13 6.12 6.08 6.06 6.03 6.00 5.97

1.54 3.57 3.59 6.01 6.15 6.13 6.09 6.06 6.03 6.00 5.97

1.54 3.56 3.60 5.09 5.95 6.16 6.15 6.15 6.15 6.15 6.15

1.54 3.57 3.59 5.61 5.56 5.38 6.59 7.21 7.76 8.01 8.06

1.54 3.57 3.59 5.65 5.67 6.06 7.40 8.21 8.29 8.15 8.10

1.54 3.57 3.59 5.66 5.71 5.58 6.85 7.57 8.59 8.21 8.12

0 1 2 3 4 5 6 7 8 9 10

f∈ 0

C C1 C2 C3 C4 C5 C6 C7 C8 C9 C10

IFD(M)

HIFD(M)

Example 1.2 (d)

EFD(M)

k=

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

0.69 1.62 2.05 2.05 2.05 2.04 2.04 2.03 2.02 2.02 2.01

0.47 1.01 3.90 4.11 3.86 3.83 3.75 3.72 3.63 3.61 3.52

0.41 1.01 3.90 4.12 3.87 3.84 3.74 3.71 3.62 3.59 3.51

0.39 1.01 3.90 4.12 3.87 3.84 3.74 3.71 3.62 3.59 3.50

0.43 1.12 2.59 3.11 4.09 4.10 4.10 4.09 4.09 4.09 4.08

0.38 1.12 2.58 2.70 4.69 5.40 5.94 5.90 5.83 5.79 5.72

0.37 1.12 2.58 2.70 4.69 5.42 5.94 5.90 5.82 5.79 5.71

0.37 1.12 2.58 2.70 4.69 5.42 5.94 5.90 5.82 5.79 5.71

0.39 1.12 2.58 2.67 3.23 4.39 6.18 6.16 6.15 6.16 6.16

0.37 1.12 2.58 2.67 3.24 2.77 6.62 7.62 8.02 7.97 7.91

0.37 1.12 2.58 2.67 3.24 2.77 6.63 7.65 8.02 7.97 7.91

0.37 1.12 2.58 2.67 3.24 2.77 6.63 7.65 8.02 7.97 7.91

0 1 2 3 4 5 6 7 8 9 10

f∈ ∗0

C C0 C∗1 C1 C∗2 C2 C∗3 C3 C∗4 C4 C∗5

IFD(M)

HIFD(M)

consider the case k = 1. The solution has a second derivative that is not only bounded, but also continuous. Thus we have f ∈ C0 and p ∈ C2 . Table 3 shows the order of accuracy computed using the EFD(M), the IFD(M), and the HIFD(M) approximations, up to M = 4, for different forcing functions f defined at (35). The order is represented by the slope of the least squares regression line of the L2 -norm of the errors using 4 different grid sizes. As expected, most schemes fail to reach the corresponding order of convergence because the solution of the Poisson problem p is not sufficiently smooth. This remark can be easily observed in Example 1.2 (a) using a discontinuous forcing function f. We notice first order convergence is obtained for the three analyzed schemes and any value M. This is a direct consequence of the poor regularity of the solution of Example 1.2 (a): the second derivative of p is not even continuous. Let us now consider the family of functions f ∈ Ck (k = 0, 1, 2, 3, . . . ) for Example 1.2 (b). Fig. 4 shows the L2 -norm error using four different grid sizes for the explicit, implicit and high-order implicit scheme for this example using M=1 and fk ∈ Ck , up to k = 5. This figure confirms that the errors of the IFD(1) and HIFD(1) schemes are smaller than for EFD(1) scheme, but the convergence rates of these schemes are no better than the second order of the EFD(1) scheme for k = 0. This is because of the discontinuity in the third derivative whereby the Taylor expansion is valid only up to the third term. However, the convergence rates of IFD(1) and HIFD(1) increase as the function gets smoother (k is bigger). Theoretically to obtain the desired order of accuracy, the EFD(1), IFD(1) and HIFD(1) schemes need that the fourth, sixth and eighth derivatives of the solution are bounded for −1 ≤ x ≤ 1, respectively. However, numerical simulations show that the desired order of accuracy of the proposed schemes are obtained even in cases where the degree of regularity of the solution p is

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

233

Fig. 4. Order of accuracy of Example 1.2 (b) using different schemes and different numbers k with f 0 (x ) = |x|.

formally insufficient to guarantee their corresponding convergence rate. Note that the EFD(1) scheme is already second-order accuracy for f ∈ C0 , the IFD(1) scheme is fourth-order accuracy for f ∈ C1 , and the HIFD(1) scheme is sixth-order accuracy for f ∈ C3 , as shown in Table 3. Even more, IFD(1) gives a machine precision error for any grid size (exact numerical solution). To demonstrate that our numerical simulations are not accidental, let us now consider Example 1.2 (c) and Example 1.2 (d). Although a second order convergence can be observed in the previous case for f ∈ C0 , the order of accuracy in example (c) is reduced to 1.5. This is directly related to the fact that the derivative of f does not exist. In general for Example 1.2 (c), all schemes need one more level of regularity to reach their corresponding convergence rate: the EFD(1) scheme is second-order for f ∈ C1 , the IFD(1) scheme is fourth-order for f ∈ C3 and the HIFD(1) scheme is sixth-order for f ∈ C4 . For Example 1.2 (d), the explicit and implicit schemes also need one more level of regularity to reach their corresponding convergence rate: the EFD(1) scheme is second-order for f ∈ C∗1 , and the IFD(1) scheme is fourth-order for f ∈ C∗3 ; however the HIFD(1) scheme is already sixth-order for f ∈ C∗3 . Notice for last example that it was only necessary that the corresponding derivatives were bounded. In much the same way one can analyze the case when M=2. For examples (b) and (c), the EFD(2) scheme has fourth-order accuracy for f ∈ C1 , the IFD(2) scheme has sixth-order accuracy for f ∈ C3 , and the HIFD(2) scheme is close to eighth-order accuracy for f ∈ C8 . We also notice that the order of accuracy of the EFD(M), IFD(M) and HIFD(M) schemes reach their maximum order in those values, even if M is increased. This behavior is directly related to the approximations at the boundaries as we discussed in Section 3.1. Thus, according to Eq. (20), the imposition of low numbers M at the boundary significantly affects the order of accuracy of the global solution. The absolute errors using the EFD(M), IFD(M) and HIFD(M) schemes for Nx = 128 and f0 (x ) = |x| is presented in Fig. 5. If the function is not sufficiently smooth, then the region where the solution has greater errors are distributed close to x = 0, where the forcing function has the discontinuity; however, once the solution is regular enough, the bigger errors are located at the boundaries for M> 1. As shown in this graph the errors for M=3 and M=4 are similar to M=2, showing the same behavior. It is important to remark that once the schemes have reached their theoretical convergence rate, they keep showing the same order of convergence notwithstanding the infinite smoothness of the right-hand side f and the solution p. In fact, no matter the regularity of the solution, the convergence rates of the implicit and high-order implicit schemes are often better or equal to the convergence rate of the explicit scheme, which means that both implicit schemes take advantage of the regularity of the solution to increase their order of accuracy. A summary of the minimum regularity required for the functions in this examples are shown in Table 4. 4.3. Example 1.3 Let us now consider an example where not only is the regularity formally insufficient to guarantee the corresponding order convergence, but also the solution of values close to the boundary has low effect in the global solution. Let us consider

234

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

EFD(1)

IFD(1)

10 -10

k=0 k=1 k=2 k=3

10 -15 -1

-0.5

0

0.5

10 -10

Error

10 -10

10 -15

10 -15

1

-1

-0.5

x

0

0.5

-1

1

10

-5

10

Error

10 -10

10 -15

10 -15

10 -15

1

10

-1

-0.5

EFD(3)

IFD(3) k=0 k=1 k=2 k=3

-5

10

-10

-1

10 -0.5

0

0.5

-1

1

-5

10

-10

10

-15

1

-1

10 -0.5

EFD(4) 10

-10

0

0.5

-1

1

k=0 k=1 k=2 k=3 k=4 k=5

10

-10

10

-1

-0.5

0

x

1

0.5

1

k=0 k=1 k=2 k=3 k=4 k=5

-10

Error

1

0.5

-5

Error 0.5

0

HIFD(4)

-5

10 -15

0

-0.5

x

10 -15

x

k=0 k=1 k=2 k=3 k=4 k=5

-15

10 -15 -0.5

1

-10

Error

10

0.5

-5

IFD(4) k=0 k=1 k=2 k=3

0

HIFD(3)

x

-5

-1

-0.5

x k=0 k=1 k=2 k=3 k=4 k=5

x

10

0.5

Error

10

-15

10

0

x

Error

10

0.5

x

Error

10

0

1

k=0 k=1 k=2 k=3 k=4 k=5

-5

Error

10 -10

-0.5

0.5

HIFD(2) k=0 k=1 k=2 k=3 k=4

Error

10 -10

0

x

IFD(2) k=0 k=1 k=2 k=3

-5

-1

-0.5

x

EFD(2) 10

k=0 k=1 k=2 k=3 k=4 k=5

10 -5

Error

10 -5

Error

10 -5

HIFD(1) k=0 k=1 k=2 k=3

-1

-0.5

0

x

Fig. 5. Absolute errors of Example 1.2 (b) using different schemes and different numbers k with f 0 (x ) = |x|.

0.5

1

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

235

Table 4 Minimum regularity of f to obtain the theoretical convergence rate for Example 1.2. The solution p always has two additional derivatives with respect to f. EFD(M)

Example 1.2 (b) Example 1.2 (c) Example 1.2 (d)

f Order f Order f Order

IFD(M)

HIFD(M)

M=1

M=2

M=1

M=2

M=1

M=2

C0 2.05 C1 2.04 C∗1 2.05

C1 4.05 C1 4.15 C∗1 3.90

C1 exact C2 3.70 C∗2 4.09

C3 6.11 C3 6.16 C∗3 5.94

C3 6.10 C4 5.95 C∗3 6.18

C8 7.70 C8 7.76 C∗4 8.02

Table 5 Convergence rate of the proposed schemes using different test functions for Example 1.3. Example 1.3 (b)

EFD(M)

IFD(M)

HIFD(M)

k=

f∈

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

0 1 2 3 4 5 6 7 8 9 10

C0 C1 C2 C3 C4 C5 C6 C7 C8 C9 C10

2.15 2.12 2.10 2.10 2.13 2.16 2.19 2.20 2.20 2.17 2.13

3.80 3.94 3.86 3.87 3.95 4.03 4.08 4.10 4.10 4.03 3.95

3.62 5.67 5.51 5.53 5.68 5.80 5.86 5.88 5.88 5.79 5.67

3.35 6.72 7.35 7.40 7.47 7.53 7.56 7.56 7.56 7.44 7.31

3.68 4.14 4.00 4.03 4.17 4.29 4.35 4.38 4.38 4.30 4.19

3.24 6.10 5.71 5.90 6.24 6.39 6.47 6.49 6.49 6.40 6.27

3.07 6.29 7.31 7.87 8.19 8.34 8.41 8.43 8.43 8.33 8.21

3.01 6.16 7.80 9.84 10.03 10.13 10.19 10.20 10.20 10.11 10.01

3.32 5.97 5.64 5.76 6.05 6.20 6.28 6.30 6.30 6.21 6.07

3.03 6.21 7.44 8.17 8.45 8.59 8.66 8.67 8.67 8.59 8.47

2.94 5.99 7.65 10.17 10.53 10.64 10.69 10.71 10.71 10.63 10.54

2.88 5.93 7.56 10.49 12.09 12.47 12.48 12.50 12.50 12.36 12.41

Example 1.3 (c)

EFD(M)

IFD(M)

HIFD(M)

k=

f∈

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

0 1 2 3 4 5 6 7 8 9 10

C0 C1 C2 C3 C4 C5 C6 C7 C8 C9 C1 0

2.14 2.12 2.10 2.10 2.13 2.16 2.19 2.20 2.19 2.17 2.13

2.46 3.93 3.85 3.87 3.95 4.03 4.08 4.10 4.08 4.03 3.95

2.14 5.46 5.50 5.54 5.69 5.80 5.86 5.88 5.85 5.79 5.67

1.98 5.54 6.94 7.41 7.48 7.54 7.56 7.56 7.52 7.44 7.30

2.20 4.13 4.00 4.04 4.18 4.29 4.36 4.38 4.36 4.30 4.19

1.90 5.24 5.63 5.94 6.25 6.40 6.47 6.49 6.47 6.40 6.26

1.76 4.95 6.15 7.90 8.20 8.34 8.41 8.43 8.40 8.33 8.20

1.67 4.91 6.40 9.22 10.02 10.14 10.19 10.20 10.18 10.11 10.01

1.95 5.32 5.60 5.79 6.07 6.21 6.28 6.30 6.27 6.20 6.07

1.74 4.89 6.22 8.19 8.46 8.59 8.66 8.67 8.65 8.58 8.47

1.67 4.82 6.21 9.21 10.46 10.64 10.70 10.71 10.68 10.62 10.54

1.61 4.81 6.13 9.20 10.92 12.50 12.50 12.52 12.47 12.42 12.25

Example 1.3 (d)

EFD(M)

IFD(M)

HIFD(M)

k=

f∈

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

0 1 2 3 4 5 6 7 8 9 10

C∗0 C0 C∗1 C1 C∗2 C2 C∗3 C3 C∗4 C4 C∗5

1.60 2.63 2.65 2.18 2.25 2.34 2.36 2.32 2.22 2.08 1.90

1.60 2.65 2.80 2.85 3.82 4.15 4.18 4.16 4.16 4.14 4.03

1.59 2.64 2.76 2.71 4.03 4.95 5.83 5.84 5.69 5.38 5.38

1.59 2.65 2.82 3.00 3.86 5.09 6.43 6.43 5.98 5.66 5.49

1.60 2.65 2.81 2.87 3.78 4.42 4.52 4.46 4.40 4.34 4.15

1.60 2.65 2.81 2.87 3.70 4.73 5.88 6.44 6.56 6.61 6.52

1.60 2.65 2.80 2.81 3.71 4.66 5.82 6.92 7.87 7.80 7.48

1.59 2.65 2.81 2.86 3.63 4.65 5.84 6.92 7.91 7.83 7.53

1.60 2.65 2.81 2.86 3.72 4.76 5.88 6.26 6.27 6.27 6.14

1.60 2.65 2.81 2.88 3.64 4.66 5.77 6.79 8.31 8.74 8.67

1.60 2.65 2.81 2.86 3.62 4.60 5.63 6.64 8.32 9.82 9.76

1.60 2.65 2.81 2.87 3.60 4.60 5.62 6.52 8.23 9.81 9.81

the same functions f as in Example 1.2 but multiplied by the Gaussian function presented in the Example 1.1. In other words, the forcing function and exact solution are selected as follows:

f (x ) = g (x )h(x ) + 2g (x )h (x ) + g(x )h (x ),

(36)

p( x ) = g ( x ) h ( x ) where g(x ) = e−4π (x+1/4 ) and functions h are given by Eq. (35) for the cases (b), (c) and (d). Notice that the regularity of f is the same as the previous examples. Table 5 shows the order of accuracy computed using the EFD(M), IFD(M), and HIFD(M) approximations for different forcing functions f, up to M = 4. As in the previous example, most schemes fail to reach the corresponding convergence rate 2

2

236

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

Table 6 Minimum regularity of f to obtain the theoretical convergence rate for Example 1.3. The solution p always has two additional derivatives with respect to f EFD(M)

Ex. 1.3 (b) Ex. 1.3 (c) Ex. 1.3 (d)

f Order f Order f Order

IFD(M)

HIFD(M)

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

M=1

M=2

M=3

M=4

C0 2.15 C0 2.14 C0 2.63

C0 3.80 C1 3.93 C∗2 3.82

C1 5.67 C2 5.50 C∗3 5.83

C2 7.35 C4 7.48 — —

C0 3.68 C1 4.13 C∗2 3.78

C1 6.10 C2 5.63 C∗3 5.88

C3 7.87 C3 7.90 C∗4 7.87

C3 9.84 C4 10.02 — —

C1 5.97 C2 5.79 C∗3 5.88

C2 7.44 C2 8.19 C∗4 8.31

C3 10.17 C4 10.46 C4 9.82

C4 12.09 C5 12.50 — —

because the solution of the Poisson problem p is not sufficiently smooth. However, in these examples, the convergence rate can be retained for bigger numbers M. For Example 1.3 (b), the second-, fourth- and sixth-order accuracy of the EFD(1), IFD(1) and HIFD(1) schemes are obtained when f is C0 , C1 , and C1 , respectively. In all cases, the degree of regularity of the solution p is formally insufficient to guarantee the corresponding order of accuracy. The effects of the discontinuity are more clearly seen in the absolute error plot. The error distribution for Nx = 128 is presented in Fig. 6. An examination of the graph of the IFD(1) scheme shows that there is an influence of the discontinuity close to zero for f ∈ C0 ; however the errors are not completely dominated by then. The resulting convergence rate is 3.68, as shown in Table 5. If we now consider a forcing function f ∈ C1 , then there is no more influence of the discontinuity in the solution and the resulting convergence rate is 4.14. The same analysis can be done for the HIFD(1) scheme, resulting in a fourth-order accuracy for a function f ∈ C1 . Now, if we consider M> 1, then more regularity of the forcing function f and solution p is generally needed. For example, the EFD(2) scheme is close to fourth-order accuracy for f ∈ C0 , the IFD(2) scheme is close to sixth-order accuracy for f ∈ C1 , and the HIFD(2) scheme is close to eighth-order accuracy for f ∈ C3 , as shown in Table 5. In much the same way one can analyze the other two examples. In general, we note that the IFD(M) scheme and the HIFD(M) scheme require more regularity of the forcing function f than the regularity of the EFD(M) scheme to obtain their corresponding convergence rate. Table 6 summarized the minimum order required for the solution of the 1D Poisson equation to reach the convergence rate. These values were chosen by examination of Table 5. It is shown that the implicit scheme needs f to have one more degree of regularity compared with the explicit scheme; meanwhile the high-order scheme needs two more degrees. We can notice that Example 1.3 (d) does not obtain twelfth-order accuracy for M = 4 because of the effect of the solution close to the boundary. 4.4. Example 2: two-dimensional problems This example examines the implicit method’s capacity for simulating 2D problems with different right-hand side forcing functions and different boundary conditions using several mesh resolutions. In all cases, the same order of the method, M = Mx = My , and the same number of sub-intervals, Nx = Ny , in the x- and y-direction are considered. Let  = (0, 1 ) × (0, 1 ) and consider the Poisson equation with forcing function f and the Dirichlet boundary conditions prescribed to satisfy the following exact solutions:

Example 2.1 :

p(x, y ) = exp(−0.5(4π )2 ((x − 1/2 )2 + (y − 1/2 )2 )),

Example 2.2 :

p(x, y ) = sin(2π x ) sin(2π y )(x3 − y4 + x2 y3 ),

Example 2.3 :

p(x, y ) =

 √ √ sin(π x )  2 sinh(π 2y ) + sinh(π 2(1 − x )) . √ sinh(π 2 )

For the first example, the 2D Poisson equation has boundary conditions which do not significantly affect the final solution of the problem. The second test problem has zero boundary values on the entire ∂  with a non-zero forcing function. For the last 2D numerical test, the Poisson equation has non-homogeneous boundary conditions and non-zero forcing function. Notice that the examples have an infinitely smooth and differentiable forcing function and exact solution of the two-dimensional Poisson problem. The exact solution of the three examples are shown at Fig. 7. The effectiveness of the proposed scheme can be demonstrated in the numerical analysis of Example 2.1. Fig. 8 and Table 7 show a comparison of the absolute errors using Example 2.1. The numerical solution has been found using three different M values. In all cases, the results show a good approximation of the exact solution. The solution errors from the explicit scheme are much higher than those from the proposed implicit schemes. The results from the HIFD(1) scheme with Nx = 16 mesh size are even more accurate than those from the EFD(1) scheme with Nx = 64 mesh size and almost the same accurate as the EFD(1) scheme with Nx = 128 mesh size. This indicates the effectiveness of the proposed scheme. The convergence analysis indicates that the EFD method is of (2M)th-order, IFD method is of (2M + 2)th-order and HIFD method is of (2M + 4)th-order of accuracy.

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

10

EFD(1)

0

10

k=0 k=1 k=2

10 -5

-10

10

-15

-1

-0.5

0

0.5

1

10

-10

10

-15

-1

-0.5

10

-10

1

10 -0.5

0

-15

-1

IFD(2)

-5

10

k=0 k=1 k=2 k=3

-10

10

0.5

-15

1

-1

10 -0.5

x

0

0.5

-10

-15

-1

k=0 k=1 k=2 k=3 k=4

-10

0

0.5

1

10

-20

-1

-0.5

x 10

k=0 k=1 k=2 k=3 k=4 k=5

10

10

-15

10 -20 -1

1

10

-0.5

0

x

10

-20

-1

-0.5

IFD(4)

-5

0.5

1

0.5

1

HIFD(4)

10 -5

k=0 k=1 k=2 k=3 k=4 k=5

-10

10

-15

10 -20 -1

0

x

Error

-10

Error

10

0.5

k=0 k=1 k=2 k=3 k=4

x

EFD(4)

-5

1

k=0 k=1 k=2 k=3 k=4 k=5

-10

Error

10

0

0.5

-10

Error

10

Error -0.5

0

HIFD(3)

10 -5

10 -15

-1

-0.5

x

10 -15

-20

1

k=0 k=1 k=2 k=3

Error

10

0.5

-10

1

IFD(3)

10 -5

k=0 k=1 k=2 k=3 k=4

0

HIFD(2)

-5

10 -15

10

-0.5

x

EFD(3)

10 -5

10

x

Error

10

-15

-1

10

0.5

-10

Error

k=0 k=1 k=2 k=3

Error 10

0

10

x

EFD(2)

-5

k=0 k=1 k=2

10 -5

x 10

HIFD(1)

0

Error

Error 10

10

10

k=0 k=1 k=2

Error

10 -5

IFD(1)

0

237

10

-0.5

0

x

0.5

1

-15

10 -20 -1

-0.5

0

x

Fig. 6. Absolute errors of Example 1.3 (b) using different schemes and different numbers k with f 0 (x ) = |x|.

0.5

1

238

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

Fig. 7. Exact solution of the two-dimensional examples.

Fig. 8. Absolute errors of Example 2.1 using different numbers M and Nx = 32.

Comparisons of the absolute errors of Example 2.2 are plotted in Fig. 9. The implicit schemes correctly simulate the Poisson equation, and are capable of resolving the homogenous boundary conditions and non-zero forcing function. Numerical solutions of the HIFD scheme is significantly higher than that of the IFD scheme under the same grid setting. It is also indicated that, even under large grid spacings, the absolute errors of the implicit schemes are still very low (6.88×10−5 for HIFD(1) and Nx = 8). Fig. 9 also shows that the maximum errors for Example 2.2 are concentrated close to the boundary resulting in similar order of accuracy for values M ≥ 2, as indicated in Table 8. As discussed in Section 3, the numerical accuracy of the EFD(1), IFD(1) and HIFD(1) scheme for the 2D Poisson equation is second-, fourth-, and sixth-order of accuracy, and the accuracy of the EFD(2), IFD(2) and HIFD(2) scheme is fourth-, and sixth- and eighth-order of accuracy. However

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

Table 7 L2 -norm and order of accuracy of the EFD, IFD and HIFD schemes for Example 2.1. EFD(1) Nx

EFD(2)

L2 -norm −2

8 16 32 64 128

7.35×10 8.54×10−3 2.01×10−3 4.95×10−4 1.23×10−4

Nx

L2 -norm

Order

4.54×10 1.67×10−3 1.17×10−4 7.62×10−6 4.82×10−7

Order

L2 -norm

8 16 32 64 128

Order

−2

3.67×10 8.40×10−4 4.78×10−5 2.92×10−6 1.81×10−7

— 5.94 4.32 4.12 4.05

L -norm −2

3.19×10 2.92×10−4 4.65×10−6 7.33×10−8 1.14×10−9

Order — 7.38 6.24 6.12 6.06

−2

4.90×10 6.11×10−4 1.33×10−5 2.34×10−7 3.76×10−9

Order

L2 -norm

Order — 6.89 5.76 5.97 6.02

IFD(3)

−2

3.01×10 2.12×10−4 2.81×10−6 4.21×10−8 6.50×10−10

— 7.79 6.52 6.20 6.08

HIFD(2)

2

L2 -norm

— 5.19 4.01 4.03 4.03

IFD(2)

HIFD(1) Nx

−2

— 3.38 2.18 2.07 2.03

IFD(1)

8 16 32 64 128

EFD(3)

L2 -norm

−2

3.12×10 9.18×10−5 3.06×10−7 1.14×10−9 4.53×10−12

Order — 9.16 8.60 8.24 8.07

HIFD(3)

2

L -norm −2

2.72×10 7.25×10−5 1.84×10−7 6.41×10−10 2.43×10−12

Order

L2 -norm

Order

— 9.32 9.01 8.35 8.13

2.73×10−2 3.90×10−5 1.95×10−8 1.60×10−11 1.75×10−14

— 10.30 11.45 10.48 9.94

Fig. 9. Absolute errors of Example 2.2 using different numbers M and Nx = 32.

239

240

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244 Table 8 L2 -norm and order of accuracy of the EFD, IFD and HIFD schemes for Example 2.2. EFD(1) Nx

EFD(2)

2

L -norm −3

8 16 32 64 128

8.03×10 2.01×10−3 5.01×10−4 1.25×10−4 3.13×10−5

Nx 8 16 32 64 128

Order

EFD(3)

2

L -norm −3

Order

L2 -norm

Order

— 3.43 3.77 3.90 3.96

6.52×10−3 7.11×10−4 5.72×10−5 4.02×10−6 2.66×10−7

— 3.48 3.80 3.92 3.96

— 2.18 2.09 2.05 2.02

6.17×10 6.94×10−4 5.71×10−5 4.06×10−6 2.70×10−7

L2 -norm

Order

L2 -norm

Order

L2 -norm

Order

3.62×10−4 2.21×10−5 1.37×10−6 8.58×10−8 5.36×10−9

— 4.39 4.19 4.09 4.05

1.91×10−4 4.80×10−6 9.92×10−8 1.78×10−9 2.99×10−11

— 5.79 5.85 5.93 5.96

2.04×10−4 5.18×10−6 1.03×10−7 1.85×10−9 3.10×10−11

— 5.80 5.87 5.93 5.97

IFD(1)

IFD(2)

HIFD(1)

IFD(3)

HIFD(2)

HIFD(3)

Nx

L2 -norm

Order

L2 -norm

Order

L2 -norm

Order

8 16 32 64 128

6.88×10−5 1.09×10−6 1.72×10−8 2.69×10−10 4.20×10−12

— 6.51 6.26 6.13 6.07

3.38×10−5 1.78×10−7 8.50×10−10 3.65×10−12 2.19×10−14

— 8.25 8.06 8.04 7.46

3.39×10−5 1.81×10−7 8.75×10−10 3.78×10−12 1.35×10−14

— 8.22 8.04 8.03 8.22

Table 9 L2 -norm and order of accuracy of the EFD, IFD and HIFD schemes for Example 2.3. EFD(1) Nx 8 16 32 64 128

EFD(2)

L2 -norm −3

4.72×10 1.20×10−3 3.02×10−4 7.57×10−5 1.89×10−5

Order — 2.15 2.08 2.04 2.02

IFD(1) Nx 8 16 32 64 128

EFD(3)

L2 -norm −3

2.35×10 2.11×10−4 1.55×10−5 1.04×10−6 6.75×10−8

Order — 3.79 3.94 3.98 3.99

IFD(2)

2

L -norm −5

3.15×10 1.90×10−6 1.18×10−7 7.37×10−9 4.60×10−10

Order — 4.41 4.19 4.09 4.05

HIFD(1) 2

L2 -norm −3

2.37×10 2.11×10−4 1.53×10−5 1.02×10−6 6.62×10−8

Order — 3.81 3.96 3.99 4.00

IFD(3)

2

L -norm −5

1.90×10 6.41×10−7 1.41×10−8 2.61×10−10 4.41×10−12

Order

L2 -norm

Order

— 5.33 5.75 5.89 5.96

2.03×10−5 6.97×10−7 1.54×10−8 2.83×10−10 4.82×10−12

— 5.31 5.75 5.90 5.94

HIFD(2) 2

HIFD(3)

Nx

L -norm

Order

L -norm

Order

L2 -norm

Order

8 16 32 64 128

8.83×10−6 1.38×10−7 2.15×10−9 3.36×10−11 5.37×10−13

— 6.54 6.27 6.14 6.04

3.90×10−6 2.14×10−8 9.78×10−11 4.17×10−13 4.50×10−14

— 8.18 8.13 8.05 3.25

3.95×10−6 2.21×10−8 1.01×10−10 4.20×10−13 3.11×10−14

— 8.16 8.12 8.09 3.80

for M > 2, the order of accuracy already reaches a maximum value (depending on the method employed) lower than the theoretical value because of the boundary approximation using Eq. (20). Fig. 10 shows the absolute errors from Example 2.3 with different M values and Nx = 32. The implicit schemes correctly simulate the general problem. However, the same problem at the boundaries can be clearly observed. Table 9 presents the convergence analysis obtained at various mesh sizes for the 2D Poisson equation using different values M. A similar analysis can be done as the previous two examples, which further enhances our conclusions. Notice that the order of accuracy for the case Nx = 128 is lost because the absolute errors should be smaller than the machine precision.

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

241

Fig. 10. Absolute errors of Example 2.3 using different numbers M and Nx = 32.

(a)

Tolerance

10

EFD(1) EFD(2) EFD(3) IFD(1) IFD(2) IFD(3) HIFD(1) HIFD(2) HIFD(3)

-2

10 -4

10

-6

10

-8

10 -10

0

500

1000

Iterations

1500

(b)

10 0

10

Tolerance

10 0

2000

EFD(1) EFD(2) EFD(3) IFD(1) IFD(2) IFD(3) HIFD(1) HIFD(2) HIFD(3)

-2

10 -4

10

-6

10

-8

10 -10

0

0.5

1

Time (sec)

1.5

Fig. 11. Efficiency analysis of Example 2.3 using Nx = 32. (a) Number of iterations and (b) CPU time (in seconds) in function of the tolerance values of the SOR method (ω = 1.3).

4.5. Efficiency of the numerical solution One of the important advantages of the proposed scheme is that in comparison with the explicit scheme implicit methods converge faster as the M value is increased. A quantitative comparison in terms of iterations and the execution time of the EFD(M), IFD(M) and HIFD(M) scheme using the SOR method is presented in Fig. 11. In this comparison, Example 2.3 was used with Nx = Ny = 32 grid size and relaxation value ω = 1.3.

242

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

EFD(1)

imaginary part

1

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

-0.5

1

imaginary part

0

0.5

1

real part IFD(1)

-1 -1

-0.5

0

0.5

1

real part IFD(2)

1

-1 -1

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

-0.5

0

0.5

1

real part HIFD(1)

1

-1 -1

-0.5

0

0.5

1

real part HIFD(2)

1

-1 -1

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

-0.5

0

real part

0.5

1

-1 -1

-0.5

0

real part

0.5

1

-1 -1

0

0.5

1

0.5

1

0.5

1

real part IFD(3)

-0.5

0

real part HIFD(3)

1

0.5

-1 -1

-0.5

1

0.5

-1 -1

EFD(3)

1

0.5

-1 -1

imaginary part

EFD(2)

1

-0.5

0

real part

Fig. 12. Eigenvalue distribution and maximum eigenvalue of the SOR (ω = 1.3) iterative matrix using Nx = 32. (For interpretation of the references to colour in the text, the reader is referred to the web version of this article.)

The results show that there is a considerable difference between the performance of the finite difference methods. The IFD(M) method produces a smaller number of iterations and consequently less CPU time, as M is increased. However the EFD(M) method has an opposite behavior, more iterations are obtained as M is increased. Thus implicit methods have not only greater precision, but also more efficient in terms of the number of iterations with the increase of M. The CPU time reaches its minimum value around M = 5. It is also shown that HIFD(2) scheme takes half of the iterations of the compact nine-point IFD(2) and one third of the explicit five-point EFD(1) scheme. The same behavior was observed for Example 2.1 and Example 2.2 (not all shown here). Numerical simulations with different grid size also give similar results. We remark that the resulting iterative matrix is always considered as a sparse matrix in the input solver. The impact of the SOR solver on the efficiency of the numerical solution of the Poisson equation is better explained by examining the spectral radius of the iterative matrix. It is known that the eigenvalue of iterative matrix with the largest amplitude, named r here, reflects the convergence rate of the method. The fastest convergence rate is achieved when the smallest value of r is achieved. Let us consider the classical stationary iterative methods: Jacobi, Gauss–Seidel (GS), and SOR

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

243

Table 10 Amplitude of the largest eigenvalue of the iterative matrix of the classical stationary iterative methods. EFD

IFD

M

Jacobi

G–S

SOR

Jacobi

G–S

SOR

Jacobi

G–S

SOR

1 2 3 4 5

0.9511 1.0638 1.1244 1.1532 1.1607

0.9045 0.9219 0.9271 0.9289 0.9294

0.8187 0.8518 0.8618 0.8653 0.8661

0.9417 0.9237 0.9099 0.9025 0.9005

0.8873 0.8542 0.8295 0.8164 0.8128

0.7873 0.7251 0.6779 0.6526 0.6458

0.9417 0.8870 0.8631 0.8565 0.8559

0.8873 0.7893 0.7488 0.7380 0.7371

0.7873 0.5996 0.5196 0.4986 0.4970

(a)

10 0

HIFD(3) 10

-2

ω ω ω ω ω

(b)

0.9 = 1.0 = 1.1 = 1.2 = 1.3 = 1.4

ω ω ω ω ω

HIFD

0.85 0.8 0.75

10 -4

= 1.0 = 1.1 = 1.2 = 1.3 = 1.4

0.7

r

Tolerance

HIFD

10

0.65

-6

0.6 10

0.55

-8

0.5 10

-10

0

200

400

Iterations

600

800

0.45

1

2

3

4

5

M

Fig. 13. Convergence analysis of Example 2.3 using Nx = 32. (a) Number of iterations of the HIFD(3) method and (b) the largest amplitude eigenvalue using the SOR method with different relaxation values ω.

iteration which are based on splittings of iterative matrix A. They lead to the following iteration matrices [15]

AJ = D−1 (L + U ), AGS = −(D + L )−1U, ASOR = −(D + ωL )−1 ((1 − ω )D − ωU ), where D is the diagonal of A and L and U are the (strict) lower and upper triangular parts. Note GS iteration overwrites the approximate solution with the new value as soon as it is computed, and the SOR iteration modifies Gauss–Seidel by adding a relaxation parameter ω to construct an iteration. It is expected that the convergence of the iterative solvers determinate by AJ , AGS and ASOR will be affected by M. This expectation is confirmed in Table 10 showing the amplitude of the largest eigenvalue of the iterative matrix of the classical stationary iterative methods. The SOR method was applied with a relaxation parameter of ω = 1.3. It indicates that the iteration for all methods are convergent (r < 1), and the convergence rate r is better for both implicit schemes IFD and HIFD, and worse for the EFD scheme as the number M is increased. We also notice that the SOR method improves considerably the convergence of the GS method. Fig. 12 shows the distribution of the eigenvalues for the EFD, IFD and HIFD scheme for the SOR iterative solver using M = 1, 2, 3, relaxation parameter ω = 1.3, and Nx = 32. We can notice that the eigenvalues of the implicit schemes are more concentrated around the origin, including r (the eigenvalue with the largest amplitude is in red). The impact of the relaxation parameter ω on the efficiency of the numerical solution of the Poisson equation is shown in Fig. 13 by examining the number of iterations and amplitude of the largest eigenvalue. In this comparison, the HIFD(M) were used. Numerical results show that the same behavior is obtained for the explicit and implicit schemes as the number M is increased using different ω values. 5. Conclusions The present paper outlines a general family of implicit schemes with any order of accuracy for the 2D Poisson equation. The newly developed schemes are accurate, efficient and have been validated by several test cases. We also proposed an alternative method to obtain compact finite difference methods of fourth- and sixth-order of accuracy. Numerical results show that eighth-order implicit scheme with SOR method converges faster than the sixth-order compact scheme, reducing its computational time considerably. We conclude that even when more points are used, the proposed implicit methods achieve not only greater precision, but also faster convergence of the solutions. Other iterative solvers, such as the Generalized Minimal Residual method and Conjugate Gradient method, were tested for efficiency. Although a similar behavior of

244

M. Uh Zapata, R. Itzá Balam / Applied Mathematics and Computation 309 (2017) 222–244

the explicit and the implicit schemes were observed, more analysis should be done and it will be reported in future work. This method can also be extended to solve the 3D case, the more general nonlinear case and time-dependent case. The research on these aspects will be also reported in our future work. Acknowledgment This work was partially supported by a Mexican Council of Science and Technology CONACYT project (CB No. 256252). References [1] H.J. Diersch, C. Fletcher, Computational techniques for fluid dynamics, in: Vol. I: Fundamental and General Techniques. Vol. II: Specific Techniques for Different Flow Categories, Springer Verlag, Berlin, 1988. [2] W.F. Spotz, High-order compact scheme for computational mechanics, in: Presented to the Faculty of the Graduate school of the university of Texas at Austin in Partial Fulfillment of the requirements for the degree of Doctor of philosophy, 1995. [3] Z.F. Tian, An high-order compact finite-difference schemes for the second dimensional Poisson equation, J. Northwest Univ. 26 (2) (1996) 109114. [4] J. Zhang., Multigrid method and fourth-order compact scheme for 2D Poisson equation with unequal mesh-size discretization, J. Comput. Phys. 179 (1) (2002) 170–179. [5] R.K. Mohanty, S. Singh, A new fourth order discretization for singularly perturbed two dimensional non-linear elliptic boundary value problems, Appl. Math. Comput. 175 (2006) 14001414. [6] M. Li, B. Fornberg, T. Tang, A compact fourth order finite difference scheme for the steady incompressible Navier–Stokes equations, Int. J. Numer. Methods Fluids 20 (1995) 11371151. [7] M. Li, T. Tang, A compact fourth-order finite difference scheme for unsteady viscous incompressible flows, J. Sci. Comput. 16 (2001) 2946. [8] P.C. Chu, C. Fan, A three-point six-order nonuniform combined compact difference scheme, J. Comput. Phys. 148 (1999) 663674. [9] H. Sun, J. Zhang, A high order finite difference discretization strategy based on extrapolation for convection diffusion equations, Numer. Methods Partial Differential Equat. 20 (1) (2004) 1832. [10] M. Nabavi, M.H.K. Siddiqui, J. Dargahi, A new 9-point sixth-order accurate compact finite difference method for the Helmholtz equation, J. Sound Vib. 307 (2007) 972982. [11] Y. Wang, J. Zhang, Sixth order compact scheme combined with multigrid method and extrapolation technique for 2D Poisson equation, J. Comput. Phys. 228 (2009) 137146. [12] S. Zhai, X. Feng, Y. He, A new method to deduce high-order compact difference schemes for two-dimensional Poisson equation., Appl. Math. Comput. 230 (2014) 9–26. [13] Y. Liu, M.K. Sen, A practical implicit finite-difference method: examples from seismic modeling, J. Geophys. Eng. 6 (3) (2009) 231. [14] J.F. Claerbout, The craft of wave-field extrapolation, in: Imaging the Earth’s Interior, Blackwell Scientific Publications, Oxford, 1985, pp. 260–265. [15] C.T. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM, 1995.