|
|
|
last edited 6 years ago by test1 |
Edit detail for Simpson's method revision 4 of 5
| 1 2 3 4 5 | ||
|
Editor: test1
Time: 2019/06/25 14:12:30 GMT+0 |
||
| Note: | ||
changed: -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. 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.
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: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,
(SEQ (|:=| (|:| |a| (|Float|)) | << | (|lo| (|segment| |sbf|)) | >> |)
(|:=| (|:| |b| (|Float|)) (|hi| (|segment| |sbf|)))
(|:=| (|:| |x| (|Expression| (|Float|)))
(|::| (|variable| |sbf|) (|Expression| (|Float|))))
(|:| |h| (|Float|)) (|:| |k| (|Integer|)) (|:| |n| (|Integer|))
(|:| |simps| (|Expression| (|Float|)))
(|:| |newsimps| (|Expression| (|Float|)))
(|:| |oe| (|Expression| (|Float|))) (|:| |ne| (|Expression| (|Float|)))
(|:| |err| (|Expression| (|Float|)))
(|:=| (|:| |sumend| (|Expression| (|Float|)))
(+ (|eval| |func| |x| (|::| |a| (|Expression| (|Float|))))
(|eval| |func| |x| (|::| |b| (|Expression| (|Float|))))))
(|:=| (|:| |sumodd| (|Expression| (|Float|)))
(|::| ((|Sel| (|Float|) |float|) 0 0 10) (|Expression| (|Float|))))
(|:=| (|:| |sumeven| (|Expression| (|Float|)))
(|::| ((|Sel| (|Float|) |float|) 0 0 10) (|Expression| (|Float|))))
(|:=| |n| 2) (|:=| |h| (/ (- |b| |a|) |n|))
(|:=| |sumeven| (+ |sumeven| |sumodd|))
(|:=| |sumodd|
(|::| ((|Sel| (|Float|) |float|) 0 0 10) (|Expression| (|Float|))))
(REPEAT (INBY |k| (SEGMENT 1 (- |n| 1)) 2)
(|:=| |sumodd|
(+ |sumodd|
(|eval| |func| |x|
(|::| (+ (* |k| |h|) |a|)
(|Expression| (|Float|)))))))
(|:=| |simps|
(*
(+ (+ |sumend| (* ((|Sel| (|Float|) |float|) 2 0 10) |sumeven|))
(* ((|Sel| (|Float|) |float|) 4 0 10) |sumodd|))
(/ |h| ((|Sel| (|Float|) |float|) 3 0 10))))
(|:=| |n| (* |n| 2)) (|:=| |h| (/ (- |b| |a|) |n|))
(|:=| |sumeven| (+ |sumeven| |sumodd|))
(|:=| |sumodd|
(|::| ((|Sel| (|Float|) |float|) 0 0 10) (|Expression| (|Float|))))
(REPEAT (INBY |k| (SEGMENT 1 (- |n| 1)) 2)
(|:=| |sumodd|
(+ |sumodd|
(|eval| |func| |x|
(|::| (+ (* |k| |h|) |a|)
(|Expression| (|Float|)))))))
(|:=| |newsimps|
(*
(+ (+ |sumend| (* ((|Sel| (|Float|) |float|) 2 0 10) |sumeven|))
(* ((|Sel| (|Float|) |float|) 4 0 10) |sumodd|))
(/ |h| ((|Sel| (|Float|) |float|) 3 0 10))))
(|:=| |oe| (|abs| (- |newsimps| |simps|))) (|:=| |simps| |newsimps|)
(REPEAT (WHILE |true|)
(SEQ (|:=| |n| (* |n| 2)) (|:=| |h| (/ (- |b| |a|) |n|))
(|:=| |sumeven| (+ |sumeven| |sumodd|))
(|:=| |sumodd|
(|::| ((|Sel| (|Float|) |float|) 0 0 10)
(|Expression| (|Float|))))
(REPEAT (INBY |k| (SEGMENT 1 (- |n| 1)) 2)
(|:=| |sumodd|
(+ |sumodd|
(|eval| |func| |x|
(|::| (+ (* |k| |h|) |a|)
(|Expression| (|Float|)))))))
(|:=| |newsimps|
(*
(+
(+ |sumend|
(* ((|Sel| (|Float|) |float|) 2 0 10) |sumeven|))
(* ((|Sel| (|Float|) |float|) 4 0 10) |sumodd|))
(/ |h| ((|Sel| (|Float|) |float|) 3 0 10))))
(|:=| |ne| (|abs| (- |newsimps| |simps|)))
(SEQ
(|:=| (|:| #1=#:G671 (|Boolean|))
(< |ne| (* |oe| ((|Sel| (|Float|) |float|) 2 0 10))))
(|exit| 1
(IF #1#
(SEQ
(|:=| (|:| #2=#:G672 (|Boolean|))
(< |oe|
(* |ne| ((|Sel| (|Float|) |float|) 165 -1 10))))
(|exit| 1
(IF #2#
(|:=| |err|
(/ |ne|
(- (/ |oe| |ne|)
(|::| ((|Sel| (|Float|) |float|) 1 0 10)
(|Expression| (|Float|))))))
(|:=| |err| |ne|))))
(|:=| |err| |ne|))))
(|:=| |oe| |ne|) (|:=| |simps| |newsimps|)
(|exit| 1
(IF (< |err| |tol|)
(|leave| 1 |$NoValue|)
|noBranch|))))
(|exit| 1 (|construct| |newsimps| |err|)))
****** level 3 ******
$x:= (lo (segment sbf))
$m:= $EmptyMode
$f:=
((((|a| #) (|tol| # #) (|sbf| # #) (|func| # #) ...)))
>> Apparent user error:
NoValueMode
is an unknown modeThis 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 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.