Simpson's 1/3 rule on the 50g

[Wes Loewer]

I dug through my old files an found a program that I wrote some 15+ years ago. Looking at it now, I might have done this differently, but here it is, warts and all.

To understand the code fully, look up the syntax of the Σ and the | (where) commands.

Code:


@ SIMPSONS(xleft,xright,fnc,xvar,steps)
@ calculates area via Simpson's Rule (using parabolas)
@ Wes Loewer
@
@ parameters:
@ xleft   - left bound
@ xright  - right bound
@ fnc     - function
@ xvar   - independent variable
@ steps   - number of steps, must be even (ie, mult. of 2)

« → xleft xright fnc xvar steps
  «
    PUSH                      @ save flags
    -105. SF                  @ approx mode ON
    xleft →NUM 'xleft' STO    @ convert to decimals
    xright →NUM 'xright' STO
                              @ make sure steps is an even integer
    steps →NUM IP 2. / CEIL 2. * 'steps' STO

    xright xleft - steps / 
    → deltax                  @ step width
    «
      @ Simpson's Rule is
      @ (d/3)*( f(x0) + 4(f(x1) + 2f(x2) + 4f(x3) + 2f(x4) + ...) + f(xn) )

      @ middle odd terms * 4
      'n' 1. steps 2. /       @ 1st 3 arg of Σ_{n=1}^{steps/2}
                              @ f(xleft + (2n-1)*deltax)
      fnc xvar xleft '2.*n-1.' deltax * + 2. →LIST |
      Σ 4. *                  @ sum up middle odd terms *4

      @ middle even terms * 2
      'n' 1. steps 2. / 1. -  @ 1st 3 arg of Σ_{n=1}^{steps/2-1}
                              @ f(xleft + (2n)*deltax)
      fnc xvar xleft '2.*n' deltax * + 2. →LIST |
      @ Σ correctly handles no terms, Σ_{n=1}^{0}
      Σ 2. *                  @ sum up middle even terms *2
      +

      @ end terms
      fnc xvar xleft  2. →LIST | →NUM @ f(xleft)
      +
      fnc xvar xright 2. →LIST | →NUM @ f(xright)
      +

      deltax * 3. /           @ * deltax / 3 to get area
    »
    POP                       @ restore flags
  »
»

Feel free to ask if there is anything that is not clear.