set.seed(123) library("numDeriv") library("tensorA") library("microbenchmark") library("papeR") ## install.packages("calculus") library("calculus") library("cubature") ## ----------------------------------------------------------------------------- ("a + b" %prod% 1i) %sum% (0 %prod% "c") %diff% (expression(d + e) %div% 3) ## ----------------------------------------------------------------------------- eval(parse(text = "a" %sum% "a"), list(a = 1)) ## ----------------------------------------------------------------------------- x <- array(letters[1:6], dim = c(2, 3)) evaluate(x, var = c(a = 1, b = 2, c = 3, d = 4, e = 5, f = 6)) ## ----------------------------------------------------------------------------- x <- array(letters[1:6], dim = c(2, 3)) var <- data.frame(a = 1:2, b = 2:3, c = 3:4, d = 4:5, e = 5:6, f = 6:7) evaluate(x, var = var) ## ----------------------------------------------------------------------------- cross(c(1, 0, 0), c(0, 1, 0)) ## ----------------------------------------------------------------------------- cross(c(1, 0, 0, 0), c(0, 1, 0, 0), c(0, 0, 0, 1)) ## ------------------------------------------------------------------------------ options(calculus.auto.wrap = FALSE) ## ----------------------------------------------------------------------------- cross(c("a", "b", "c"), c("d", "e", "f")) ## ----------------------------------------------------------------------------- options(calculus.auto.wrap = TRUE) ## ----------------------------------------------------------------------------- mxdet(matrix(1:4, nrow = 2)) ## ----------------------------------------------------------------------------- options(calculus.auto.wrap = FALSE) ## ----------------------------------------------------------------------------- mxdet(matrix(letters[1:4], nrow = 2)) ## ----------------------------------------------------------------------------- options(calculus.auto.wrap = TRUE) ## ---- warning=FALSE----------------------------------------------------------- n <- 4 e <- letters[1:n^2] grid <- expand.grid(lapply(1:n^2, function(e) runif(2))) colnames(grid) <- e microbenchmark( numeric = { x <- apply(grid, 1, function(e) det(matrix(e, nrow = n))) }, symbolic = { x <- evaluate(mxdet(matrix(e, nrow = n)), grid) }) ## ----------------------------------------------------------------------------- mxinv(matrix(1:4, byrow = TRUE, nrow = 2)) ## ----------------------------------------------------------------------------- options(calculus.auto.wrap = FALSE) ## ----------------------------------------------------------------------------- mxinv(matrix(letters[1:4], byrow = TRUE, nrow = 2)) ## ----------------------------------------------------------------------------- options(calculus.auto.wrap = TRUE) ## ---- warning=FALSE----------------------------------------------------------- n <- 4 e <- letters[1:n^2] grid <- expand.grid(lapply(1:n^2, function(e) runif(2))) colnames(grid) <- e microbenchmark( numeric = { x <- apply(grid, 1, function(e) solve(matrix(e, nrow = n))) }, symbolic = { x <- evaluate(mxinv(matrix(e, nrow = n)), grid) }) ## ----------------------------------------------------------------------------- a <- matrix(1:4, nrow = 2, byrow = TRUE) b <- matrix(letters[1:4], nrow = 2, byrow = TRUE) a %mx% b ## ----------------------------------------------------------------------------- A <- array(1:24, dim = c(2, 3, 4)) attributes(A) ## ----------------------------------------------------------------------------- index(A) <- c("i", "j", "k") attributes(A) ## ----------------------------------------------------------------------------- epsilon(2) ## ----------------------------------------------------------------------------- epsilon(3) ## ----------------------------------------------------------------------------- delta(n = 3, p = 1) ## ----------------------------------------------------------------------------- x <- array(1:8, dim = c(2, 2, 2)) print(x) ## ----------------------------------------------------------------------------- contraction(x) ## ----------------------------------------------------------------------------- index(x) <- c("i", "j", "i") contraction(x) ## ----------------------------------------------------------------------------- index(x) <- c("i", "j", "i") contraction(x, drop = FALSE) ## ----------------------------------------------------------------------------- A <- array(1:6, dim = c(i = 2, j = 3)) print(A) ## ----------------------------------------------------------------------------- B <- array(1:4, dim = c(k = 2, i = 2)) print(B) ## ----------------------------------------------------------------------------- C <- array(letters[1:4], dim = c(i = 2, i = 2)) print(C) ## ----------------------------------------------------------------------------- einstein(A, B, C) ## ---- warning=FALSE---------------------------------------------------------------- a <- array(1:1e+06, dim = c(a = 2, i = 5, j = 100, k = 50, d = 20)) b <- array(1:1e+05, dim = c(a = 2, j = 100, i = 5, l = 100)) Ta <- tensorA::to.tensor(a) Tb <- tensorA::to.tensor(b) microbenchmark( calculus = calculus::einstein(a, b), tensorA = tensorA::einstein.tensor(Ta, Tb) ) ## ----------------------------------------------------------------------------- 1:3 %inner% letters[1:3] ## ----------------------------------------------------------------------------- matrix(1:6, byrow = TRUE, nrow = 2, ncol = 3) %dot% letters[1:3] ## ----------------------------------------------------------------------------- 1:3 %outer% letters[1:3] ## ----------------------------------------------------------------------------- 1:3 %kronecker% letters[1:3] ## ----------------------------------------------------------------------------- derivative(f = "sin(x)", var = "x") ## ----------------------------------------------------------------------------- sym <- derivative(f = "sin(x)", var = c(x = 0)) num <- derivative(f = function(x) sin(x), var = c(x = 0)) ## ---- echo=FALSE-------------------------------------------------------------- print(c(Symbolic = sym, Numeric = num)) ## ----------------------------------------------------------------------------- sym <- derivative(f = "sin(x)", var = c(x = 0), order = 4) num <- derivative(f = function(x) sin(x), var = c(x = 0), order = 4) ## ---- echo=FALSE-------------------------------------------------------------- print(c(Symbolic = sym, Numeric = num)) ## ----------------------------------------------------------------------------- derivative(f = "y^2 * sin(x)", var = c("x", "y"), order = c(1, 2)) ## ----------------------------------------------------------------------------- f <- function(x, y) y^2 * sin(x) derivative(f, var = c(x = 0, y = 0), order = c(1, 2), accuracy = 6) ## ----------------------------------------------------------------------------- derivative("x^2 * y^2", var = c("x", "y")) ## ----------------------------------------------------------------------------- derivative("x^6 * y^6", var = c("x", "y"), order = 6) ## ----------------------------------------------------------------------------- f <- function(x, y) x^2 * y^2 derivative(f, var = c(x = 1, y = 2)) ## ----------------------------------------------------------------------------- f <- function(x, y) c(x * y, x^2 * y^2) derivative(f, var = c(x = 1, y = 2)) ## ----------------------------------------------------------------------------- f <- function(x) c(sum(x), prod(x)) derivative(f, var = c(0, 0, 0)) ## ---- echo=FALSE, results='asis', warning=FALSE, message=FALSE---------------- f.sym <- c("x^2 * exp(x)", "x * sin(x^2)", "x * log(x^2)", "exp(sin(x))") f.num <- function(x) { c(x^2 * exp(x), x * sin(x^2), x * log(x^2), exp(sin(x))) } x <- log(seq(1.1, 1000, by = 0.1)) s <- evaluate(derivative(f.sym, "x"), data.frame(x = x)) c <- sapply(x, function(x) calculus::derivative(f.num, x, accuracy = 4)) n <- sapply(x, function(x) c(numDeriv::jacobian(f.num, x))) var <- c("x^2 exp(x)", "x sin(x^2)", "x log(x^2)", "exp(sin(x))") df1 <- as.data.frame(s/t(c) - 1) colnames(df1) <- var df1$pkg <- factor("calculus") df2 <- as.data.frame(s/t(n) - 1) colnames(df2) <- var df2$pkg <- factor("numDeriv") x <- summarise(rbind(df1, df2), group = "pkg", test = FALSE, digits = 1e+16, sanitize = TRUE, quantiles = FALSE, count = TRUE) x # xtable(x, digits = -1) ## ---- warning=FALSE----------------------------------------------------------- x <- rep(0, 1000) f <- function(x) sum(x) microbenchmark(`calculus O(2)` = calculus::derivative(f, x, accuracy = 2), `calculus O(4)` = calculus::derivative(f, x, accuracy = 4), `calculus O(6)` = calculus::derivative(f, x, accuracy = 6), `calculus O(8)` = calculus::derivative(f, x, accuracy = 8), numDeriv = numDeriv::grad(f,x)) ## ----------------------------------------------------------------------------- partitions(n = 2, length = 3, fill = TRUE, perm = TRUE, equal = FALSE) ## ----------------------------------------------------------------------------- taylor("exp(x)", var = "x", order = 2) ## ----------------------------------------------------------------------------- f <- function(x, y) log(y) * sin(x) taylor(f, var = c(x = 0, y = 1), order = 2) ## ----------------------------------------------------------------------------- D(D(expression(exp(-x^2/2)), "x"), "x") ## ----------------------------------------------------------------------------- f <- expression(exp(-x^2/2)) for (i in 1:14) f <- D(f, "x") object.size(f) ## ----------------------------------------------------------------------------- sim <- ode(f = c("x", "x * (1 + cos(10*t))"), var = c(x = 1, y = 1), times = seq(0, 2 * pi, by = 0.001), timevar = "t") ## ---- fig.cap='Solution to the system of ordinary differential equations in Section \\ref{sec:ode} obtained with the function \\code{ode} using the Runge-Kutta method.', fig.pos='h', fig.height=3---- par(mar = c(4, 4, 1, 4)) times <- as.numeric(rownames(sim)) matplot(times, sim, type = "l", lty = 1, col = 1:2, xlab = "Time", ylab = "Value") legend("topleft", legend = c("x", "y"), col = 1:2, cex = 0.8, fill = 1:2) gradient("x * y * z", var = c("x", "y", "z")) ## ----------------------------------------------------------------------------- gradient("x * y * z", var = c("x", "y", "z"), coordinates = "spherical") ## ----------------------------------------------------------------------------- gradient("x * y * z", var = c("x", "y", "z"), coordinates = c(1, "x", "x*sin(y)")) ## ----------------------------------------------------------------------------- f <- function(x, y, z) x * y * z gradient(f, var = c(x = 1, y = pi/2, z = 0), coordinates = "spherical") ## ----------------------------------------------------------------------------- f <- function(x) x[1] * x[2] * x[3] gradient(f, var = c(1, pi/2, 0), coordinates = "spherical") ## ----------------------------------------------------------------------------- f <- function(x) c(prod(x), sum(x)) gradient(f, var = c(3, 2, 1)) ## ----------------------------------------------------------------------------- f <- function(x) c(prod(x), sum(x)) gradient(f, var = c(3, 2, 1), coordinates = "cylindrical") ## ----------------------------------------------------------------------------- f <- function(x) x^2 jacobian(f, var = c(1)) ## ----------------------------------------------------------------------------- f <- function(x, y, z) x * y * z hessian(f, var = c(x = 3, y = 2, z = 1)) ## ----------------------------------------------------------------------------- f <- function(x, y, z) c(x * y * z, x + y + z) h <- hessian(f, var = c(x = 3, y = 2, z = 1)) ## ----------------------------------------------------------------------------- h[1, , ] h[2, , ] ## ----------------------------------------------------------------------------- f <- c("x^2", "y^2", "z^2") divergence(f, var = c("x", "y", "z")) ## ----------------------------------------------------------------------------- f <- c("sqrt(r)/10", "sqrt(r)") divergence(f, var = c("r", "phi"), coordinates = "polar") ## ----------------------------------------------------------------------------- f <- matrix(c("x^2", "y^2", "z^2", "x", "y", "z"), nrow = 2, byrow = TRUE) divergence(f, var = c("x", "y", "z")) ## ----------------------------------------------------------------------------- f <- function(x, y, z) { matrix(c(x^2, y^2, z^2, x, y, z), nrow = 2, byrow = TRUE) } divergence(f, var = c(x = 0, y = 0, z = 0)) ## ----------------------------------------------------------------------------- f <- c("x^3 * y^2", "x") curl(f, var = c("x", "y")) ## ----------------------------------------------------------------------------- f <- c("x", "-y", "z") curl(f, var = c("x", "y", "z")) ## ----------------------------------------------------------------------------- f <- matrix(c("x", "-y", "z", "x^3 * y^2", "x", "0"), nrow = 2, byrow = TRUE) curl(f, var = c("x", "y", "z")) ## ----------------------------------------------------------------------------- f <- "x^3 + y^3 + z^3" laplacian(f, var = c("x", "y", "z")) ## ----------------------------------------------------------------------------- f <- array(c("x^3 + y^3 + z^3", "x^2 + y^2 + z^2", "y^2", "z * x^2"), dim = c(2, 2)) laplacian(f, var = c("x", "y", "z")) ## ----------------------------------------------------------------------------- i <- integral(f = "x", bounds = list(x = c(0, 1))) i$value ## ----------------------------------------------------------------------------- i <- integral(f = function(x) x, bounds = list(x = c(0, 1))) i$value ## ----------------------------------------------------------------------------- i <- integral(f = "y * x", bounds = list(x = c(0, 1), y = 2)) i$value ## ----------------------------------------------------------------------------- i <- integral(f = "y * x", bounds = list(x = c(0, 1), y = c(0, 1))) i$value ## ----------------------------------------------------------------------------- i <- integral(f = 1, bounds = list(r = c(0, 1), theta = c(0, 2 * pi)), coordinates = "polar") i$value ## ----------------------------------------------------------------------------- i <- integral(f = 1, bounds = list(r = c(0, 1), theta = c(0, pi), phi = c(0, 2 * pi)), coordinates = "spherical") i$value ## ----------------------------------------------------------------------------- V <- "1/(4 * pi * r)" ## ----------------------------------------------------------------------------- var <- c("r", "theta", "phi") E <- -1 %prod% gradient(V, var = var, coordinates = "spherical") ## ----------------------------------------------------------------------------- i <- integral(E[1], bounds = list(r = 1, theta = c(0, pi), phi = c(0, 2 * pi)), coordinates = "spherical") i$value