Updated 2012-01-25 23:38:56 by dkf

Arjen Markus (15 june 2004) While looking for some background information on special mathematical functions I found a very convenient formula to approximate zeroth and first order Bessel functions (for small enough arguments). It is based on Gaussian quadrature, a well-known approximation of definite integrals. So, here are the humble beginnings of a new mathematical module ...
``` # bessel.tcl --
#    Evaluate the most common Bessel functions via quadrature formulas
#

# namespace special
#    Create a convenient namespace for the "special" mathematical functions
#
namespace eval ::math::special {
#
# Define a number of common mathematical constants
#
variable pi 3.1415926
}

# J0 --
#    Zeroth-order Bessel function
#
# Arguments:
#    x         Value of the x-coordinate
# Result:
#    Value of J0(x)
#
proc ::math::special::J0 {x} {
variable pi
#
# Determine the number of components (very crude heuristic)
# (We use a quadrature formula here and use the symmetry of
# the cos function to reduce the number of actually computed
# components by a gfactor 2.
#

set ncomps [expr {2*int(\$x/2.0+0.5)}]
if { \$ncomps < 6 } {
set ncomps 6
}

set result 0.0
for { set i 1 } { \$i < \$ncomps } { incr i 2 } {
set result [expr {\$result + cos( \$x * cos((2*\$i-1)*\$pi/(2.0*\$ncomps)) ) }]
}

expr {2.0*\$result/\$ncomps}
}

# J1 --
#    First-order Bessel function
#
# Arguments:
#    x         Value of the x-coordinate
# Result:
#    Value of J1(x)
#
proc ::math::special::J1 {x} {
variable pi
#
# Determine the number of components (very crude heuristic)
# (We use a quadrature formula here and use the symmetry of
# the cos function to reduce the number of actually computed
# components by a gfactor 2.
#

set ncomps [expr {2*int(\$x/2.0+0.5)}]
if { \$ncomps < 6 } {
set ncomps 6
}

set result 0.0
for { set i 1 } { \$i < \$ncomps } { incr i 2 } {
set factor [expr {cos((2*\$i-1)*\$pi/(2.0*\$ncomps))}]
set result [expr {\$result + \$factor * sin( \$x * \$factor) }]
}

expr {2.0*\$result/\$ncomps}
}

# J1/2 --
#    Half-order Bessel function
#
# Arguments:
#    x         Value of the x-coordinate
# Result:
#    Value of J1/2(x)
#
proc ::math::special::J1/2 {x} {
variable pi
#
# This Bessel function can be expressed in terms of elementary
# functions. Therefore use the explicit formula
#
if { \$x != 0.0 } {
expr {sqrt(2.0/\$pi/\$x)*sin(\$x)}
} else {
return 0
}
}

# some tests --
#
foreach x {0.0 2.0 4.4 6.0} {
puts "J0(\$x) = [::math::special::J0 \$x] - J1(\$x) = [::math::special::J1 \$x] \
- J1/2(\$x) = [::math::special::J1/2 \$x]"
}```