Green's functions in rectangular domains with Dirichlet or Neumann boundary conditions

Journal o f Electrostatics, 13 (1982) 227--248 Elsevier Scientific Publishing Company, Amsterdam -- Printed in The Netherlands

227

G R E E N ' S F U N C T I O N S IN R E C T A N G U L A R DOMAINS WITH D I R I C H L E T OR NEUMANN BOUNDARY CONDITIONS

ANDRI~ MORAND and JACK ROBERT Laboratoire de Genie Electrique de Paris, Ecole Sup@rieure d'Electricit@, Universit@s de Paris VI et Paris X I (Laboratoire associ@ au CNRS, No. 127), Plateau du Moulon, 91190 Gif-sur- Yvette (France)

(Received September 25, 1981 ; accepted in revised form April 27, 1982)

Summary The Green functions relative to rectangular boundaries are expressed by means of a transformation which maps a circle onto a rectangle. The transformation, established in the first part of the paper, allows the determination of Green's functions relative to the potential and to the flux function for both Dirichlet and Neumann boundary conditions. The existence of multiple Green's function expressions corresponding to a given set of boundary values is studied in the last kind of conditions. The general expressions of Green's functions within infinite flat strips, as well as semi-infinite ones, are derived from the preceding results. Several applications imply the use of Green's functions by means of either an integration or a numerical process. The integration is applied to a unitary potential while the computation resorts to the numerical Green functions. In spite of the small number of iterations these functions provide satisfactory graphical results for both Dirichlet and Neumann problems. A comparison with other methods used in the determination of potentials considers both the theoretical and the numerical aspects of the process.

1. I n t r o d u c t i o n G r e e n ' s f u n c t i o n s are a p o w e r f u l t o o l in the d e t e r m i n a t i o n o f p o t e n t i a l and field lines within two- or t h r e e - d i m e n s i o n a l d o m a i n s [1 ]. In this paper, owing t o t h e use o f c o n f o r m a l mappings, the s t u d y is restricted t o plane domains. T h e first p o i n t to be dealt w i t h is a t r a n s f o r m a t i o n w h i c h maps a circular d o m a i n o n t o a rectangle. Then, u n d e r b o t h Dirichlet and N e u m a n n b o u n d a r y c o n d i t i o n s , this t r a n s f o r m a t i o n leads t o the G r e e n ' s f u n c t i o n expressions w i t h i n a n y r e c t a n g u l a r d o m a i n . A s t u d y o f infinite strips is included because these d o m a i n s m a y be c o n s i d e r e d as d e g e n e r a t e rectangles in w h i c h t h e general expressions apply. In t h e last part p o t e n t i a l f u n c t i o n s are o b t a i n e d b y m e a n s o f integrations involving t h e G r e e n f u n c t i o n s a n d b y c o m p u t a t i o n s resorting t o t h e n u m e r i c a l G r e e n f u n c t i o n s . T h e s t u d y ends with c o m p a r i s o n t o o t h e r processes o f p o t e n t i a l d e t e r m i n a t i o n . Part o f this c o m p a r i s o n constitutes A p p e n d i x 2, while A p p e n d i x 1 c o n t a i n s e x p l a n a t i o n s a b o u t the multiplicity o f G r e e n ' s f u n c t i o n s w i t h i n N e u m a n n b o u n d a r y c o n d i t i o n s . 0304-3886/82/0000--0000/$02.75

© 1982 Elsevier Scientific Publishing Company

228 2. T r a n s f o r m a t i o n (T) In a rectangular d o m a i n t h e Green's f u n c t i o n expressions are o b t a i n e d by means o f a c o n f o r m a l mapping, t h e r e a f t e r d e n o t e d as (T), which associates the rectangle to a circular domain. This m a p p i n g is the p r o d u c t o f t h r e e transf o r m a t i o n s : a bilinear t r a n s f o r m a t i o n , a S c h w a r z - - C h r i s t o f f e l t r a n s f o r m a t i o n and a L a n d e n t r a n s f o r m a t i o n . T h e first one maps the inside o f the circle in Fig. l ( a ) o n t o the u p p e r half plane in Fig. l(b}. T h e n the half plane is m a p p e d

,¢

Px--i ]

,'///~'//,

(F) >

/////

~Y

D "/////,/,

2K'(M)//C

@ ////×

~(M) >x

Fig. 1. Domains related by conformal mappings. An inversion transforms the circle (a) into the half plane (b). A Schwarz--Christoffel transformation, merged with a Landen transformation, maps (b) onto the rectangle (c). The mapping of (a) onto (c) is completed by the transformation (T). o n t o a rectangle b y the s e c o n d t r a n s f o r m a t i o n [2, 3]. T h e affixes o f a c o m m o n p o i n t P within the half plane and within the rectangle are respectively n o t e d as and Z. So, b y use o f the J a c o b i a n elliptic f u n c t i o n sn, the Schwarz--Christoffel t r a n s f o r m a t i o n is: ~* = s n ( Z I m ) = s n ( Z , k )

(1)

The p a r a m e t e r rn and the m o d u l u s k (which is the square r o o t o f m) can be evaluated f r o m t h e ratio o f the rectangle sides [4, 5]. Lastly, a L a n d e n ascending t r a n s f o r m a t i o n [6] yields simplifications in t h e f o r m u l a e [7], while it replaces t h e p a r a m e t e r m b y a new one [8] t h e r e a f t e r d e n o t e d M. T h e c o m p l e x variable Z has a r e d u c e d value such t h a t the c o o r d i n a t e s o f the rectangle vertices (Fig. 1(c)) are either equal t o zero or to the c o m p l e t e elliptic integrals o f t h e first kind, K(M) and K'(M) [9]. T h e overall t r a n s f o r m a t i o n (T) maps the circle in Fig. l ( a ) o n t o the rectangle in Fig. l ( c ) . T h e relations b e t w e e n t h e Cartesian c o o r d i n a t e s (x, y), or t h e polar ones (p, 0), o f a c o m m o n p o i n t P inside t h e circle and t h e

229 Cartesian coordinates (X, Y) of the homologous point in the rectangle are [10] x/R

= p cos0/R

= -x/-MsnX/D

y/R

= psinO/R

= dnX cnY/D

(p/R) 2 = (dnY- x/~cnX snY)/D

(2)

with D = dnY + x/-M-cnX snY R is the circle radius. The Jacobian elliptic functions sn, cn and dn with the argument X are characterized by the parameter M. The parameter turns into M1, a number complementing M to unity, when this argument is Y. In Section 3 the relations (2) are used in determining Green's functions but they actually have a wider field of applications with some examples described in Morand

[10]. 3. Green's functions In this paper Green's functions correspond to boundary conditions of either the Dirichlet or the Neumann kind. Though the mixed boundary conditions or the Poisson equation conditions are not considered below it is worth noting the present m e t h o d also applies to such cases. With Dirichlet or Neumann boundary conditions there exist two Green's functions, denoted K and G after Durand [11]. They are the two components of a complex Green function [12], hereafter noted as Wg W~(z) = g ( x , y ) + j V ( x , y )

(3)

This expression should be compared to the definition of a complex potential [131 W(z) = g(x,y) + jV(x,y)

(4)

U and V are respectively the flux function and the potential function. We can notice that the complex potentials (3) and (4) differ, by a factor j, from their definition in Durand [13] but this does not affect the determinations of G and K. It results, from (3) and (4), that the Green function G, which is involved in potential determinations in Section 4.1, is a potential function too and K is the corresponding flux function. These functions G and K are expressed within a circular domain and are then transposed inside a rectangle by means of (2). Boundary conditions of Dirichlet and Neumann kind are characterized by the respective subscripts D and N to the Green functions G, K and Wg. This notation avoids any confusion with the complete elliptic integral of the first kind K (M). 3.1. D i r i c h l e t c o n d i t i o n s

Let P~ be a point upon the boundary (F), with the affix zi in the circular

230

geometry. The Green complex potential corresponding to ordinary point P ( z ) is [14] : WgD (z ;z;) =

Pi and to any

j (z + z/)

(5)

2~R(z - z i )

The separation of WgD into its real and imaginary parts leads to expressions for the GD and KD Green's functions within the circle. Then by use of (2) we get: (1 - Mll/2) cnX snY d n Y i

G D (X, Y;Xi, Yi) =

2~(dnY d n Y i - M snX snX i - dnX cnY dnX i cnYi)

(6)

(1 - Mll/: ) (dnX cnY snX~ - snX dnX i cnY~)

KD(X,Y;Xi, Yi)=

2~(dnY dnYi - M snX snX i - dnX cnY dnX i cnYi)

(7)

X and Y are the coordinates of P, while Xi and Yi are those of Pi. In order to simplify the writing, the elliptic function parameters have been omitted. As in (2) this parameter is M when the argument is X, or X~, and it is M1 when the argument is Y or Y~. According to the location of Pi on the boundary, the coordinate Xi, or Y~, has a special value giving simple values to the corresponding functions, and so there only remains in (6) and (7) the functions with the respective argument Y~ or Xi. As an example, upon the horizontal side AB (see Fig. 1(c)) the functions cnY; and dnY~ are replaced by the unit value. So, in applications, the Green's functions GD and KD have actual expressions simpler than (6) and (7). Fig. 2 shows, within a square, a set of level curves for GD and KD, or, stated in another way, it shows equipotential and field lines pertaining to an electric dipolar straight line, normal to the figure plane, when its trace Pi KD=O -.01 -.02:03

=05

-,2

-1

-.5

-,3

Fig. 2. G r e e n ' s f u n c t i o n s G D a n d K D w i t h i n a s q u a r e .

231 is located at the origin of coordinates [15]. In this case dnXi, cnYi and dnYi equal one and snX~ is zero. A last point must be specified before we apply Green's functions to potential determinations. It concerns the evaluation of the boundary element dl in the vicinity of Pi. This element is evaluated in the circular geometry and then it is transferred into the rectangle plane by means of (2). Let ¢ be the polar angle for Pi on the circle (F) the radius of which is R. In the elementary arc expression, Rd0, the quantity dO is evaluated through the differentiation of (2), while R is eliminated by means of R

=

M1/2/(1

-

Mll/2)

(S)

This condition characterizes the preservation of central symmetries from the circular domain to the rectangle [16]. Like we did for GD and KD, we have to distinguish several cases according to the location of Pi on the boundary. On the horizontal side AB (see Fig. 1(c)) the differential element is dl = (1 + M1/2) cnXi dXi

(9)

On the vertical side BC it becomes

dl = M1~/~(1 + M~I/2) (snYJdnYi) d Y i

(10)

On the two other sides, CD and AD, we have to change the signs of (9) and (10) respectively.

3.2. Neumann conditions In Neumann's problems the complex potential W, at any ordinary point P(z), is defined as [17] :

W(z) = Cst - (j/n) f

(dV/dn) ln(z - zi) dl

(11)

F In this expression dV/dn is the normal derivative of the potential V at a point Pi on the boundary (F). In such conditions the complex Green function is:

WgN(Z;Zi) = - ( j / u ) ln(z - zi)

(12)

After the use of (2), the imaginary and the real part of this function are:

1 ( GN(X, Y;Xi, Yi) = - - In 2.

KN(X, Y;Xi, Yi) .

1

c n Y d n Y i - dnXi cnY i ( d n Y + M 1 / 2 c n X s n Y ) h . . . . . . . arc~an[ M 1/2 snXi (dnY +M 1:2 cnX snY) - M 1/2 snX dnY i )

.

.

( d n Y + M l / 2 c n X s n Y ) dnYi ) dnY d n V ; - ~ M s ~ [ d - ~ c - - ~ d n X , cnYi (13)

.

[ d~ .

(14) The above remarks concerning the elliptic function parameters and the simplifications due to special values of X i and Yi apply to (13) and (14). The

232 KN=O

-.I

~2

=3

-.4

Fig. 3. G r e e n ' s f u n c t i o n s G N and K N. The drawings are derived f r o m eqns. (13) and (14). The reversal o f the G N e q u i p o t e n t i a l curvature is e m p h a s i z e d b y t h e drawing o f the curve G N = 0.084.

curves in Fig. 3 r e p r e s e n t G N a n d K N inside a square w h e n Pi is l o c a t e d at t h e origin o f axes. In t h e case o f K N t h e angle given b y t h e a r c t a n g e n t f u n c t i o n was c h o s e n within the interval [0,~ ], a n d a c o n s t a n t p o t e n t i a l was a d d e d in every p o i n t so as t o get zero u p o n t h e Y axis. A set of e q u i p o t e n t i a l curves, c o m p u t e d f r o m t h e series e x p a n s i o n of t h e G N f u n c t i o n , is s h o w n in D u r a n d ' s treatise [ 1 8 ] . T h o u g h these curves bear no g r a d u a t i o n , it is o b v i o u s t h e y do n o t r e p r e s e n t the s a m e f u n c t i o n as in Fig. 3 b e c a u s e G N has a c o n s t a n t value along t h e u p p e r side of t h e square. T h e r e f o r e , as e x p l a i n e d in A p p e n d i x 1, d i f f e r e n t G r e e n ' s f u n c t i o n s m a y exist as s o l u t i o n s t o a given N e u m a n n p r o b l e m . T h e n , w i t h t h e n o t a t i o n s o f Fig. l ( a ) , t h e c o m p l e x G r e e n f u n c t i o n (12) can be r e p l a c e d b y [17] :

J WgN(Z;Zi) = ---

In

(2jR(z-zi)

)

(15)

T h e c o r r e s p o n d i n g e x p r e s s i o n s f o r GN a n d KN are:

GN(X, Y;Xi,Yi) KN(X, Y;Xi, Yi)

1 2n

= --

in

(

(dnY-M'/2snX)(dnYi-M'/2snXi) ) dnYdnY i - - M snX snX i - - d n X c n Y d n X i cnY i

1 { M 1/2 c n X s a Y (dnYi Z M1/2 snXi) = -- arctan 7r ~ d n X c n Y (dnYi - M 1/2 snXi) - dnXi cnYi ( d n Y -

(16) ]

M 1/2snX)/ (17)

The integration of (13) and (16) along the boundary must give two potential functions only differing from each other by a constant value [17 ]. The same remark applies to (14) and (17). These statements will be verified in Section 4.2.

233

3.3. Infinite strips An infinite flat strip (or a half one) can be derived from a rectangular domain. Thus, the general expressions of Green's functions within these strips are readily obtained from the results of Sections 3.1 and 3.2. They will be compared to the expressions which were known before in special cases. An infinite flat strip is obtained when the parameter M, which characterizes the elliptic functions, is set equal to one. Then K'(M) equals n/2, while K(M) is infinite, and so the rectangle is transformed into an infinite horizontal strip the width of which is ~ (Fig, 4). It has been stated above that any ~Y

/////////////,

"if, / / / / / / / / / / / / / P J,

Pi

(F) /////////////~

0 / / " ^'~i " / / / / / / / / /

>×

Fig. 4. Infinite flat strip. elliptic function occurring in the Green's function expression is dependent upon the parameter M when its argument is X or X~, and upon the complementary parameter M, when the argument is Y or Yi. Inside the domains of Fig. 4 these functions respectively degenerate into hyperbolic and circular functions. Therefore, the V D function, as given by (6), becomes:

GD(X,Y;Xi,Yi) .

1 ( 2u

.

.

.

.

sinYcoshXi

.

)

coshX coshX~ - sinhX sinhXi - cosY cosYi

(18)

The cosY i term has the respective value + 1 or - 1 according to the location of P~ on the lower or the upper side of the strip. Furthermore, we must, as above, replace the differential element dl b y dXi, i.e. in the present case we divide (18) b y coshX i, in order to enable a comparison with a known expression for G D. The same remarks hold for the K D function as deduced from (7) 1

(

cos___yYsinhX___A-- sin__~_~c_osYA

KD(X,Y;Xi,Yi) = 2--~- coshX coshX i - sinhX sinhXi - cosY cosY/

)

(19)

When Pi(Xi,Yi) is at the origin of coordinates we find the known expressions [19] of these functions again, without any change for G D and with an additional constant for KD. As for the GN and K N Green functions we must use (16) and (17) because the comparison rests on expressions issued, by conformal mappings, from the half plane of Fig. 9 (b) when the point Po is rejected to infinity (see Appendix

234

1). The previous t r a n s f o r m a t i o n s o f elliptic f u n c t i o n s into circular and hyperbolic ones lead to

GN(X,Y;Xi,Yi

) - - - - In [ e x p ( X + 2~ -

Xi)

( c o s h X coshXi - sinhX sinhX i -

- c o s Y cosY/)]

K N(X, Y; Xi, Yi)

1

= ~ arctan 7r

(

(20)

sinYe Xi

- - -- - c o s Y e xi - c o s Y / e - x

)

(21 )

The p r e c e d i n g r e m a r k s a b o u t the c o s Y i t e r m and the differential e l e m e n t are still valid. When Pi is at t h e origin o f c o o r d i n a t e s , the G N and K s f u n c t i o n s are after D u r a n d [20] and with t h e p r o p e r n o t a t i o n s -

GN(X,Y;O,O)

- -~

KN(X,Y; 0,0)

= --

-

2~

{In2 + In[eX(coshX- cosY)] }

1 arctan (

e -X

sinY ----

cosY

'

)

(22)

(23)

When t h e a r c t a n g e n t f u n c t i o n is d e t e r m i n e d w i t h i n the interval [0, ~ ] we can verify that, e x c e p t f o r an a d d i t i o n a l c o n s t a n t , (20) a n d (21) are t h e same as t h e k n o w n expressions (22) a n d (23).

3.4. Semi-infinite strips When the p a r a m e t e r M is zero t h e n u m b e r s K(M) and K'(M) are respectively equal to ~/2 and to infinity, a n d the rectangle degenerates into a semiinfinite strip (Fig. 5 (a)). U n d e r such c o n d i t i o n s , G r e e n ' s f u n c t i o n s involve the ratio o f t w o null n u m b e r s a n d are, t h e r e f o r e , u n d e t e r m i n e d . This situation arises f r o m t h e fact t h a t the L a n d e n t r a n s f o r m a t i o n , w h i c h leads t o (2), is no longer valid w h e n M is zero. H o w e v e r , the calculations m a y be achieved b y m e a n s o f a s y m p t o t i c expansions. As a first step, we settle the p a r a m e t e r M Y

~'//////////// (F}

p X

(F)

Pi

Pi =~12 ~//' 3 ~" ~12

x

o 7/~////////

® Fig. 5. S e m i - i n f i n i t e flat

® strip.

~×

235 to a value close to zero. Then we refer to formulas relating elliptic functions to circular functions, such as [21,22] sn(XlM) = s i n X - '

M cosX (X - sinX cosX)

(24)

With M in the vicinity of zero, the complementary parameter, idl, is close to unity and the corresponding elliptic functions are replaced by hyperbolic ones, for example [23,24] sn(Y]M,) = tanhY + (M/4 cosh2Y) (sinhY coshY- Y)

(25)

These approximations lead to expressions of Green's functions which include a principal part not dependent on M and other terms tending to zero with M. When applied to (6) this process leads to

GD(X,Y, Xi,Yi) = _~(

cosX sinhY )(26) sin2X + sin2Xi + sinh2Y + sinh2Yi - 2 sinX sinXi coshY coshYi

This function has a known expression in the case of Fig. 5 (b). In order to complete the comparison (26) is subjected to the relevant change of axes 1 sinhX sinY GD(X,Y;Xi,Yi) = - 2n (coshX c o s Y - coshX~ cosYi) 2 + sinh2X sin2Y

(27)

This last expression covers any location of P~(Xi,Yi) but simplifications occur whether this point is situated upon the horizontal sides or upon the vertical one. The same process is applied to (9) and (10) and it leads to 2 sinYi dYi

dl =

1+ 2 sinhX/dXi

for Xi = 0

(28)

for Yi = 0 or

The minus sign is relative to the upper side of the strip. When P; lies on the vertical side AB (Fig. 5 (b)), the GD function has the following known expression [15] :

GD(X,Y;O, Yi) =

sinhXsinYsinYi

@(

)

(coshX cosY - cosYi) 2 + sinh2X sin2y

(29)

The insertion of (28) into (27) leads to the same expression. When P~ is located upon the horizontal sides the GD function can be related to the equivalent function expressed within an infinite strip [25]. This statement has been verified by use of (18). The same process applies to KD Green's function

KD(X,Y;Xi,Yi)

_1_1( 2n

---~ coshX - - _cosY ~ - - ------~, coshX, cosY, (coshX c o s l - c o s , ~ i

cosY O -

+ sinh2Z sin2Y

)(30)

236 Like in (27), simplifications occur in the functions which depend on Xi or Yi according to the location of Pi on the boundary. In the special case with P~ on the vertical side AB (Fig. 5 {b)) we have expressed K D by means of the relations quoted in Durand [20]. This process gives the same result as the combination of (28) and (30) does. ~ h e n Pi is located upon the horizontal sides we have verified that the expression of K D could be found again by use of (19). As regards Neumann's conditions, the two couples of relations (13) and (14) or (16) and (17) lead to the same Green functions because, when the rectangle turns into a semi-infinite strip, the second electric line (see Appendix 1) intervening in the determination of (16) and (17) is rejected to infinity, as it was already the case for (13) and (14). So, with the notations of Fig. 5 (b), we have 1 GN(X,Y;Xi,Yi) ln[(coshX c o s Y - coshXi c o s Y i ) 2 + sinh2X sin2Y] (31) 2n 1 ( sinY coshX ) KN (X,Y;Xi,Yi) = -- arctan . . . . r coshXi cosY/ - coshX cosY

(32)

Like in (27) and (30), any location of Pi upon the strip boundary is included in these expressions. In the special case where Pi is on the vertical side AB (31) and (32) give two already known expressions again [26]. 4. Determinations of potentials One of the most important applications of Green's functions we are dealing with, is the determination of the potential V, as well as the flux function U, inside a rectangular domain with a given set of boundary conditions. The first application is based upon the usual process of integration while the two others appeal to the concept of numerical Green's functions.

4.1. Unitary potential We call the potential function obtained inside a domain " u n i t a r y potential" when over the whole boundary the potential has no other reduced values then zero and unity. Basically, with Dirichlet boundary conditions, the potential function V at any c o m m o n point P(X,Y) is related to the numerical values of potential along the boundary by [27]: V(X,Y) = f

V(Xi,Yi) GD(X,Y;Xi,Yi) dl

(33)

F

Inside a rectangle the G D Green function is given by (6) while the differential element dl corresponds to either (9) or (10). V~ith such expressions for G D and dl, the integration should likely be completed but with particular distributions V(Xi, Yi). This is just the case with a unitary potential distribution (Fig. 6).

237

Y 2 K'(M}

D I

C

V=O V=O

A } V=l -K{M) X1

I

X2

,

+K(M) X

Fig. 6. Unitary potential distribution.

The reduced potential value equals unity between P, and P2 and it is zero anywhere else on the boundary (F). Choosing both points P, and P2 on the horizontal side AB yields simplifications but the m e t h o d is valid through any location of these points on (F). V/ith Y/as zero in (6) and by use of the differential element (9), the potential (33) becomes: 1

x~

V(X,Y) = - - f n Jx,

M cnX snY cnX i dX~ 2 ( d n Y - M snX snXi - - dnX cnY dnXi)

(34)

In order to complete this integration we consider the following function:

F(X,Y, Xi) = a r c t a n [ d n Y

McnXsnYf(Xi__._))

]

- d n X a n y - M snX f ( X , )

(35)

with snXi f(x~)

=

1 + dnXi

~'e can verify the differential of (35) equals the expression under the sum sign in (34). Thus we have: 1 x~ 1 V(X,Y) =--~ fx dF =--n [F(X,Y, X2) - F(X,Y,X,)] 1

1

V(X,Y) = -7T

arctan

M cnX snY [f{X2) - f(X,)] dn Y - dnX cnY - M snX [f(X2) +

f(X,)] + (dnY + dnX cnY) f(X2) f(Xl) (36)

So, with somewhat different notations, we find again a result we had established by another process [28]. In the same way, the replacement of G n by the K D Green function in (34) leads to the flux function U(X,Y). Like for V, this result is the same as in Morand [28]. It can be noted that this reference exhibits general expres-

238 sions, for U and V, applicable to any location of P, or P2 on the rectangle sides. Through the permutation of U and V the previous unitary potential distribution is transformed into a set of two electric lines, of inverse polarities, with their traces in P~ and P2 (Fig. 6). Now the rectangle sides coincide with field lines. These boundary conditions are the same as those obtained in Section 3.2 for Green's functions relative to Neumann's conditions. Although (17) and (36) exhibit analogies they are not quite identical because P~ and P2 have not the same location in both cases. On the contrary, starting with the potential V as given by Morand [28] we can exactly find again the K N Green function expressed by (17). Likewise, the flux function from the last reference gives the expression (16) for the G N Green function again. So, we can immediately generalize the G N and K s expression to every location of Pi and P0 (see Fig. 9 in Appendix 1) through the selection of the corresponding expressions for U and V in l~orand [28].

4.2. Numerical Green's functions The above analytical determination of potential is restricted to those cases in which the integration (33) is feasible. On the contrary, by means of c o m p u t a tions this determination can be made whatever the distribution V(Xi, Yi) may be. Among the set of numerical processes we shall consider the use of numerical Green's functions. These functions have been introduced in order to replace integrations using Green's functions, such as (33), by sums resting on a finite number of terms [29]. Thus, under Dirichlet's conditions, the potential V is given by

V(X,Y) = ~ V(Xi,Yi) ~D(X,Y;Xi,Yi) i

(37)

The numerical Green function, ~ D, is connected to the G D Green function by [30]:

~D(X,Y;Xi,Yi) = lim [h GD(X,Y;Xi,Yi)]

(38)

h -90

In this expression h is a boundary element in the neighbourhood of Pi. Vve notice from (38) that ~ D has a finite value everywhere in the domain, even in Pi where this function tends to unity. A determination of potential by use of (37) is worked out, as an example, within a square whose boundary distribution V(Xi,Yi) is a cosine function over half the horizontal side AB (Fig. 7(b)) and is zero elsewhere. The length h, in (38), is chosen equal to K/8. We selected this value in order to show that quite satisfactory results could be obtained with very few iterations in (37). These iterations are completed at the points with the abscissae X i = n(K/8)

(39)

239 (e)

U= 0 Ix

-.05

V= 0

/J

]

\\

I

//

..:o.£~.Ji___i I / //#x ,i 7"-

\\\

(b) 2

Ill

t##'5"4

-.3

A

B

-K

// X tV 1

-K/2

0

+ K/2

B +K

X

Fig. 7. (a) Equipotential and field lines computed by means of the numerical Green functions for Dirichlet boundary conditions. The arrows signal the points used in the computations. (b) Potential distribution along side AB. T h e integer n is c o m p r i s e d b e t w e e n - 4 and +4 inclusively. T h e p o t e n t i a l distribution is t h e r e f o r e V(Xi,O) = cos(n X J K )

= c o s ( n n/8)

(40)

As regards t h e G D f u n c t i o n , in (38), we m u s t m u l t i p l y (6) b y t h e ratio d l i d X i , as derived f r o m (9), because, whereas the differential e l e m e n t dl is relative t o G r e e n ' s f u n c t i o n s in S e c t i o n 3.2, it is p r e s e n t l y the differential e l e m e n t dX~ which c o n c e r n s t h e sum (38), as well as the integration (33). Thus, we have

M.

V(X, Y) = 1 6 ,

n = --4

t\ dn Y - M snX snXi - d n X

c n Y dnXi

)

,:4,)

Taking into a c c o u n t t h e simplifications d u e t o the s y m m e t r y and t o the zeros which o c c u r w h e n n equals -+4, we get V(X,Y) = 3 n=O*

MK 8.

X

c o s ( n , / 8 ) cnX snY c n ( n K / 8 ) [ d n X c n Y d n ( n K / 8 ) - d n Y ]

(42)

[ d n Y - d n X c n Y d n ( n K / 8 ) ] 2 - M 2 sn2X sn 2 (n K / 8 )

T h e star indicates the c o r r e s p o n d i n g t e r m m u s t be divided b y two. Likewise, in the case o f the flux f u n c t i o n U ( X , Y ) , the process applies to the K D f u n c t i o n

240

(7) and it gives

MK U(X,Y) = 81r

×

cos(n 7r/8) snX cra¥ c n ( n K / 8 ) [dnX c n Y - d n Y d n ( n K / 8 ) ]

3

[dnY - dnX cnY dn(n K/8)] 2 - MZsnZX sn2(n K/8)

n=0*

(43),

A set of level curves for U and V are drawn in Fig. 7 (a). They exhibit distortions in the vicinity of the points (39) at which the G D and KD functions have infinite values. On the contrary, and in spite of the large intervals between these points, in the remainder of the domain the curves have rather regular shapes which may be enough for applications in which a high precision is not required. Obviously this precision is improved when the length h, in (38), is shortened, i.e. when the number of iterations in (42) and (43) is increased. An expression similar to (38) defines a numerical Green function, N (X,Y;Xi,Yi), for Neumann boundary conditions. The potential is obtained by means of an expression similar to (37)

V(X,Y) = ~ i

d V ( X i ' Y i ) ~ N(X,Y;Xi,Yi) dn

(44)

As an example, within a square, this expression is applied to the determination of the potential whose normal derivative is -+1 upon the horizontal side AB (Fig. 8(b)) and is zero anywhere else on the boundary. The space between two

(a)

V=O

.05 ]

i U=O

I

/

/

I

\\ .1 I

""

"\ .2 dV/dY .3

(b)

+1 ,5 -K(M) A

.7B

AI

0

IB

+K(M)

-1

Fig. 8. (a) E q u i p o t e n t i a l a n d field lines c o m p u t e d b y m e a n s of t h e n u m e r i c a l G r e e n f u n c t i o n s for N e u m a n n b o u n d a r y c o n d i t i o n s . T h e a r r o w s s h o w t h e p o i n t s used in t h e c o m p u t a t i o n s . (b) P l o t of t h e n o r m a l derivative o f p o t e n t i a l a l o n g side AB.

X

241

successive points Pi is chosen equal to K/4. Equation (13) is merged with (9) to take into account the proper differential element. As the normal derivative of the potential is zero at the middle point of AB and as the numerical Green function is zero in A and B, we only have six effective points to consider. Arrows signal these points in Fig. 8. In such conditions the potential V(X,Y) is

V(X,Y) -

K(1 + MI(2)

4n

3 I n Kcn

,=1

× (MsnXsn(nK/4)

arctanh dnY - dnX cnY dn(n

4

)1 K/4)

(45)

As regards the flux function a similar process, starting with (14), gives

U(X,Y) = ~tcn nK

~,=1! ~

K(1 + MI'~)

4~

X

l McnXsnYsn(nK/4)

[arctan\dnYd~nK~ --~-nY

.)+arctantM~2sn(nK/4) 1 \ dn(nK/4) ) ] (46)

Equipotential curves, for V and U, are represented in Fig. 8 (a). Although there are very few iterations in (45) and (46), these curves exhib{t a regular shape in the major part of the square. A comparison, merely qualitative and restricted to V, may be established with the curves represented in Durand's treatise [31]. In this last case the potential was computed from its series expansion as given by the m e t h o d of separation of variables. Excepting a small area at the lower part of the square, we can verify the two sets of curves are quite similar. As the various expressions of G N , o r K s , give the same function V, or U, save for an additive constant, the same property must characterize the numerical Green functions. In order to assess this statement, the process was resumed with G s and K N given by (16) and {17). The resulting expressions are

V(X,Y) 3

Elan --

.=1{

K(1 + M1/2)

4~

X

sdnXsn'n '4' cnY dn(n g/4) )

arctanh[Ml/2sn(nK/4) ]]

f

(47)

242 K(1 + M~/2) u(x,Y)

n~=l[

-

47r

x

n K arctan I.d M cnX snYsn(nK/4)

cn ~

n Y d n ( n K / 4 ) - dnX cnY

)]

(48)

J

We readily verify these expressions and the previous ones, (45) and (46), only differ from each other through a constant value.

4.3. Comparison to other methods of determination of potential In order to justify the practical applications of the results in Section 3 we must a t t e m p t to characterize their overall accuracy, their handling easiness in computations and their possible advantages over the usual processes of determination of potential. This third point is obviously the most important one but it is the last to be considered as it depends upon the two others. First, according to their usefulness in the discussion, some details about the elliptic functions computations are brought back to mind. The simplest way to do such computations seems to use numerical tables [32] but it is quite impractical, and time consuming, when it comes to repetitive computations. A second process resorts to the series expansions of the elliptic functions [ 33]. Any ordinary programmable pocket calculator can be used for this purpose. For instance, by means of a HP 67 calculator, we have got five digit exact values in times ranging up to tens of seconds (the longer times occur when the elliptic function parameter is close to one). Five digit results were chosen in order to enable comparisons with the numerical tables, but such an accuracy is not necessary when the results are used to characterize practical devices. In such cases, the comparison between the device and its theoretical model may exhibit voltage errors readily reaching the per cent range [34]. So, we can shorten the number of iterations and, therefore, the computing time without any significant loss of precision upon the overall results. However, the most important consequence of these practical limitations is that a third process of elliptic function computations can be used. It combines the approximative expressions quoted at the beginning of Section 3.4 with successive applications of Landen's transformations, of either the ascending or the descending kind [35]. Although the result accuracy rises with the number of these successive transformations, usually one or two such transformations are sufficient for rectangles with side ratios ranging from 1/1 to about 1/4. With long-shaped rectangles these Landen's transformations become useless and the approximative formulas apply readily. So, within any kind of rectangle, the computations just imply the use of circular or hyperbolic functions. They could be performed by means of everybody's pocket calculator, but programmable computing aids avoid excessive waste of time. For example, the third computing process was applied to several determinations of the sets of curves included in this paper.

243 Each original set of curves was drawn inside a square the side of which was 200 m m (8 in) long. The voltage and field lines values were computed using 1 mm increments completed by linear interpolations, and the curves were plotted point by point. When using a HP 41 programmable calculator we got a new point within a delay matching the elapse of time required to plot the last computed point. So, even with very cheap equipment, the results in this paper are easy to use. As for the accuracy it is only a matter of computing duration which is closely related to the computing capabilities. The comparison with the usual methods giving the potential values within Dirichlet or Neumann boundary conditions is not straightforward for many factors such as programming easiness, computing m e m o r y capabilities, time consumption, cost . . . . , are involved. Besides, even the practical usefulness of Green's functions is contested [36]. Thus, in order to clarify the discussion we shall halve the comparison. First, the above Green's function expressions are compared to the previously known equivalent results. {As these results are exposed at length in Durand's treatise many references to this author are given.) The comparison has already been outlined in several places throughout the paper and it also appears in Appendix 2. There, the comparison, restricted to one kind of Green's function, starts with the most general expressions and, by successive simplifications, leads to the simplest ones. At each step of the process it is clear the already known results are one order of complexity higher than the above ones. In the comparison with the usual numerical methods, the analytical character of the above results can be an important advantage (see for ex; mple Section 4.1). In our opinion this advantage holds for numerical apphcations because, when using Green's functions the potential is known everywhere inside the domain with a special potential distribution over the boundary. Then, the user has but to sum up the contributions from each point on the boundary. This task resorts either to the numerical Green functions, as in Section 4.2, or to the classical computing processes. In these cases, the absence of relation between any point inside the domain and its neighbours may be important for parallel organizations of computations. The possibility to do computations for points chosen at will is also of great interest. On the contrary the usual numerical methods, such as finite difference elements, involve at the same time a large number of regularly spaced points inside the domain. Problems may arise about convergence towards steady values, the accuracy of results and the increase of point density in any useful area. The very great number of m e m o r y elements used in such computations implies the rejection of less powerful computing tools, like those we used. As for the determination of the field lines associated to a just-computed set of equipotential curves we only have, when using Green's functions, to replace the G function by the equivalent K function. The other methods seem to lack such simplicity. So, in potential and field lines determinations within rectangular boundaries, we think the above Green's function expressions represent a powerful tool,

244 the c o m p l e x i t y o f which is an o r d e r o f m a g n i t u d e below t h a t of the already k n o w n m e t h o d s b o t h in t h e o r e t i c a l and numerical respects. 5. C o n c l u s i o n The main p o i n t in this w o r k is the d e t e r m i n a t i o n o f the general expressions of G r e e n ' s f u n c t i o n s within r e c t a n g u l a r domains. As regards the N e u m a n n b o u n d a r y c o n d i t i o n s , the existence o f an infinity o f G r e e n ' s f u n c t i o n s giving the same p o t e n t i a l f u n c t i o n has been emphasized. The general expressions of G r e e n ' s f u n c t i o n s within infinite strips p r o c e e d from the rectangle d e g e n e r a t i o n and k n o w n results are o b t a i n e d as special cases o f these expressions. O w i n g t o its practical i m p o r t a n c e the d e t e r m i n a t i o n o f b o t h p o t e n t i a l and flux f u n c t i o n is the subject o f three applications. In the first one the potential and the flux f u n c t i o n c o r r e s p o n d i n g to a u n i t a r y d i s t r i b u t i o n arc o b t a i n e d b y means o f an integration. The t w o o t h e r d e t e r m i n a t i o n s resort to the numerical Green f u n c t i o n s u n d e r Dirichlet or N e u m a n n b o u n d a r y conditions. S a t i s f a c t o r y graphical results are o b t a i n e d f r o m very few iterative steps. The c o m p a r i s o n to the usual p o t e n t i a l d e t e r m i n a t i o n processes enhances t h e interest o f t h e evinced results b o t h in t h e o r e t i c a l and numerical respects. Lastly, it m u s t be n o t e d the f u n c t i o n s we have c o n s i d e r e d t h r o u g h o u t this p a p e r pertain to electrostatics but, as G r e e n ' s f u n c t i o n s a p p l y to every q u a n t i t y satisfying Laplace's e q u a t i o n , the a p p l i c a t i o n range o f the above results widens to o t h e r i m p o r t a n t d o m a i n s o f physics such as t h e r m o statics or h y d r o d y n a m i c s . References 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

E. Durand, Electrostatique, Vol. 3, Masson, Paris, 1966, pp. 115--203. E. Durand, Electrostatique, Vol. 2, Masson, Paris, 1966, pp. 309--11. E. Weber, Electromagnetic Theory, Dover, New York, 1965, pp. 329--33. E. Weber, Electromagnetic Theory, Dover, New York, 1965, p. 355. E.T. Whittaker and G.N. Watson, Modern Analysis, Cambridge University Press, London, 1963, p. 479. E.T. Whittaker and G.N. Watson, Modern Analysis, Cambridge University Press, London, 1963, p. 507. A. Morand, Ph.D. Thesis, Universit~ Paris XI - Orsay, 1980, pp. 32--4. L.M. Milne-Thomson, Jacobian Elliptic Function Tables, Dover, New York, 1950, p. 19. I.S. Gradshteyn and I.M. Ryzhik, Table of Integrals, Series and Products, Academic Press, New York, 1965, p. 905. A. Morand, Ph.D. Thesis, Universit~ Paris XI - Orsay, 1980, pp. 287--94. E. Dur~ ld, Electrostatique, Vol. 2, Masson, Paris, 1966, p. 120. E. Durand, Electrostatique, Vol. 3, Masson, Paris, 1966, p. 118. E. Durand, Electrostatique, Vol. 2, Masson, Paris, 1966, p. 238. E. Durand, Electrostatique, Vol. 3, Masson, Paris, 1966, p. 141. E. Durand, Electrostatique, Vol. 3, Masson, Paris, 1966, p. 116. A. Morand, Ph.D. Thesis, Universit~ Paris XI - Orsay, 1980, p. 69. E. Durand, Electrostatique, Vol. 3, Masson, Paris, 1966, p. 143. E. Durand, Electrostatique, Vol. 3, Masson, Paris, 1966, p. 135.

245 19 E. Durand, Electrostatique, Vol. 3, Masson, Paris, 1966, p. 131. 20 E. Durand, Electrostatique, Vol. 3, Masson, Paris, 1966, p. 132. 21 M. Abramowitz and I.S. Stegun, Handbook of Mathematical Functions, Dover, New York, 1970, p. 573. 22 L.M. Milne-Thomson, Jacobian Elliptic Function Tables, Dover, New York, 1950, p. 20. 23 M. Abramowitz and I.S. Stegun, Handbook of Mathematical Functions, Dover, New York, 1970, p. 574. 24 L.M. Milne-Thomson, Jacobian Elliptic Function Tables, Dover, New York, 1950, p. 21. 25 E. Durand, Electrostatique, Vol. 2, Masson, Paris, 1966, p. 116. 26 E. Durand, Electrostatique, Vol. 3, Masson, Paris, 1966, p. 133. 27 E. Durand, Electrostatique, Vol. 2, Masson, Paris, 1966, p. 110. 28 A. Morand, Ph.D. Thesis, Universit~ Paris XI - Orsay, 1980, p. 35. 29 E. Durand, Electrostatique, Vol. 3, Masson, Paris, 1966, p. 167. 30 E. Durand, Electrostatique, Vol. 3, Masson, Paris, 1966, p. 168. 31 E. Durand, Electrostatique, Vol. 2, Masson, Paris, 1966, p. 95. 32 L.M. Milne-Thomson, Jacobian Elliptic Function Tables, Dover, New York, 1950, pp. 41--110. 33 E.T. Whittaker and G.N. Watson, Modern Analysis, Cambridge University Press, London, 1963, pp. 508--12. 34 J.C. Gibbings, ASME J. Basic Engng., 88 (1966) 139--42. 35 H. Hancock, Theory of Elliptic Functions, Dover, New York, ]958, pp. 250--1. 36 E. Weber, Electromagnetic Theory, Dover, New York, 1965, p. 526. 37 E. Durand, Electrostatique, Vol. 3, Masson, Paris, 1966, p. 120. 38 J.B. Diaz and R.B. Ram, Appl. Anal., 7 (1978) 255--8. 39 E. Durand, Electrostatique, Vol. 3, Masson, Paris, 1966, p. 178. 40 E.T. Whittaker and G.N. Watson, Modern Analysis, Cambridge University Press, London, 1963, p. 524.

Appendix 1 Multiplicity of Green function expressions with Neumann boundary conditions In S e c t i o n 3.2 it has b e e n p o i n t e d o u t t h a t t w o d i f f e r e n t sets o f curves, and, therefore, two different Green's functions, correspond to the same Neumann p r o b l e m . When we r e f e r t o t h e physical i n t e r p r e t a t i o n o f G r e e n ' s f u n c t i o n s these t w o results can be c o n s i d e r e d as t w o special cases f r o m a general expression. In t h a t r e s p e c t t h e c o m p l e x G r e e n f u n c t i o n s (12}, in t h e circular d o m a i n , c o r r e s p o n d s t o an infinite electric line n o r m a l to t h e figure plane, w i t h its t r a c e in Pi o n t h e b o u n d a r y [ 1 7 ] . I t c a n easily be verified t h a t t h e curves o f Fig. 3 are, w i t h i n t h e circular g e o m e t r y , t r a n s f o r m e d i n t o a set o f c o n c e n t r i c circles c e n t r e d in Pi f o r GN and into t h e f a m i l y o f o r t h o g o n a l straight lines f o r K N . Likewise, (15) is issued, b y m e a n s o f a bilinear t r a n s f o r m a t i o n [ 1 4 ] , f r o m an initial d o m a i n w h i c h is a h a l f p l a n e w i t h an infinite electric line in Pi [ 3 7 ] , This c o n c e p t o f a u n i q u e infinite line is n o t realistic a n d it seems b e t t e r t o c o n s i d e r t w o parallel straight lines w i t h o p p o s e d charge densities. A c t u a l l y such a c o u p l e o f lines exists in t h e a b o v e h a l f p l a n e g e o m e t r y since a s e c o n d

247

the figure quoted after Durand [18] corresponds, when completed by symmetry, to the case in which the points Pi and PO are placed at the middles of two opposite sides. We can easily verify on (16) that the location of Pi at the point (-K, +K’) sets an equipotential curve of GN in coincidence with the Y axis. Appendix

2

Sums of series Commonly, Green’s functions are expressed by means of series expansions, and useful comparisons can be made with the expressions demonstrated in this paper. As an example the Gn Green function for Dirichlet’s conditions is [26] : G,(X,Y,Xi,O)

=

$

X

m C(

Sinh[(2~+l)(n/2K)(2K’-Y)]~0~[(2~+l)(nX/2K)]~0~[(2~+l)(nXi/2K)] sinh[(2p

p=o

+ l)(nK’/K)]

sinh[p (n/K) (2K’ - Y)] sin[p (n/K)X]sin[p sinh[2p(nK’/K)]

(n/K)X,‘] t

(A.1)

Here, the notations of Durand [ 261 have been modified in order to apply to Fig. 10. The location of Pi is restricted to the horizontal side AB. Due to the uniqueness of the solution of Laplace’s equation [ 381 the relation (A.l) must be identical to (6), for they are two different mathematical expressions of the same potential. Actually the differential element relative to G, in (A.l) is equal to dXi, while the element applying to (6) is the quantity dl as given

Fig. 10. Rectangular system.

247

the figure quoted after Durand [18] corresponds, when completed by symmetry, to the case in which the points Pi and PO are placed at the middles of two opposite sides. We can easily verify on (16) that the location of Pi at the point (-K, +K’) sets an equipotential curve of GN in coincidence with the Y axis. Appendix

2

Sums of series Commonly, Green’s functions are expressed by means of series expansions, and useful comparisons can be made with the expressions demonstrated in this paper. As an example the Gn Green function for Dirichlet’s conditions is [26] : G,(X,Y,Xi,O)

=

$

X

m C(

Sinh[(2~+l)(n/2K)(2K’-Y)]~0~[(2~+l)(nX/2K)]~0~[(2~+l)(nXi/2K)] sinh[(2p

p=o

+ l)(nK’/K)]

sinh[p (n/K) (2K’ - Y)] sin[p (n/K)X]sin[p sinh[2p(nK’/K)]

(n/K)X,‘] t

(A.1)

Here, the notations of Durand [ 261 have been modified in order to apply to Fig. 10. The location of Pi is restricted to the horizontal side AB. Due to the uniqueness of the solution of Laplace’s equation [ 381 the relation (A.l) must be identical to (6), for they are two different mathematical expressions of the same potential. Actually the differential element relative to G, in (A.l) is equal to dXi, while the element applying to (6) is the quantity dl as given

Fig. 10. Rectangular system.

248

by (9). By means of this last equation, and taking into account the null value for the ordinate .of P/, (6) becomes V D (X,

M ( cnXsnYcnX i ) Y ;Xi, O ) = -2~ d n Y - M snX snXi - d n X c n Y d n X i

(A.2)

This expression is the sum of the series (A.1). ~,hen special cases are considered several sums of series can be expressed. So, when Pi coincides with the origin the last sum in (A.1) disappears and the equality between (A.2) and (A.1) leads to

1 K

,_:, sinh[(2p + 1 ) ( ~ / 2 K ) ( 2 K ' - Y)]cos[(2p + 1)OrX/2K)] sinh[(2p + 1)(~K'/K)]

p=O

M cnX shY

(A.3)

27r(dnY - dnX cnY) At the centre of the rectangle (X = 0, Y = K') eqn. (A.3) becomes

1

K M 1/2 :

(A.4)

p=0 cosh[(2p + 1)(~K'/2K)] A quick numerical verification of (A.4) can be made when the rectangle is a square. Then the parameter M is 0.5 and the elliptic integrals K(M) and K'(M) equal each other. In Durand's treatise there is a table [39] giving the numerical value of the left-hand part in (A.4) up to twelve figures. Taking for K(0.5) the value [40] K(0.5) = 1.85 407 468

(A.5)

the calculation of the right-hand term in (A.4), by means of a pocket calculator, allows the eight first figures to be found again. In a similar fashion, the comparison between the above formulas of K D , G N o r K N and the corresponding series expansions can lead to various new sums of series.