Implicit locally one-dimensional methods for two-dimensional diffusion with a non-local boundary condition

Mathematics and Computers in Simulation 49 (1999) 331±349

Implicit locally one-dimensional methods for two-dimensional diffusion with a non-local boundary condition Mehdi Dehghan* Department of Applied Mathematics, Faculty of Mathematics and Computer Science, Amirkabir University of Technology, Hafez Avenue No. 424, Tehran 15914, Iran Received 1 April 1998; accepted 1 April 1999

Abstract Two new second-order finite difference techniques based upon the classical 3-point backward time centered space (BTCS) method and the Crank±Nicolson scheme, and also a fourth-order finite difference scheme based on Crandall's method for onedimensional diffusion, are used to solve the two-dimensional time dependent diffusion equation with non-local boundary conditions. In these cases locally one-dimensional (LOD) techniques are used to extend the one-dimensional techniques to solve the two-dimensional problem. The stability properties and truncation error of these methods are discussed and the results of a numerical experiment for these three methods are presented. Error estimates are also tabulated. The results of numerical testing shows that these schemes uses less central processor (CPU) time than the fully implicit schemes. # 1999 IMACS/ Elsevier Science B.V. All rights reserved. Keywords: Finite differences schemes; Heat equation; LOD techniques; Non-local boundary value problems; Numerical integration techniques; Implicit methods; Partial differential equations; Stability CPU time

1. Introduction The problem to which the three finite difference methods are applied is the two-dimensional timedependent diffusion equation: @u @ 2 u @ 2 u ˆ ‡ @t @x2 @y2

(1)

with initial condition given by: u…x; y; 0† ˆ f …x; y†;

0  x; y  1;

ÐÐÐÐ * Corresponding author. 0378-4754/99/$20.00 # 1999 IMACS/Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 4 7 5 4 ( 9 9 ) 0 0 0 5 6 - 7

(2)

332

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

and boundary conditions: u…0; y; t† ˆ g0 …y; t†;

0  t  T; 0  y  1

(3)

u…1; y; t† ˆ g1 …y; t†;

0  t  T; 0  y  1

(4)

u…x; 1; t† ˆ h1 …x; t†;

0  t  T; 0  x  1

(5)

u…x; 0; t† ˆ h0 …x†…t†;

0  t  T; 0  x  1

(6)

with the non-local boundary condition: Z 1 Z d…x† u…x; y; t† dx dy ˆ m…t†; 0  x; y  1 0

(7)

0

where f, g0, g1, h0, h1, d and m are known functions, while the functions u and  are unknown. This kind of problem arises in many important applications in heat transfer, control theory, thermoelasticity and medical science [1±7]. The schemes on which the solution of (1) is based are described in Section 2. The method of incorporating (7) with  unknown is described in Section 3. The results produced by the three numerical methods in solving the two-dimensional diffusion equation for a test problem are described in Section 4. A discussion on their computational efficiency is given in Section 5. Section 6 summarizes the findings of this article.

2. Finite difference methods The numerical methods developed in this article have three common elements. Firstly, the LOD time splitting method is used. Secondly, the methods under consideration (the unconditionally stable secondorder 3-point BTCS scheme), the Crank±Nicolson method or the unconditionally stable fourth-order Crandall's method are used to solve the resulting one-dimensional diffusion equations. Thirdly, a numerical integration procedure is used to approximate the non-local boundary condition (7). Commencing with the initial condition (2) for each nˆ0, 1, 2, . . ., N, the process of stepping from time tn to tn‡1 is carried out in two stages. In the first stage, in advancing from tn ˆ nk to the time tn‡1=2 ˆ …tn ‡ k=2†, the partial differential equation: 1 @u @2u ˆ x 2 2 @t @x

(8)

is solved numerically at the spatial points (xi, yj), i ˆ 1; 2; . . . ; M ÿ 1, for each j ˆ 0; 1; . . . ; M. Commencing with previously computed values uni; j , i; j ˆ 1; 2; . . . ; M ÿ 1, and boundary values: un0; j ˆ g0 …yj ; tn †;

j ˆ 0; 1; 2; . . . ; M

(9)

unM; j ˆ g1 …yj ; tn †;

j ˆ 0; 1; 2; . . . ; M

(10)

n‡1=2

results in the set of approximate values ui; j intermediate time tn‡1=2 .

, i ˆ 1; 2; . . . ; M ÿ 1, j ˆ 0; 1; . . . ; M, being found at the

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

333

Then in advancing from the time tn‡1=2 to tn‡1 ˆ …tn ‡ k† the equation: 1 @u @2u ˆ y 2 2 @t @y

(11) n‡1=2

is solved numerically at spatial points (xi, yj), commencing with initial values ui; j , n‡1=2 n‡1=2 and ui; M , i ˆ 1; 2; . . . ; M ÿ 1. Note that i; j ˆ 1; 2; . . . ; M ÿ 1 and using as boundary values ui;0 the boundary conditions (3)±(6) are not used at the intermediate time tn‡1=2 . This is because in the time interval tn to tn‡1=2 , the process of diffusion in the x-direction has been applied with a diffusion coefficient which is twice that in the original Eq. (1) as can be seen by rearranging (8) in the form: @u @2u ˆ 2x 2 : @t @x

(12)

This gives results not compatible with the boundary conditions (3) and (4) being applied at the intermediate time; the given boundary values g0, g1, h0, and h1 at the intermediate time in (1) are based on x and y diffusion coefficients in (1) being applied over the half-time step [9]. Note that the values of n‡1=2 ui; j , i; j ˆ 1; 2; . . . ; M, are not approximate solutions to the original problem. 2.1. The classical 3-point BTCS method The use of the 3-point BTCS method [8], in the x and y sweeps of the LOD method for solving the two one-dimensional problems (8) and (11), is described in the following. In the y-sweep of the LOD procedure, we have for each i ˆ 1; . . . ; M ÿ 1: n‡1=2

n‡1=2

ÿsui; jÿ1 ‡ …1 ‡ 2s†ui; j

n‡1=2

ÿ sui; j‡1 ˆ uni; j

(13) n‡1=2

refers to values of ui; j which is applied with j ˆ 2; 3; . . . ; M ÿ 2, where s ˆ k=h2 . The notation ui; j computed at the intermediate stage, that is, at the end of time …tn ‡ k=2†. This procedure is unconditionally von Neumann stable and solvable [8]. To avoid any calculation using the intermediate values the following fourth-order implicit formulae [17] for i ˆ 1; . . . ; M ÿ 1 will be used to provide two more equations to the original linear system: n‡1=2

ÿ 12…30s ÿ 17†ui;1

n‡1=2

‡ 24…30s ‡ 19†ui;2

ˆ 12s…6s ‡ 11†uni;0

‡ …180s3 ÿ 420s2 ÿ 7s ‡ 204†uni;1 ÿ 8…90s3 ÿ 165s2 ‡ 61s ÿ 57†uni;2 ‡ 18s…60s2 ÿ 100s ‡ 47†uni;3 ÿ 20s…6s ÿ 7†…6s ÿ 1†uni;4 ‡ s…6s ÿ 1†…30s ÿ 17†uni;5 n‡1=2

(14)

n‡1=2

ÿ 12…30s ÿ 17†ui;Mÿ1 ‡ 24…30s ‡ 19†ui;Mÿ2 ˆ 12s…6s ‡ 11†uni;M ‡ …180s3 ÿ 420s2 ÿ 7s ‡ 204†uni;Mÿ1 ÿ 8…90s3 ÿ 165s2 ‡ 61s ÿ 57†uni;Mÿ2 ‡ 18s…60s2 ÿ100s ‡ 47†uni;Mÿ3 ÿ20s…6sÿ7†…6s ÿ 1†uni;Mÿ4 ‡ s…6s ÿ 1†…30s ÿ17†uni;Mÿ5

(15)

The resulting system of linear algebraic equations is tridiagonal and can be solved by using the very fast Thomas algorithm [10]. This system is not diagonally dominant, but the numerical results showed that it is solvable up to s ˆ 2.

334

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

Fig. 1. Computational stencil for the x-sweep of the 3-point BTCS method.

In the above procedure we avoid the use of any boundary values of the original problem at the intermediate time levels. However, the formulae (14) and (15) are stable [9] only for s  1:4986. This affects the unconditional stability of (13). In the x-sweep the following formula is used with i ˆ 1; 2; . . . ; M ÿ 1, for each j ˆ 1; 2; . . . ; M ÿ 1: n‡1=2

n‡1 n‡1 n‡1 ÿsuiÿ1; j ‡ …1 ‡ 2s†ui; j ÿ sui‡1; j ˆ ui; j

:

(16)

The resulting system is diagonally dominant, which guarantees that it is solvable. Values of un‡1 i; j on the boundaries x ˆ 0; 1 and y ˆ 0; 1 for the local problem are provided by the boundary conditions (3)±(6). The computational molecule used for the x-sweep for the 3-point BTCS formula is given in Fig. 1. The modified equivalent partial differential equation of this method in the x-sweep is in the following form [12]: @u @ 2 u x …
x†2 @4u ÿ x 2 ÿ …1 ‡ 6s† 4 ‡ O……
x†4 † ˆ 0: @t @x @x 12

(17)

There is no value of s for which the above formula is fourth-order accurate, so the truncation error is always proportional to (
x)2. The computational molecule used for the x-sweep for the …2; 5† method used in formula (14) and (15) is given in Fig. 2 and the modified equivalent partial differential equation of this method in the x-sweep is in the following form [17]:  @u @2u …
x†4 60228 5 5028 4 24233 3 217679 2 ÿ x 2 ÿ x s ‡ s ÿ s ‡ s 1296s7 ‡ 7344s6 ‡ 5 @t @x 5 5 5 15 …11 ‡ 6s†  877129 424589 @ 6 u ‡ O……
x†6 † ˆ 0: (18) ‡ sÿ 6 60 180 @x

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

335

Fig. 2. Computational stencil for the x-sweep of the implicit (2, 5) method.

2.2. Crank±Nicolson method The Crank±Nicolson method [11] is used in the following. Eq. (8) is solved numerically over the half-time step tn to tn‡1=2 using: n‡1=2

n‡1=2

ÿsui; jÿ1 ‡ 2…1 ‡ s†ui; j

n‡1=2

ÿ sui; j‡1 ˆ suni; jÿ1 ‡ 2…1 ÿ s†uni; j ‡ suni; j‡1

(19)

applied for j ˆ 2; 3; . . . ; M ÿ 2, for each i ˆ 0; . . . ; M. This procedure is unconditionally von Neumann stable and solvable [11] for all s > 0, which is the same as for the 3-point BTCS method in (13). Then for each i ˆ 0; 1; . . . ; M, values at points adjacent to the boundary y ˆ 0 are calculated using the classical FTCS formula in the following form: n‡1=2

ui;1

ˆ suni;0 ‡ …1 ÿ 2s†uni;1 ‡ suni;2

(20)

while values at points adjacent to the boundary y ˆ 1 are calculated using: n‡1=2

ui;Mÿ1 ˆ suni;M ‡ …1 ÿ 2s†uni;Mÿ1 ‡ suni;Mÿ2 :

(21)

Both (20) and (21) are stable only for 0 < s  1=2. The Crank±Nicolson formula is then used to solve Eq. (11) over the time interval tn‡1/2 to tn‡1 as follows. For i ˆ 1; 2; . . . ; M ÿ 1, and each j ˆ 1; 2; . . . ; M ÿ 1, use: n‡1=2

n‡1=2

n‡1 n‡1 ‡ 2…1 ‡ s†un‡1 ÿsuiÿ1;j i; j ÿ sui‡1;j ˆ suiÿ1;j ‡ 2…1 ÿ s†ui; j

n‡1=2

‡ sui‡1;j :

(22)

Note that the above procedure avoids the use of any boundary values of the original problem at the intermediate time levels. However, it is stable only for 0 < s  1=2, so it is only useful in this range when the 3-point FTCS scheme is less accurate, that is when j1 ÿ 6sj < 1, i.e., s > 1=3. Hence, it is more accurate in the range 1=3 < s  1=2.

336

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349 n‡1=2

n‡1=2

The formulae (19), (20) and (21) are also used for i ˆ 0 and i ˆ M, so the values of u0;j and uM;j , j ˆ 0; 1; . . . ; M, which are required in using formulae (22), have already been found. Values of un‡1 i; j on the boundaries x ˆ 0; 1 and y ˆ 0; 1 for the local problem are provided by the boundary conditions (3)±(6). This procedure is unconditionally von Neumann stable, and solvable for all s > 0, which is the same as for the 3-point BTCS method in (13). The formulae (20) and (21) are used to avoid the use of any intermediate values from boundary values of the original problem at the intermediate time levels [8]. Unfortunately, the limited stability of (20) and (21) affects the stability of (19) and (22) when applied to the interior grid points which makes the method have the same range of stability as the LOD …1; 3† FTCS. The best way to overcome this problem is to develop an implicit …2; 4† finite-difference scheme which is unconditionally stable and solvable and use it for i ˆ 1 and i ˆ M ÿ 1 to provide two more equations to the linear system (19). Some efforts were made to develop finite difference schemes in the form: n‡1=2

c1 ui

n‡1=2

‡ c2 ui‡1

ˆ c3 uniÿ1 ‡ c4 uni ‡ c5 uni‡1 ‡ c5 uni‡2 :

(23)

But all were either not diagonally dominant or they had a very limited range of stability. However, one would use (14) and (15) and get better results than with (20) and (21). The computational molecule of the x-sweep of this method is given in Fig. 3 and the modified equivalent partial differential equation for the Crank±Nicolson formula in the x-sweep is in the following form [12]: @u @ 2 u x …
x†2 @ 4 u ÿ x 2 ÿ ‡ O……
x†4 † ˆ 0: @t @x 12 @x4

(24)

This equation has a truncation error which is always O……
x†2 † so is second-order accurate for all s > 0 in the range of stability. Also, note that it has a leading error which is independent of s.

Fig. 3. Computational stencil for the x-sweep of the Crank±Nicolson method and the x-sweep of the Crandall's scheme.

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

337

2.3. The fourth-order Crandall's method Crandall's method [13] is used in the following: Eq. (8) is solved numerically over the half-time step tn to tn‡1/2, by using: n‡1=2

n‡1=2

…1ÿ6s†ui; jÿ1 ‡ 2…5 ‡ 6s†ui; j

n‡1=2

‡ …1ÿ 6s†ui; j‡1 ˆ …1‡ 6s†uni; jÿ1 ‡ 2…5 ÿ 6s†uni; j ‡ …1‡ 6s†uni; j‡1

(25) for j ˆ 2; 3; . . . ; M ÿ 2, with each i ˆ 0; . . . ; M. This procedure is unconditionally von Neumann stable [13] for all s > 0. To avoid any calculation using the intermediate values, the fourth-order implicit formulae (14) and (15) for i ˆ 0; 1; . . . ; M will be used to provide two more equations to the original linear system. Crandall's method is then used to solve Eq. (11) over the time interval tn‡1/2 to tn‡1, as follows: For i ˆ 1; 2; . . . ; M ÿ 1, and each j ˆ 1; 2; . . . ; M ÿ 1, use: n‡1=2

n‡1=2

n‡1 n‡1 …1ÿ 6s†un‡1 iÿ1;j ‡2…5 ‡ 6s†ui; j ‡…1ÿ 6s†ui‡1;j ˆ …1‡6s†uiÿ1;j ‡ 2…5 ÿ 6s†ui; j

n‡1=2

‡ …1 ‡ 6s†ui‡1;j : (26)

In the above procedure we avoid the use of any boundary values of the original problem at the intermediate time levels. However, formulae (14) and (15) are stable only for s  1:4986. This affects the unconditional stability of (25) and (26). This formulae is not diagonally dominant, but the numerical results showed that it is solvable only if s  0:5. Values of un‡1 i; j on the boundaries x ˆ 0; 1 and y ˆ 0; 1 for the local problem are provided by the boundary conditions (3)±(6) and the computational molecule of the x-sweep of this method is the same as the Crank±Nicolson formula shown in Fig. 3. The modified equivalent partial differential equation for the Crandall's formula in the x-sweep has the following form [12]: @u @2u …
x†4 @6u ÿ x 2 ÿ x …20s2 ÿ 1† 6 ‡ O……
x†6 † ˆ 0: @t @x @x 240

(27)

p This equation has a truncation error which is generally O……
x†4 †, except for s ˆ 1= 20, when it is sixth-order p for the local problem. This would only be useful if (14) and (15) was also sixth-order for s ˆ 1= 20. 3. The non-local boundary treatment The presence of an integral term in a boundary condition can greatly complicate the application of standard numerical techniques. The accuracy of the quadrature must be compatible with that of the discretization of the differential equation. In this case it means that a higher order of accuracy than the second-order composite trapezoidal rule [16] is required. In the following, the fourth-order Simpson's composite ``one-third'' rule [14] is used. Consider the integral: Z d…x† u…x; y; tn‡1 † dy: (28) H…x; tn‡1 † ˆ 0

338

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

Application of Simpson's composite ``one-third'' rule gives: Z

1

0

…M=2†ÿ1 M=2 X X h n‡1 n‡1 H0n‡1 ‡ 4 H…x; tn‡1 † dx ' H2iÿ1 ‡2 H2in‡1 ‡ HM 3 iˆ1 iˆ1

Hin‡1

Z ˆ

d…ih†

0

!

u…xi ; y; tn‡1 † dy

(29)

(30)

in which: Hin‡1 ˆ

Z

2li h

0

u…xi ; y; tn‡1 † dy ‡

Z

d…ih† 2li h

u…xi ; y; tn‡1 † dy

(31)

where li ˆ ‰d…ih†=2hŠ;

(32)

and ‰
Š represents the integer part of the argument. Substituting in the second integral of (31): zi ˆ y=h ÿ 2li

(33)

yields: Z

d…ih†

2li h

u…xi ; y; t

n‡1

Z † dy ˆ h

i 0

u…xi ; zi ; tn‡1 † dzi

(34)

where i ˆ …d…ih†=h† ÿ 2li :

(35)

Replacement of u in the integral with a quadratic interpolating polynomial (the Newton's forwarddifference formula) through the grid values concerned gives [15] Z 0

i

u…xi ; zi ; tn‡1 † dzi ˆ

Z

i 0

‰u…xi ; y2li ; tn‡1 †

‡ zi
u…xi ; y2li ; tn‡1 † ‡ 12zi …zi ÿ 1†
2 u…xi ; y2li ; tn‡1 †Š dzi ‡ O…h4 †:

(36)

where
u…xi ; y2li ; tn‡1 † ˆ u…xi ; y2li‡1 ; tn‡1 † ÿ u…xi ; y2li ; tn‡1 †

(37)


2 u…xi ; y2li ; tn‡1 † ˆ u…xi ; y2li‡2 ; tn‡1 † ÿ 2u…xi ; y2li‡1 ; tn‡1 † ‡ u…xi ; y2li ; tn‡1 †:

(38)

and

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

339

Integrating (37) and collecting like terms then substituting in (31) means that: li li ÿ1 X X h n‡1 n‡1 ˆ ‰u…xi ; 0; t † ‡ 4 u…xi ; y2jÿ1 ; t † ‡ 2 u…xi ; y2j ; tn‡1 † ‡ u…xi ; y2li ; tn‡1 † 3 jˆ1 jˆ1

Hin‡1

‡ 3i …1 ÿ 3i =4 ‡ i2 =6†u…xi ; y2li ; tn‡1 † ‡ 3i2 …1 ÿ i =3†u…xi ; y2li ‡1 ; tn‡1 † ‡ …i2 =4†…2i ÿ 3†u…xi ; y2li ‡2 ; tn‡1 †Š ‡ O…h4 †:

(39)

However, at all interior points we have computed u…xi ; yj ; tn‡1 † to say, rth-order, where: r u…xi ; y2li ; tn‡1 † ˆ ui;n‡1 j ‡ O…h †:

(40)

Substituting these approximations in (40) gives: " li li ÿ1 X X h n‡1 n‡1 n‡1 n‡1 2 n‡1 u ‡4 ui;2jÿ1 ‡ 2 un‡1 Hi ˆ i;2j ‡ ui;2li ‡ 3i …1 ÿ 3i =4 ‡ i =6†ui;2li 3 i;0 jˆ1 jˆ1 i 2 n‡1 q ‡3i2 …1 ÿ i =3†un‡1 ‡ … =4†…2 ÿ 3†u i i;2li ‡1 i i;2li ‡2 ‡ O…h †

(41)

where q ˆ minfr; 4g: Putting: V…t

n‡1

Z †ˆ

1

0

(42)

H…x; tn‡1 † dx

(43)

and using the approximation: vn‡1

…M=2†ÿ1 M=2 X X h n‡1 n‡1 n‡1 H0 ‡ 4 ˆ H2iÿ1 ‡ 2 H2in‡1 ‡ HM 3 iˆ1 iˆ1

then gives: n‡1

v

…M=2†ÿ1 M X X h2 n‡1 n‡1 n‡1 ˆ u2iÿ1;0 ‡2 un‡1 u0;0 ‡ 4 2i;0 ‡ uM;0 9 iˆ1 iˆ1

Note that: n‡1

v

h ˆ n‡1 3

Z

1 0

! ‡ O…hq †

(44)

‡ Rn‡1 ‡ O…hq †:

(45)

!

h0 …x† dx ‡ Rn‡1 ‡ O…hq †;

(46)

where Rn‡1 is the summation in vn‡1 excluding the values at the boundary y ˆ 0. Since vn‡1 is an approximation to the left hand-side of (7) it follows that: n‡1



mn‡1 ÿ Rn‡1 ‡ O…hqÿ1 †; ˆ R1 …h=3† 0 h0 …x† dx

Z

1 0

h0 …x† dx 6ˆ 0:

(47)

340

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

As seen above, the order of convergence of  depends on two things: firstly, the order of the finitedifference formula used at interior grid points, and secondly, the order of the numerical quadrature used to approximately evaluate (7). For example, if un‡1 is evaluated using a fourth-order formula, when q ˆ 4; n‡1 is only third-order convergent. If un‡1 is found at interior points by a second formula when q ˆ 2 then n‡1 is only first-order convergent. In order to handle the non-local boundary condition (7) the numerical integration procedure which was described is used: n‡1 ˆ

mn‡1 ÿ Rn‡1 ‡ O…hqÿ1 †; R1 …h=3† 0 h0 …x† dx

q ˆ minfr; 4g:

(48)

The order of convergence of  is 1 if un‡1 is evaluated using the second-order 3-point BTCS formula or the Crank±Nicolson method, and it is 3 if un‡1 is evaluated using the fourth-order Crandall's formula. Note that in solving the linear systems (16), (22), or (26), the values of un‡1 i;0 are not needed and the . So in the non-local problem the values at the interior grid points can be found before evaluating un‡1 i;0 numerical integration described can be employed without any difficulty and hence the main advantage of first solving the two-dimensional equation in the y-sweep can be seen. If the reverse order is taken it means that the much slower iteration technique described in [16] must be used to handle the non-local boundary condition. The same is true for the other schemes developed here. This procedure is, therefore, much more efficient than using fully implicit methods. 4. Numerical test Consider (1)±(7) with x ˆ y the functions f, h0, h1, g0, g1, m,  and d as defined by following: f …x; y† ˆ exp…x ‡ y†

(49)

g0 …y; t† ˆ exp…y ‡ 2t†

(50)

g1 …y; t† ˆ exp…1 ‡ y ‡ 2t†

(51)

h0 …x† ˆ exp…x†

(52)

h1 …x† ˆ exp…1 ‡ x ‡ 2t†

(53)

…t† ˆ exp…2t†

(54)

m…t† ˆ …exp…1† ÿ 1†2 exp…2t†

(55)

d…x† ˆ x…1 ÿ x†

(56)

for which the exact solution is: u…x; y; t† ˆ exp…x ‡ y ‡ 2t†:

(57)

In Table 1 the results are shown for uNi; j with h ˆ 0:05 and s ˆ 1=2 at T ˆ 1:0, using the 3-point BTCS and the Crank±Nicolson method and Crandall's method, with given boundary values

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

341

Table 1 Results for u with T ˆ 1:0; h ˆ 0:05; s ˆ 1=2 x 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

y 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Exact u 9.025013 11.023176 13.463738 16.444647 20.905243 24.532530 29.964100 36.598234 40.447304

Error LOD BTCS

LOD Crank±Nicolson

ÿ0.210ÿ3 ÿ0.710ÿ3 ÿ0.110ÿ2 ÿ0.20ÿ2 ÿ0.2 10ÿ2 ÿ0.210ÿ2 ÿ0.210ÿ2 ÿ0.110ÿ2 ÿ0.510ÿ3

ÿ0.610ÿ4 ÿ0.210ÿ3 ÿ0.310ÿ3 ÿ0.510ÿ3 ÿ0.610ÿ3 ÿ0.610ÿ3 ÿ0.510ÿ3 ÿ0.40ÿ3 ÿ0.110ÿ3

LOD Crandall ÿ0.210ÿ ÿ0.810ÿ7 ÿ0.210ÿ6 0.210ÿ6 ÿ0.310ÿ6 ÿ0.310ÿ6 ÿ0.210ÿ6 ÿ0.110ÿ6 ÿ0.310ÿ7

everywhere. Note that the errors obtained when using the 3-point BTCS formula are generally about 10 000 times larger, and when using the Crank±Nicolson scheme are generally a 1000 times larger than those produced using Crandall's method. When the absolute value of the error: eni; j ˆ u…ih; jh; nk† ÿ uni; j

(58)

at the point (0.5, 0.5) at time T ˆ 1:0 was graphed against h on a logarithmic scale for various values of s, it was found that the slope of lines was always close to 2 for the 3-point BTCS formula and the Crank±Nicolson formula, and close to 4 for Crandall's formula (see Figs. 4±6). These results illustrate the theoretical orders of accuracy evident from the modified equivalent equation discussed earlier in Section 2 for the local problem. Fig. 5 shows that the accuracy of the Crank±Nicolson method is not altered as s increases. This is because of the fact that, for the same value of h, the leading error term in (24) does not depend on the

Fig. 4. Relation between error in u and grid spacing for the 3-point BTCS method.

342

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

Fig. 5. Relation between error in u and grid spacing for the Crank±Nicolson method.

Fig. 6. Relation between error in u and grid spacing for the Crandall's scheme.

value of s. It is clear from Fig. 6 that the greatest error obtained when using the Crandall's method is less than the smallest error obtained when using the 3-point BTCS scheme or the Crank±Nicolson method. The results for  with h ˆ 0:05; s ˆ 1=2 using three methods to solve the non-local boundary value problem are shown in Table 2. This table shows that the results obtained when using the 3-point BTCS method or the Crank±Nicolson scheme have errors about 100 times larger than the results when the Crandall's formula is applied. When the absolute value of the error: en ˆ …nk† ÿ n

(59)

at the point (0.5, 0.5) at time T ˆ 1:0 was graphed against h on a logarithmic scale for various values of s, it was found that the slope of lines was always close to 1 for the 3-point BTCS formula, and the

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

343

Table 2 Results for  with h ˆ 0:05; s ˆ 1=2 t 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Exact  1.221403 1.491825 1.822119 2.225541 2.718282 3.320117 4.055200 4.953032 6.049647 7.389056

Error LOD BTCS

LOD Crank±Nicolson

LOD Crandall

0.210ÿ2 0.310ÿ2 0.310ÿ2 0.410ÿ2 0.510ÿ2 0.610ÿ2 0.710ÿ2 0.910ÿ2 0.110ÿ1 0.110ÿ1

0.510ÿ3 0.710ÿ3 0.810ÿ3 0.110ÿ2 0.110ÿ2 0.210ÿ2 0.210ÿ2 0.210ÿ2 0.310ÿ2 0.310ÿ2

0.110ÿ4 0.110ÿ4 0.210ÿ4 0.210ÿ4 0.210ÿ4 0.310ÿ4 0.410ÿ4 0.510ÿ4 0.610ÿ4 0.710ÿ4

Crank±Nicolson formula, but near 3 for the Crandall's formula (see Figs. 7±9). These results reflects the orders of convergence referred to earlier in Section 2. The interesting feature of these figures is that the maximum discretization error produced by Crandall's scheme is smaller than the minimum discretization error obtained when using the 3-point BTCS method or the Crank±Nicolson scheme.

Fig. 7. Relation between error in  and grid spacing for the 3-point BTCS method.

Fig. 8. Relation between error in  and grid spacing for the Crank±Nicolson method.

344

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

Fig. 9. Relation between error in  and grid spacing for the Crandall's scheme.

5. Central processor time The theoretical estimate of the CPU time taken for one time step of each method follows. Tests showed that, for the Sun IPX SPARC station, the operations of addition, subtraction and multiplication requires practically the same amount of CPU time, about 0.86 ms. In the following, this is considered as one unit of time. It can be seen that the CPU time needed for the local problem when using each of the methods described here is proportional to: TBTCS ˆ …M ÿ 1†…6…M ÿ 1† ‡ 6  2 ‡ 5  2† ‡ …M ÿ 1†…6…M ÿ 1† ‡ 2  2† ˆ 2…6M 2 ‡ M ÿ 7† (60) TCN ˆ …M ‡ 1†…6…M ÿ 3† ‡ 5…M ÿ 5† ‡ 7  2 ‡ 5  2† ‡ …M ÿ 1†…5…M ÿ 3† ‡ 6…M ÿ 1† ‡ 7  2† ˆ 2…11M 2 ÿ 13M ÿ 6†

(61)

and TCR ˆ …M ‡ 1†…6…M ÿ 1† ‡ 5…M ÿ 3† ‡ 6  2 ‡ 5  2† ‡ …M ÿ 1†…5…M ÿ 3† ‡ 6…M ÿ 1† ‡ 7  2† ˆ 2…11M 2 ÿ 3M ‡ 4†:

(62)

The CPU time used for the 3-point BTCS, the Crank±Nicolson scheme and the Crandall's method, respectively, are reported in the Tables 3±5. The number of operations needed at each time step for the three methods developed in this paper with s ˆ 1=6 are given in Table 6, and the CPU time required for s ˆ 1=6 is reported in Table 7. The absolute value of the discretization error at the point …0:5; 0:5† at time T ˆ 1:0 is graphed against the CPU time on a logarithmic scale for various values of s in Figs. 10±12. The theoretical ratio of the CPU time required for the BTCS procedure relative to the Crank± Nicolson procedure is therefore: R1c ˆ TBTCS =TCN ˆ

2…6M 2 ‡ M ÿ 7† : 2…11M ÿ 13M ÿ 6†

(64)

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

345

Table 3 CPU time in seconds for the BTCS scheme M

s ˆ 0:25

s ˆ 0:33

s ˆ 0:50

20 30 40 50 60 70

15.2 69.7 225.8 511.6 1060.6 1955.2

13.3 50.6 181.8 397.1 801.2 1470.5

10.2 34.6 119.8 260.5 603.6 989.4

Table 4 CPU time in seconds for the Crank±Nicolson scheme M

s ˆ 0:25

s ˆ 0:33

s ˆ 0:50

20 30 40 50 60 70

19.3 97.8 324.1 734.7 1539.7 2887.3

14.6 72.8 249.9 576.5 1168.6 2146.2

12.7 49.4 159.5 373.4 781.5 1430.7

Table 5 CPU time in seconds for the Crandall's scheme M

s ˆ 0:25

s ˆ :33

s ˆ 0:50

20 30 40 50 60 70

21.2 103.2 340.5 770.4 1608.6 3032.5

15.9 77.1 261.7 605.3 1217.8 2250.6

12.8 52.5 168.6 384.0 817.1 1503.4

Table 6 CPU time in seconds for s ˆ 1=6 M 20 30 40 50 60 70

CPU time LOD BTCS

LOD Crank±Nicoloson

LOD Crandall

16.0 75.5 253.0 594.9 1135.5 2111.0

25.7 123.8 421.6 924.3 1947.7 3652.2

27.2 130.2 441.5 962.8 1997.2 3724.1

346

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

Table 7 Number of operations in each step M

TBTCS

TCN

TCR

20 30 40 50 60 70

4826 10846 19266 30086 43306 58926

8268 19008 34148 53288 77628 105968

8688 19628 34968 54708 78848 107388

Similarly, the theoretical ratio of the CPU time required for the Crank±Nicolson method relative to Crandall's method is: R2c ˆ TCN =TCR ˆ

2…11M 2 ÿ 13M ÿ 6† : 2…11M ÿ 3M ‡ 4†

(65)

As M increases R1c and R2c tends to the limiting value of 6=11 ˆ 0:5454 . . . and 1.0, respectively. The theoretical and experimental ratio of the CPU time required for the schemes discussed in this article are, respectively, reported in Tables 8 and 9.

Fig. 10. Relation between the CPU times and the error for the 3-point BTCS method used in the LOD procedure.

Fig. 11. Relation between the CPU times and the error for the 3-point Crank±Nicolson method used in the LOD procedure.

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

347

Fig. 12. Relation between the CPU times and the error for the 3-point Crandall's scheme used in the LOD procedure.

Table 8 The theoretical and experimental ratio of the CPU time required for the BTCS scheme relative to the Crank±Nicolson scheme M 20 30 40 50 60 70

R1c Theoretical

Experimental

0.584 0.571 0.564 0.564 0.557 0.556

0.621 0.610 0.601 0.595 0.583 0.578

Table 9 The theoretical and experimental ratio of the CPU time required for the Crank±Nicolson scheme relative to the Crandall's scheme M 20 30 40 50 60 70

R2c Theoretical

Experimental

0.951 0.968 0.974 0.977 0.985 0.987

0.944 0.951 0.955 0.960 0.974 0.975

348

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

6. Conclusion In this paper three methods, namely the second-order 3-point BTCS scheme, the second-order Crank±Nicolson method and the fourth-order Crandall's method are used to solve the two-dimensional diffusion equation with both local and non-local boundary conditions through a LOD procedure which employed those one-dimensional schemes to apply them in each direction. In order to prevent losing the accuracy a suitable treatment on the boundaries at the intermediate time level was suggested. So, near the boundaries in the first half-time interval the 3-point FTCS formula was used for the Crank± Nicolson method and a fourth-order …2; 5† implicit finite difference formula was used for the 3-point BTCS scheme and the fourth-order Crandall's method. The latter worked very well for twodimensional non-local diffusion problem because of its low error. These conditionally stable methods seems particularly suited in the sense of implementing easily for parabolic partial differential equations with continuous boundary conditions. A comparison with the fully implicit schemes for the model problem clearly demonstrates that the new techniques use less CPU time. The only disadvantage of these methods was their limited range of stability. This was because of avoiding the use of the boundary values at the intermediate time levels, as this made these procedures to be dependent to some other conditional schemes to evaluate the values near the boundaries. In the LOD procedure first the equation was solved in the y-sweep. With this appropriate choice, using the much slower iteration technique described in [16] was avoided. References [1] S. Wang, Y. Lin, A finite difference solution to an inverse problem determining a control function in a parabolic partial differential equation, Inverse Problems 5 (1989) 631±640. [2] J.R. Cannon, J. van der Hoek, Diffusion subject to the specification of mass, J. Math. Anal. Appl. 115 (1986) 517±529. [3] J.R. Cannon, S. Prez-Esteva, J. van der Hoek, A Galerkin procedure for the diffusion equation subject to the specification of mass, Siam. J. Numer. Anal. 24 (1987) 499±515. [4] J.R. Cannon, J. van der Hoek, Implicit finite difference schemes for the diffusion of mass in porous media, in: B.J. Noye, (Ed.), Numerical Solutions of Partial Differential Equations, North-Holland, Amsterdam, 1982, pp. 527±539. [5] J.R. Cannon, A. Matheson, On the numerical solution of diffusion subject to the specification of mass, Int. J. Eng. Sci. 31 (1993) 347±355. [6] J.R. Cannon, Y. Lin, S. Wang, An implicit finite difference scheme for the diffusion equation subject to mass specification, Int. J. Eng. Sci. 28 (1990) 573±578. [7] W.A. Day, A decreasing property of solutions of a parabolic equation with applications in thermoelasticity and other theories, Q. Appl. Math. 41 (1983) 475±486. [8] A.R. Mitchell, D.F. Griffiths, The Finite Difference Methods in Partial Differential Equations, Wiley, New York, 1980. [9] B.J. Noye, K.J. Hayman, LOD and ADI methods for the two-dimensional diffusion equation, Int. J. Comput. Math. 51 (1994) 215±228. [10] L.H. Thomas, Elliptic problems in linear difference equations over a network, Technical Report, Watson Scientific Computing Laboratory, Columbia University, New York, 1949. [11] J. Crank, P. Nicolson, A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type, Proc. Cambridge Phil. Soc. 43(1) (1947) 50±67. [12] R.F. Warming, B.J. Hyett, The modified equation approach to the stability and accuracy analysis of finite difference methods, J. Comput. Phys. 14(2) (1974) 159±179. [13] S.H. Crandall, An optimum implicit recurrence formulae for the heat conduction equation, Q. Appl. Math. 13 (1955) 318±320.

M. Dehghan / Mathematics and Computers in Simulation 49 (1999) 331±349

349

[14] C.F. Gerald, Applied Numerical Analysis, 5th ed., Addison-Wesley, California, 1995. [15] S. Wang, A numerical method for the heat conduction subject to moving boundary energy specification, Numer. Heat Transfer 130 (1990) 35±38. [16] J.R. Cannon, Y. Lin, A. Matheson, The solution of the diffusion equation subject to specification of mass, Appl. Anal. 50(2) (1993) 1±11. [17] B.J. Noye, Numerical Solution of Partial Differential Equations. (lecture notes) (1992).