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')
# Load 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")
f %gradient% varArguments
| 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
gradient or jacobian array.
Examples
# gradient with respect to x
gradient(f = "sin(x)", var = "x")
"sin(x)" %gradient% "x"
# gradient with respect to x and evaluate in x = 0
gradient(f = "sin(x)", var = c("x" = 0))
"sin(x)" %gradient% c(x=0)
# gradient with respect to (x,y)
gradient(f = "y*sin(x)", var = c("x","y"))
"y*sin(x)" %gradient% c("x","y")
# jacobian with respect to (x,y)
f <- c("y*sin(x)", "x*cos(y)")
gradient(f = f, var = c("x","y"))
f %gradient% 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))
f %gradient% c(x=0,y=0)
# gradient in spherical coordinates
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% varArguments
| 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% varArguments
| 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% varArguments
| 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% varArguments
| 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[1], 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)