The numerical solution of the dirichlet problem for the Helmholtz equation

Appl. Math. Lett. Vol. 9, No. 2, pp. 85-89, 1996 Pergamon

Copyright@1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0893-9659/96 $15.00 + 0.00 S0893-9659(96)00018-3

T h e N u m e r i c a l S o l u t i o n of t h e D i r i c h l e t P r o b l e m for t h e H e l m h o l t z E q u a t i o n J. B. R. DO VAL AND M. G. ANDRADE FO.* Departamento de TelemAtica, Faculdade de Engenharia Eldtrica Universidade Estadual de Campinas - UNICAMP C.P. 6101, 13081-970 - Campinas - SP, Brazil jbosco©dt, fee. unicamp, br

(Received and accepted June 1995) Abstract--A numerical treatment for the boundary value problem involving the Helmholtz equation Au -- A2u = 0 is presented. The precision error for square grids is of order O(h 2) and a stricter evaluation for the method is obtained when compared with the precision error for the standard finite difference scheme. K e y w o r d s - - E l l i p t i c differential equations, Numerical solutions for partial differential equations, Helmholtz equation, Nonstandard difference approximations.

1. I N T R O D U C T I O N T h e finite difference m e t h o d is t h e s t a n d a r d n u m e r i c a l m e t h o d for b o u n d a r y value p r o b l e m s . I t is o b t a i n e d by t h e d i s c r e t i z a t i o n of t h e differential o p e r a t o r , a n d the t r u n c a t i o n of t e r m s of the c o r r e s p o n d e n t T a y l o r series yields t h e b o u n d for t h e precision error. T h e m e t h o d here is b a s e d on a n a p p r o x i m a t i o n for an integral on t h e circle a s s o c i a t e d w i t h t h e s o l u t i o n for t h e H e h n h o l t z e q u a t i o n in t h e interior of this p a r t i c u l a r d o m a i n , a n d it s t a n d s as an a l t e r n a t i v e to t h e t r a d i t i o n a l finite difference m e t h o d . T h e circle c o n t a i n s t h e four n e i g h b o r i n g p o i n t s to a p o i n t in t h e mesh grid, a n d a p a r t i a l covering for t h e d o m a i n defined b y t h e grid p o i n t s a n d circles c e n t e r e d at each i n t e r i o r p o i n t is i n t r o d u c e d . As a result, an i m p r o v e m e n t in t h e precision is o b t a i n e d , a n d t h e usual finite difference m e t h o d for t h e H e l m h o l t z e q u a t i o n is seen as an a p p r o x i m a t i o n of t h e p r o p o s e d m e t h o d . T h i s conclusion is s t r a i g h t f o r w a r d , since t h e modified Bessel function of zero order, which arises in t h e p r o p o s e d m e t h o d , possesses t h e first two t e r m s in its series r e p r e s e n t a t i o n i d e n t i c a l to t h e t e r m s a p p e a r i n g in t h e s t a n d a r d finite difference m e t h o d . T h e precision error o b t a i n e d is of o r d e r O(h2), which is p r o b a b l y t h e precision error limit for t h e e q u a t i o n s t u d i e d . In [1], it was shown t h a t t h e r e is no d i s c r e t i z a t i o n involving five p o i n t s in a s q u a r e grid for t h e L a p l a c e e q u a t i o n which a t t a i n s a t r u n c a t i o n error inferior to O(h4), t h u s a p r e c i s i o n error of o r d e r O(h2). One should e x p e c t a similar t y p e of error limit for t h e H e h n h o l t z equation.

2. P R O B L E M

FORMULATION

AND

BASIC DEFINITIONS

L e t ~t be an o p e n a n d l i m i t e d s u b s e t of t h e E u c l i d e a n 7~2 space, w i t h a b o u n d a r y given by 0 ~ . L e t ft = ~ t2 0 ~ d e n o t e t h e closure of ~ , a n d (x, y) a p o i n t in Ft. Let u : ft ~-* 7~ be a real function *This work was partially supported by the National Research Council-CNPQ, Brazil, under Grant 300.721/86-2 (NV) EE and FAEP - UNICAMP, Brazil, under Grant 0252/93. Typeset by AA4S-TEX 85

86

J.B.R.

DO VAL AND M. G. ANDRADE FO.

defined on ~. We denote by C°(lt) (or C°(Oit)) the class of bounded and continuous functions on it (or 0it) and by Ck(it) the class of functions that are continuous and possess continuous partial derivatives in ft of order k. Let us consider a function u(x, y) E C4 (it) N C°(~) satisfying the following standard Dirichlet problem:

z ~ ( x , y) - ~2u(z, y) = o,

u(x, y) = ~(x, y),

(x, y) e it,

Problem (E)

(x, y) e sit.

We assume that the basic conditions for the existence of a unique solution are satisfied, e.g., [2]. In order to solve problem (E) numerically for a general domain it, it is necessary to introduce a discrete grid on ~ and evaluate u approximately at each point (node). We adopt the uniform rectangular grid or the square grid according to the following definition. DEFINITION 1. Let Sh denote a uniform rectangular or square grid with step h covering a rectangular region ~ such that for some N and M,

Sh := {(xi,yj) I xi = xo + i h , yj = Yo +Jh; i = 0 , 1 , 2 , . . . , N and j = 0, 1 , 2 , . . . , M } . A node Sh is denoted by zi,j = (xi, yj), and a discrete function ¢ : Sh ~-~ ~ , called a grid function, can be associated to any function ¢ : it ~-~ 7~ such that ¢i,j = ¢(zi,j). A node zi,j Of Sh is internal if the distance from zi,j to Oft is larger than h. The total number of interior points is denoted by NoMo, with No = N - 1 and Mo = M - 1 . A

For v E 7~iv and A E T~NxN, vi and aid denote a generic element o f v and A, respectively, and we adopt the usual Euclidean norm for v and the corresponding induced norm for A [3]. For a function u E Ck(it), let Dk be an upper bound for the k th order differentials of u, namely Dk:=max

sup

Oxk_Jy j

:(x,y) Eit

:0_
.

(1)

Also, we denote by I0(.) the modified Bessel function of zero order.

3. M A I N R E S U L T THEOREM 1. Let u(x, y) e C4(ft) be the solution for the boundary value problem (E), and let Sh be the square grid given in Definition 1, and ui,j be the corresponding grid function. Let Ui,j be a grid function that satisfies the following difference equation for i = 1, 2 , . . . , No, j -- 1 , 2 , . . . , M 0 : 1 U~d - 4Io(Ah) (Ui-l,j + U~+ld + U~,j+I + Uid-1). (2)

It follows for p = 2 or p = oo that

[[U-u[]p < a(h,A) (\ h 44! D 4 + O ( h 6) ) ,

(3)

where

a(h, A) --

x/NoMo I0(Ah) + 2 sin 2 (_~h) _ 1'

for p = 2,

1 I0(Ah) - 1'

for p = co.

COROLLARY 1. Let u(x, y) 6 C4(it) be the solution for the boundary value problem (E), and Sh be the square grid according to Definition 1. If Ui,j is a grid function satisfying (2), then

Itu - ~1too = O( h2).

(4)

Helmholtz Equation

87

REMARK 1. The difference equation in Theorem 1 can be compared with the standard difference equation obtained for the Helmholtz equation when the grid Sh is adopted, e.g., [4]. Using the Taylor series, the difference equation (2) has a standard similar counterpart given by Ui,j --

1

4To(Ah)

(Ui-l,j

~- U i + l , j ~ V i , j + l -}- U i , j - 1 ) ,

where

To(Ah) = 1 +

(5)

()~h) 2

The truncation error for (5) is obtained by direct development of the Taylor series for u(x, y) in a point in Sh, and by applying the basic relation Au = A2u into the sum of the four series involving the four neighboring points. Denoting the truncation error by r{ui,j}, we obtain 7" { Ui,j } = Ui_ l, j ~- Ui + l, j -~- Ui,j + 1 -}- U i , j _ 1 --

4T0(Ah) u~,j ,

and the analysis in [4] yields the following bound, which is justified for a function u E C4(~): 2h 4 (04ui,j 04ui,j ~ 4h4D4 r{ui,j}=~ \ Ox4 + Oy4 j + O ( h 5) <_ 4----(--.+O(hS)" The precision error can now be evaluated using an upper bound suggested by [5]: 1 Ilu - uII2 ~ ][A-Xl[2 II~-{u}ll2 = m i n l < k < N o M o [Ak] [IT{u}H2' where A is the resulting matrix of the linear set of equations yielded by the application of (5) at each node of Sh, and Ak, k = 1 , 2 , . . . , NoMo denote the eigenvalues of A. Now, we bring to bear the following evaluation for the dominant eigenvalue of A [6]: min

l
IAkl = 4 (T0(Ah) - 1) + 8 sin2 ( ~ )

.

Therefore, the precision error for the Taylor scheme possesses the following bound:

IJu - vll~ _< T0(~h) + 2sin 2 ( ~ )

- 1 \

4!

+ O (h 5)

(6)

REMARK 2. The basic difference between expressions (2) and (5) is the modified Bessel function in the denominator of (2), and the polynomial T0(Ah) in the denominator of (5). Distinct evaluations for the precision error are obtained as seen in (3) and (6), with an advantage for the method presented in Theorem 1. This can be verified first by noticing that the denominator of (3) presents the infinite series I0(Ah) in the place of the polynomial T0(Ah) arising in the standard finite difference method. In addition, T0(Ah) corresponds to the first two terms of the series representation of I0(Ah), which clearly shows the consistency of the two schemes, and demonstrates that the Taylor series development represents an approximation for the method proposed here. The ratio between the two coefficients appearing in (6) and (3) yields the following expression: R = I0(Ah)+2sin 2 ( ~ - ~ ) - 1 = 1 + T0(Ah) + 2sin 2 (_~h) _ 1 To(~h) + 2sin ~ ( ~ )

1

~=~

88

J.B.R.

DO VAL AND M. G. ANDRADE F o .

r i0.

0.i

a

b

c

d

0.00

1

2

5

i0.

20.

50.

i00.

~%

Figure 1. A plot for }-~k=2 /(To(Ah) + 2sin2(-~) - 1), N = 500; (a) h = 0.1; (b) h ----0.05; (c) h--- 0.025; (d) h-- 0.01. In Figure 1, a plot is shown for the second t e r m on the right-hand side of the expression for the ratio r picturing the improvement provided by the method for different values of parameters A and grid size h. The justification behind Theorem 1 relies on the solution of the Dirichlet problem defined locally on an intertwined set of balls covering the domain. This procedure conducts to the solution of problem (E) as assured by the Alternating Method of Schwarz, e.g., [7,8], applied to the set of balls. The Dirichlet problem for the Helmholtz equation defined in the interior of a single ball possesses an analytical form; see [9]. For our purposes, the interest is in the evaluation of u at the center point z; in this situation, the formula in [9] yields 1

ui,j - 2;rI0(Ah)

2"

~(h, 8) dO.

(7)

A numerical approximation for the integral above is proposed and the numerical error is evaluated as follows. THEOREM 2. Let u be a function of class C 4 in a neighborhood o[the circle 01~h(Z ) . The following

ewluation holds: Jo

u(h, 0) dO -

-~

k=l

<

h4D4 4! + o (hs).

T h e grid function ui,j could be evaluated at a point z as in (7), provided t h a t the solution for (E) on the circle OB is known. Recognizing t h a t we can only known function u(x, y) in four regularly spaced points on the circle OB, as yielded by Sh, we ought to use these point values for solving the integral in (7) numerically. Theorem 2 provides the error evaluation for the integral in (7) using the trapezoidal method for regularly spaced points. This simple method attains great precision thanks to the fact t h a t 8 --* u(h, 8) is a periodical function [10]. The proofs of Theorems 1 and 2 are available upon request and will be published elsewhere. Preliminary numerical evaluations have been confirming, so far, the gain in precision with respect to the standard finite difference method. When the product A × h is a large number, basically larger t h a n one, the gain obtained with the method is more evident, as expected from Theorem 1 and the comparisons in Figure 1. In fact, the method has been showing an even greater advantage in precision than t h a t indicated by the theoretical evaluation presented here, but further experiments should be necessary to strengthen this conclusion.

Helmholtz Equation

89

REFERENCES 1. G. Birkhoff and S. Gulati, Optimal few-point discretizations of linear source problems, SIAM Y. Numer. Anal. 11 (4), 700-728 (1974). 2. D. Gilbrag and N.S. Trundinger, Elliptic Partial Differential Equations of Second Order, Springer-Verlag, New York, (1983). 3. G.H. Golub and C.F.V. Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, (1984). 4. L.V. Kantorovich and V.I. Krylov, Approximate Methods of Higher Analysis, Intersc. Publ., New York, (1964). 5. E. Isaacson and H.B. Keller, Analysis of Numerical Methods, John Wiley, New York, (1966). 6. P. DuChateau and D.W. Zachmann, Partial Differential Equations, McGraw-Hill, New York, (1986). 7. 1%. Courant and D. Hilbert, Methods of Mathematical Physics (Partial Differential Equations), Vol. II, John Wiley, New York, (1962). 8. T.P. Mathew, Schwarz alternating and iterative refinement methods for mixed formulations of elliptic problems, Part I and Part II: Algorithms and numerical results, Numer. Math. 65 (4), 445-492 (1993). 9. J. G6rowski, On some properties of the solution of the Dirichlet problem for the Helmholtz equation in the interior and exterior of a circle, Demonstration Math. 19 (2), 303-315 (1986). 10. P.J. Davis and P. Rabinowitz, Methods of Numerical Integration, Academic Press, New York, (1975).