login  home  contents  what's new  discussion  bug reports help  links  subscribe  changes  refresh  edit

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
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
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:

fricas
)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:

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

Our approximations:

fricas
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
or "\$" to specify which version of the function you need.