Bation.1 Bur.au or ..»^-^"-

Library, ».W- ^^^^ JUN 2 9 1964 *

^ecUnlcciL fiote 92o. 2/6

COMPUTATION OF HANKEL FUNCTIONS

LESLIE A. BERRY

U. S. DEPARTMENT OF COMMERCE

NATIONAL BUREAU OF STANDARDS
THE NATIONAL BUREAU OF STANDARDS

The National Bureau of Standards is a principal focal point in the Federal Government for assuring

maximum application of the physical and engineering sciences to the advancement of technology in industry and commerce. Its responsibilities include development and maintenance of the national stand- ards of measurement, and the provisions of means for making measurements consistent with those standards ; determination of physical constants and properties of materials ; development of methods

for testing materials, mechanisms, and structures, and making such tests as may be necessary, particularly for government agencies; cooperation in the establishment of standard practices for incorporation in codes and specifications; advisory service to government agencies on scientific and technical problems; invention and development of devices to serve special needs of the Government; assistance

to industry, business, and consumers in the development and acceptance of commercial standards and

simplified trade practice recommendations; administration of programs in cooperation with United

States business groups and standards organizations for the development of international standards of

practice; and maintenance of a clearinghouse for the collection and dissemination of scientific, technical, and engineering information. The scope of the Bureau's activities is suggested in the following

listing of its four Institutes and their organizational units. Institute for Basic Standards. Electricity. Metrology. Heat. Radiation Physics. Mechanics. Applied Mathematics. Atomic Physics. Physical Chemistry. Laboratory Astrophysics.* Radio Standards Laboratory: Radio Standards Physics; Radio Standards Engineering.** Office of Standard Ref- erence Data. Institute for Materials Research. Analytical Chemistry. Polymers. Metallurgy. Inorganic Materials. Reactor Radiations. Cryogenics.** Office of Standard Reference Materials.

Central Radio Propagation Laboratory.** Ionosphere Research and Propagation. Troposphere

and Space Telecommunications. Radio Systems. Upper Atmosphere and Space Physics.

Institute for Applied Technology. Textiles and Apparel Technology Center. Building Research.

Industrial Equipment. Information Technology. Performance Test Development. Instrumentation.

Transport Systems. Office of Technical Services. Office of Weights and Measures. Office of Engineering Standards. Office of Industrial Services.

* NBS Group, Joint Institute for Laboratory Astrophysics at the University of Colorado.

* * Located at Boulder, Colorado.
NATIONAL BUREAU OF STANDARDS

Issued June 12, 1964

COMPUTATION OF HANKEL FUNCTIONS

Leslie A. Berry

Central Radio Propagation Laboratory

National Bureau of Standards

Boulder, Colorado

NBS Technical Notes are designed to supplement the Bureau's regular publications program. They provide a means for making available scientific data that are of transient or limited interest. Technical Notes may be

listed or referred to in the open literature.

For sale by the Superintendent of Documents, U. S. Government Printing Office Washlnaton, D.C. 20402

Price: 2 cents
Foreword

This work was sponsored in part by Rome Air Development

Center, Griffiss Air Force Base, New York, in connection with Tasks

4 and 5 of D. O. No. AF-30(602) -2488, and in part by NBS Project

No. 85111.

11
Contents

Page

Foreword ii

Abstract V

1

Introduction 1

2. Complex Order, v or z Large 2

2.1. Determination of T] 4

3. Real Order 5

3. 1. Small |z I 5

3. 2. Large |z | 6

3. 3, The Convergence Factor 7

4. Discussion of Results 9

5. Conclusions 10

6. References 11

7. List of Symbols 13

Table 1 14

Figures 1 through 5 15

111

Computation of Hankel Functions

Leslie A. BerryMethods of computing solutions to Bessel's equation,

z^ —5- + z —^ + (z^ - v^) y = 0, for complex v and z, and

dz dz

for small, real v, are given in detail. The methods are

most suitable for use on an electronic digital computer.

With a computer which carries eight decimal digits, the

methods presented for small, real v will return at least

six correct significant digits for all z. When v and z are

both complex, the accuracy depends on |z | (since asymptotic

forms are used) and the ratio v/z. Relative error curves

are shown as a function of |v/z | parametric in |z |. For

V sufficiently far from z, seven correct digits are obtained,

which for v — z, the relative error is of the order 10 ^

/ \z\.

1. Introduction

Solutions to electromagnetic or sound-wave problems in spherical

or cylindrical co-ordinates usually involve solutions of Bessel's equation, z^ —X- + z -r- + (z^ - v^)y = 0. The Hankel functions [Watson,

dz dz

1958] Hy (z) and H^ (z) are one set of linearly independent solutions.

Other solutions can be written in terms of these solutions [see (19)].

The superscript indicates the "kind, " v, the "order, " and z, the

"argument" of the functions, and either v, or z, or both may be complex. The purpose of this paper is to present in one place detailed

instructions for computing the Hankel functions on an electronic digital

computer, and to investigate quantitatively the relative errors incurred.

No new mathematical development is involved, but since considerable

algebraic reduction and systemization of scattered results were necessary,

publication of the results seems desirable.

Section 2 details methods of conciputing Hy(z) when v and z are

complex. If V — z, Hi(z) must be computed accurately, so techniques

are given in section 3 for computing Hv(z) when v is real and small

and z is unrestricted. The accuracy of the methods is discussed in

section 4,
2, Complex Order, v or z Large

In this section, it is assumed only that either v or z is large.

What is meant by "large" will be discussed quantitatively in section 4,

but the asymptotic expressions used here [(1) and table 1 ] accurately

describe the qualitative behavior of the Hankel functions over raost of

the V and z planes, since the minimum relative error is about one per

cent for |v| — |z| ^l (figure 5).

Two approximations are used. In the first, which will be called

the Hankel approximation [Fok, 1934],

H^/^{z) s^exp [(-1)1^-1 in/6] Vl - ri coth T] Hi^^(-i?), k=l, 2, (1)

where

^ = V (tanh r\ - r]), (2)

cosh r\ = v/z, and |arg v| < tt/2 [Olver, 1954]. For |arg v| > n/2, (3)

(9) is used.

The principal square root is used; i. e. ,

arg vw = g- arg w,

and - TT < arg w = TT.

for v/z ^1 [Fok, 1946],

-^ = 1 0-Sr [-10 -f^) -0-1^)1- <^)

The first term on the right has been used extensively as an

approximation to - i§ in radio propagation theory [Fok, 1946; Bremmer,

1949].
The second approximation, called the Debye approximation, will

be written in terras of [Watson, 1958]

^-— ivtanhri ^^^

^'" (^j tanh T] J

The combination of S\ (z) and S^ (z) to be used depends on the

location of v/z in the complex plane. The first four A„ are [British

Association Mathematical Tables, 1952]:

Ao = 1,

15 ,2 ^^ = 8 - 24 ^°'^ '^'

3 77 ^,2 385 ^, 4 ^^ = 128 - 576 "°*^ ^+3456 "°* ^'

5 1521 , 2 17017 , 4 17017 , e >3 - l024 - 256^ ^°^^ ^+ I3i240

^°* ^

" 248832 ^°^^ ^- ^^^

Forms for A4 -- A7 can also be found in the above reference. In

(5), use the principal square root. The series in (5) can be written

00 Y

r(m+|) A„ A^coth 11 3 Ag coth^ ri 1 5 A3 coth^ T] L r(^) (±ivtanhTl)»

" V

^ v^ v^

m =0

(7)

This is an asymptotic series and should be truncated at the smallest

term; i. e. , all terms preceding the smallest term should be summed.

The regions of the plane to be considered are shown in figure 1,

and the expressions for Hy(z) are given in table 1. These expressions

are correct only if Re(z) ^ 0, but if Re(z) < 0, the formulas [Watson,

1958]
and

sinTTv smTTv

where k is any integer, can be used to evaluate the desired function.

Although expressions can be derived for Hv(z) in regions 6b, 7b, and

9, it is simpler to compute H_y(z) [so that -v/z is in 6a, 7a, or 8]

and use [Watson, 1958]

h1^^(z) = exp [(-1)^ vni] h1.V(z). k = 1, 2. (9)

2. 1. Determination of T]

Directions for computing T] and for determining the correct region

of the plane, figure 1, will now be given. Let T\ = a +i3, where a and

3 are real. Then a + i3 = cosh"''" v/z. The infinitely many-valued

function, cosh"'', is made single valued by requiring that ^ 3 < TT [Watson, 1958]. 3 = if v/z is real and v/z ^ 1. Specifically

a=-log(|^±y^ -1 I). (10)

The sign is chosen so that a has the same sign as Im(v/z). If

Im(v/z) = 0, and |v/z j > 1, the positive sign is chosen; if |v/z | ^ 1,

11 = Cos"^ v/z. Then

3 = Cos"^ [Re(v/z)/cosh a] (U)

These directions are exact for use in (5), and are valid for (1) over

most of the plane. (For a discussion of the region of validity for the

Hankel approximation, consult Olver [1954].) For (1), when v = z, arg (-i§) can be determined by examining the first term ^f (4); when v

is not near z, the Debye approximation is superior.

The m in table 1 is the largest integer less than

{{1 - a tanh a) tan 3 + 3}/TT.
Now, let

f(a, P) = 1 - e cot 3 - a tanh a . (12)

The signs of f(a, 0), f(a, 0) + n cot 0, Re(v/z), and Im(v/z) uniquely

determine in which of the regions shown in figure 1 v/z is located.

Coliamn 2 of table 1 shows how the region can be determined.

3. Real Order

To use (1), a routine for computing Hi(z) must be available. In

this section, methods are given for computing Hv(z) for small, real v. The restriction that v be real is probably not necessary if z^ is carefully defined, but it is retained in this section,

3. 1. Small |z |

If V is not an integer [Watson, 1958]

h1^^(z) = {-1)^ [exp [(-1)'^ vHi] J^(z) - J_^(z)3/i sin vTT, k = 1, 2,

(13)

where Jv(z) is the Bessel function of the first kind, defined by [Watson,

1958]

v(z) = y

(-l)''(z/2)^-"^''

L m ! r (v+ m+ 1)

m =0

(v/2)" r (z/2)^ (z/2)^ {z/2)^

r(v+l) V 1 ! (v+1) 2 ! (v+l)(v+2) 3 ! (v+l)(v+2)(v+3)

(14)

Jy(z) is made single-valued by specifying arg (z/2)^ = v arg (z/2) for

A - 2tt < arg (z/2) ^ A. Ordinarily, A = TT, but may have different

values for some applications.

If V is an integer, (13) is indeterminate. To compute Hn(z) for

integral n and small z, compute Jo(z) and Ji(z) using (14) and then find

the Bessel function of the second kind, Yo(z), from [Watson, 1958]
Y,(.)=|[<Y+log(f))J.(z)+f '^(^'•^,'^' ]. (15)

m=i

where y is Euler's constant, y = . 5772156649 • • •. The Jn(z), m > 1 are found with the recurrence formula [Watson, 1958],

J.+ i(2) =^ J„(z) - J„-i(z). (16)

Then [Watson, 1958]

y,(z) = {j,(z) Yo(z) +:^] /Jo(z) (17)

and

Finally,

Yn,.-i(z) =^ Y„(z) - Y„-,(z). (18)

H^^(z) . J„(z) + (-l)^-n Yn(z). (19)

3. 2. Large |z | As |z I gets large, (14) converges slowly and loses precision.

More serious is the precision lost in (13) and (19), since in half of the

z -plane the terms on the right are nearly negatives of each other when

Im(z) is large. For example, four significant digits are lost in (19)

when n = 1, z = 5 exp (i 85°), It is therefore necessary to use the

asymptotic series [Watson, 1958]-

jj(l)|^j ^ exp(iz-i(v+ i)T./ 2) T^,.2i^)_ ,.„< ^^g , < 2„)_ (20a)

y n z/2

jj(^3)|^j __ exp(-iz-H(v+|)n£2) ^^j2i.), (-2n<argz<n), (20b)

^J n z/2
where

T,(z) = wy

'^-'-^'"^-'-^''---'^-'-'^"-"''

. (21) ^ n! (4zr

n = l

Unless V = n + -g, for integral n, Ty(z) is an asynaptotic series and

should be stopped at the smallest term. If v = n + ^, Tv(z) terminates,

and (20) is the finite series for Hn + i(z).

Using a computer which carries eight decimal digits in floating

point operations, about four correct significant digits in Hy(z) are

obtained over all of the z-plane, using the methods presented above.

This accuracy is increased to at least six digits if Dingle's [1958,

1959 ] convergence factor is used in (21), as described below.

3. 3. The Convergence Factor

In a series of papers, R. B, Dingle has developed a general

theory [Dingle, 1958] of convergence factors and has applied this theory

to many special functions [Dingle, 1958, 1959]. When the product of a convergence factor, Cb(z), and the nth term of an asymptotic series is

added to the previous n - 1 terms, the sum is the exact sum of the

entire series. In practice, of course, the sum is not "exact," but the

use of a convergence factor does significantly extend the range of

application of an asymptotic series.

In the particxolar case of the Hankel functions, (21) is written:

N-l

(4v^ - l^)(4vn! (4z)"

Y

(4v"-l^)(4v"-3")---(4v^-(2n-l)^)

^y^^) = ^ + i n!(4zr

nz:l

(4vS.l3)...(4v^_(2N-l)^)

+ N! (4z)N ^n(v, z), (22)

where N is chosen so that the term corresponding to N-l is the smallest

term of (21). Then, [ Dingle, 1959] for s = N - v - |,
r f^ .\ ~\ r(v + |)r(N+v-t + |) t .(t), .

t=o

' ^ ' (N+v-^) ' ^ ' (N+v-i)(N+v-|-) ' ^ '

where [Dingle, 1958]

(23)

Ai°^z) r{s + l) J 1+x/z '

2 A ^z) . (s + z+l) A^ '(z) - z,

.(2) (i)/.x , a(°), 2zA^, '(z) = (s + z) A^, '(z) + A^^ >{z) - 1,

and

tzAl*)(z) = (s + z + 2-t) a[*"')(z) + A['~^h^), t>2. (24)

When |s I ^ 1, or |z | ^ 1, and |arg z | < 3TT/4 [Dingle, 1958],

a'/N^)-

^

s + z

1 - z(s - 2z) z (s^ - 8sz + 6z^)

(s + z)*^ (s + z) (s + z)'

•)• (25)

(0) Otherwise, A^ '(z) can be found by numerical integration. For the

general routine used with this paper, A^ '(z) was found with an eight

point Gaussian quadrature formula:

x^ e-'^ f(x) dx ^ ^ W, f(x,), (26)

^=1

since s ^ 5 for use in (24). The W^ and x^ are taken from Rabinowitz

and Weiss [1959]. Then Cfg(v, z) is approximated by the first four

terms of (23).

8
4, Discussion of Results

Fortran subroutines were written to compute Hj^(z) and Ha(z) and

3 3 Hn(z) using the methods of section 3. The values for Hj^(z) and H2(z)

3, 3 were checked against the Tables of Modified Hankel Functions of Order

One -Third and Their Derivatives [The Staff of the Computation

Laboratory, 1945]. The computed values agree with the tabular values

to at least six significant digits where the tables give six digits. (The

tables give eight decimal places. ) The values of Hn(z) computed by the

routine were checked against NBS tables for Jo(z) and J3^(z), [Mathematical Tables Project, 1947] and Yo(z) and Y (z) [Computation

Laboratory, 1950] and gave comparable results.

The Debye and Hankel approximations were programmed separately and checked for real order by comparing the values they returned

with the correct (to at least six significant digits) answers given by the

methods of section 3. A subroutine valid for small, complex v and z was obtained from SHARE [Goldstein, Kresge', and Chen, I960], and

the asymptotic routines described here returned answers correct to

the expected accuracy, providing an independent check of the routines'

logic.

Defining

„ , . I answer - approximation i ,_^, Relative error = ^^ , (27)

answer

figures 2, 3, and 4 show relative error curves for both approximations

as functions of Iv/z|, parametric in |z | for |arg z| = 0°, 5 , and 15°,

respectively. It is apparent from these curves that the Debye approximation is the better except when v — z. Even when |v/z | = 1, if

|arg v/z I and |z | are not small, the Debye approximation is superior,

as shown in figure 5. Nevertheless, any general routine for computing

Hv(z) must include the Hankel approximation. The relative -error

curves all show that, for any z, there is an € > such that, if jv-z | < e, the Debye approximation is useless.

If series (7) is truncated at the smallest term, the relative error

in the Debye approximation is of the order of the magnitude of that

term. Examination of figures 2 through 5 yields an empirical formula.
[relative error I Sf 10 ^

I \z\, (28)

for the Hankel approximation (1).

A general routine for Hy(z) was written, incorporating both the

Debye and Hankel approximations. The routine uses the Debye approximation unless the smallest term of series (7) is larger than 10"^/ (z |,

in which case the Hankel approximation is used. This routine and

other computer routines described in this paper are available from

SHARE [Berry, 1963],

5. Conclusions

Using an electronic computer, solutions to Bessel's equation of

real order and complex argiiment can be computed very accurately.

Even when the order is complex, computations can be made sufficiently

accurate for most applications. The subroutines incorporating the

methods should be valuable in radio and sound wave computations.

10
6. References

Berry, L. A. (1963), SHARE subroutines HNKL and HFl 3, SHARE

Distribution No s. 1488 and 1489. Also order SHARE Distribution

Nos. 1490--1493 for necessary subroutines,

Bremmer, H. (1949), Terrestrial radio waves; theory of propagation

(Elsevier Publ. Co., New York, N. Y. ).

Computation Laboratory, NBS (1950), Tables of the Bessel functions

Yo(z) and Yi(z) for connplex arguments (Columbia University

Press, N. Y. ).

British Association Mathematical Tables (1952), Bessel functions of

positive integer order, 10 (Cambridge at the University Press).

Dingle, R. B. (Apr. 1958), Asymptotic expansions and converging

factors, Parts I-IU, Proc. Roy. Soc. (London) Ser. A: 244 , 456-491; (Jan. 1959), Parts IV, V, Proc. Roy. Soc. (London)

Ser. A: 249, 270-292.

Fok, V. A. (1934), New asymptotic expression for the Bessel functions,

Doklady, A. N. USSR ]_, 3, 97-99.

Fok, V. A. (1946), Diffraction of radio waves around the earth's

surface (Publishers of the Academy of Science, USSR, Moscow).

Goldstein, M. , M. Kresge, and S. Chen (I960), NUBES 3, SHARE

Distribution No. 469.

Mathematical Tables Project, NBS (1947), Tables of the Bessel functions

Jj,(z) and Ji(z) for complex arguments, second ed. (Columbia

University Press, N. Y. ).

Olver, F. W. J. (1954), The asymptotic expansion of Bessel functions

of large order, Phil. Trans. Roy. Soc. (London), Ser. A: 247,

328-368.

Rabinowitz, P., and G. Weiss (1959), Tables of abscissa and weights

for numerical evaluation of integrals of the form | e~''x"f(x) dx,

MTAC, Xni, 285-294. o

11
The Staff of the Computation Laboratory (1945), Tables of the modified

Hankel functions of order one -third and of their derivatives

(Harvard University Press, Cambridge, Massachusetts).

Watson, G. N. (1958), A treatise on the theory of Bessel functions,

second ed. (Cambridge at the University Press).

12
7. List of Symbols

H\ (z) Hankel fxinction of argument z, order v, and kind k,

k = 1, 2.

Jy(z) Be s sal function of the first kind, of argixment z and order

V.

Yy{z) Bessel function of the second kind, of argument z and

order v, sometimes called Neumann's function and denoted

N,(z).

Tv(z) The series defined in (21).

T] = a + ip = cosh"^ —

z

(k)

S^ (z) The series defined in (5).

(k)

Aj, The coefficient of S^ (z), see (6).

f(a, 3) = 1 - 3 cot 3 - a tanh a

m = Largest interger less than [(1 - a tanh a) tan + 3}/Tr.

C|s|(z) Convergence factor for asymptotic series, see section 3. 3.

13
Region

Determined

by „;»,z, „;^',.,

1 Ua, p)>

Ua, P)+TT cot P> s;»,., sj^'u,

2 f(a.P)<0

f(ff, P)+TT cot (3>0 Im V

^ <

z

s;»,., - s;^',., s;^',.,

3 f(«.P)<o

f(Qr, P)+TTCOtp>

Im .Vj <

z

s;",., s;^',., - s;" ,.,

4 f(<», P)>

f{a, P)+Trcotp<

Im (v/z) >

S^ (z) „ (2) 2vTri (1) S (z)-e S (z) V V

5 Ha. P)>

f(a, p)+-iTcotp<

Im (v/z)<

„ (1), , -2vTTi„ (2), , S^ (z)-e S^ (z) s;^',.,

6a f(a,P)<

f(a,P)+Tr cotp<0

Im (v/z) > Re (v/z) >

P (m+l)vTii sin mvir

1 sin VTT s^%)J%) V V

S^^'(z) + e<"^+^>^^^ ^^" "^^" S<Hz) V sin VTT V

6b f(a,p)<0

f(Q', p)+Tr cotp < Im (v/z) >

Re (v/z) <

Compute H (z) and

-V use (9)

Compute H (z) and

use (9)

7a Ha,?,) <0

f(a,P)+Trcotp < Im (v/z) < Re (v/z) <

<i''\) g-(m+l)vTri sinmvTT^(2)

V sin VTT v

r ^-(m+DvTTi sin mvul (2|(1)

L sin VTT 1 V * ' v * '

7b f(a,P) <0

f(Q',P)+TTCOtP <0 Im (v/z) < Re (v/z) <

Compute H (z) and

use (9)

Compute H '(z) and use (9)

8 Im (v/z) = Re (v/z) > 1

1 c -

(1), > e (2), V S^ (z) - S^ (z) ^S<^\z).s(2\z)

2 V v

9 Im (v/z) = Re (v/z) < - 1 Compute H (z) and

-V use (9)

Compute H (z) and

use (9)

Table 1

Ik
Im(v/z)

v/z plane

'8— Re(V'z)

Figure 1. Regions of the v/z plane (table 1).

15
10 -I

10"

10'

gio--

u

LlI >

10-6

10"

10-8

RELATIVE ERROR

I org Vz| = 0"

DEBYE

HANKEL

|Z|=IO,

0.1

Figure 2.

0.2 10

Relative errors of the Debye and Hankel approximations for

|arg v/z I = 0°, for various Iz 1.

16
10"

10-2

RELATIVE ERROR

I org Vz I = 5°

DEBYE

HANKEL

10"

gio-^

a:

lij

LJ>

iij 10"^

10-6

10"

lO"®"-

0.1

Figure 3.

0.2 5 10

E.elative errors of the Debye and Hankel approximations for

(arg v/z I = 5° for various |z |.

IT
10-I

RELATIVE ERROR

|arg Vz\ = IS'*

10"

10-3

^ IZl =100

1,0-

UJUJ>

<

^10-5

|Z| =1000

10-6

101-7

10-8

|Z| = I

0.1

Figure 4,

0.2 0.5 I 2 5 10

Relative errors of the Debye and Hankel approximations for

|arg v/z I = 15°.

18
RELATIVE ERROR

|v/2| = I DEBYE HANKEL

20 30 40

larg v/z|, DEGREES

50 60

Figure 5. Relative errors of the Debye and Hankel approximations

for |v/z 1 = 1, as a function of |arg v/z |, for various |z

19


1