12711


FACTA UNIVERSITATIS  
Series: Electronics and Energetics Vol. 37, No 3, September 2024, pp. 517 – 529 

https://doi.org/10.2298/FUEE2403517V 

 

© 2024 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND 

Original scientific paper 

A HIGHLY ACCURATE ALGORITHM FOR COMPUTATION OF 

COMPLEX-VALUED BESSEL, NEUMANN AND 

HANKEL FUNCTIONS OF INTEGER ORDER  

Slavko Vujević, Tonći Modrić 

University of Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval 

Architecture, Split, Croatia  

ORCID iDs:  Slavko Vujević  https://orcid.org/0000-0002-1898-1746 
 Tonći Modrić  https://orcid.org/0000-0002-1411-1591 

 

Abstract. In this paper, a highly accurate algorithm for computation of complex-

valued Bessel, Neumann and Hankel functions of integer order is given. The algorithm 

enables the computation of these functions in the entire complex plane with quadruple 

precision, which can be reduced to double precision. The complex values of the Bessel 

and Neumann functions of the zeroth and first order can be computed in a special way 

for small, medium-sized and large arguments in the first quadrant of the complex plane. 

The mapping of functions from the first quadrant to the other quadrants is described by 

simple formulas. Bessel and Neumann functions of higher positive integer order can be 

computed using forward and backward recurrence relations. Two types of Hankel 

functions are linear combinations of the Bessel and Neumann functions. Bessel, 

Neumann and Hankel functions of negative integer order are equal to positive order 

functions up to the sign. 

Key words: Bessel function, Neumann function, Hankel function, complex-valued 

function, Gauss-Legendre quadrature 

1. INTRODUCTION 

Bessel functions of the first kind, Bessel functions of the second kind (Neumann 

functions) and Bessel functions of the third kind (Hankel functions) of the integer order 

occur frequently in mathematics, but also in various applied scientific fields, such as 

electromagnetic problems [1–10] and earthquake engineering problems [11]. Important 

formulas and various numerical methods for computation of these functions can be found 

in references [12–21]. 

 
Received May 13, 2024; revised June 25, 2024; accepted July 03, 2024 

Corresponding author: Tonći Modrić 

University of Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture, Ruđera 
Boškovića 32, 21000 Split, Croatia 

E-mail: tmodric@fesb.hr   



518 S. VUJEVIĆ, T. MODRIĆ 

 

In this paper, a new highly accurate algorithm for the computation of complex-valued 

Bessel, Neumann and Hankel functions of integer order is developed. The computation of 

Bessel, Neumann and Hankel functions of integer order in all quadrants of the complex 

plane is based on the computation of Bessel and Neumann functions of the zeroth and 

first order in the first quadrant of the complex plane. For numerical reasons, it is 

advisable to use the scaling of Bessel, Neumann and Hankel functions to avoid the 

numerical instability that occurs with large magnitudes of complex variables. 

Section 2 introduces scaled and unscaled Bessel, Neumann and Hankel functions. 

Section 3 contains equations for computation of Bessel and Neumann functions of the 

zeroth and first order in the first quadrant of the complex plane for small, medium-sized 

and large magnitudes of the complex argument. Equations for the translation of these 

functions from the other quadrants into the first quadrant are given in Sections 4–6. The 

computation of these functions of higher positive integer orders is presented in Section 7, 

and the validation of the proposed model is descriptively presented in Section 8. 

2. SCALED AND UNSCALED BESSEL, NEUMANN AND HANKEL FUNCTIONS 

The unscaled functions of integer order ... 2, 1, ,0 =n considered in this paper are 

Bessel functions ),(zJn Neumann functions )(zYn and Hankel functions, which are 

defined by: 

 )(j )()()1( zYzJzH nnn +=  (1) 

 )(j )()()2( zYzJzH nnn −=  (2) 

where yxz += j is a complex independent variable, which consists of real part x and 

imaginary part y, whereas j is the imaginary unit.
)1(

nH and
)2(

nH denote Hankel functions 

of the first and second kind, respectively. 

To avoid numerical stability problems, computation of scaled Bessel functions is 

introduced in the proposed model. In the first quadrant of the complex plane, as the 

imaginary part of the complex variable increases, the values of Bessel functions also 

increase.  

Scaled and unscaled functions are related by: 

 ( ) e ( )
ys

n nB z B z
−

=   (3) 

where: 

  (1) (2)( ) ( ), ( ), ( ), ( )n n n n nB z J z Y z H z H z  (4) 

  (1)s (2)s( ) ( ), ( ), ( ), ( )s s s

n n n n nB z J z Y z H z H z  (5) 

The following equations are valid for negative integer order: 

 ( ) ( 1) ( )n

n nB z B z− = −   (6) 

 ( ) ( 1) ( )s n s

n nB z B z− = −   (7) 



 Highly Accurate Computation of Complex-Valued Bessel, Neumann and Hankel Functions of Integer Order 519 

 

3. COMPUTATION OF BESSEL AND NEUMANN FUNCTIONS  

OF THE ZEROTH AND FIRST ORDER IN THE FIRST QUADRANT 

3.1. Small magnitude of the complex argument 

Computation of Bessel and Neumann functions of the zeroth and first order in the first 

quadrant of the complex plane ( 0, 0)x y  for small magnitude of the complex argument: 

 






=
precisionquadruplefor,5

precisiondoublefor,15
1bz  (8) 

is based on the computation of truncated power series [12, 14], which can be rewritten as: 

 

2

0
0

( )
20

im

i
i

z
J z a

=

  
      

  (9) 

 

2

1
0

( )
2 1 20

im
i

i

az z
J z

i=

  
     +   

  (10) 

 ( )
2

0 0
1

2 2
( ) ln ( )

2 20

im

i
i

z z
Y z C J z a i

=

    
  +  −             

  (11) 

 

2

1 1
0

( )2 2
( ) ln ( )

2 2 1 20

im
i

i

a iz z z
Y z C J z

z i=

     
 − +  +  −           +    

  (12) 

where: 

 

2

0 1

10
1 ; for 1, 2, ..., 40i ia a a i m

i
−

 
= = −  =  

 
 (13) 

 
1

1
( ) for 1

i

k

i i
k=

 =   (14) 

 

1
2 ( ) ,      for  1

( ) 1

1, for 0

i i
i i

i


 + 

 = +
 =

 (15) 

whereas C = 0.57721566490153286060651209008240243 is known as Euler constant. 

The parameter b1, introduced by Equation (8), as well as parameter m, which is equal 

to or less than 40, introduced by Equations (9) – (13), are determined based on our numerous 

numerical tests. The values of functions computed by proposed algorithm were compared 

with the values computed by the free online software package Wolfram Alpha. 

For Bessel and Neumann functions of the zeroth and first order, the summation is 

terminated after 40=m or after the following conditions are met: 

 

2 2

0
0

10 for  ( )
20 20

m im
q

m i
i

z z
a a J z−

=

   
      
   

  (16) 



520 S. VUJEVIĆ, T. MODRIĆ 

 

 

2 2

1
0

10 for  ( )
1 20 1 20

m im
qm i

i

a az z
J z

m i

−

=

   
      

+ +   
  (17) 

 

2 2

0
1

( ) ( ) 10 for  ( )
20 20

m im
q

m i
i

z z
a m a i Y z−

=

    
              

  (18) 

 

2 2

1
0

( ) ( )
10 for  ( )

1 20 1 20

m im
qm i

i

a m a iz z
Y z

m i

−

=

     
        + +    

  (19) 

The parameter q in Equations (16) – (19), determined by our numerical tests depending 

on the desired accuracy, have the following values:  

 






=
precisionquadruplefor,34

precisiondoublefor,17
q  (20) 

3.2. Medium-sized magnitude of the complex argument 

Computation of Bessel and Neumann functions of the zeroth and first order in the first 

quadrant of the complex plane ( 0, 0)x y  for medium-sized magnitude of the complex 

argument, ,21 bzb  where: 

 






=
precisionquadruplefor,50

precisiondoublefor,40
2b  (21) 

is based on the application of Gauss-Legendre quadrature to integral representations of 

Bessel and Neumann functions: 

 

/ 2

0

0

2
( ) cos( sin )J z z d



 


=     (22) 

 

/ 2

1

0

2
( ) sin( sin ) sinJ z z d



  


=      (23) 

 
/2

sinh

0

0 0

2 2
( ) sin ( sin ) e z tY z z d dt

 
− =      −  

 
   (24) 

 
/2

sinh

1

0 0

2 2
( ) cos ( sin ) sin e sinhz tY z z d t dt

 
− = −       −   

 
   (25) 

where ϑ and t are independent variables. 

The parameter b2, introduced by Equation (21), is determined on the basis of our numerous 

numerical tests. 

Due to the simpler numerical computation of the improper integral in Eq. (25), a new 

integral representation of 1( )Y z can be obtained by substituting sinh cosh e tt t −= − : 



 Highly Accurate Computation of Complex-Valued Bessel, Neumann and Hankel Functions of Integer Order 521 

 

 
/2

sinh

1

0 0

2 2 2
( ) cos( sin ) sin e t z tY z z d dt

z

 
− − = − −       +  

   
   (26) 

Computation of the scaled Bessel functions of the zeroth and first order for medium-sized 

magnitude of complex argument, 1 2 ,b z b  using the Gauss-Legendre quadrature in 

the first quadrant of the complex plain ( 0;  0)x y  is based on the following integral 

representations: 

 

/ 2

j sin j sin

0 0

0

1
( ) e ( ) (e e )s y y z y zJ z J z d



  


− − +   − −  =  =  +   (27) 

 

/ 2

j sin j sin

1 1

0

( ) e ( ) (e e ) sins y y z y zj
J z J z d



   


− − −   − +  =  =  −    (28) 

Computation of the scaled Neumann functions of the zeroth and first order for medium-

sized magnitude of complex argument, ,21 bzb  using the Gauss-Legendre quadrature in 

the first octant of the complex plane, where ,0  yx is based on the following integral 

representations: 

 

/ 2

j sin j sin sinh

0 0

0 0

j 2
( ) e ( ) (e e ) es y y z y z y z tY z Y z d dt


  

 



− − −   − +   − − =  =  −  −     (29) 

 

/ 2

j sin j sin

1 1

0

sinh

0

2 e 1
( ) e ( ) (e e ) sin

2
e

y
s y y z y z

t y z t

Y z Y z d
z

dt



   
 



−
− − +   − −  



− − − 


=  = − −  +  



+  





 (30) 

Our numerous numerical tests have shown that the infinite upper limits of the improper 

integrals in Eqs. (29) and (30) can be replaced by finite limits of the integrals without loss 

of accuracy: 

 ( )
  +−−


−−

0

0

sinhsinhj

0

sinh eee
mt

txytytzy dtdt  (31) 

 ( )
  ++−−


−−−

1

0

sinhsinhj

0

sinh eee
mt

txyttytzyt dtdt  (32) 

where finite upper limits of the integrals can be defined by equations: 

 






 −
==+ −

x

y
ttxy mm

75
sinh          75sinh 1

00  (33) 

 75sinh 11 =++ mm txyt  (34) 

Upper integral limit tm1, described by non-linear Equation (34), can be computed by 

Gauss-Newton iterative method using the following equation: 

 
( )

( )
1 1

1 1 1

1

( ) sinh ( ) 75
( ) ( ) ; 1, 2, ...

1 cosh ( )

m k m k

m k m k

m k

y t x t
t t k

x t
+

+ +  −
= − =

+ 
 (35) 

with an initial value: 



522 S. VUJEVIĆ, T. MODRIĆ 

 

 
1 0

1 0

75
( ) sinh m

m

y t
t

x

− − − 
=  

 
 (36) 

Computation of the scaled Neumann functions for medium-sized magnitude of complex 

argument, ,21 bzb  in the second octant of the complex plane, where ,0 xy is 

also carried out using Gauss-Legendre quadrature. In this case, due to the greater 

attenuation of the integrands of improper integrals, the following expressions are used: 

 0 0 0

2
( ) j ( ) ( j )Y z J z K z=  −  − 


 (37) 

 1 1 1

2
( ) j ( ) j ( j )Y z J z K z


=  +   −   (38) 

where K0 and K1 are modified Bessel functions of the second kind of the zeroth and first 

order, respectively. 

From the integral representation of modified Bessel functions [12, 13]: 

 
j cosh

0

0

( j ) e z tK z dt



 −  =   (39) 

 
j cosh

1

0

( j ) e coshz tK z t dt



 −  =    (40) 

where, due to numerical reasons, by substituting ,esinhcosh ttt −+= a new integral 

representation for K1 can be obtained: 

 
j cosh

1

0

j
( j ) e ej z t z tK z dt

z



 − +  −  =  +   (41) 

After substituting Eqs. (22), (23), (39) and (41) into Eqs. (37) and (38), a new integral 

representation of unscaled and scaled Neumann functions in the second octant of the 

complex plane can be written as: 

 
j cosh

0 0

0

2
( ) j ( ) e z tY z J z dt





 =  −    (42) 

 
j j cosh

1 1

0

2 2
( ) e j ( ) j ez t z tY z J z dt

z 



 − +  = −  +  +   
   (43) 

 
j cosh

0 0 0

0

2
( ) e ( ) j ( ) es y s y z tY z Y z J z dt





− − +  =  =  −    (44) 

 
j j cosh

1 1 1

0

2 2
( ) e ( ) e j ( ) j es y y z s t y z tY z Y z J z dt

z 



− − +  − − +  =  = −  +  +   
   (45) 

Infinite upper limits of improper integrals in Eqs. (44) and (45) can be replaced by finite 

limits of the integrals without loss of accuracy: 

 ( )
  +−


+−

0

0

cosh1coshj

0

coshj eee
mt

tytxtzy dtdt  (46) 



 Highly Accurate Computation of Complex-Valued Bessel, Neumann and Hankel Functions of Integer Order 523 

 

 ( )
  +−−


+−−

1

0

cosh1coshj

0

coshj ee e
mt

tyttxtzyt dtdt  (47) 

where finite upper limits of the integrals can be defined for y ≤ 37 by equations: 

 






 −
==+ −

y

y
ttosyy mm

75
cosh          75hc 1

00  (48) 

 75cosh 11 =++ mm tyyt  (49) 

Upper integral limit tm1, described by non-linear Equation (49), can be computed by 

Gauss-Newton iterative method using the following equation: 

 
( )

( )
1 1

1 1 1

1

( ) cosh ( ) 75
( ) ( ) ; 1, 2, ...

1 sinh ( )

m k m k

m k m k

m k

t y y t
t t k

y t
+

+ +  −
= − =

+ 
 (50) 

with an initial value: 

 
1 0

1 0

75
( ) cosh m

m

y t
t

y

− − − 
=  

 
 (51) 

For y > 37 it can be taken that ,010 == mm tt and therefore the truncated improper 

integrals given by Eqs. (46) and (47) are equal to zero. 

In the proposed algorithm, all finite integrals and truncated improper integrals are 

computed using Gauss-Legendre quadrature, where the total number of integration points 

is 100 for quadruple precision and 40 for double precision. In both cases, the algorithm is 

performed in quadruple precision. 

 

3.3. Large magnitude of the complex argument 

For large magnitude of the complex argument, ,2bz  asymptotic approximations of 

Bessel and Neumann functions of the zeroth and first order in the first quadrant of the 

complex plain ( 0, 0)x y  can be written as [12]: 

 
0 0 0

2
( ) ~ ( ) cos ( ) sin 

4 4
J z P z z Q z z

z

 



    
  − +  −    

     
 (52) 

 
1 1 1

2
( ) ~ ( ) cos ( ) sin 

4 4
J z Q z z P z z

z

 



    
  − +  −    

     
 (53) 

 
0 0 0

2
( ) ~ ( ) sin ( ) cos 

4 4
Y z P z z Q z z

z

 



    
  − −  −    

     
 (54) 

 
1 1 1

2
( ) ~ ( ) sin ( ) cos 

4 4
Y z Q z z P z z

z

 



    
  − −  −    

     
 (55) 

where auxiliary functions 0 0 1( ),  ( ),  ( )P z Q z P z and ( )zQ1 can be written as: 



524 S. VUJEVIĆ, T. MODRIĆ 

 

 

2
2

1
0 2 2

1

(2 1)

( ) 1 ;      ( 1)
(2 ) 8

k

m
kk i

kk k
k

i
p

P z p
z k

=

=

 −

 + = − 
 




!
 (56) 

 

2
2

1
0 2 2

1

(2 1)
1

( ) 1 ;      ( 1)
8 (2 1) 8

k

m
kk i

kk k
k

i
q

Q z q
z z k

=

=

 +
 

  + = −  
  +  




!
 (57) 

 

2
2

1
1 2 2

1

4 (2 1)

( ) 1 ;  ( 1)
(2 ) 8

k

m
kk i

kk k
k

i
f

P z f
z k

=

=

 −  − 
 + = − 

 




!
 (58) 

 

2
2

1
1 2 2

1

4 (2 1)
3

( ) 1 ;     ( 1)
8 (2 1) 8

k

m
kk i

kk k
k

i
g

Q z g
z z k

=

=

 −  +  
  + = −  

  +  




!
 (59) 

The parameter m in Equations (56) – (59) is based on our numerous numerical tests equal 

to or less than 25. The convergence is controlled in the same way as in the computation 

of Bessel and Neumann functions of the zeroth and first order given by Eqs. (16) – (19). 

Moreover, the unique convergence condition can be written as follows: 

 
2 2

1

1 10 ;     { ,  , ,  }
m

qm k
k k k k km k

k

c c
c p q f g

z z

−

=

 +    (60) 

Using the following mathematical identities: 

 
2

 cossin 

4
sin 

zz
z

−
=







 
−  (61) 

 
2

 cossin 

4
 cos

zz
z

+
=







 
−  (62) 

 
2 2(e 1) sin j (e 1) cos 

( , ) e sin 
2

y y
y x x

F x y z
− −

− +  −  − 
=  =  (63) 

 
2 2(e 1) cos j (e 1) sin

( , ) e cos 
2

y y
y x x

G x y z
− −

− −  +  − 
=  =  (64) 

asymptotic approximations of the scaled Bessel functions can be expressed by the 

following equations: 

 0 0 0 0

0

( , ) ( ( ) ( )) ( , ) ( ( ) ( ))
( ) ~s F x y P z Q z G x y P z Q z

J z
z

 + +  −


 (65) 

 1 1 1 1

1

( , ) ( ( ) ( )) ( , ) ( ( ) ( ))
( ) ~s F x y P z Q z G x y P z Q z

J z
z

 + −  −


 (66) 

 0 0 0 0

0

( , ) ( ( ) ( )) ( , ) ( ( ) ( ))
( ) ~s F x y P z Q z G x y P z Q z

Y z
z

 − −  +


 (67) 



 Highly Accurate Computation of Complex-Valued Bessel, Neumann and Hankel Functions of Integer Order 525 

 

 1 1 1 1

1

( , ) ( ( ) ( )) ( , ) ( ( ) ( ))
( ) ~s F x y P z Q z G x y P z Q z

Y z
z

−  − −  +


 (68) 

3.4. Computation of Neumann functions by Wronskians 

In the case when all four Bessel and Neumann functions of the zeroth and first order 

are computed, as well as if: 

 0 ( ) 0sJ z   (69) 

the following Wronskians can be used for computation of the unscaled and scaled 

Neumann function of the first order in the first quadrant of the complex plane: 

 
1 0

1

0

( ) ( ) 2
( )

( )

J z Y z z
Y z

z J z





   −
=

 
 (70) 

 
2

1 0
1

0

( ) ( ) 2 e
( )

( )

s s y
s

s

J z Y z z
Y z

z J z





−   − 
=

 
 (71) 

4. COMPUTATION OF FUNCTIONS OF THE ZEROTH AND FIRST ORDER  

IN THE SECOND QUADRANT 

In the second quadrant of the complex plain ( 0, 0)x y  computation of the unscaled 

and scaled Bessel, Neumann and Hankel functions of integer order n can be translated in the 

first quadrant using the following expressions: 

 * *( ) ( 1) ( ( ))n

n nJ z J z= −  −  (72) 

 
* * * *( ) ( 1) ( ( )) j 2 ( ( ))n

n n nY z Y z J z = −  − +   −   (73) 

 
(1) * * * *( ) ( 1) j ( ( )) ( ( ))n

n n nH z Y z J z = −   − − −   (74) 

 
(2) * * * *( ) ( 1) 3 ( ( )) j ( ( ))n

n n nH z J z Y z = −   − −  −   (75) 

 * *( ) ( 1) ( ( ))s n s

n nJ z J z= −  −  (76) 

 
* * * *( ) ( 1) ( ( )) j 2 ( ( ))s n s s

n n nY z Y z J z = −  − +   −   (77) 

 
(1) * * * *( ) ( 1) j ( ( )) ( ( ))s n s s

n n nH z Y z J z = −   − − −   (78) 

 
(2) * * * *( ) ( 1) 3 ( ( )) j ( ( ))s n s s

n n nH z J z Y z = −   − −  −   (79) 

where asterisk * denotes the complex conjugation. 



526 S. VUJEVIĆ, T. MODRIĆ 

 

5. COMPUTATION OF FUNCTIONS OF THE ZEROTH AND FIRST ORDER  

IN THE THIRD QUADRANT 

In the third quadrant of the complex plain ( 0, 0)x y  computation of the unscaled 

and scaled Bessel and Neumann functions of integer order n can be translated in the first 

quadrant using the following expressions: 

 ( ) ( 1) ( )n

n nJ z J z= −  −  (80) 

 ( ) ( 1) ( ) j 2 ( )n

n n nY z Y z J z = −  − −   −   (81) 

 
(1) ( ) ( 1) 3 ( ) j ( )n

n n nH z J z Y z = −   − +  −   (82) 

 
(2) 1( ) ( 1) ( ) j ( )n

n n nH z J z Y z+  = −  − +  −   (83) 

 ( ) ( 1) ( )s n s

n nJ z J z= −  −  (84) 

 ( )( ) ( 1) ( ) j 2 ( )s n s s

n n nY z Y z J z= −  − −   −  (85) 

 
(1) ( ) ( 1) 3 ( ) j ( )s n s s

n n nH z J z Y z = −   − +  −   (86) 

 
(2) 1( ) ( 1) ( ) j ( )s n s s

n n nH z J z Y z+  = −  − +  −   (87) 

6. COMPUTATION OF FUNCTIONS OF THE ZEROTH AND FIRST ORDER  

IN THE FOURTH QUADRANT 

In the fourth quadrant of the complex plain ( )0,0  yx computation of the unscaled 

and scaled Bessel and Neumann functions of integer order n can be translated in the first 

quadrant using the following expressions: 

 * *( ) ( ( ))n nJ z J z=  (88) 

 * *( ) ( ( ))n nY z Y z=  (89) 

 
(1) * * * * (2) * *( ) ( ( )) j ( ( )) ( ( ))n n n nH z J z Y z H z= +  =  (90) 

 
(2) * * * * (1) * *( ) ( ( )) j ( ( ))  ( ( ))n n n nH z J z Y z H z= −  =  (91) 

 * *( ) ( ( ))s s

n nJ z J z=  (92) 

 * *( ) ( ( ))s s

n nY z Y z=  (93) 

 
(1) * * * * (2) * *( ) ( ( )) j ( ( ))  ( ( ))s s s s

n n n nH z J z Y z H z= +  =  (94) 

 
(2) * * * * (1) * *( ) ( ( )) j ( ( ))  ( ( ))s s s s

n n n nH z J z Y z H z= −  =  (95) 



 Highly Accurate Computation of Complex-Valued Bessel, Neumann and Hankel Functions of Integer Order 527 

 

7. COMPUTATION OF FUNCTIONS OF HIGHER POSITIVE INTEGER ORDER 

For computation of Bessel, Neumann and Hankel functions of higher positive integer 

orders, forward and backward recurrence relations can be used [14], [19–24]. Forward 

recurrence relations can be written as: 

 1 2

2 ( 1)
( ) ( ) ( ); 2, 3, ...n n n

n
B z B z B z n

z
− −

 −
=  − =  (96) 

 1 2

2 ( 1)
( ) ( ) ( ); 2, 3, ...s s s

n n n

n
B z B z B z n

z
− −

 −
=  − =  (97) 

The Miller backward recurrence algorithm [14], [19–24] for the unscaled and scaled 

functions can be based on the assumptions that: 

 ( ) ( )n nB z S F z=   (98) 

 ( ) ( )s

n nB z S F z=   (99) 

where: 

 0 1

0 1

( ) ( )
or

( ) ( )

B z B z
S S

F z F z
= =  (100) 

 0 1

0 1

( ) ( )
or

( ) ( )

s sB z B z
S S

F z F z
= =  (101) 

where expressions with a larger denominator are chosen. 

The Miller backward recurrence relation can be written as: 

 1 2

2 ( 1)
( ) ( ) ( ); , 1, ..., 0n n n

n
F z F z F z n N N

z
+ +

 +
=  − = −  (102) 

where: 

 1 2( ) 1 ; ( ) 0N NF z F z+ += =  (103) 

Algorithms for estimating the integer parameter N can be found in [14]. 

8. MODEL VALIDATION 

The proposed model for highly accurate computation of Bessel, Neumann and Hankel 

functions of integer order in the entire complex plane for complex variables of arbitrary 

magnitude was implemented in a FORTRAN program. To determine the accuracy of the 

obtained results and the numerical stability of the model, a comparison was made with the 

publicly available program package Wolfram Alpha. While the developed FORTRAN 

program uses double-double precision computing, Wolfram Alpha can employ numbers with 

an arbitrary number of decimal places. 

The minimum number of matching digits, when comparing the results obtained in the 

program written based on the proposed model with the results obtained in the Wolfram 

Alpha program package, for complex variables of arbitrary magnitude and order of 

functions, is 30, which is a satisfactory level of accuracy of this model. 



528 S. VUJEVIĆ, T. MODRIĆ 

 

9. CONCLUSION 

The proposed algorithm for computation of complex-valued Bessel, Neumann and 

Hankel functions of integer order provides highly accurate results in the entire complex plane 

with quadruple precision, which can be reduced to double precision. The computation of these 

functions in all quadrants of the complex plane is based on the computation of Bessel and 

Neumann functions of the zeroth and first order in the first quadrant of the complex plane, 

using truncated power series for small arguments, Gauss-Legendre quadrature for medium-

sized arguments and asymptotic approximations of the Bessel functions for large arguments. 

Bessel, Neumann and Hankel functions of higher positive integer order can be computed 

using forward and backward recurrence relations, whereas these functions of negative integer 

order are equal to positive order functions up to the sign. For numerical reasons, the functions 

can be scaled. The obtained results were compared with the results provided by the Wolfram 

Alpha software package. The methodology presented in this paper is also applicable to the 

computation of modified Bessel functions, which is an area of future research. 

REFERENCES 

[1] S. Vujević and D. Lovrić, "High-accurate numerical computation of internal impedance of cylindrical 

conductors for complex arguments of arbitrary magnitude", IEEE Trans. Electromagn. Compat., vol. 56, 
no. 6, pp. 1431–1438, December 2014. 

[2] D. Lovrić and S. Vujević, "Accurate computation of internal impedance of two-layer cylindrical 

conductors for arguments of arbitrary magnitude", IEEE Trans. Electromagn. Compat., vol. 60, no. 2, pp. 

347–353, April 2018. 

[3] J. A. B. Faria, "A Matrix Approach for the Evaluation of the Internal Impedance of Multilayered 

Cylindrical Structures", PIER B, vol. 28, pp. 351–367, 2011. 
[4] K. Kubiczek and M. Kampik, "Highly Accurate and Numerically Stable Matrix Computations of the 

Internal Impedance of Multilayer Cylindrical Conductors", IEEE Trans. Electromagn. Compat., vol. 62, 

no. 1, pp. 204–211, February 2020. 
[5] K. Kubiczek, "Computation of the Characteristic Parameters of Coaxial Waveguides Used in Precision 

Sensors", Sensors, vol. 23, 2324, February 2023. 

[6] J. Acero, C. Carretero, I. Lope, R. Alonso, J. M. Burdío, "Analytical solution of the induced currents in 
multilayer cylindrical conductors under external electromagnetic sources", Applied Mathematical 

Modelling, vol. 40, Issues 23–24, pp. 10667–10678, December 2016. 

[7] K. Kubiczek and M. Kampik, "Fast and Numerically Stable Analytical Computations for the Power 
Induced in Cylindrical Multilayered Conductors Under External Magnetic Fields", IEEE Trans. 

Electromagn. Compat., vol. 65, no. 1, pp. 292–299, February 2023. 

[8] G. S. Liodakis, T. N. Kapetanakis, M. P. Ioannidou, A. T. Baklezos, N. S. Petrakis, C. D. Nikolopoulos, 
and I. O. Vardiambasis, "Electromagnetic Wave Scattering by a Multiple Core Model of Composite 

Cylindrical Wires at Oblique Incidence", Appl. Sci., vol. 12, 10172, October 2022. 

[9] R. Gordon, A. Choudhury, and T. Lu, "Gap plasmon mode of eccentric coaxial metal waveguide", Opt. 
Express, vol. 17, no. 17, pp. 5311–5320, March 2009. 

[10] S. Vujević and I. Krolo, "Computation of spectral-domain Green's functions of the infinitesimal current 

source in a planar multilayer medium", PIER B, vol. 100, pp. 55–71, 2023. 
[11] E. Zlatanović, V. Šešov, D. Lukić and Z. Bonić, "Mathematical interpretation of seismic wave scattering 

and refraction on tunnel structures of circular cross-section", Facta Universitatis, Series: Architecture 

and Civil Engineering, vol. 18, no. 3, pp. 241‒260, 2020. 
[12] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and 

Mathematical Tables. National Bureau of Standards, Washington, D.C., 1964, pp. 355–389. 

[13] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products. Elsevier Academic Press, 
Amsterdam, 2007. 

[14] S. Zhang and J. Jin, Computation of Special Functions. John Wiley & Sons, New York, 1996, pp. 126–201. 

[15] D. E. Amos, A Subroutine Package for Bessel Functions of a Complex Argument and Nonnegative 
Order. Sandia National Laboratories, Albuquerque, New Mexico, 1985. 



 Highly Accurate Computation of Complex-Valued Bessel, Neumann and Hankel Functions of Integer Order 529 

 

[16] D. E. Amos, "Algorithm 644: A portable package for Bessel functions of a complex argument and 
nonnegative order", ACM Trans. Math. Software, vol. 12, no. 3, pp. 265–273, September 1986. 

[17] J. P. Coleman and A. J. Monaghan, "Chebyshev expansions for the Bessel function Jn(z) in the complex 

plane", Mathematics of Computation, vol. 40, no. 161, pp.343–366, January 1983. 
[18] J. P. Coleman, "A Fortran subroutine for the Bessel function Jn(x) of order 0 to 10", Computer Physics 

Communications, vol. 21, pp. 109–118, 1980. 

[19] C. F. du Toit, "The numerical computation of Bessel functions of the first and second kind for integer 
orders and complex arguments", IEEE Trans. Antennas Propagat., vol. 38, no. 9, pp. 1341–1349, September 

1990. 

[20] C. F. du Toit, "Evaluation of some algorithms and programs for the computation of integer-order Bessel 

functions of the first and second kind with complex arguments", IEEE Antennas Propagat. Magaz., vol. 

35, no. 3, pp. 19–25, June 1993. 

[21] C. F. du Toit, "Bessel functions Jn(z) and Yn(z) of integer order and complex argument", Computer 
Physics Communications, vol. 78, pp. 181–189, 1993. 

[22] M. Goldstein and R. M. Thaler, "Recurrence techniques for the calculation of Bessel functions", 

Mathematics of Computation, vol. 13, no. 66, pp. 102–108, 1959. 
[23] F. W. J. Olver, "Error analysis of Miller's recurrence algorithm", Mathematics of Computation, vol. 18, 

no. 85, pp. 65–74, 1964. 

[24] W. Gautschi, "Computational aspects of three-term recurrence relations", SIAM Review, vol. 9, no. 1, pp. 
24–82, 1967.