While working on this Rosetta Code task: "http://rosettacode.org/wiki/Numerical_Integration"; I came across the following question regarding the base library implementation of Simpson's Rule.
I compared 3 J implementations of Simpson's rule: * simpsonBL adapted from the integrat.ijs script in the Base Library (now migrated to JAL as 'math/misc/integrat') * simpsonJR by John Randall (JR) in the forum "http://www.jsoftware.com/pipermail/programming/2008-August/011747.html"; * simpsonWK implemented as per the Wikipedia (WK) page for Simpson's Rule The issue is I need to specify twice the number of steps for simpsonBL and simpsonJR in order to get the same answer as for simpsonWK. E.g. integrating sin function from 0 to 1p1 (actual answer is 2): sin simpsonBL 0 1p1 10 2.00010951732 sin simpsonJR 0 1p1 10 2.00010951732 sin simpsonWK 0 1p1 10 2.00000678444 sin simpsonBL 0 1p1 20 2.00000678444 sin simpsonJR 0 1p1 20 2.00000678444 This makes sense as simpsonWK calculates the function of interest for the mid-point of each sub-interval (similar to doubling steps) whereas the others don't. I'm wondering whether the simpson conjunction in 'math/misc/integrat' is "correct" as is or whether it should be changed to include "steps=. +:steps" so that it gives the same answer as simpsonWK given the same inputs? As a side note, applying the same change to simpsonJR solves the problem it currently has with odd numbers of steps. I've adapted them below to use similar formats and variable names to help illustrate the differences. simpsonBL=: 1 : 0 NB. adapted from ~addons/math/misc/integrat.ijs 'a b steps'=. 3{.y,128 size=. (b-a) % steps vals=. u a + size * i.>:steps size * +/vals * 3%~ 1,((<:steps)$4 2),1 ) simpsonJR=: 1 : 0 NB. adapted from John Randall's post 'a b steps'=. 3{.y,128 size=. (b-a) % steps coeffs=. 1 4 1 +//.@(*/) 0=2 | i.<:steps (size%3) * coeffs +/ .* u size * i.>:steps ) simpsonWK=: 1 : 0 NB. adapted from Wikipedia formula (less efficient) 'a b steps'=. 3{.y,128 size=. (b-a) % steps pts=. a + size * i.>:steps size * +/ 2 (6 %~ [: +/ 4&*...@u@-:@(+/) , u)\ pts ) ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
