On a numerical algorithm for computing fields in laminated media

0041-5553/86 $lO.GG+G.OG U.S.S.R. Comput.Maaths.Math.Phys.,Vo1.26,No.4,pp.178-183,1986 Printed in Great Britain G 1988 pergamon .JournalsLtd.

ON A NUMERICAL ALGORITHM FOR COMPUTING FIELDS IN ~MINATED MEDIA*

S.I. SMAGIN

The Fourier-Bessel integrals that occur when solving field distribution problems in plane-parallel laminated media are considered. An algorithm is given for the approximate evaluation of these integrals; it is based on deformation of the contour of integration in the complex plane.

1. An important class of problems of mathematical physics consists of problems of computing the fundamental solutions of second-order elliptic equations whose coefficients depend on one coordinate. These are encountered e.g., in the theory of soundings of laminated media by different types of.physical fields. Moreover, without a reasonably accurate and rapid method of solving such problems, it is impossible to study the fields in laminated media with inclusions by the method of integral equations. The general approach to their solution is as follows. The initial problem is reduced by a Fourier or Bessel transformation to the solution of the ordinary differential equation

Knowing ~(2. z’, I.),

we find by an inverse transformation the required fundamental solution G(M,M')=

(2)

6 is the where JO is the Bessel function of the 1st kind of zero order, r-[(2-rr)'+(~--y')*J~, Dirac delta function, k'(s) is a piecewise continuous, and p(z) a twice piecewise differentiable, complex function with discontinuities at the points z-z1,l=i,2, . . ..n-i. k’(z)-x(z)-ti0(2),

x, 0, q, t;ao,

P(+m+ed,

Ik’l’O,

b-t+t%

(3)

lPI~IP*l"4 Re [A*--k”(z)]“>O.

With z=+ we understand the integrals (2) as the limits of the integrals as are piecewise constant, then k*(z) -k+=x&iCir, A,==(h2-kt2)‘h, p(z) -pr==q&itr,

z-+2,. If

k*(z), p(z)

(3')

2,..., n. Z‘~z
regardless of what a jump of the coefficients or of their first derivatives, occurs at the point. The efficiency of a numerical.algorithm for computing the fundamental solutions of laminated media is largely determined by the quality of the algorithm for evaluating integrals since (2) approximately. The greatest difficulties arise in the case of large I‘and 1{/--~/'1, we need to be able to evaluate with high accuracy the integrals of rapidly oscillating functions. In the present papar the algorithm for the numerical integration of (21 is based on shifting the path of integration to a complex domain in which no such difficulties arise. We also prove some claims about the domains of analyticity and the asymptotic behaviour of the integrands, whereby the contour can be optimally deformed, and numerical results are quoted. This approach was earlier used in /1/ for the Helmholts equation in a medium with any number of homogeneous layers. It was shown to be applicable to the system of Maxwell's equations in /2/. In /3/thediscussion was confined to media with pure imaginary squares of wave numbers, while the domain where the integrands are analytic was not fully described. The esscence of tha algorithm is that, after standard transformations,the integrals (2) are reduced to the form

*Zh,vychhisl.l~at.mat,Fiz.,26,8,1234-1242,1986 178

179

(4)

where H,"' is the Hankel function of 1st kind of zero order, ~(2,I',h)SLu(ZIZ'.-h). It can not more be shown that u(z,z',A)is bounded at regular points, while it tends to 0 as h-a slowly than IA\_',and for H,"' we have the asymptotic form /4/ (5) Hence there is; the possibility of deforming the path of integration in (4) in such a way that the domain bounded by the new contour and the line Im(A)-0 consists of only regular points. In the indices of the exponential functions in (4) and (5) there then appear negative real parts, which ensure rapid decrease of the integrands. In short, the numerical integration of rapidly oscillating functions can be reduced to the much simpler problem of integrating slowly varying expressions with a known rapidly decreasing weighting function of the exponential type. From all the possible paths we have to choose in each specific case the one most convenient fox computations. Since, in order to choose the contour, we need to know the analytic properties and asymptotic form of the integrands, we shall study these questions next. 2. Since the domain in which the solution of problem (1) depends analytically on the parameter 1 is the same as the set of regular points, we will examine the solvability of the problem for complex values of h. Theorem 1. Let qfz)#O, t(z)*o. Then, for problem (If to have a unique solution, it that the parametex h belong to the domain

suffices

where

Proof. Notice that Eq.(l) is a symbolic form of the relation

lim u=O, *-*m which holds for all tp of space NR') of finite infinitely differentiable functions. The solution of problem (7) exists for any complex a and can be written as the sum of two terms, one of which takes account of the delta function on the right-hand side and can be constructed explicitly, while the other is the solution of an operator equation of the 2nd kind with completely continuous operator, mapping the Sobolev space W,'(R') with weight (cr+t)M' into itself. Assume that the solution of (7) is not unique and that u,, ua are two distinct solutions. Then, V-U,-u,EW,'(R') satisfies the homogeneous Eq.(7) and the relations r,(n)= j~P,cl.~z+[n(~z-~a-u)+t(~~~-~)~~~,')d~=O, -cQ

(84 (8b)

which are obtained from (7) by replacing q=R by the complex conjugate v' to u and separating real and imaginary parts with the aid of (3). It is easily verified directly that Let X+(z). q(3)r%“-~‘-x~z~l+t(a)~2%?-e~a~l~0,

(9)

whence it follows that

the equality being possible if u-0.

If X+=Qs(z), we have the inequality

t(z)~~~-fl~--X(Z)]-_4(z)~z~~--D(zf]~O. Consequently, I%(v)>0 for all

MO

with

&E

fl Q*(a). z&t’

(IO)

180

The theorem is proved.

Fig.2

Fig.1

Our sufficient condition (6) for problem (1) to be uniquely solvable, or what amounts to the same thing, for its solution to be analytic with respect to h is obtained under very general assumptions about the coefficients of the equation. If stronger constraints are imposed on the coefficients, our condition can be refined. Theorem

2.

Let

p(r)=Re(p(z))=q(z)>O, Xmnr= 2,

emin=

X(Z).

infCl(Z),

8 max -=

ER'

Then,

in

the

domain

the solution of problem (1) is an analytic function of Proof.

then I,(v)>0 ad

Sup 8 (2). IER'

h.

The condition n'<6*-xnuris equivalent to S'-~*-xmu>O. Since

with hE(hlt)'<$,*--~,.). Similarly, I,(u)>0 for h~(h~~<&,,./(2~)), theequalitybeingpossibleonlyif c-0. The theoremisproved. The domain Q is shown in Fig.1 (not hatched).

Iz(v)GO

for all v*O for

~~@h=%d(2&)),

Note 1. The case considered in /3/, when Re(k*(z))~O, Im(p(z))=O, is reducible to the above. The domain of analyticity is here the entire complex plane, except for the part bounded by and the line n--E. It is essentially wider than the domain hyperbolas n-Ck./(ZE), n-:8ms./(2E) Inl
3.

Let ~(z)=ak*(z), a>&

IP(z)

x(z)@, e(z)fO.Then, in the unhatched domainof

Fig.2,

U

Q=(~,~l["z~~~~~~~]z
n h w~+d45 ( ZER’ (I I k’ (4 I


of (1) is an analytic function of 1.. The proof is just the smae as for Theorem 1.

the solution

Theorem

all 1=1,2 ,...,

4. Let the coefficients of (1) be piecewise constant functions of z. 12, we have at least one of the four inequalities: *qr[Ik~*C,*‘+Re(hZ)]*~r ft,[ Ihtl’C,*‘+Re(hr2)]

Im(h,‘)>O,

(114

Fq, Im(h,') >O,

(11b)

where

c,=q-+

aI

=

sin(ImPd 4) Re(ad

sh(Re(4) 4 Im Qd ’

d-z,-r-z,, then parameter h Proof.

belongs to the domain where the solution of (1) is analytic.

We write integrals (8) in a more convenient form for future study: 11(u) =

2 [&‘f I=1

1u,’ I*dz + ‘I

”

Re(bi2;p-l~1m (at) x

If, for

“I-1 j

1I.J I’

dz] =

0,

181

and estimate the terms which contain the integrals of The general solution of Eq.(?) is

juS'l'in terms of integrals of

/vi'.

~~A~~-rexpfh~(z-zr-,)IfAllexp[--hr(z--zr)], A---A,.-0, where ZE(Z,, z,_,), 1=1,Z,...,n. Substituting (13) into we obtain

exp(- n) - exp (- VI”) a2 -

atI’

(13)

(12) and omitting the intermediate working,

.

As&-J Q

(Aat-14 +

(1&t-l I* + I 4 I")exp I- Re04 41 x sin(Im(Qd,) s”I!$!$~’ ’ 1 Im(h,) I) 'I-1 d,

i

’

exp'- ") - exp'- 'I*)(A,,_,A; I, -

+ AS~A:I_J) >

at*

I 1.~I*(I AI-ZI’ +

I -411’)expI- Re (&)dJ x

sh(~e(~J)d~) _ ( Rs(%)d, I where p-k&.

Consequently,

We similarly obtain the upper bound:

Now let the first inequality (lla) hold for aI.; I==%,&..., n. arrive at the inequality

for all uf0.

Then, using (14), we

If the second inequality flla) holds, we obtain, using (15),

for &O.

me two inequalities (lib) are considered similarly. They respectively ensure that I*(u) is positive and negative for all v*O. The theorem is proved. Note 2. As IAl- we have CPO, and inequalities (11) transform into relations (9), then (10)* If la/--, Cz(X)-+I , and instead of the first inequalities (lla) and (lib) we write a(2&+31)W, qr(2+e,f8mu/(2~)}.andis the same as the similardomain in the case Pf*)=RefP@))-gWW. When choosing the contours for evaluating integrals (2), an important role is played by the behaviour of the solution of problem (1) as IhI-=. Theorem 5. As we have the following asymptotic formsofthesolutionofproblem (1): bl-4"2 ii_ x u(z,z',~)=- fP(Z) “a PmrlexpI- a@(z-d)]

E 1 Pm

j4

%

P&

Pii1 - P1+ few[2& (2 - a.l)l + 0 (I h I-9 + 1 + 0 (I h Pii1 + p*+

1 z>z*,

u(z,z’,h)=-

7 [

1 + I 2 <

where

2’9

0 (1 h I-‘)

P (2) Pmtl +

‘1sPm+exP~g&-z’)l

1

f p~-xp~Pi-1 -I- p*-

‘6

0 Iexp

I-

j_

2%

12 -

ZP;+, Pj+ f

~dl+

+ Pj’

l-9), x

(16b)

Pj;l

0 fl h

I-91 , I

182 Proof.

The solution of (1) can be written as u(z,z',h)-_[p(z)]"H(z, z',h),

where u" satisfies in each of the intervals (2,-&l),l=f, 2,...,n, the equation

C

ii,,” -

P’i)

the matching conditions [p’“Z] -0, with

z=z~,I=l, 2,...,n-I,

rp-'(P"H).7=&

(Q3)

and the conditions at infinity lim ti=O;

w

Z-*03

is the jump of function W on passing through the point z=zl. It was shown in /5/ that homogeneous Eq.(17) has two linearly independent solutions with asymptotic behaviours WI

n,(z,h)--e~"[l+O(lAl-~)], a,(~,h)-e-~'[1+0(111-')]. Hence, in each of intervals (21, ZM). 14, sought in the form

2,...,

n,

the asymptotic form of the solution can be

u(z,h)-=P'%i(Z, ~)-p"(AII-*exp[ho(2-2,_,)] [I+

0(lJ4-‘)l+Aa~expi-b(z-a)]

Substituting (20) into (18) and into (19) and neglecting quantities of order obtain a system of 2n linear algebraic equations in As+,, A,,: (pl+)“* [AH-I sx~(- ~1) + -&I All+%exp (- w+d = 0, (pl+)l/n[Am

exp (-

A nl+Bexp (AI= A*,,=().

19 -

(piiJ’*[&+I +

I=&2

,...,

O(lhl-I),

we

(214

.421] -(pLJ"~[&+I -

n+dl= ~;‘~,h

cm

[l+O(lhl-a)]}.

(21b) n--i,

The solution of system (21) can be found by the method given in /6/. Omitting the intermediate algebra involved in solving system (21), and substituting the result into (20), we obtain the required asymptotic forms (16). The theorem is proved. 3. Analysis of (5) and (16) shows that the optimal contours (with respect to rate of decrease as ]A]+-) for integrals (4) are the straight lines

On them, the integrands decrease as lhl-'"exp(-IhlXIM-M'I) and IhI-'exp(-IhI IM-n-i']), and there are no oscillations, since the ratios of the imaginary to the real parts of the indices of the exponential functions are respectively

‘“,,;“=-;I -

and

(E-lEEI)IY-Y'/Iz--l

1

and have moduli not exceeding 1. When computing the integrals we can use the quadrature formulae for functions specified in the interval (0,-) with exponential weight, see /7/. Another method is to transform the integrals with infinite limits into integrals with finite limits. In either case we can use quadrature formulae of high algebraic order of accuracy, since analytic functions are integrated. Table 1

Different considerations can be used for choosing the contour. For instance, for evaluating a large number of integrals of type (4), corresponding to different values of M contour for which Im(h)-_IRe(h)], A-+-. Here again, there and M', it is best to use a c-on are no oscillations, while the integrated functions decrease like IL]-"Xexp[-]h] (Iz-z’I+r)/y2] Moreover, the condition Im(l)=-]Re(A)(as ]A.]+~ ensures and lnl-'exp[-_lhl(lz-z'lfly-_y'l)/Y2]. that h belongs to the domain where the solution of (1) is analytic in all the cases considered above.

183 z~O enables the required sheet of Riemann surface to be chosen. It will be satisfied if the necessary cuts are made. We find the lines of the cuts from the relations Re(h)-0, l-1, 2,...,n. They define pieces of hyperbolas n==x,l (ZE) from branching points h-fe’*“Ik,I to A==*tim. Integration was performed over the contour p={h(n=l%(},which was used in /l/ with M, M'=R*. By the replacement b=V(l--(VI)-'Xexp[insign (V)[4],the integrals over r are written as

c(M.M’)=~iu(z,Z.h(V))x s=P 4 h(V)IY-Y'lf msigl(")]](l-_IVl)-ZdV I[ and computed by the Gauss quadrature formulae. The results for initial data

are shown in Table 1. Numerical experiments shows that, by means of our algorithm, we can efficiently solve problems of computing the fields in laminated media with different values of M and M’ which are commonly encountered in practice. REFERENCES 1. SMAGIN S.I., Calculation of Green's function for Helmholtz's equation with one-dimensional piecewise constant wave number, In: Conditionally well-posed problems of mathematical Physics in the interpretation of geophysical observations (Matem. probl. geofiz.), VTs SO AN SSSR, Novosibirsk, 1978. 2. MAZALOV V.N. and SMAGIN S.I., The electromagnetic field of a point dipole in a horizontally laminated medium, In: Communications on applied mathematics. Models and equations, VTs AN SSSR, MOSCOW, 1980. 3. V'YURKOV V.V. and DREIZIN YU.A., On the computation of the fields of dipoles in conducting laminated media, Dokl Akad. Nauk SSSR, 262, 5, 1108-1112, 1982. 4. GRADSHTEIN I.S. and RYZHIK I.M., Tables of integrals, sums, series, and products (Tablitsy integralov, summ, ryadov i prozvedenii), Nauka, Moscow, 1962. 5. NAIMARK M.A., Linear differential operators (Lineinye differentsial'nyieoperatory), Nauka, Moscow, 1969. 6. DMITRIEV V.I., A General method of computing the electromagnetic field in laminated medium, In: Computational methods and programming (Vychisl.metody i programmirovanie), 10, Izd-vo MGU, Moscow, 1968. 7. KRYLOV V.I., Approximate evaluation of integrals (Priblizhennoevyzhislenie integralov), Nauka, MOSCOW, 1967.

Translated by D.E.B.