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.