# High Dimensional Numerical and Symbolic Calculus in R   Efficient C++ optimized functions for numerical and symbolic calculus. It includes basic symbolic arithmetic, tensor calculus, Einstein summing convention, fast computation of the Levi-Civita symbol and generalized Kronecker delta, Taylor series expansion, multivariate Hermite polynomials, accurate high-order derivatives, differential operators (Gradient, Jacobian, Hessian, Divergence, Curl, Laplacian) and Monte Carlo integration in arbitrary orthogonal coordinate systems: cartesian, polar, spherical, cylindrical, parabolic or user defined by custom scale factors.

## Quickstart

``````# Install calculus
install.packages('calculus')

require('calculus')``````

## Usage

### `derivative`: Numerical and Symbolic Derivatives

Description

Compute symbolic derivatives based on the `D` function, or accurate and reliable numerical derivatives based on finite differences.

Usage

``derivative(f, var = "x", order = 1, accuracy = 2, stepsize = NULL, deparse = TRUE)``

Arguments

Argument Description
`f` function, expression or character array.
`var` character vector, giving the variable names with respect to which derivatives will be computed. If a named vector is provided, derivatives will be computed at that point. See examples.
`order` integer vector, giving the differentiation order for each variable. See details.
`accuracy` accuracy degree for numerical derivatives.
`stepsize` finite differences stepsize for numerical derivatives. Auto-optimized by default.
`deparse` logical. Return character instead of expression or call?

Details

The function behaves differently depending on the length of the `order` argument.

If `order` is of length 1, then the n-th order derivative is computed for each function with respect to each variable: where F is the tensor of functions and is the tensor of variable names with respect to which the n-th order derivatives will be computed.

If `order` matches the length of `var` , then it is assumed that the differentiation order is provided for each variable. In this case, each function will be derived ni times with respect to the i-th variable, for each of the j variables: where F is the tensor of functions to differentiate.

If `var` is a named vector, e.g. `c(x = 0, y = 0)` , derivatives will be computed at that point. Note that if `f` is a function, then `var` must be a named vector giving the point at which the numerical derivatives will be computed.

Value

array of derivatives.

Examples

``````# derive f with respect to x
derivative(f = "sin(x)", var = "x")

# derive f with respect to x and evaluate in x = 0
derivative(f = "sin(x)", var = c("x" = 0))

# derive f twice with respect to x
derivative(f = "sin(x)", var = "x", order = 2)

# derive f once with respect to x, and twice with respect to y
derivative(f = "y^2*sin(x)", var = c("x","y"), order = c(1,2))

# compute the gradient of f with respect to (x,y)
derivative(f = "y*sin(x)", var = c("x","y"))

# compute the Jacobian of f with respect to (x,y)
f <- c("y*sin(x)", "x*cos(y)")
derivative(f = f, var = c("x","y"))

# compute the Hessian of f with respect to (x,y)
g <- derivative(f = "y^2*sin(x)", var = c("x","y"))
derivative(f = g, var = c("x","y"))

# compute the Jacobian of f with respect to (x,y) and evaluate in (0,0)
f1 <- function(x, y) y*sin(x)
f2 <- function(x, y) x*cos(y)
derivative(f = c(f1, f2), var = c("x"=0,"y"=0))``````

### `gradient`: Numerical and Symbolic Gradient

Description

Compute the gradient or jacobian of functions, expressions and characters.

Usage

``````gradient(f, var, accuracy = 2, stepsize = NULL, coordinates = "cartesian")

Arguments

Argument Description
`f` function, expression or character array.
`var` character vector, giving the variable names with respect to which derivatives will be computed. If a named vector is provided, derivatives will be computed at that point.
`accuracy` accuracy degree for numerical derivatives.
`stepsize` finite differences stepsize for numerical derivatives. Auto-optimized by default.
`coordinates` coordinate system to use. One of: `cartesian` , `polar` , `spherical` , `cylindrical` , `parabolic` , `parabolic-cylindrical` or a character vector of scale factors for each varibale.

Value

Examples

``````# gradient with respect to x
gradient(f = "sin(x)", var = "x")

# gradient with respect to x and evaluate in x = 0
gradient(f = "sin(x)", var = c("x" = 0))

# gradient with respect to (x,y)
gradient(f = "y*sin(x)", var = c("x","y"))

# jacobian with respect to (x,y)
f <- c("y*sin(x)", "x*cos(y)")
gradient(f = f, var = c("x","y"))

# jacobian with respect to (x,y) and evaluate in (x = 0, y = 0)
f <- c(function(x, y) y*sin(x), function(x, y) x*cos(y))
gradient(f = f, var = c(x=0,y=0))

gradient('r*theta*phi', var = c('r','theta','phi'), coordinates = 'spherical')

# numerical gradient in spherical coordinates
f <- function(r, theta, phi) r*theta*phi
gradient(f, var = c('r'=1, 'theta'=pi/4, 'phi'=pi/4), coordinates = 'spherical')``````

### `hessian`: Numerical and Symbolic Hessian

Description

Compute the hessian matrix of functions, expressions and characters.

Usage

``````hessian(f, var, accuracy = 2, stepsize = NULL, coordinates = "cartesian")
f %hessian% var``````

Arguments

Argument Description
`f` function, expression or character.
`var` character vector, giving the variable names with respect to which derivatives will be computed. If a named vector is provided, derivatives will be computed at that point.
`accuracy` accuracy degree for numerical derivatives.
`stepsize` finite differences stepsize for numerical derivatives. Auto-optimized by default.
`coordinates` coordinate system to use. One of: `cartesian` , `polar` , `spherical` , `cylindrical` , `parabolic` , `parabolic-cylindrical` or a character vector of scale factors for each varibale.

Value

hessian matrix.

Examples

``````# hessian with respect to x
hessian(f = "sin(x)", var = "x")
"sin(x)" %hessian% "x"

# hessian with respect to x and evaluate in x = 0
hessian(f = "sin(x)", var = c("x" = 0))
"sin(x)" %hessian% c(x=0)

# hessian with respect to (x,y)
hessian(f = "y*sin(x)", var = c("x","y"))
"y*sin(x)" %hessian% c("x","y")

# hessian in spherical coordinates
hessian('r*theta*phi', var = c('r','theta','phi'), coordinates = 'spherical')

# numerical hessian in spherical coordinates
f <- function(r, theta, phi) r*theta*phi
hessian(f, var = c('r'=1, 'theta'=pi/4, 'phi'=pi/4), coordinates = 'spherical')``````

### `divergence`: Numerical and Symbolic Divergence

Description

Compute the divergence of functions, expressions and characters.

Usage

``````divergence(f, var, accuracy = 2, stepsize = NULL, coordinates = "cartesian")
f %divergence% var``````

Arguments

Argument Description
`f` function, expression or character array.
`var` character vector, giving the variable names with respect to which derivatives will be computed. If a named vector is provided, derivatives will be computed at that point.
`accuracy` accuracy degree for numerical derivatives.
`stepsize` finite differences stepsize for numerical derivatives. Auto-optimized by default.
`coordinates` coordinate system to use. One of: `cartesian` , `polar` , `spherical` , `cylindrical` , `parabolic` , `parabolic-cylindrical` or a character vector of scale factors for each varibale.

Value

divergence array.

Examples

``````# divergence of a vector field
f <- c('x^2','y^3','z^4')
divergence(f, var = c('x','y','z'))
f %divergence% c('x','y','z')

# numerical divergence of a vector field
f <- c(function(x,y,z) x^2, function(x,y,z) y^3, function(x,y,z) z^4)
divergence(f, var = c('x'=1,'y'=1,'z'=1))
f %divergence% c('x'=1,'y'=1,'z'=1)

# divergence of array of vector fields
f1 <- c('x^2','y^3','z^4')
f2 <- c('x','y','z')
a <- matrix(c(f1,f2), nrow = 2, byrow = TRUE)
divergence(a, var = c('x','y','z'))
a %divergence% c('x','y','z')

# divergence in polar coordinates
f <- c('sqrt(r)/10','sqrt(r)')
divergence(f, var = c('r','phi'), coordinates = 'polar')``````

### `curl`: Numerical and Symbolic Curl

Description

Compute the curl of functions, expressions and characters.

Usage

``````curl(f, var, accuracy = 2, stepsize = NULL, coordinates = "cartesian")
f %curl% var``````

Arguments

Argument Description
`f` function, expression or character array.
`var` character vector, giving the variable names with respect to which derivatives will be computed. If a named vector is provided, derivatives will be computed at that point.
`accuracy` accuracy degree for numerical derivatives.
`stepsize` finite differences stepsize for numerical derivatives. Auto-optimized by default.
`coordinates` coordinate system to use. One of: `cartesian` , `polar` , `spherical` , `cylindrical` , `parabolic` , `parabolic-cylindrical` or a character vector of scale factors for each varibale.

Value

curl array.

Examples

``````# curl of a vector field
f <- c('x*y','y*z','x*z')
curl(f, var = c('x','y','z'))
f %curl% c('x','y','z')

# irrotational vector field
f <- c('x','-y','z')
curl(f, var = c('x','y','z'))
f %curl% c('x','y','z')

# numerical curl of a vector field
f <- c(function(x,y,z) x*y, function(x,y,z) y*z, function(x,y,z) x*z)
curl(f, var = c('x'=1,'y'=1,'z'=1))
f %curl% c('x'=1,'y'=1,'z'=1)

# curl of array of vector fields
f1 <- c('x*y','y*z','z*x')
f2 <- c('x','-y','z')
a <- matrix(c(f1,f2), nrow = 2, byrow = TRUE)
curl(a, var = c('x','y','z'))
a %curl% c('x','y','z')

# curl in polar coordinates
f <- c('sqrt(r)/10','sqrt(r)')
curl(f, var = c('r','phi'), coordinates = 'polar')``````

### `laplacian`: Numerical and Symbolic Laplacian

Description

Compute the laplacian of functions, expressions and characters.

Usage

``````laplacian(f, var, accuracy = 2, stepsize = NULL, coordinates = "cartesian")
f %laplacian% var``````

Arguments

Argument Description
`f` function, expression or character array.
`var` character vector, giving the variable names with respect to which derivatives will be computed. If a named vector is provided, derivatives will be computed at that point.
`accuracy` accuracy degree for numerical derivatives.
`stepsize` finite differences stepsize for numerical derivatives. Auto-optimized by default.
`coordinates` coordinate system to use. One of: `cartesian` , `polar` , `spherical` , `cylindrical` , `parabolic` , `parabolic-cylindrical` or a character vector of scale factors for each varibale.

Value

laplacian array.

Examples

``````# laplacian of a scalar field
f <- 'x^2+y^2+z^2'
laplacian(f, var = c('x','y','z'))
f %laplacian% c('x','y','z')

# laplacian of scalar fields
f <- c('x^2','y^3','z^4')
laplacian(f, var = c('x','y','z'))
f %laplacian% c('x','y','z')

# numerical laplacian of scalar fields
f <- c(function(x,y,z) x^2, function(x,y,z) y^3, function(x,y,z) z^4)
laplacian(f, var = c('x'=1,'y'=1,'z'=1))
f %laplacian% c('x'=1,'y'=1,'z'=1)

# laplacian of array of scalar fields
f1 <- c('x^2','y^3','z^4')
f2 <- c('x','y','z')
a <- matrix(c(f1,f2), nrow = 2, byrow = TRUE)
laplacian(a, var = c('x','y','z'))
a %laplacian% c('x','y','z')

# laplacian in polar coordinates
f <- c('sqrt(r)/10','sqrt(r)')
laplacian(f, var = c('r','phi'), coordinates = 'polar')``````

### `integral`: Monte Carlo Integration

Description

Integrate multidimensional functions, expressions, and characters in arbitrary orthogonal coordinate systems .

Usage

``integral(f, ..., err = 0.01, rel = TRUE, coordinates = "cartesian", verbose = TRUE)``

Arguments

Argument Description
`f` function, expression or characters.
`...` integration bounds.
`err` acuracy requested.
`rel` logical. Relative accuracy? If `FALSE` , use absolute accuracy.
`coordinates` coordinate system to use. One of: `cartesian` , `polar` , `spherical` , `cylindrical` , `parabolic` , `parabolic-cylindrical` or a character vector of scale factors for each varibale.
`verbose` logical. Print on progress?

Value

`list` with components

Name Description
`value` the final estimate of the integral.
`abs.error` estimate of the modulus of the absolute error.

Examples

``````# integrate character
integral('sin(x)', x = c(0,2*pi), rel = FALSE, verbose = FALSE)

# integrate expression
integral(parse(text = 'x'), x = c(0,1), verbose = FALSE)

# integrate function
integral(function(x) exp(x), x = c(0,1), verbose = FALSE)

# multivariate integral
integral(function(x,y) x*y, x = c(0,1), y = c(0,1), verbose = FALSE)

# surface of a sphere
integral('1', r = 1, theta = c(0,pi), phi = c(0,2*pi), coordinates = 'spherical', verbose = FALSE)

# volume of a sphere
integral('1', r = c(0,1), theta = c(0,pi), phi = c(0,2*pi), coordinates = 'spherical', verbose = FALSE)``````
``````
# Electric charge contained in a region of space
# Based on the divergence theorem and Maxwell's equations

# electric potential of unitary point charge
V <- '1/(4*pi*r)'

# electric field
E <- -1 %prod% gradient(V, c('r', 'theta', 'phi'), coordinates = 'spherical')

# electric charge
integral(E, r = 1, theta = c(0,pi), phi = c(0,2*pi), coordinates = 'spherical', verbose = FALSE)
``````

### `partitions`: Partitions of an Integer

Description

Fast algorithm for generating integer partitions.

Usage

``partitions(n, max = 0, length = 0, perm = FALSE, fill = FALSE, equal = TRUE)``

Arguments

Argument Description
`n` positive integer.
`max` maximum integer in the partitions.
`length` maximum number of elements in the partitions.
`perm` logical. Permute partitions?
`fill` logical. Fill partitions with zeros to match `length` ?
`equal` logical. Return only partition of `n` ? If `FALSE` , partitions of all integers less or equal to `n` are returned.

Value

`list` of partitions, or `data.frame` if `length>0` and `fill=TRUE` .

Examples

``````# partitions of 4
partitions(4)

# partitions of 4 and permute
partitions(4, perm = TRUE)

# partitions of 4 with max element 2
partitions(4, max = 2)

# partitions of 4 with 2 elements
partitions(4, length = 2)

# partitions of 4 with 3 elements, fill with zeros
partitions(4, length = 3, fill = TRUE)

# partitions of 4 with 3 elements, fill with zeros and permute
partitions(4, length = 3, fill = TRUE, perm = TRUE)

# partitions of all integers less or equal to 3
partitions(3, equal = FALSE)

# partitions of all integers less or equal to 3, fill to 2 elements and permute
partitions(3, equal = FALSE, length = 2, fill = TRUE, perm = TRUE)``````

### `taylor`: Taylor Series

Description

Compute the Taylor series approximation of functions, expressions or characters.

Usage

``taylor(f, var = "x", order = 1, accuracy = 2, stepsize = NULL)``

Arguments

Argument Description
`f` function, expression or character
`var` character. The variables of `f` .
`order` the order of the Taylor approximation.
`accuracy` accuracy degree for numerical derivatives.
`stepsize` finite differences stepsize for numerical derivatives. Auto-optimized by default.

Value

`list` with components

Name Description
`f` the Taylor series.
`order` the approximation order.
`terms` `data.frame` containing the variables, coefficients and degrees of each term in the Taylor series.

Examples

``````# univariate taylor series
taylor('exp(x)', var = 'x', order = 3)

# univariate taylor series of arbitrary functions
taylor(function(x) exp(x), var = 'x', order = 3)

# multivariate taylor series
taylor('sin(x*y)', var = c('x','y'), order = 6)

# multivariate taylor series of arbitrary functions
taylor(function(x,y) sin(x*y), var = c('x','y'), order = 6)``````

### `hermite`: Hermite Polynomials

Description

Compute univariate and multivariate Hermite polynomials.

Usage

``hermite(order, sigma = 1, var = "x")``

Arguments

Argument Description
`order` integer. The order of the Hermite polynomial.
`sigma` the covariance matrix of the Gaussian kernel.
`var` character. The variables of the polynomial.

Details

Hermite polynomials are obtained by successive differentiation of the Gaussian kernel where is a d-dimensional square matrix and is the vector representing the order of differentiation for each variable.

Value

`list` of Hermite polynomials with components

Name Description
`f` the Hermite polynomial.
`order` the order of the Hermite polynomial.
`terms` `data.frame` containing the variables, coefficients and degrees of each term in the Hermite polynomial.

Examples

``````# univariate Hermite polynomials up to order 3
hermite(3)

# univariate Hermite polynomials with variable z
hermite(3, var = 'z')

# multivariate Hermite polynomials up to order 2
hermite(order = 2, sigma = matrix(c(1,0,0,1), nrow = 2), var = c('z1', 'z2'))``````

### `kronecker`: Generalized Kronecker Delta

Description

Compute the Generalized Kronecker Delta.

Usage

``kronecker(n, p = 1)``

Arguments

Argument Description
`n` number of elements for each dimension.
`p` order of the generalized Kronecker delta, `p=1` for the standard Kronecker delta.

Value

array representing the generalized Kronecker delta tensor.

Examples

``````# Kronecker delta 3x3
kronecker(3)

# generalized Kronecker delta 3x3 of order 2 -> 3x3 x 3x3
kronecker(3, p = 2)``````

### `levicivita`: Levi-Civita Symbol

Description

Compute the Levi-Civita totally antisymmetric tensor.

Usage

``levicivita(n)``

Arguments

Argument Description
`n` dimension

Value

array representing the Levi-Civita tensor.

Examples

``````# Levi-Civita tensor in 2-d
levicivita(2)

# Levi-Civita tensor in 3-d
levicivita(3)``````

### `trace`: Tensor Contraction

Description

Sum over repeated indices in a tensor. Can be seen as a generalization of the trace.

Usage

``trace(x, i = NULL, drop = TRUE)``

Arguments

Argument Description
`x` array.
`i` subset of repeated indices to sum up. If `NULL` , the tensor contraction takes place on all repeated indices of `x` .
`drop` logical. Drop summation indices? If `FALSE` , keep dummy dimensions.

Value

array.

Examples

``````# trace of numeric matrix
x <- matrix(1:4, nrow = 2)
trace(x)

# trace of character matrix
x <- matrix(letters[1:4], nrow = 2)
trace(x)

# trace of a tensor (sum over diagonals)
x <- array(1:27, dim = c(3,3,3))
trace(x)

# tensor contraction over repeated indices
x <- array(1:27, dim = c(3,3,3))
index(x) <- c('i','i','j')
trace(x)

# tensor contraction over specific indices only
x <- array(1:16, dim = c(2,2,2,2))
index(x) <- c('i','i','k','k')
trace(x, i = 'k')

# tensor contraction keeping dummy dimensions
x <- array(letters[1:16], dim = c(2,2,2,2))
index(x) <- c('i','i','k','k')
trace(x, drop = FALSE)``````

### `einstein`: Numerical and Symbolic Einstein Summation

Description

Implement the Einstein notation for summation over repeated indices.

Usage

``einstein(..., drop = TRUE)``

Arguments

Argument Description
`...` arbitrary number of indexed arrays.
`drop` logical. Drop summation indices? If `FALSE` , keep dummy dimensions.

Value

array.

Examples ``````a <- array(1:10, dim = c(2,5))
b <- array(1:45, dim = c(5,3,3))
c <- array(1:12, dim = c(3,4))
d <- array(1:15, dim = c(5,3))

index(a) <- c('i','j')
index(b) <- c('j','k','k')
index(c) <- c('k', 'l')
index(d) <- c('j', 'k')

einstein(a,b,c,d)``````

## Documentation

https://cran.r-project.org/package=calculus