This routine provides Simpson's method for numerical integration. Although Axiom already provides a Simpson's method, this version has a syntax that will be intuitive to anyone who has used the integrate() function.

)abbrev package SIMPINT SimpsonIntegration
SimpsonIntegration(): Exports == Implementation where
  F        ==> Float
  SF       ==> Segment F
  EF       ==> Expression F
  SBF      ==> SegmentBinding F
  Ans      ==> Record(value:EF, error:EF)

  Exports ==> with
   simpson : (EF,SBF,EF) -> Ans
   simpson : (EF,SBF) -> Ans

  Implementation ==> add
   simpson(func:EF, sbf:SBF, tol:EF) ==
     a : F := lo(segment(sbf))
     b : F := hi(segment(sbf))
     x : EF := variable(sbf) :: EF

     h : F
     k : Integer
     n : Integer

     simps : EF
     newsimps : EF

     oe : EF
     ne : EF
     err : EF

     sumend : EF := eval(func, x, a::EF) + eval(func, x, b::EF)
     sumodd : EF := 0.0 :: EF
     sumeven : EF := 0.0 :: EF

     -- First base case -- 2 intervals ----------------
     n := 2
     h := (b-a)/n
     sumeven := sumeven + sumodd
     sumodd := 0.0 :: EF

     for k in 1..(n-1) by 2 repeat
       sumodd := sumodd + eval( func, x, (k*h+a)::EF )

     simps := ( sumend + 2.0*sumeven + 4.0*sumodd )*(h/3.0)

     -- Second base case -- 4 intervals ---------------
     n := n*2
     h := (b-a)/n
     sumeven := sumeven + sumodd
     sumodd := 0.0 :: EF

     for k in 1..(n-1) by 2 repeat
       sumodd := sumodd + eval( func, x, (k*h+a)::EF )

     newsimps := ( sumend + 2.0*sumeven + 4.0*sumodd )*(h/3.0)

     oe := abs(newsimps-simps)                  -- old error
     simps := newsimps

     -- general case -----------------------------------
     while true repeat
       n := n*2
       h := (b-a)/n

       sumeven := sumeven + sumodd
       sumodd := 0.0 :: EF

       for k in 1..(n-1) by 2 repeat
         sumodd := sumodd + eval( func, x, (k*h+a)::EF )

       newsimps := ( sumend + 2.0*sumeven + 4.0*sumodd )*(h/3.0)

       -- This is a check of Richardson's error estimate.
       -- Usually p is approximately 4 for Simpson's rule, but
       -- occasionally convergence is slower

       ne := abs( newsimps - simps )            -- new error

       if ( (ne<oe*2.0) and (oe<ne*16.5) ) then -- Richardson should be ok
         -- p := log(oe/ne)/log(2.0)
         err := ne/(oe/ne-1.0::EF)              -- ne/(2^p-1)
       else
         err := ne                              -- otherwise estimate crudely

       oe := ne
       simps := newsimps

       if( err < tol ) then
         break

     [ newsimps, err ]

   simpson(func:EF, sbf:SBF) ==
     simpson( func, sbf, 1.e-6::EF )

   Compiling FriCAS source code from file 
      /var/zope2/var/LatexWiki/1602341404191315960-25px001.spad using 
      old system compiler.
   SIMPINT abbreviates package SimpsonIntegration 
   processing macro definition F ==> Float 
   processing macro definition SF ==> Segment Float 
   processing macro definition EF ==> Expression Float 
   processing macro definition SBF ==> SegmentBinding Float 
   processing macro definition Ans ==> Record(value: Expression Float,error: Expression Float) 
   processing macro definition Exports ==> -- the constructor category 
   processing macro definition Implementation ==> -- the constructor capsule 
------------------------------------------------------------------------
   initializing NRLIB SIMPINT for SimpsonIntegration 
   compiling into NRLIB SIMPINT 
   compiling exported simpson : (Expression Float,SegmentBinding Float,Expression Float) -> Record(value: Expression Float,error: Expression Float)
Time: 4.01 SEC.

   compiling exported simpson : (Expression Float,SegmentBinding Float) -> Record(value: Expression Float,error: Expression Float)
Time: 0.03 SEC.

(time taken in buildFunctor:  0)

;;;     ***       |SimpsonIntegration| REDEFINED

;;;     ***       |SimpsonIntegration| REDEFINED
Time: 0.08 SEC.

   Cumulative Statistics for Constructor SimpsonIntegration
      Time: 4.12 seconds

   finalizing NRLIB SIMPINT 
   Processing SimpsonIntegration for Browser database:
--->-->SimpsonIntegration((simpson (Ans EF SBF EF))): Not documented!!!!
--->-->SimpsonIntegration((simpson (Ans EF SBF))): Not documented!!!!
--->-->SimpsonIntegration(constructor): Not documented!!!!
--->-->SimpsonIntegration(): Missing Description
------------------------------------------------------------------------
   SimpsonIntegration is now explicitly exposed in frame initial 
   SimpsonIntegration will be automatically loaded when needed from 
      /var/zope2/var/LatexWiki/SIMPINT.NRLIB/code

This simpson() function overloads the already existing function and either may be used. To see available simpson() functions, do:

)display op simpson

There are 3 exposed functions called simpson :
   [1] (Expression Float,SegmentBinding Float,Expression Float) ->
            Record(value: Expression Float,error: Expression Float)
            from SimpsonIntegration
   [2] (Expression Float,SegmentBinding Float) -> Record(value:
            Expression Float,error: Expression Float)
            from SimpsonIntegration
   [3] ((Float -> Float),Float,Float,Float,Float,Integer,Integer) ->
            Record(value: Float,error: Float,totalpts: Integer,success:
            Boolean)
            from NumericalQuadrature

To compute an integral using Simpson's rule, pass an expression and a BindingSegment? with the limits. Optionally, you may include a third argument to specify the acceptable error.

The exact integral:

integrate( sin(x), x=0..1 ) :: Expression Float

(1)

Our approximations:

simpson( sin(x), x=0..1 )

(2)

simpson( sin(x), x=0..1, 1.e-10 )

(3)