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.

spad
)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 )
spad
   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| (|::| ((|elt| (|Float|) |float|) 0 0 10) (|Expression| (|Float|)))) (REPEAT (STEP |k| 1 2 (- |n| 1)) (LET |sumodd| (+ |sumodd| (|eval| |func| |x| (|::| (+ (* |k| |h|) |a|) (|Expression| (|Float|))))))) (LET |newsimps| (* (+ (+ |sumend| (* ((|elt| (|Float|) |float|) 2 0 10) |sumeven|)) (* ((|elt| (|Float|) |float|) 4 0 10) |sumodd|)) (/ |h| ((|elt| (|Float|) |float|) 3 0 10)))) (LET |ne| (|abs| (- |newsimps| |simps|))) (SEQ (LET #1=#:G691 (< | << ne >> | (* |oe| ((|elt| (|Float|) |float|) 2 0 10)))) (|exit| 1 (IF #1# (SEQ (LET #2=#:G692 (< |oe| (* |ne| ((|elt| (|Float|) |float|) 165 -1 10)))) (|exit| 1 (IF #2# (LET |err| (/ |ne| (- (/ |oe| |ne|) (|::| ((|elt| (|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:= ((((|ne| # #) (|newsimps| # # #) (|sumodd| # # # ...) (|k| # #) ...)))
>> 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] ((Float -> Float),Float,Float,Float,Float,Integer,Integer) -> Record(value: Float,error: Float,totalpts: Integer,success: Boolean) from NumericalQuadrature

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_318602826(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 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. simpson( sin(x), x=0..1, 1.e-10 )
There are no library operations named simpson having 3 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) Float
Perhaps you should use "@" to indicate the required return type, or "$" to specify which version of the function you need.




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