Artifact 30a503c2f0968d25a2836d0c839aef7ecac24709:
# bessel.tcl --
# Evaluate the most common Bessel functions
#
# TODO:
# Yn - finding decent approximations seems tough
# Jnu - for arbitrary values of the parameter
# J'n - first derivative (from recurrence relation)
# Kn - forward application of recurrence relation?
#
# namespace special
# Create a convenient namespace for the "special" mathematical functions
#
namespace eval ::math::special {
#
# Define a number of common mathematical constants
#
::math::constants::constants pi
#
# Export the functions
#
namespace export J0 J1 Jn J1/2 J-1/2 I_n
}
# J0 --
# Zeroth-order Bessel function
#
# Arguments:
# x Value of the x-coordinate
# Result:
# Value of J0(x)
#
proc ::math::special::J0 {x} {
Jn 0 $x
}
# J1 --
# First-order Bessel function
#
# Arguments:
# x Value of the x-coordinate
# Result:
# Value of J1(x)
#
proc ::math::special::J1 {x} {
Jn 1 $x
}
# Jn --
# Compute the Bessel function of the first kind of order n
# Arguments:
# n Order of the function (must be integer)
# x Value of the argument
# Result:
# Jn(x)
# Note:
# This relies on the integral representation for
# the Bessel functions of integer order:
# 1 I pi
# Jn(x) = -- I cos(x sin t - nt) dt
# pi 0 I
#
# For this kind of integrands the trapezoidal rule is
# very efficient according to Davis and Rabinowitz
# (Methods of numerical integration, 1984).
#
proc ::math::special::Jn {n x} {
variable pi
if { ![string is integer -strict $n] } {
return -code error "Order argument must be integer"
}
#
# Integrate over the interval [0,pi] using a small
# enough step - 40 points should do a good job
# with |x| < 20, n < 20 (an accuracy of 1.0e-8
# is reported by Davis and Rabinowitz)
#
set number 40
set step [expr {$pi/double($number)}]
set result 0.0
for { set i 0 } { $i <= $number } { incr i } {
set t [expr {double($i)*$step}]
set f [expr {cos($x * sin($t) - $n * $t)}]
if { $i == 0 || $i == $number } {
set f [expr {$f/2.0}]
}
set result [expr {$result+$f}]
}
expr {$result*$step/$pi}
}
# 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.0
}
}
# J-1/2 --
# Compute the Bessel function of the first kind of order -1/2
# Arguments:
# x Value of the argument (!= 0.0)
# Result:
# J-1/2(x)
#
proc ::math::special::J-1/2 {x} {
variable pi
if { $x == 0.0 } {
return -code error "Argument must not be zero (singularity)"
} else {
return [expr {-cos($x)/sqrt($pi*$x/2.0)}]
}
}
# I_n --
# Compute the modified Bessel function of the first kind
#
# Arguments:
# n Order of the function (must be positive integer or zero)
# x Abscissa at which to compute it
# Result:
# Value of In(x)
# Note:
# This relies on Miller's algorithm for finding minimal solutions
#
namespace eval ::math::special {}
proc ::math::special::I_n {n x} {
if { ! [string is integer $n] || $n < 0 } {
error "Wrong order: must be positive integer or zero"
}
set n2 [expr {$n+8}] ;# Note: just a guess that this will be enough
set ynp1 0.0
set yn 1.0
set sum 1.0
while { $n2 > 0 } {
set ynm1 [expr {$ynp1+2.0*$n2*$yn/$x}]
set sum [expr {$sum+$ynm1}]
if { $n2 == $n+1 } {
set result $ynm1
}
set ynp1 $yn
set yn $ynm1
incr n2 -1
}
set quotient [expr {(2.0*$sum-$ynm1)/exp($x)}]
expr {$result/$quotient}
}
#
# some tests --
#
if { 0 } {
set tcl_precision 17
foreach x {0.0 2.0 4.4 6.0 10.0 11.0 12.0 13.0 14.0} {
puts "J0($x) = [::math::special::J0 $x] - J1($x) = [::math::special::J1 $x] \
- J1/2($x) = [::math::special::J1/2 $x]"
}
foreach n {0 1 2 3 4 5} {
puts [::math::special::I_n $n 1.0]
}
}