https://github.com/cran/pracma
Raw File
Tip revision: a7001ff1805634d18a10fa371b38dbe3e48f8c9e authored by HwB on 30 October 2011, 00:00:00 UTC
version 0.8.1
Tip revision: a7001ff
.Rapp.history
library(pracma)
quadgk <- function(f, a, b, tol = .Machine$double.eps^0.5, ...) {#
    stopifnot(is.numeric(a), length(a) == 1,#
              is.numeric(b), length(b) == 1)#
    eps <- .Machine$double.eps#
#
    fun <- match.fun(f)#
    f <- function(x) fun(x, ...)#
    #
    if (a == b)     return(0)#
    else if (a > b) return(-1 * quadgk(f, b, a, tol = tol))#
#
    # Nodes and weights for Gauss-Kronrod (7, 15)#
    n15 <- c(-0.9914553711208126, -0.9491079123427585, -0.8648644233597691,#
             -0.7415311855993944, -0.5860872354676911, -0.4058451513773972,#
             -0.2077849550078985,  0.0,                 0.2077849550078985,#
              0.4058451513773972,  0.5860872354676911,  0.7415311855993944,#
              0.8648644233597691,  0.9491079123427585,  0.9914553711208126)#
    n7  <- c(-0.9491079123427585, -0.7415311855993944, -0.4058451513773972,#
               0.0,#
               0.4058451513773972, 0.7415311855993944,  0.9491079123427585)#
    #
    w15 <- c(0.02293532201052922, 0.06309209262997855,  0.1047900103222502,#
             0.1406532597155259,  0.1690047266392679,   0.1903505780647854,#
             0.2044329400752989,  0.2094821410847278,   0.2044329400752989,#
             0.1903505780647854,  0.1690047266392679,   0.1406532597155259,#
             0.1047900103222502,  0.06309209262997855,  0.02293532201052922)#
    w7  <- c(0.1294849661688697,  0.2797053914892767,   0.3818300505051189,#
             0.4179591836734694,#
             0.3818300505051189,  0.2797053914892767,   0.1294849661688697)#
#
    .gkadpt <- function(f, a, b, tol = tol) {#
        # use nodes and weights from the environment#
        x15 <- 0.5 * ((b - a) * n15 + b + a)#
        x7  <- 0.5 * ((b - a) * n7  + b + a)#
        Q7  <- sum(w7  * f(x7))  * (b-a)/2#
        Q15 <- sum(w15 * f(x15)) * (b-a)/2#
        cat(a, b, G7, K15, "\n")#
#
        if (!is.finite(Q7) || !is.finite(Q15)) {#
            warning("Infinite or NA function value encountered.")#
            return(Q15)#
        } else if (abs(Q15 - Q7) < tol) {#
            return(Q15)#
        } else if (abs(b-a) < 16*eps) {#
            warning("Minimum step size reached; singularity possible.")#
            return(Q2)#
        } # else#
#
        Q2 <- gaussgk(f, (a+b)/2, b, tol = tol)#
        Q1 <- gaussgk(f, a, (a+b)/2, tol = tol)#
#
        return(Q1 + Q2)#
    }#
#
    # start the recursive procedure#
    .gkadpt(f, a, b, tol = tol)#
}
quadgk(sin, 0, 1)
quadgk <- function(f, a, b, tol = .Machine$double.eps^0.5, ...) {#
    stopifnot(is.numeric(a), length(a) == 1,#
              is.numeric(b), length(b) == 1)#
    eps <- .Machine$double.eps#
#
    fun <- match.fun(f)#
    f <- function(x) fun(x, ...)#
    #
    if (a == b)     return(0)#
    else if (a > b) return(-1 * quadgk(f, b, a, tol = tol))#
#
    # Nodes and weights for Gauss-Kronrod (7, 15)#
    n15 <- c(-0.9914553711208126, -0.9491079123427585, -0.8648644233597691,#
             -0.7415311855993944, -0.5860872354676911, -0.4058451513773972,#
             -0.2077849550078985,  0.0,                 0.2077849550078985,#
              0.4058451513773972,  0.5860872354676911,  0.7415311855993944,#
              0.8648644233597691,  0.9491079123427585,  0.9914553711208126)#
    n7  <- c(-0.9491079123427585, -0.7415311855993944, -0.4058451513773972,#
               0.0,#
               0.4058451513773972, 0.7415311855993944,  0.9491079123427585)#
    #
    w15 <- c(0.02293532201052922, 0.06309209262997855,  0.1047900103222502,#
             0.1406532597155259,  0.1690047266392679,   0.1903505780647854,#
             0.2044329400752989,  0.2094821410847278,   0.2044329400752989,#
             0.1903505780647854,  0.1690047266392679,   0.1406532597155259,#
             0.1047900103222502,  0.06309209262997855,  0.02293532201052922)#
    w7  <- c(0.1294849661688697,  0.2797053914892767,   0.3818300505051189,#
             0.4179591836734694,#
             0.3818300505051189,  0.2797053914892767,   0.1294849661688697)#
#
    .gkadpt <- function(f, a, b, tol = tol) {#
        # use nodes and weights from the environment#
        x15 <- 0.5 * ((b - a) * n15 + b + a)#
        x7  <- 0.5 * ((b - a) * n7  + b + a)#
        Q7  <- sum(w7  * f(x7))  * (b-a)/2#
        Q15 <- sum(w15 * f(x15)) * (b-a)/2#
        cat(a, b, Q7, Q15, "\n")#
#
        if (!is.finite(Q7) || !is.finite(Q15)) {#
            warning("Infinite or NA function value encountered.")#
            return(Q15)#
        } else if (abs(Q15 - Q7) < tol) {#
            return(Q15)#
        } else if (abs(b-a) < 16*eps) {#
            warning("Minimum step size reached; singularity possible.")#
            return(Q2)#
        } # else#
#
        Q2 <- gaussgk(f, (a+b)/2, b, tol = tol)#
        Q1 <- gaussgk(f, a, (a+b)/2, tol = tol)#
#
        return(Q1 + Q2)#
    }#
#
    # start the recursive procedure#
    .gkadpt(f, a, b, tol = tol)#
}
quadgk(sin, 0, 1)
quadgk(sin, 0, pi)
quadgk(sin, 0, pi, tol=1e-15)
quadgk <- function(f, a, b, tol = .Machine$double.eps^0.5, ...) {#
    stopifnot(is.numeric(a), length(a) == 1,#
              is.numeric(b), length(b) == 1)#
    eps <- .Machine$double.eps#
#
    fun <- match.fun(f)#
    f <- function(x) fun(x, ...)#
    #
    if (a == b)     return(0)#
    else if (a > b) return(-1 * quadgk(f, b, a, tol = tol))#
#
    # Nodes and weights for Gauss-Kronrod (7, 15)#
    n15 <- c(-0.9914553711208126, -0.9491079123427585, -0.8648644233597691,#
             -0.7415311855993944, -0.5860872354676911, -0.4058451513773972,#
             -0.2077849550078985,  0.0,                 0.2077849550078985,#
              0.4058451513773972,  0.5860872354676911,  0.7415311855993944,#
              0.8648644233597691,  0.9491079123427585,  0.9914553711208126)#
    n7  <- c(-0.9491079123427585, -0.7415311855993944, -0.4058451513773972,#
               0.0,#
               0.4058451513773972, 0.7415311855993944,  0.9491079123427585)#
    #
    w15 <- c(0.02293532201052922, 0.06309209262997855,  0.1047900103222502,#
             0.1406532597155259,  0.1690047266392679,   0.1903505780647854,#
             0.2044329400752989,  0.2094821410847278,   0.2044329400752989,#
             0.1903505780647854,  0.1690047266392679,   0.1406532597155259,#
             0.1047900103222502,  0.06309209262997855,  0.02293532201052922)#
    w7  <- c(0.1294849661688697,  0.2797053914892767,   0.3818300505051189,#
             0.4179591836734694,#
             0.3818300505051189,  0.2797053914892767,   0.1294849661688697)#
#
    .gkadpt <- function(f, a, b, tol = tol) {#
        # use nodes and weights from the environment#
        x15 <- 0.5 * ((b - a) * n15 + b + a)#
        x7  <- 0.5 * ((b - a) * n7  + b + a)#
        Q7  <- sum(w7  * f(x7))  * (b-a)/2#
        Q15 <- sum(w15 * f(x15)) * (b-a)/2#
        cat(a, b, Q7, Q15, "\n")#
#
        if (!is.finite(Q7) || !is.finite(Q15)) {#
            warning("Infinite or NA function value encountered.")#
            return(Q15)#
        } else if (abs(Q15 - Q7) < tol) {#
            return(Q15)#
        } else if (abs(b-a) < 16*eps) {#
            warning("Minimum step size reached; singularity possible.")#
            return(Q2)#
        } # else#
#
        Q2 <- .gkadpt(f, (a+b)/2, b, tol = tol)#
        Q1 <- .gkadpt(f, a, (a+b)/2, tol = tol)#
#
        return(Q1 + Q2)#
    }#
#
    # start the recursive procedure#
    .gkadpt(f, a, b, tol = tol)#
}
quadgk(sin, 0, pi, tol=1e-15)
quadgk(sin, 0, pi, tol=1e-15) - 2
quadgk(sin, pi, 0, tol=1e-15) - 2
quadgk(sin, pi, 0, tol=1e-15) +2
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
quadgk(f1, 0, pi) - 12.0703463163896
quadgk(f2, -1, 1)- 0.549360306778006
quadgk(f3, 0, 2) -  -0.460501533846733
quadgk(f4, 0, pi)-  2.0
quadgk(f4, 0, pi, tol=1e-10)-  2.0
quadgk(f4, 0, pi, tol=1e-12)-  2.0
quadgk <- function(f, a, b, tol = .Machine$double.eps^0.5, ...) {#
    stopifnot(is.numeric(a), length(a) == 1,#
              is.numeric(b), length(b) == 1)#
    eps <- .Machine$double.eps#
#
    fun <- match.fun(f)#
    f <- function(x) fun(x, ...)#
    #
    if (a == b)     return(0)#
    else if (a > b) return(-1 * quadgk(f, b, a, tol = tol))#
#
    # Nodes and weights for Gauss-Kronrod (7, 15)#
    n15 <- c(-0.9914553711208126, -0.9491079123427585, -0.8648644233597691,#
             -0.7415311855993944, -0.5860872354676911, -0.4058451513773972,#
             -0.2077849550078985,  0.0,                 0.2077849550078985,#
              0.4058451513773972,  0.5860872354676911,  0.7415311855993944,#
              0.8648644233597691,  0.9491079123427585,  0.9914553711208126)#
    n7  <- c(-0.9491079123427585, -0.7415311855993944, -0.4058451513773972,#
               0.0,#
               0.4058451513773972, 0.7415311855993944,  0.9491079123427585)#
    #
    w15 <- c(0.02293532201052922, 0.06309209262997855,  0.1047900103222502,#
             0.1406532597155259,  0.1690047266392679,   0.1903505780647854,#
             0.2044329400752989,  0.2094821410847278,   0.2044329400752989,#
             0.1903505780647854,  0.1690047266392679,   0.1406532597155259,#
             0.1047900103222502,  0.06309209262997855,  0.02293532201052922)#
    w7  <- c(0.1294849661688697,  0.2797053914892767,   0.3818300505051189,#
             0.4179591836734694,#
             0.3818300505051189,  0.2797053914892767,   0.1294849661688697)#
#
    .gkadpt <- function(f, a, b, tol = tol) {#
        # use nodes and weights from the environment#
        x15 <- 0.5 * ((b - a) * n15 + b + a)#
        x7  <- 0.5 * ((b - a) * n7  + b + a)#
        Q7  <- sum(w7  * f(x7))  * (b-a)/2#
        Q15 <- sum(w15 * f(x15)) * (b-a)/2#
#
        if (!is.finite(Q7) || !is.finite(Q15)) {#
            warning("Infinite or NA function value encountered.")#
            return(Q15)#
        } else if (abs(Q15 - Q7) < tol) {#
            return(Q15)#
        } else if (abs(b-a) < 16*eps) {#
            warning("Minimum step size reached; singularity possible.")#
            return(Q2)#
        } # else#
#
        Q2 <- .gkadpt(f, (a+b)/2, b, tol = tol)#
        Q1 <- .gkadpt(f, a, (a+b)/2, tol = tol)#
#
        return(Q1 + Q2)#
    }#
#
    # start the recursive procedure#
    .gkadpt(f, a, b, tol = tol)#
}
all.equal(quadgk(f1, 0, pi), 12.0703463163896,#
          tolerance = 1e-13)#
all.equal(quadgk(f2, -1, 1),  0.549360306778006,#
          tolerance = 1e-15)#
all.equal(quadgk(f3, 0, 2),  -0.460501533846733,#
          tolerance = 1e-14)#
all.equal(quadgk(f4, 0, pi, tol = 1e-12),  2.0,#
          tolerance = 1e-12)
setwd('/Users/HwB/Work/myRforge/optimist/pkg/pracma/tests')
all.equal(quadgk(f1, 0, pi), 12.0703463163896,#
          tolerance = 1e-13)#
all.equal(quadgk(f2, -1, 1),  0.549360306778006,#
          tolerance = 1e-15)#
all.equal(quadgk(f3, 0, 2),  -0.460501533846733,#
          tolerance = 1e-13)#
all.equal(quadgk(f4, 0, pi, tol = 1e-12),  2.0,#
          tolerance = 1e-12)
back to top