Adaptive Simpson Quadrature

Description

Numerically evaluate an integral using adaptive Simpson's rule.

Usage

1
simpadpt(f, a, b, tol = 1e-6, ...)

Arguments

f

univariate function, the integrand.

a, b

lower limits of integration; must be finite.

tol

relative tolerance

...

additional arguments to be passed to f.

Details

Approximates the integral of the function f from a to b to within an error of tol using recursive adaptive Simpson quadrature.

Value

A numerical value or vector, the computed integral.

Note

Based on code from the book by Quarteroni et al., with some tricks borrowed from Matlab and Octave.

References

Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, Springer-Verlag, Berlin Heidelberg.

See Also

quad, simpson2d

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
myf <- function(x, n) 1/(x+n)  # 0.0953101798043249 , log((n+1)/n) for n=10
simpadpt(myf, 0, 1, n = 10)    # 0.095310179804535

##  Dilogarithm function
flog  <- function(t) log(1-t) / t  # singularity at t=1, almost at t=0
dilog <- function(x) simpadpt(flog, x, 0, tol = 1e-12)
dilog(1)  # 1.64493406685615
          # 1.64493406684823 = pi^2/6

## Not run: 
N <- 51
xs <- seq(-5, 1, length.out = N)
ys <- numeric(N)
for (i in 1:N) ys[i] <- dilog(xs[i])
plot(xs, ys, type = "l", col = "blue",
             main = "Dilogarithm function")
grid()
## End(Not run)

Questions? Problems? Suggestions? or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.