JPhrases/SimpsonIntegraton

Here's a little function to do numeric integration based on Simpson's method. I found this in a 1995 article about J in Byte magazine online; the article is here. I've updated the syntax slightly to reflect a simplification in J since that time.

NB.* Simpson.ijs: numeric integration using Simpson's method: from
NB. http://www.byte.com/art/9509/img/505079a2.htm (Byte magazine article, 9/1995)

NB. form: verb simpson int
NB.   verb is the monadic function to be integrated.
NB.   int has 2 or 3 elements:
NB.     [0] lower bound of interval
NB.     [1] upper bound of interval
NB.     [2] number of subintervals (default 128)
NB. result is the integral
NB. e.g. 43.75 = ^&3 simpson ] 3 4
Simpson=: 1 : 0
   'lower upper int'=. 3{.y,128
   size=. int%~upper-lower
   val=. u lower+size*i.>:int
   size*+/val*3%~1,1,~4 2$~<:int
)

Some Examples of Use

First, a few that are easy to check.

   ] Simpson 0 1    NB. straight line: x=y (half the unit square)
0.5
   1: Simpson 0 1   NB. straight line: y=1 (unit square)
1
   *: Simpson 0 1   NB. parabola from 0 to 1
0.33333333

The integral of x^2 is 3%~x^3, so the one-third above is correct. Here's a picture so we can see that this also looks about right.

Some slightly more complex examples:

   1&o. Simpson o.0 1  NB. Sine from zero to pi
2

This also looks plausible.

Another one that's easy to check due to symmetry of sine (the positive wave up in the first half of the cycle should be the same as the negative wave down in the second half):

   1&o. Simpson o.0 2  NB. Sine from zero to 2*pi
_1.7234975e_16
   NB. i.e. equals zero ex numerical imprecision of floating point

We see that the imprecision grows for an increasingly wide range but it is probably acceptable for most practical applications. Here we plot what should be all zeroes as we integrate from plus or minus pi to plus or minus 1000 pi and it's off by a fraction of a billionth at the extreme, without increasing the default resolution of only 128 points (segments of integration).

   plot (1&o.) Simpson"(1) (o._1 1)*"(1 0)>:i.1000
   pd 'save png c:\amisc\j\simpsonIncreasingInaccuracy.png'
19778

One last one using sine. I added one to it because the area plot doesn't handle negative areas.

   ([: >: 1 o. ]) Simpson o. 0 2   NB. 1+sine from zero to 2*pi
6.2831853
   o.2                             NB. Makes sense: takes up half the 2 by 2 pi rectangle.
6.2831853
   'type area' plot (o. 0 2);'[: >: 1 o. ]'
   pd 'save png c:\amisc\j\plot1+Sin02pi.png'
14219