EDN -- 11.09.95 algorithm yields precise bessel functio

-November 09, 1995

Design Ideas:November 9, 1995

Algorithm yields precise Bessel function

Brian Geelen,
Satellite Microwave & Communications Ltd
Weert, The Netherlands

You frequently need to determine the Bessel function of the first kind, Jn(x), for the analysis of modulated and complex waveforms. Such analyses could include, for example, routine verification of transmitter phase- or frequency-modulation indices, by the measurement of the signal-suppression terms relative to an unmodulated-carrier reference level. Another example is the evaluation of communication-channel performance through sophisticated computer simulations.

Conventional solutions for Jn(x) using polynomial approximations are of limited use because of low precision and constraints on order and range arguments. Practical computing considerations generally restrict series expressions, although capable of any desired precision, to arguments of small values. By using series solutions for J0(x) and J1(x), which are easy to evaluate to more than eight decimal places for practical useful values of x&20, however, you can accurately determine higher orders of n through recursion. In addition, you avoid the need to use asymptotic expansions at a higher crossover point while considering these practical values. This approach results in a precise and efficient procedure that you can implement with a simple algorithm. To ensure a stable calculation, solve for Jn(x) with an upward recursion for x>n and a downward recursion for x

The following equation (Ref 1) gives the general power-series form for Jn(x):

 Equation 1

For J0(x), you can rewrite Eq 1 in the following summation form, which is more suitable for programming implementation:

 Equation 2

You can evaluate Eq 2 with high efficiency by using calculation to a length, where the last series term decreases to a value reflecting the required overall precision. In this way, you minimize computational effort and maintain the needed result accuracy. For example, with J0(5), you need evaluate only the loop to m=31 for the summation increment to be smaller than 1 x 10–41; you thus obtain a result accurate to more than 10 decimal places. Listing 1 gives an example implementation for J0(x) using Eq 2. The program is written in Fortran using double-precision declarations. You can obtain a solution for J1(x) in a similar way using the following equation:

 Equation 3

 LISTING 1—SERIES SOLUTION FOR J0(x) C DETERMINATION OF Jo(x) THROUGH A SERIES SOLUTI0N C AM=-1.0D0 FACT=1.0D0 AJ0=0.0D0 C DO 100 IJ@-(0)=1,40 C AM=AM+1.0D0 FACT=FACT+AM IF(AM.LT.0.5)FACT=1.0D0 A1=(-1.0D0)**AM A2=(0.5D0*AX)**(2.0D0*AM) A3=FACT*FACT AJ0=((A1*A2)/A3)+AJ0 DEL=(A1*A2)/A3 DEL2=DABS(DEL) IF(DEL2.LT.0.5D-40)GOTO 110 100 CONTINUE 110 CONTINUE

You can obtain higher order solutions for Jn(x), where n[greater than or equal to]2, by simply using recursion. This process allows you to determine higher orders to the same precision you obtain with the lower order input-parameter variables. To maintain stability against accumulation of round-off error, however, use upward recursion for x>n and downward recursion for xn by using the following equation, as shown in the example in Listing 2:

 Equation 4

 LISTING 2—FORTRAN ROUTINE FOR UPWARD RECURSION C DETERMINATION OF Jn(x) THROUGH UPWARD RECURSION C STABLE WHILE X IS > N C N1=INT(AN) IF(AX.GT.FLOAT(N1)GOTO 200 C C IF X

You obtain downward recursion with Eq 4 using a process of trial values and high arbitrary starting-seed values, as shown in the example in Listing 3. By combining the separate solutions for J0(x) and J1(x) with upward and downward recursion, you can determine Jn(x) from a concise computer program. This approach offers a precise and efficient algorithm for the Bessel function of the first kind, useful over a broad range of arguments and orders. EDN BBS/DI_SIG#1788 (dial (617) 558-4241 with modem settings 300/1200/2400 8,N,1) contains a Fortran listing and a menu-driven, executable Fortran routine. The program can solve for Jn(x) for either a constant n with x as the variable or a constant x with n as the variable. (DI #1788)

 LISTING 3—FORTRAN ROUTINE FOR DOWNWARD RECURSION C DETERMINATION OF Jn(x) THROUGH DOWNWARD RECURSION C STABLE WHILE X IS

References

1. Abramowitz, M and Stegun, IA, editors, Handbook of Mathematical Functions,
Dover Publications, 1965.