# ::math::integrate --
#
# calculate the area under a curve defined by a
# set of (x,y) data pairs. the x data must increase
# monotonically throughout the data set for the
# calculation to be meaningful, therefore the
# monotonic condition is tested, and an error is thrown
# if the x value is found to be decreasing.
#
# Arguments:
# (x,y) pairs in list format
#
# Results:
# Area and error bound
# (calculation is "Simpson's rule")
#
# at least 5 data pairs are required, and if the number
# of data pairs is even, a padding value of (x0,0) will
# be added.
#
# 04-13-2000 -- by: Phil Ehrens at The LIGO Lab
namespace eval math {}
proc ::math::integrate { xy_pairs } {
set length [ llength $xy_pairs ]
if { $length < 10 } {
return -code error "at least 5 x,y pairs must be given"
}
;## are we dealing with x,y pairs?
if { [ expr $length % 2 ] } {
return -code error "unmatched xy pair in input"
}
;## are there an even number of pairs? Augment.
if { ! [ expr $length % 4 ] } {
set xy_pairs [ concat [ lindex $xy_pairs 0 ] 0 $xy_pairs ]
}
set x0 [ lindex $xy_pairs 0 ]
set x1 [ lindex $xy_pairs 2 ]
set xn [ lindex $xy_pairs end-1 ]
set xnminus1 [ lindex $xy_pairs end-3 ]
if { $x1 < $x0 } {
return -code error "monotonicity broken by x1"
}
if { $xn < $xnminus1 } {
return -code error "monotonicity broken by xn"
}
;## handle the assymetrical elements 0, n, and n-1.
set sum [ expr {[ lindex $xy_pairs 1 ] + [ lindex $xy_pairs end ]} ]
set sum [ expr {$sum + (4*[ lindex $xy_pairs end-2 ])} ]
set data [ lrange $xy_pairs 2 end-4 ]
set xmax $x1
set i 1
foreach {x1 y1 x2 y2} $data {
incr i
if { $x1 < $xmax } {
return -code error "monotonicity broken by x$i"
}
set xmax $x1
incr i
if { $x2 < $xmax } {
return -code error "monotonicity broken by x$i"
}
set xmax $x2
set sum [ expr {$sum + (4*$y1) + (2*$y2)} ]
}
if { $xmax > $xnminus1 } {
return -code error "monotonicity broken by xn-1"
}
set h [ expr { ( $xn - $x0 ) / $i } ]
set area [ expr { ( $h / 3.0 ) * $sum } ]
set err_bound [ expr { ( ( $xn - $x0 ) / 180.0 ) * pow($h,4) * $xn } ]
return [ list $area $err_bound ]
}Arjen Markus I have a set of procedures ready that will allow one to integrate functions and sets of ordinary differential equations. I just proposed adding this to the Tcllib.
