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 FriCAS already provides a Simpson's method, this version has a syntax that will be intuitive to anyone who has used the integrate() function.

spad
)abbrev package SIMPINT SimpsonIntegration
SimpsonIntegration(): Exports == Implementation where
F        ==> Float
SF       ==> Segment F
EF       ==> Expression F
SBF      ==> SegmentBinding F
Ans      ==> Record(value:F, error:F)
Exports ==> with
simpson : (EF,SBF,F) -> Ans
simpson : (EF,SBF) -> Ans
Implementation ==> add
simpson(func:EF, sbf:SBF, tol:F) ==
a : F := low(segment(sbf))
b : F := high(segment(sbf))
x : EF := variable(sbf) :: EF
h : F
k : Integer
n : Integer
simps : F
newsimps : F
oe : F
ne : F
err : F
sumend : F := retract(eval(func, x, a::EF) + eval(func, x, b::EF))
sumodd : F := 0.0 :: F
sumeven : F := 0.0 :: F
-- First base case -- 2 intervals ----------------
n := 2
h := (b-a)/n
sumeven := sumeven + sumodd
sumodd := 0.0 :: F
for k in 1..(n-1) by 2 repeat
sumodd := sumodd + retract(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 :: F
for k in 1..(n-1) by 2 repeat
sumodd := sumodd + retract(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 :: F
for k in 1..(n-1) by 2 repeat
sumodd := sumodd + retract(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::F)              -- 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::F )
spad
   Compiling FriCAS source code from file
/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/6914547373735257695-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,Float) -> Record(value: Float,error: Float)
Time: 0.66 SEC.
compiling exported simpson : (Expression Float,SegmentBinding Float) -> Record(value: Float,error: Float)
Time: 0.01 SEC.
(time taken in buildFunctor:  0)
;;;     ***       |SimpsonIntegration| REDEFINED
;;;     ***       |SimpsonIntegration| REDEFINED
Time: 0 SEC.
Cumulative Statistics for Constructor SimpsonIntegration
Time: 0.67 seconds
finalizing NRLIB SIMPINT
Processing SimpsonIntegration for Browser database:
--->-->SimpsonIntegration(constructor): Not documented!!!!
--->-->SimpsonIntegration((simpson ((Record (: value (Float)) (: error (Float))) (Expression (Float)) (SegmentBinding (Float)) (Float)))): Not documented!!!!
--->-->SimpsonIntegration((simpson ((Record (: value (Float)) (: error (Float))) (Expression (Float)) (SegmentBinding (Float))))): Not documented!!!!
--->-->SimpsonIntegration(): Missing Description
; compiling file "/var/aw/var/LatexWiki/SIMPINT.NRLIB/SIMPINT.lsp" (written 25 JUN 2019 02:33:51 PM):
; /var/aw/var/LatexWiki/SIMPINT.NRLIB/SIMPINT.fasl written
; compilation finished in 0:00:00.128
------------------------------------------------------------------------
SimpsonIntegration is now explicitly exposed in frame initial
SimpsonIntegration will be automatically loaded when needed from
/var/aw/var/LatexWiki/SIMPINT.NRLIB/SIMPINT

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

fricas
)display op simpson
There are 3 exposed functions called simpson :
[1] (Expression(Float),SegmentBinding(Float),Float) -> Record(value
: Float,error: Float)
from SimpsonIntegration
[2] (Expression(Float),SegmentBinding(Float)) -> Record(value: Float
,error: Float)
from SimpsonIntegration
[3] ((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 )
 (2)
Type: Record(value: Float,error: Float)
fricas
simpson( sin(x), x=0..1, 1.e-10 )
 (3)
Type: Record(value: Float,error: Float)

 Subject:   Be Bold !! ( 14 subscribers )
Please rate this page: