[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]" } ====== <> Mathematics | Numerical Analysis