JPhrases/SimpsonIntegration
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. 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
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