Changes http://wiki.axiom-developer.org/SimpsonsMethod/diff
--
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.
\begin{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 )
\end{spad}
This simpson() function overloads the already existing function and either may
be used. To see available simpson() functions, do:
\begin{axiom}
)display op simpson
\end{axiom}
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.
\begin{axiom}
integrate( sin(x), x=0..1 ) :: Expression Float
simpson( sin(x), x=0..1 )
simpson( sin(x), x=0..1, 1.e-10 )
\end{axiom}
--
forwarded from http://wiki.axiom-developer.org/[EMAIL PROTECTED]
