Fundamental solutions of nine-point discrete Laplacians

Applied Numerical Mathematics 10 (1992) 325-334 North-Holland

325

APNUM341

Robert E. Lynch Departments of Computer Sciences and Mathematics, Purdue University, West Lafayette, IN 47907, USA

Abstract

Lynch, R.E., Fundamental solutions of nine-point discrete Laplacians, Applied Numerical Mathematics 10 (1992) 325-334. Fundamental solutions of a family of nine-point discretizations of the Laplace operator and their asymptotic expansions are presented.

We analyze the “fundamental

solutions” G of the discrete Laplace operator:

(1a ) subject to

first differences of Gj,k( Cu)tend to zero as j2 + k 2 + 00,

(lb)

where La is the nine-point difference operator defined by

L&

= (2a - 4)& + (1 - qq-l,k +

+a["j-I,k-l

+

U;.+l,k-1

+ q+l,k f q&-l + &+,I +

uj+l,k+*

+

LS.-l,k+i]*

In stencil form this is

Correspondence to: R.E. Lynch, Departnent of Computer Sciences, Ptudue LJniversity, West Lafayette, IN 47907, USA. Telephone: (317) 494-1997. Fw: (3j 7) 444-0739. E-mall: [email protected]. 0168-9274/92/$05.00

0 1992 - Elsevier Science Publishers B.V. All rights reserved

R.E. Lynch / Nine-point discrete Laplacians

326

For smooth u with Uj,k = U(jh, kh), hs2L,uj ,k = V2Uj 5k + 0(h2).

When cy= 0, L, is the “standard” (“five-point star”) operator of discrete potential theory. For cy= $, it is Polya’s bilinear finite element approximation of the Laplacian (see Birkhoff and Lynch [2, pp. 190~1.911).The only value of (Ywhich yields a higher-than-order-2 accuracy for the Laplace or the Poisson equation is a! = $; this gives the O(h6) “optimal” nine-point discrete Laplacian for which V2Uj1k + ~h2V4uj ,k + ~h4 h-'L,/3Ujk= ,

a4

V6Uj ,k % 2 -V2Uj,k

ax2ay2

1+

0(h6)

(for application to the Poisson equation, see [2, p. 921). With A$$a) denoting the asymptotic expansion through O( R-“) terms, we obtain the result 2vGj 9k((Y)N 27~Ayl(a) , =logR+y+; +-

+ 0(Rm6)

8 log / 1 \ -

+(3a-l)-

cos 40 12R2

(90 0L2- i8) cos40 - (22%~~- 150~~+ 25) cos80 + 0( R-6), 240R4

where R2 = j2 + k2, y = 0.5772156649..

.

(2)

is Euler’s constant, and o = arctam k/j).

revious results

We are unaware of published results about G(a) with a! different from zero. The solution of [la>--(lb) with Q = 0 is well known (see McCrea and Whipple [9]-or Stohr [13,III], Sobolev [10,11], Duffin [4], Duffin and Shelly [6], van der Pol [14]). De Boor, Hiillig and Riemenschneidcr [3] do discuss more general difference operators, but use (la) with a! = 0 as an example. They, and some of the others sources, give the “diagonal” values of G(0): j+ 1, &-2,....

(3)

From these and G,-,(O) = 0, all the other vaiues can be computed by use of the difference equation and symmetry. Duffin and Shaffer [5] derived the complete asymptotic expansion of the doubac Fourier transforms of modifications of functions of the form A myn + F,(x, y ), where r2 = x 2 + Y2, where m and II are nonnegative integers, and where F, has partial derivatives of all orders. They give [5, Theorem 41 27FGj,k(O)Y A:2i(0) + 0( Rm4) 5 cos 4U = log R + y + ; log 2 - 12R2 + 0(R-4).

(4)

R.E. Lynch / Nine-point discrete Laplacians

327

Earlier, McCrea and Whipple 19, (9.6)] gave A$(O), the first three terms on the rightmost side of (4).

3. Analysis As can be verified by direct substitution, a solution of (la) can be written as

(eijX+ikY -

1 Gj,k(a)

/‘“I’”2eu-4+2((1-~)[cos~+co~y]

s

=

0

1) dx

0

dy

+acosxcosy)’

(5>

where the denominator in the integrand is equal to (Lcueijx+iky)/eijx+i’y. Clearly G,,(a) = 0 and it is a consequence of the asymptotic formula (2) that (5) satisfies the second boundary condition (lb). Like the case a! = 0 treated by Duffin and Shelley [6] and Duffin and Shaffer [S], this double integral is absolutely convergent, because in the integrand the numerator is O(r) as r2 = x2 + y2 &0 and the denominator is O( r2). Thus it is permissible to evaluate the double integral as an integrated integral. The integration with respect to x in (5) can be carried out explicitly in terms of a contour integral. With the substitutions w = eiX, dx = dw/iw, that integral is equal to 9/

=

Ceiix - 1

2v

0

cw’- 1 w2- (2A/B)w

1 dx = _si

2A-2Bcosx

where the contour is the unit circle:

w =

eiX_The denominator

-1 1 dw’

(6)

in the integrand has zeros

--

A-@-B2 W_=

and

B

-

w+= l/w_.

-

The simple pole at w _ is the only singularity inside the circle, whence

I

C{[ A

3=7T

1

v%%+B)I

LFiF

(7)

We now determine G,,,(a) in (5) from (6) in which A=a-2+(1-cw)cosy, B=a-l-acosy, C = eiky. Set t = u + iv = eiy; then A2 - B2 = [3 - 2~ - (1 - 2~) cos y](l -- cos y)

cy- 1) = (2 t2 L;f 2 (

-bat

+ l)(t - 1)2

(8)

R.E. Lynch / Nine-point discrete Laplacians

328

Fig. 1.

with and

6-4~~ bcr= --_ I-2cu

D

(Y=

3-2a-&i3 l-2CY

l

Specifically, for a < & 0, is positive; 01/z = 0; for i < a! < 1, Da is negative; and for 1 < cy, 0, is complex with 1Da I= 1. For cy< $ and 4 < cyc 1, equation (5) becomes (see !6), (7), and (8)) 1 GO,k(a)

=

Because ( t k - 1) =

-

(t -

z

k-l I’==,;;;-

l)(l -t

1

-dt. t - 1 ,/(~cY- l>(t - DJ(t - l/D,)

!----

+ t* +

l

l

l

+tk-I)

and G&cr) = 0, it follows that

k-=0, l,....

(9) Figure l(a) illustrates the situation for 0 <(Y < : (when 0.17.. . >, 0, > 0) and Fig. l(b) illustrates the situation for $ < LY< 1 (when 0 > 0, > - 1); the path of integration in (9) is along the arc abc of the unit circ!e. The integrand has branch points at 0, and l/D,, indicated by “bullets” in the figure; the dashed line indicates the branch cut we take between these points, beginning at Da and going to the right along the real axis to l/D,. The integral is equal to zero along the closed contour abcdefa which does not enclose any singularities. With t = 0, + pe”, dt = pie” d6, and p c 1, the integrand is O(G); hence the contribution along the arc def tends to zero as p 10. With t = Da + peie, dt = eie dp, the integrals along the intervals cd and fa are -

and -

I-D, / p=o

(Da + pei2T)k fi e’“(l/D,

-Da -pe’*“)

ei2Tdp

H.E. Lynch / Nine-point dkcrete Laplacians

respectively. Consequently, (9) becomes Go,k+,(a) = G, 3,(d Jam

=

1

1

Uk -&u

-

Tr/ u=&

329

+J;“#(d, where

with U= (2cy - l)(u*-&UU + 1).

Given the values of Y0 and Z., the other values of ,a;Ccan be found by a recurrence relation. Because Ukdi7 = (2a - 1)

(k + l)uk+’ - (k + $)uk + kuk-’

m

we have

For 0 < Q!< $ and for $ < QI< 1, the value of the integral Z&a) is given by Dwight [7, p. 75, X380.001]. We can use either to evaluate the limit as a! tends to $. Thus with a! = f + p, 1

\/2p+l log [ /--

-%o(i+P) = Tm

log[l+m+p/4+

I ***I

n/S ‘1-k o(@). Tr In summary, we have

GO,lW=4w

=

I‘9

if QI=&

T

[&-T-+1

1

(11)

if $
Finally, from Dwight [7, p* 75, #380.011] we obtain 1 sr(a)

= ,i-

LtDU$dU

2~~3 = =-@)

2 + Tr(2LY- 1)

(12) l

The values Gj,k(~) of the fundamental solution of the “optimal” nine-point written as (see Table 2) Gjk(t)=A,.+Bjk/7F+Cik~, . _ 9 ,

1

bPk&n

Cm be

R.E. Lynch / Nine-point dkcrete Laplacians

330

where the coefficients A, 3, and C, are rational numbers. Garrett Birkhoff ’ pointed out that since 1, l/v, and fi are rationally independent, the coefficients Bi,k and Cj,k must satisfy = 0 exactly. Similarly, Gj,k(O)= Dj,k i- Ej,k/T, where D and E are rational. L1/3”j,k

4. Asymptotic expansions The asymptotic expansion for Gi_i(0) in (3) is easily obtained from the expansion of the .P Ia logarithmic derivative of the gamma function. Combine (6.3.2) and (6.3.18) of Abramowitz and Stegun [l, pp. 258-2591 to get 1 l+z+***+:

1 J-1

mlogj+y-22--

1‘

1

1 + --’ 12j2 120j4

t. -

B - 2m - . . . 2mj2m

(13)

7

where B, is the mth Bernoulh number. Rewrite this relationship with 2j in place of j and subtract from it one-half of (13), then set j = Rfi and simplify to obtain 1 7 2nGj i(O) w log 2 + 12~2 - -240~4 +

B2m(22m-1- 1) l

l

l

+

22m+lmR2m

+

... .

From the analysis of Duffin and Shaffer [5], a straightforward calculation (see below) yields the asymptotic expansion of Gj,k(a) to as many terms as one wants for general a! < 1; in particular one gets A;‘&!) through O(l/R4) terms as displayed in (2). Note that for the ‘“optimal” nine-point discretization of the Laplacian, L ,,3, the coefficient of R-* in (2) is zero. We outline how (2) was obtained (more details are given in Lynch [8]). The function F used in the derivation of Duffin and Shaffer [5, p. 5941 is replaced with F(x, y; r~) = - 12~ - 4 + 2((1 - cy)[cos X + cos y] + fXcos X cos y]] -*

(14)

and from (7) above, their double integral Y2 reduces to Y2 = /“dypF(x, & /t

=

dY

?r

=

y; cw)dx

~[3-2a!-(1-2c~)cosy](1-cosy)

2/ T

0

’ Personal communication.

-log&+3

log

(15)

R.E. Lynch / Nine-point discrete Laplaciam

331

Table 1 Fundamental solution Gj,k(O) k=4

-

-

-

k=3

-

-

-

k=2

-

k=l

2

k=O

1

-

7F

0

1

j=O

4 j=l

12 --5v

2

1

r-z 1 --

1

-+3v

z

I

105T

23 4

-

176

-

2

3-z

23 -2 3rr

40 --7r

17

12

IT j=2

T

4 118

4

--4

1

20-z

49 4 184

j=4

j=3

where the change of variable V(y;

1 + cos y

a) = 3-2a-(l-2ar)cosy

l/2

1

(16)

is a slight modification of the ingenious change V( y; O>used by Duffin and Shaffer. The rest of their analysis [5, pp. 594-5951 remains unchanged to produce 8 2vGj,k(a)

m 2r~I$“l(a) = log R + 7 + $ log 1 (

-

1

.

The other terms of the asymptotic expansion are obtained from the Taylor series representation, F,(rcos&

rsin8;

CY)= -F(rcos&

rsin8;

(Y)-F2

= n,(e) + A1(8)r2 + A2(8)r4 +

.. l

?

and the use of [5, Theorem 21 (see also [5, pp. 595-5961).

Lists of values of Gj ,(O) for small values of j and kcare given in most of the references given in Section 2 above. A few values are given in Tables 1 and 2, Tables 3 and 4 list some values of Gj,,(1/3). Relative errors,

(17)

-

0

k=l

k=O

j=O

-

1

I

k=2

k=3

k=4

k-5

Table 3 Fundamental

k=4 k=3 k=2 k=l k=O

Table 2 Fundamental

j=l

6

&-

3

j=2

46 --3 7T

6

12 3& J-+2 --6

26

3

r

3 --2

886

78--_-

84

j=l

j=2

0.42441 0.38662 0.36338

j=3

276 --2

36 --7F

516 -+?r

3

586

T

72

6

-522 + 45 2

12416

11475 5976 - -22140 2 T

-

0.31831 0.25000

in decimal notation

Gj,,(1/3)

I

GJO)

-

solution

solution

6

102385fi

5412806 3

j=4

2646 --P 3

T

840

816 607fi -_------_84 w 6

3438 2820 - - 9966 r

46824 -+ T

468888-

-

- 44460

491280 - v

j=3

0.48808 0.46221 0.44035 0.43028

w

2

3

464444986

2

IT

j=5

6

1122496

8844 --1806fi+y r

v

--

6

?r

3

627

- 13782

10200

7668 + 39289fi

----

4064088 + 2990043fi 2 - 3883982 v 571437 294960 332286

42117flO

80442003

j=4

0.53355 0.5 1394 0.49596 0.48240 0.47699

3

s d am.

$ c, z R

$ -. ;3

$

c S * \

?p PJ

R.E. Lynch / Nine-point discrete Laplacians Table 4 Fundamental solution G&/3) _-

k=6 k=S k=4 k=3 k=2 k=l k=O

I

333

in decimal notation

0

0.34530 0.28868

0.455 15 0.41779 0.39954

0.5 1963 0.49374 0.47284 0.46437

0.56541 0.54577 0.52801 0.5 1506 0.5 1022

0.60092 0.58513 0.57023 0.55757 0.54887 0.54575

0.62994 0.61674 0.60404 0.59253 0.58316 0.57695 0.57477

j=O

j=l

j=2

j=3

j=4

j=5

j=6

Table 5 Rel&ative error SCO; 4) k=4 k=3 k=2 k=l k=O

I

-

-

-

j=O

j-1

l.O6E-03 1.38E-01

-.,

-

1.58E-03 1.32E-03 2.25E-03

1.37E-03 1.72E-03 8.07E-04 - 2.00E-04

9.65E-04 l.l9E-03 9.66E-Cl4 3.01E-04 - 9.45E-05

j=2

j=3

j=4

Table 6 Relatative error S%l/3; 4) k=6 k=5 k=4 k=3 k=2 k=l k=O

I

-

-

j=O

j=l

- 2.30E-03 1.51E-02

-

-

-

-

- 5.06E-05 1.22E-04 - 1.35E-04

- 3.77E-06 - 5.94E-07 2.52E-05 - 4.16E-05

- 5.82E-07 - 5.79E-07 2.14E-06 2.69E-06 - 7.22E-06

- 1.50E-07 - 1.61E-07 1.73E-07 8.34E-07 1.35E-07 - 1.55E-06

j=2

j=3

j=4

j=5

5.02E-07 - 7.65E-08 7.8OE-09 1.54E-07 2.47E-07 - 6.49E-08 - 4.43E-07 j=6

are given in Tables 5 and 6. Some entries in these tables were obtained with the help of MACSYMA developed by the MathLab Group at Massachusetts Institute of Technology and made available to the public by the Department of Energy.

I thank Garrett Birkhoff for several helpful suggestions.

R.E. Lynch / Nine-point discrete Laplacians

334

ferences [I] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functivq

National Bureau of Standards, ,4pplied Mathematics Series 55 (National Bureau of Standards, Gaithersburgh, MD, 1964). [2] G. Birkhoff and R.E. Lynch, Numerical Solution of Elliptic Problems (SIAM, Philadelphia, PA, 1984). [3] C. de Boor, K. H ti11ig and S. Riemenschneider, Fundamental solutions for multivariate difference equations, J. Amer. Math. 111 (1989) 403-415. [4] R.J. Duffin, Basic properties of discrete analytic functions, Duke Math. J. 23 (1956) 335-363. [5] R.J. Duffin and D.H. Shaffer, Asymptotic expansion of double Fourier transforms, Duke Math. J. 27 (1960) 581-596. [6] R.J. Duffin and E.P. Shelly, Difference equations of polyharmonic type, Duke Math. J. 25 (1958) 209-2:;,8. [7] H.B. Dwight, Tables of Integrals and Other Mathemutical Data (Maxmillan, New York, 4th ed., 1961).

[8] R.E. Lynch, Fundamental solutions of g-point discrete Laplacians: derivation and tables, Report CSD-7‘R-92004, Department of Computer Science, Purdue University, West Lafayette, IN (1992). [9] VV’. I 1. McCrea and F.J.W. Whipple, Random paths in two and three dimensions, Proc. Royal Sot. Edinbugh 60 (194(i) 281-298.

[lo] S.L. Sobolev, On the uniqueness of solution of difference equations of elliptic type, Akad. Nauk SSSR q:‘N.S.) Dokl. 87 (1952) 179-182 (in Russian). El11 S.L. Sobolev, On a difference equation, Akad. Nauk SSSR (N.S.) Dokl 87 (1952) 341-343 (in Russian). [12] S.L. Sobolev, Letter to Editor, Akad Nauk SSSR (N.S.) Dokl. 88 (1953) 740 (in Russian). [13] A. Stijhr, Uber einige lineare partielle Differenzengleichungen mit konstanten Koefizienten, I, II, and III, Math. Nachr. 3 (1950) 281-298. 1141 B. van der Pol, The finite difference analogy of the periodic wave equation and the potential equation, in: M. Kac, Probability and Related Topics in Physical Sciences (Interscience Publishers, New York, 1959) ApFn,ndix IV.