Mathematics
f(z)    Complex Numbers and Functions     Tech Note
Complex Bessel Functions:  J0 and I0

This page shows how to compute the Bessel Functions J0 and I0.

CJ0:  Complex Bessel function J0

This series expansion was implemented directly, with a limit imposed on the maximum number of terms computed, or on the size of the last term added to the series.

CJ0 function

// z = J0(a)
FUNCTION CJ0 (CONST a:  TComplex):  TComplex;
  CONST
    MaxTerm    : BYTE  = 35;
    EpsilonSqr : TReal = 1.0E-20;
 
  VAR
    addflag:  BOOLEAN;
    aTemp  :  TComplex;
    i      :  BYTE;
    SizeSqr:  TReal;
    term   :  TComplex;
    zSQR25 :  TComplex;
BEGIN
  aTemp := CConvert (a,cfRectangular);
  RESULT := ComplexOne;                // term 0
  zSQR25 := Cmult(aTemp, aTemp);
  zSQR25.x := 0.25 * zSQR25.x;
  zSQR25.y := 0.25 * zSQR25.y;
  term := zSQR25;
  RESULT := CSub(RESULT,zSQR25);       // term 1
  addflag := FALSE;
  i := 1;
  REPEAT
    term := CMult(zSQR25, term);
    INC (i);
    addflag := NOT addflag;
    term.x := term.x / SQR(i);
    term.y := term.y / SQR(i);
    IF   addflag
    THEN RESULT := CAdd(RESULT, term)        // sum := sum + term
    ELSE RESULT := CSub(RESULT, term);       // sum := sum - term
    SizeSqr := SQR(term.x) + SQR(term.y)
  UNTIL (i > MaxTerm) OR (SizeSqr < EpsilonSqr)
END {CJ0};

Special z values were used to test the above CJ0 function.  An array, b, of test values was defined as follows:

Test values “z”

// "Polar" definitions
b[ 1] := CSet( 0.0, 0.0, cfPolar);
b[ 2] := CSet( 0.5, 0.0, cfPolar);
b[ 3] := CSet(10.0, 0.0, cfPolar);
 
b[ 4] := CSet( 0.5, PI/6, cfPolar);
b[ 5] := CSet(10.0, PI/6, cfPolar);
 
b[ 6] := CSet( 0.5, PI/4, cfPolar);
b[ 7] := CSet(10.0, PI/4, cfPolar);
 
b[ 8] := CSet( 0.5, 5*PI/12, cfPolar);
b[ 9] := CSet(10.0, 5*PI/12, cfPolar);
 
b[10] := CSet( 0.5, PI/2, cfPolar);
b[11] := CSet(10.0, PI/2, cfPolar);

The test results were computed using this code:

Code to compute CJ0 test values

Memo.Lines.Add('Complex Bessel function CJ0 = J0(z)');
Memo.Lines.Add('');
Memo.Lines.Add('                                ' +
               '              J0(z)');
Memo.Lines.Add('              z(Polar)          ' +
               '           rectangular');
Memo.Lines.Add('     -------------------------  ' +
               '----------------------------------');
 
FOR k := Low(b) TO High(b) DO
BEGIN
  z := CJ0(b[k]);
 
  Memo.Lines.Add(   Format('%2d  %-s %-s',
                            [k, CToPolarStr(b[k],  9,4),
                                CToRectStr(z, 16,10)]))
END;

Here are the test results:

Sample computed values of CJ0

Complex Bessel function CJ0 = J0(z)
 
                                              J0(z)
              z(Polar)                     rectangular
     -------------------------  ----------------------------------
 1     0.0000 * CIS(   0.0000)     1.0000000000 +    0.0000000000i
 2     0.5000 * CIS(   0.0000)     0.9384698072 +    0.0000000000i
 3    10.0000 * CIS(   0.0000)    -0.2459357645 +    0.0000000000i
 4     0.5000 * CIS(   0.5236)     0.9682684872 -    0.0532808827i
 5    10.0000 * CIS(   0.5236)    -5.0471892329 -   18.1437389326i
 6     0.5000 * CIS(   0.7854)     0.9990234640 -    0.0624932184i
 7    10.0000 * CIS(   0.7854)   138.8404659416 -   56.3704585539i
 8     0.5000 * CIS(   1.3090)     1.0546148557 -    0.0321025326i
 9    10.0000 * CIS(   1.3090) -1546.3765101169 - 1270.8491488690i
10     0.5000 * CIS(   1.5708)     1.0634833707 +    0.0000000000i
11    10.0000 * CIS(   1.5708)  2815.7166284663 +    0.0000000000i

The above results match values in the Table of the Bessel Functions J0(z) and J1(z) for Complex Arguments.    The “z” values above were defined in polar coordinates to be consistent with the table of values.

 

CI0:  Complex Hyperbolic Bessel function I0

The terms in the CI0 series are the same as the CJ0 series but they are all additive (instead of the alternating signs with the CJ0 function), so CI0 is somewhat simpler than CJ0:

CI0 function

// z = I0(a)
FUNCTION CI0 (CONST a:  TComplex):  TComplex;
  CONST
    MaxTerm    : BYTE  = 35;
    EpsilonSqr : TReal = 1.0E-20;
 
  VAR
    aTemp  :  TComplex;
    i      :  BYTE;
    SizeSqr:  TReal;
    term   :  TComplex;
    zSQR25 :  TComplex;
BEGIN
  aTemp := CConvert (a,cfRectangular);
  RESULT := ComplexOne;                // term 0
  zSQR25 := Cmult(aTemp,aTemp);
  zSQR25.x := 0.25 * zSQR25.x;
  zSQR25.y := 0.25 * zSQR25.y;
  term := zSQR25;
  RESULT := CAdd(RESULT, zSQR25);      // term 1
  i := 1;
  REPEAT
    term := CMult(zSQR25, term);
    INC (i);
    term.x := term.x / SQR(i);
    term.y := term.y / SQR(i);
    RESULT := CAdd(RESULT, term);          // sum := sum + term
    SizeSqr := SQR(term.x) + SQR(term.y)
  UNTIL (i > MaxTerm) OR (SizeSqr < EpsilonSqr)
END {CI0};

I am not aware of any published tables for complex I0 values, so only a minimal set of  test cases are included here.

For purely imaginary arguments:

The first three test values in these results are correct given the formula above and the J0 results above. 

Sample computed values of CI0

Complex Bessel function CI0 = I0(z)
 
                              I0(z)
          z                rectangular
     ------------  ----------------------------
 1    0.0 +  0.0i   1.000000000 +  0.000000000i
 2    0.0 +  0.5i   0.938469807 +  0.000000000i
 3    0.0 + 10.0i  -0.245935764 +  0.000000000i
 4    1.0 +  0.0i   1.266065878 +  0.000000000i
 5    2.0 +  0.0i   2.279585302 +  0.000000000i
 6    3.0 +  0.0i   4.880792586 +  0.000000000i
 7    1.0 +  1.0i   0.937608477 +  0.496529948i
 8    5.0 +  3.0i -22.771493068 + 10.300893255i

Test values 4 through 6 above  were verified to be correct using Table 9.8 in Abramowitz and Stegun. 

I have no independent verification for test values 7 and 8 above, but I believe they are correct.

References

Abramowitz and Stegun, Handbook of Mathematical Functions, With Formulas, Graphs, and Mathematical Tables, Dover Publications, 1972.

Chapter 9, “Bessel Functions of Integer Order”:

-         Equation 9.1.12, p. 360, J0

-         Equation 9.6.12, p. 375, I0

Spanier and Oldham, An Atlas of Functions, Hemisphere Publishing Corp, 1987.

Chapter 49, The Hyperbolic Bessel Functions I0(x) and I1(x)

Chapter 52, The Bessel Coefficients J0(x) and J1(x)

Bowman, Frank, Introduction to Bessel Functions, Dover Publications, 1958.

Table of the Bessel Functions J0(z) and J1(z) for Complex Arguments, Columbia University Press, New York, 1947.  

Download ComplexMathLibrary

Updated 18 Feb 2002


since 24 June 2001