## integrate

```# ::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
#
# 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.

 Category Mathematics