Edit detail for Simpson's method revision 3 of 5

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/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/1602341404191315960-25px001.spad
      using old system compiler.
   SIMPINT abbreviates package SimpsonIntegration 
------------------------------------------------------------------------
   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)
****** comp fails at level 5 with expression: ******
error in function simpson 

(SEQ
 (LET |n|
   (* |n| 2))
 (LET |h|
   (/ (- |b| |a|) |n|))
 (LET |sumeven|
   (+ |sumeven| |sumodd|))
 (LET |sumodd|
   (|::| ((|Sel| (|Float|) |float|) 0 0 10) (|Expression| (|Float|))))
 (REPEAT (INBY |k| (SEGMENT 1 (- |n| 1)) 2)
         (LET |sumodd|
           (+ |sumodd|
              (|eval| |func| |x|
                      (|::| (+ (* |k| |h|) |a|) (|Expression| (|Float|)))))))
 (LET |newsimps|
   (*
    (+ (+ |sumend| (* ((|Sel| (|Float|) |float|) 2 0 10) |sumeven|))
       (* ((|Sel| (|Float|) |float|) 4 0 10) |sumodd|))
    (/ |h| ((|Sel| (|Float|) |float|) 3 0 10))))
 (LET |ne|
   (|abs| (- |newsimps| |simps|)))
 (SEQ
  (LET (|:| #1=#:G671 (|Boolean|))
    (< | << ne >> | (* |oe| ((|Sel| (|Float|) |float|) 2 0 10))))
  (|exit| 1
   (IF #1#
       (SEQ
        (LET (|:| #2=#:G672 (|Boolean|))
          (< |oe| (* |ne| ((|Sel| (|Float|) |float|) 165 -1 10))))
        (|exit| 1
         (IF #2#
             (LET |err|
               (/ |ne|
                  (- (/ |oe| |ne|)
                     (|::| ((|Sel| (|Float|) |float|) 1 0 10)
                      (|Expression| (|Float|))))))
             (LET |err|
               |ne|))))
       (LET |err|
         |ne|))))
 (LET |oe|
   |ne|)
 (LET |simps|
   |newsimps|)
 (|exit| 1
  (IF (< |err| |tol|)
      (|leave| 1 |$NoValue|)
      |noBranch|)))
****** level 5  ******
$x:= ne
$m:= (Float)
$f:=
((((#:G671 #) (|ne| # #) (|newsimps| # # #) (|sumodd| # # # ...) ...)))

   >> Apparent user error:
   Cannot coerce ne 
      of mode (Expression (Float)) 
      to mode (Boolean)

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

)display op simpson

There is one exposed function called simpson :
   [1] ((D3 -> D3),D3,D3,D3,D3,Integer,Integer) -> Record(value: D3,
            error: D3,totalpts: Integer,success: Boolean)
            from NumericalQuadrature(D3) if D3 has FPS

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 )

   There are no library operations named simpson having 2 argument(s) 
      though there are 1 exposed operation(s) and 0 unexposed 
      operation(s) having a different number of arguments. Use HyperDoc
      Browse, or issue
                              )what op simpson
      to learn what operations contain " simpson " in their names, or 
      issue
                             )display op simpson
      to learn more about the available operations.

   Cannot find a definition or applicable library operation named 
      simpson with argument type(s) 
                             Expression(Integer)
                     SegmentBinding(NonNegativeInteger)

      Perhaps you should use "@" to indicate the required return type, 
      or "$" to specify which version of the function you need.