-- Simpson integration for discrete values
-- by Georg Gollmann <gollmann@edvz.tuwien.ac.at>

Simpson(Y,h) = sum(Y[i] * simpCoeff(h,i),i,1,count(Y)) /
                                   3 when count(h) > 0,
  h/3 * (Y[1] + Y[count(Y)] +
  4 * sum(Y[2*i], i, 1, (count(Y)-1)/2) +
  2 * sum(Y[2*i+1], i, 1, (count(Y)-3)/2))

simpCoeff(x,i) = (2 when i/2 = trunc(i/2), 1 otherwise) *
 (x[i+1] - x[i-1], -- general case
  x[i+1] - x[i], -- first section
  x[i] - x[i-1]) -- last section


-- Examples
-- Curve defined by equidistant points
Yh = {0, 1.3, 2, 2.3, 2.5, 2.4, 2.2, 1.8, 0}
h = 1.5
Simpson(Yh,h):22.300

-- Same curve with unequal spacing of points
x = {0, 1.5, 3.5,   6, 8.5, 10.5, 12}
y = {0, 1.3, 2.2, 2.5, 2.3,  1.8,  0}
Simpson(y, x):22.317