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.

fricas
(1) -> <spad>
fricas
)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>
fricas
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: 1.47 SEC.
compiling exported simpson : (Expression Float,SegmentBinding Float) -> Record(value: Float,error: Float) Time: 0 SEC.
(time taken in buildFunctor: 128)
;;; *** |SimpsonIntegration| REDEFINED
;;; *** |SimpsonIntegration| REDEFINED Time: 0 SEC.
Cumulative Statistics for Constructor SimpsonIntegration Time: 1.47 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 08 SEP 2024 01:02:32 AM):
; wrote /var/aw/var/LatexWiki/SIMPINT.NRLIB/SIMPINT.fasl ; compilation finished in 0:00:00.164 ------------------------------------------------------------------------ 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

\label{eq1}0.4596976941 \<u> 318602826(1)
Type: Expression(Float)

Our approximations:

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

\label{eq2}\begin{array}{@{}l}
\displaystyle
\left[{value ={0.4596983187 \_ 9846143679}}, \: \right.
\
\
\displaystyle
\left.{error ={0.6126922587 \_ 0430639838 E - 6}}\right] 
(2)
Type: Record(value: Float,error: Float)
fricas
simpson( sin(x), x=0..1, 1.e-10 )

\label{eq3}\begin{array}{@{}l}
\displaystyle
\left[{value ={0.4596976941 \_ 4137428151}}, \: \right.
\
\
\displaystyle
\left.{error ={0.9513291004 \_ 0963440174 E - 11}}\right] 
(3)
Type: Record(value: Float,error: Float)




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