## `besselI, besselJ, besselK, besselY` -- the Bessel functions

#### Introduction

`besselJ(v, z)`, `besselI(v, z)`, `besselY(v, z)`, and `besselK(v, z)` represent the Bessel functions:

    J(v,z) = (z/2)^v/PI^(1/2)/gamma(v+1/2)*int(cos(z*cos(t))*sin(t)^(2v),t=0..PI),

    I(v,z) = (z/2)^v/PI^(1/2)/gamma(v+1/2)*int(exp(z*cos(t))*sin(t)^(2v),t=0..PI),

    Y(v,z) = (J(v,z)*cos(v*PI)-J(-v,z))/sin(v*PI),

    K(v,z) = PI/2*(I(-v,z)-I(v,z))/sin(v*PI)

`besselJ(v,z)` and `besselY(v,z)` are the Bessel functions of the first and second kinds, respectively; `besselI(v,z)` and `besselK(v,z)` are the corresponding modified Bessel functions.

#### Call(s)

`besselI(v, z)`
`besselJ(v, z)`
`besselK(v, z)`
`besselY(v, z)`

#### Parameters

`v, z` | - | arithmetical expressions

#### Returns

an arithmetical expression.

### Overloadable:

`z`

#### Side Effects

When called with floating point arguments, these functions are sensitive to the environment variable `DIGITS` which determines the numerical working precision.

#### Details

- The Bessel functions are defined for complex arguments v and z.
- A floating point value is returned if either of the arguments is a floating point number and the other argument is numerical. For most exact arguments the Bessel functions return an unevaluated function call. Special values at index v=0 and/or argument z=0 are implemented. Explicit symbolic expressions are returned, when the index v is a half integer. Cf. example 2.
- For nonnegative integer indices v some of the Bessel functions have a branch cut along the negative real axis. A jump occurs when crossing this cut. Cf. example 3.
- If floating point approximations are desired for arguments that are exact numerical expressions, then we recommend to use `besselJ(v, float(x))` rather than `float(besselJ(v, x))`. In particular, for half integer indices the symbolic result `besselJ(v,x)` is costly to compute. Further, floating point evaluation of the resulting symbolic expression may be numerically unstable. Cf. example 4.

#### Example 1

Unevaluated calls are returned for exact or symbolic arguments:

    >> besselJ(2, 1 + I), besselK(0, x), besselY(v, x)

                  besselJ(2, 1 + I), besselK(0, x), besselY(v, x)

Floating point values are returned for floating point arguments:

    >> besselI(2, 5.0), besselK(3.2 + I, 10000.0)

               17.50561497, 1.423757712e-4345 + 4.555796986e-4349 I

#### Example 2

Bessel functions can be expressed in terms of elementary functions if the index is an odd integer multiple of 1/2:

    >> besselJ(1/2, x), besselY(3/2, x)

                                    1/2 /            cos(x) \
                              1/2  2    | - sin(x) - ------ |
                      sin(x) 2          \              x    /
                      -----------, --------------------------
                       1/2   1/2            1/2   1/2
                      x    PI              x    PI

    >> besselI(7/2, x), besselK(-7/2, x)

           1/2 /  1   \1/2 /         / 15     \           / 6   15 \ \
          2    | ---- |    | cosh(x) | -- + 1 | - sinh(x) | - + -- | |,
               \ x PI /    |         |  2     |           | x    3 | |
                           \         \ x      /           \     x  / /

                        2    3       / PI \1/2
             (15 x + 6 x  + x  + 15) | -- |    exp(-x)
                                     \ 2  /
             -----------------------------------------
                               3  1/2
                              x  x

#### Example 3

The negative real axis is a branch cut of the Bessel functions for non-integer indices v. A jump occurs when crossing this cut:

    >> besselI(-3/4, -1.2), besselI(-3/4, -1.2 + I/10^10),
       besselI(-3/4, -1.2 - I/10^10)

          - 0.7606149199 - 0.7606149199 I,

             - 0.76061492 - 0.7606149199 I, - 0.76061492 + 0.7606149199 I

#### Example 4

The symbolic expressions returned by Bessel functions with half integer indices may be unsuitable for floating point evaluation:

    >> y := besselJ(51/2, PI)

          /  1/2 / 450675225   52650   1466947857375           \ \
          | 2    | --------- - ----- - ------------- + ... + 1 | | / PI
          |      |      4         2           6                | |
          \      \    PI        PI          PI                 / /

Floating point evaluation of this exact result is subject to numerical cancellation. The following result is dominated by round-off:

    >> float(y)

                                    8.862488737

The numerical working precision has to be increased to obtain a more accurate result:

    >> DIGITS:= 39: float(y)

           0.00000000000000000000116013005751977784273169237677941647977

Direct floating point evaluation via the Bessel function yields a correct result within working precision:

    >> DIGITS := 5: besselJ(51/2, float(PI))

                                    1.1601e-21

    >> delete y, DIGITS:

#### Example 5

The functions `diff`, `float`, `limit`, and `series` handle expressions involving the Bessel functions:

    >> diff(besselJ(0, x), x, x), float(ln(3 + besselI(17, sqrt(PI))))

                    besselJ(1, x)
                    ------------- - besselJ(0, x), 1.098612289
                          x

    >> limit(besselJ(2, x^2 + 1)*sqrt(x), x = infinity)

                                         0

    >> series(besselY(3, x)/x, x = infinity, 3)

           1/2 / 1 \3/2    /     7 PI \
          2    | - |    sin| x - ---- |
               \ x /       \      4   /
          ----------------------------- +
                        1/2
                      PI

                 1/2 / 1 \5/2    /     7 PI \
             35 2    | - |    cos| x - ---- |
                     \ x /       \      4   /    / / 1 \7/2 \
             -------------------------------- + O| | - |    |
                             1/2                 \ \ x /    /
                         8 PI

#### Background

- The Bessel functions are regular (holomorphic) functions of z throughout the z-plane cut along the negative real axis, and for fixed z<>0, each is an entire (integral) function of v.
- besselJ(v,z) and besselY(v,z) satisfy Bessel's equation in w(v,z):

      z^2*diff(w,z,z)+z*diff(w,z)+(z^2-v^2)*w = 0.

  Correspondingly, besselI(v,z) and besselK(v,z) satisfy the modified Bessel equation:

      z^2*diff(w,z,z)+z*diff(w,z)-(z^2+v^2)*w = 0.

- When the index v is an integer, the Bessel functions are governed by reflection formulas:

      besselI(-v,z) = besselI(v,z),
      besselJ(-v,z) = (-1)^v*besselJ(v,z),
      besselK(-v,z) = besselK(v,z),
      besselY(-v,z) = (-1)^v*besselJ(v,z).

#### Changes

- `besselI` and `besselJ` are new functions, `besselJ` and `besselY` have been rewritten and enhanced.
- Reference: M. Abramowitz and I. Stegun, ``Handbook of Mathematical Functions'', Dover Publications Inc., New York (1965).

Do you have questions or comments?
Copyright © SciFace Software GmbH & Co. KG 2000