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

# Edit detail for Simpson's method revision 5 of 5

 1 2 3 4 5 Editor: test1 Time: 2019/06/25 14:33:52 GMT+0 Note:

changed:
-  Ans      ==> Record(value:EF, error:EF)
Ans      ==> Record(value:F, error:F)

changed:
-   simpson : (EF,SBF,EF) -> Ans
simpson : (EF,SBF,F) -> Ans

changed:
-   simpson(func:EF, sbf:SBF, tol:EF) ==
-     a : F := lo(segment(sbf))
-     b : F := hi(segment(sbf))
simpson(func:EF, sbf:SBF, tol:F) ==
a : F := low(segment(sbf))
b : F := high(segment(sbf))

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

changed:
-     sumodd := 0.0 :: EF
sumodd := 0.0 :: F

changed:
-       sumodd := sumodd + eval( func, x, (k*h+a)::EF )
sumodd := sumodd + retract(eval( func, x, (k*h+a)::EF ))

changed:
-     sumodd := 0.0 :: EF
sumodd := 0.0 :: F

changed:
-       sumodd := sumodd + eval( func, x, (k*h+a)::EF )
sumodd := sumodd + retract(eval( func, x, (k*h+a)::EF ))

changed:
-       sumodd := 0.0 :: EF
sumodd := 0.0 :: F

changed:
-         sumodd := sumodd + eval( func, x, (k*h+a)::EF )
sumodd := sumodd + retract(eval( func, x, (k*h+a)::EF ))

changed:
-         err := ne/(oe/ne-1.0::EF)              -- ne/(2^p-1)
err := ne/(oe/ne-1.0::F)              -- ne/(2^p-1)

changed:
-     simpson( func, sbf, 1.e-6::EF )
simpson( func, sbf, 1.e-6::F )


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.

)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
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 )
   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,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)