R Under development (unstable) (2023-11-08 r85496 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ## > ## q u a d r a t u r e . R Test suite > ## > > > quad <- pracma::quad > quadl <- pracma::quadl > quadgk <- pracma::quadgk > quadgr <- pracma::quadgr > quadinf <- pracma::quadinf > quad2d <- pracma::quad2d > dblquad <- pracma::dblquad > > simpson2d <- pracma::simpson2d > simpadpt <- pracma::simpadpt > gauss_kronrod <- pracma::gauss_kronrod > clenshaw_curtis <- pracma::clenshaw_curtis > romberg <- pracma::romberg > > gaussLegendre <- pracma::gaussLegendre > gaussHermite <- pracma::gaussHermite > gaussLaguerre <- pracma::gaussLaguerre > > ## F i n i t e I n t e r v a l s > > f1 <- function(x) exp(x)*sin(x) # [0, pi] 12.0703463163896 = 1/2*(1+e^pi) > f2 <- pracma::runge # [-1, 1] 0.549360306778006 > f3 <- function(x) 1/(x^3 - 2*x - 5) # [0, 2] -0.460501533846733 > f4 <- function(x) abs(sin(10*x)) # [0, pi] 2.0 > > # quad (Adaptive Simpson) > all.equal(quad(f1, 0, pi, tol=1e-12), 12.0703463163896, + tolerance = 1e-12) [1] TRUE > all.equal(quad(f2, -1, 1, tol=1e-12), 0.549360306778006, + tolerance = 1e-12) [1] TRUE > all.equal(quad(f3, 0, 2, tol=1e-12), -0.460501533846733, + tolerance = 1e-12) [1] TRUE > all.equal(quad(f4, 0, pi, tol=1e-12), 2.0, + tolerance = 1e-12) [1] TRUE > > # quadl (Adaptive Lobatto) > all.equal(quadl(f1, 0, pi, tol=1e-9), 12.0703463163896, + tolerance = 1e-12) [1] TRUE > all.equal(quadl(f2, -1, 1, tol=1e-9), 0.549360306778006, + tolerance = 1e-12) [1] TRUE > all.equal(quadl(f3, 0, 2, tol=1e-9), -0.460501533846733, + tolerance = 1e-12) [1] TRUE > all.equal(quadl(f4, 0, pi, tol=1e-12), 2.0, + tolerance = 1e-12) [1] TRUE > > # quadgr (Gauss-Richardson) > all.equal(quadgr(f1, 0, pi, tol=1e-12)$value, 12.0703463163896, + tolerance = 1e-13) [1] TRUE > all.equal(quadgr(f2, -1, 1, tol=1e-12)$value, 0.549360306778006, + tolerance = 1e-15) [1] TRUE > all.equal(quadgr(f3, 0, 2, tol=1e-12)$value, -0.460501533846733, + tolerance = 1e-15) [1] TRUE > all.equal(quadgr(f4, 0, pi, tol=1e-12)$value, 2.0, + tolerance = 1e-15) [1] TRUE > > # quadgk (Adaptive Gauss-Kronrod) > all.equal(quadgk(f1, 0, pi), 12.0703463163896, + tolerance = 1e-13) [1] TRUE > all.equal(quadgk(f2, -1, 1), 0.549360306778006, + tolerance = 1e-15) [1] TRUE > all.equal(quadgk(f3, 0, 2), -0.460501533846733, + tolerance = 1e-13) [1] TRUE > all.equal(quadgk(f4, 0, pi, tol = 1e-12), 2.0, + tolerance = 1e-12) [1] TRUE > > # Adaptive Simpson (simpadpt) > all.equal(simpadpt(f1, 0, pi, tol=1e-12), 12.0703463163896, + tolerance = 1e-13) [1] TRUE > all.equal(simpadpt(f2, -1, 1, tol=1e-12), 0.549360306778006, + tolerance = 1e-12) [1] TRUE > all.equal(simpadpt(f3, 0, 2, tol=1e-12), -0.460501533846733, + tolerance = 1e-13) [1] TRUE > all.equal(simpadpt(f4, 0, pi, tol=1e-12), 2.0, + tolerance = 1e-14) [1] TRUE > > # Gauss-Kronrod > all.equal(gauss_kronrod(f1, 0, pi)$value, 12.0703463163896, + tolerance = 1e-13) [1] TRUE > all.equal(gauss_kronrod(f2, -1, 1)$value, 0.549360306778006, # BAD + tolerance = 1e-2) [1] TRUE > all.equal(gauss_kronrod(f3, 0, 2)$value, -0.460501533846733, # Bad + tolerance = 1e-5) [1] TRUE > all.equal(gauss_kronrod(f4, 0, pi)$value, 2.0, # BAD + tolerance = 1e-0) [1] TRUE > > # Clenshaw-Curtis > all.equal(clenshaw_curtis(f1, 0, pi, n = 128), 12.0703463163896, + tolerance = 1e-12) [1] TRUE > all.equal(clenshaw_curtis(f2, -1, 1, n = 128), 0.549360306778006, + tolerance = 1e-12) [1] TRUE > all.equal(clenshaw_curtis(f3, 0, 2, n = 128), -0.460501533846733, + tolerance = 1e-12) [1] TRUE > all.equal(clenshaw_curtis(f4, 0, pi, n = 1024), 2.0, # Bad + tolerance = 2e-5) [1] TRUE > > # romberg > all.equal(romberg(f1, 0, pi, tol=1e-12)$value, 12.0703463163896, + tolerance = 1e-12) [1] TRUE > all.equal(romberg(f2, -1, 1, tol=1e-12)$value, 0.549360306778006, + tolerance = 1e-12) [1] "Mean relative difference: 1.660053e-10" > all.equal(romberg(f3, 0, 2, tol=1e-12)$value, -0.460501533846733, # BAD + tolerance = 1e-3) [1] TRUE > all.equal(romberg(f4, 0, pi, tol=1e-12)$value, 2.0, + tolerance = 1e-12) [1] "Mean absolute difference: 2" > > f5 <- function(x) log(x)*sin(x)/x # pi/2 * gamma , cannot be computed ! > f6 <- function(x) sin(x)^2 * exp(-x) # [0, Inf] , 0.4 > f7 <- function(x) sin(x)^2 * exp(-x^2) # [-Inf, Inf] , (e-1)*sqrt(pi)/(4*e) > x7 <- (exp(1)-1) * sqrt(pi) / (2*exp(1)) > > # quadinf > all.equal(quadinf(f6, 0, Inf), 0.4, tolerance = 1e-15) [1] "Modes: list, numeric" [2] "names for target but not for current" [3] "Length mismatch: comparison on first 1 components" > all.equal(quadinf(f7, -Inf, Inf), x7, tolerance = 1e-15) [1] "Modes: list, numeric" [2] "names for target but not for current" [3] "Length mismatch: comparison on first 1 components" [4] "Component 1: Mean relative difference: 5.529293e-14" > > all.equal(quadgr(f6, 0, Inf)$value, 0.4, tolerance = 1e-11) [1] TRUE Warning messages: 1: In sin(x) : NaNs produced 2: In sin(x) : NaNs produced 3: In sin(x) : NaNs produced 4: In sin(x) : NaNs produced > all.equal(quadgr(f7, -Inf, Inf)$value, x7, tolerance = 1e-9) [1] TRUE Warning messages: 1: In sin(x) : NaNs produced 2: In sin(x) : NaNs produced 3: In sin(x) : NaNs produced 4: In sin(x) : NaNs produced 5: In sin(x) : NaNs produced 6: In sin(x) : NaNs produced 7: In sin(x) : NaNs produced > > gL <- gaussLaguerre(64) > all.equal(sum(gL$w * sin(gL$x)^2), 0.4, tolerance = 1e-15) [1] TRUE > gH <- gaussHermite(64) > all.equal(sum(gH$w * sin(gH$x)^2), x7, tolerance = 1e-14) [1] TRUE > > f8 <- function(x, y) y * sin(x) # [0, pi/2]x[0, 1] , 1/2 > f9 <- function(x, y) ifelse(x^2 + y^2 <= 1, 1-x^2-y^2, 0) > > # quad2d > all.equal(quad2d(f8, 0, pi/2, 0, 1), 0.5, tolerance = 1e-15) [1] TRUE > all.equal(quad2d(f9, -1, 1, 0, 1, n = 128), pi/4, tolerance = 1e-6) [1] TRUE > > # dblquad > all.equal(dblquad(f8, 0, pi/2, 0, 1), 0.5, tolerance = 1e-15) [1] TRUE > #all.equal(dblquad(f9, -1, 1, 0, 1), pi/4, tolerance = 1e-6) > # disabled because of problems with Fedora and Solaris > > # simpson2d > all.equal(simpson2d(f8, 0, pi/2, 0, 1), 0.5, tolerance = 1e-9) [1] TRUE > all.equal(simpson2d(f9, -1, 1, 0, 1), pi/4, tolerance = 1e-5) [1] TRUE > > # Integrals with singularities at boundaries: > f11 <- function(t) log(1-t) / t # [1, 0] pi^2/6 , dilogarithm > f12 <- function(t) log(-log(t)) # [0, 1] gamma = 0.57721 56649 01532 ... > f13 <- function(t) 1 / sqrt(t) # [0, 1] 2.0 > > all.equal(quad(f11, 1, 0, tol = 1e-12), 1.64493406684823, + tolerance = 1e-10) [1] TRUE > all.equal(quad(f12, 0, 1, tol = 1e-12), -0.577215664901533, + tolerance = 5e-10) [1] TRUE > all.equal(quad(f13, 0, 1, tol = 1e-12), 2.0, + tolerance = 1e-4) # Bad [1] TRUE > > all.equal(quadl(f11, 1, 0, tol = 1e-12), 1.64493406684823, + tolerance = 1e-12) [1] TRUE > all.equal(quadl(f12, 0, 1, tol = 1e-12), -0.577215664901533, + tolerance = 5e-12) [1] TRUE > all.equal(quadl(f13, 0, 1, tol = 1e-12), 2.0, + tolerance = 1e-7) # Bad [1] TRUE > > all.equal(quadgr(f11, 1, 0, tol = 1e-12)$value, 1.64493406684823, + tolerance = 1e-12) [1] TRUE > all.equal(quadgr(f12, 0, 1, tol = 1e-12)$value, -0.577215664901533, + tolerance = 5e-12) [1] TRUE > all.equal(quadgr(f13, 0, 1, tol = 1e-12)$value, 2.0, + tolerance = 1e-12) [1] TRUE > > all.equal(simpadpt(f11, 1, 0, tol = 1e-12), 1.64493406684823, + tolerance = 1e-11) [1] TRUE > all.equal(simpadpt(f12, 0, 1, tol = 1e-12), -0.577215664901533, + tolerance = 5e-11) [1] TRUE > all.equal(simpadpt(f13, 0, 1, tol = 1e-10), 2.0, + tolerance = 1e-7) # Bad [1] TRUE > > ## E o F > > proc.time() user system elapsed 0.35 0.06 0.43