5078
Analyzing the Modified Bessel Function of the Second Kind with IDL
This Help Article discusses "Analyzing the Modified Bessel Function of the Second Kind."
Determining Bessel Function Accuracy. A recurrence relationship between Bessel functions of differing order can be used to determine how accurately IDL is computing Bessel functions.
For Related information, see the following Help Articles:
Analyzing the Modified Bessel Function of the Second Kind. In the following example, the recurrence relationship

where K(x) is the modified Bessel function of the second kind of order n - 1, n, or n + 1 is used. Subtracting the Bessel function of order n from both sides results in the following equation.

Evaluating the left side of this equation will reveal how accurately IDL computes the modified Bessel function of the second kind.
The resulting plots are for n equal to 1 through 6. All of these plots show that this Bessel function is calculated within machine tolerance.
Code example:
PRO analyzingBESELK
; Derive x values.
x = (DINDGEN(1000) + 1.)/200. + 5.
; Initialize display window.
WINDOW, 0, TITLE = 'Modified Bessel Functions'
; Display the first 8 orders of the modified Bessel
; function of the second kind.
PLOT, x, BESELK(x, 0), /XSTYLE, /YSTYLE, $
XTITLE = 'x', YTITLE = 'f(x)', $
TITLE = 'Modified Bessel Functions of the Second Kind'
OPLOT, x, BESELK(x, 1), LINESTYLE = 1
OPLOT, x, BESELK(x, 2), LINESTYLE = 2
OPLOT, x, BESELK(x, 3), LINESTYLE = 3
OPLOT, x, BESELK(x, 4), LINESTYLE = 4
OPLOT, x, BESELK(x, 5), LINESTYLE = 5
OPLOT, x, BESELK(x, 6), LINESTYLE = 0
OPLOT, x, BESELK(x, 7), LINESTYLE = 1
; Initialize display window for recurrence relations.
WINDOW, 1, XSIZE = 896, YSIZE = 512, $
TITLE = 'Testing the Recurrence Relations'
!P.MULTI = [0, 2, 3, 0, 0] ; for multiple displays
; Initialize title variable.
nString = ['0', '1', '2', '3', '4', '5', '6', '7']
; Display recurrence relationships for order 1 to 6.
; NOTE: the results of these relationships should be
; very close to zero.
FOR n = 1, 6 DO BEGIN
equation = x*(BESELK(x, (n - 1)) - $
BESELK(x, (n + 1))) + 2.*FLOAT(n)*BESELK(x, n)
PLOT, x, equation, /XSTYLE, /YSTYLE, CHARSIZE = 1.5, $
TITLE = 'n = ' + nString[n] + ': Orders of ' + $
nString[n - 1] + ', ' + nString[n] + ', and ' + $
nString[n + 1]
PRINT, 'n = ' + nString[n] + ': '
PRINT, 'minimum = ', MIN(equation)
PRINT, 'maximum = ', MAX(equation)
ENDFOR
; Return display window back to its default setting, one
; display per window.
!P.MULTI = 0
END
_____________________________________________________
Reviewed by BC on 09/05/2014