`library(calculus)`

The package integrates seamlessly with cubature for
efficient numerical integration in `C`

. The function `integral`

provides the interface for multidimensional integrals of
`functions`

, `expressions`

, and
`characters`

in arbitrary orthogonal
coordinate systems. If the package cubature is not
installed, the package implements a naive Monte Carlo integration by
default. The function returns a `list`

containing the
`value`

of the integral as well as other information on the
estimation uncertainty. The integration bounds are specified via the
argument `bounds`

: a list containing the lower and upper
bound for each variable. If the two bounds coincide, or if a single
number is specified, the corresponding variable is not integrated and
its value is fixed. For arbitrary orthogonal coordinates \(q_1\dots q_n\) the integral is computed
as:

\[ \int J\cdot f(q_1\dots q_n) dq_1\dots dq_n \]

where \(J=\prod_i h_i\) is the Jacobian determinant of the transformation and is equal to the product of the scale factors \(h_1\dots h_n\).

Univariate integral \(\int_0^1xdx\):

```
<- integral(f = "x", bounds = list(x = c(0,1)))
i $value
i#> [1] 0.5
```

that is equivalent to:

```
<- integral(f = function(x) x, bounds = list(x = c(0,1)))
i $value
i#> [1] 0.5
```

Univariate integral \(\int_0^1yxdx|_{y=2}\):

```
<- integral(f = "y*x", bounds = list(x = c(0,1), y = 2))
i $value
i#> [1] 1
```

Multivariate integral \(\int_0^1\int_o^1yxdxdy\):

```
<- integral(f = "y*x", bounds = list(x = c(0,1), y = c(0,1)))
i $value
i#> [1] 0.25
```

Area of a circle \(\int_0^{2\pi}\int_0^1dA(r,\theta)\)

```
<- integral(f = 1,
i bounds = list(r = c(0,1), theta = c(0,2*pi)),
coordinates = "polar")
$value
i#> [1] 3.141593
```

Volume of a sphere \(\int_0^\pi\int_0^{2\pi}\int_0^1dV(r,\theta,\phi)\)

```
<- integral(f = 1,
i bounds = list(r = c(0,1), theta = c(0,pi), phi = c(0,2*pi)),
coordinates = "spherical")
$value
i#> [1] 4.188794
```

As a final example consider the electric potential in spherical coordinates \(V=\frac{1}{4\pi r}\) arising from a unitary point charge:

`<- "1/(4*pi*r)" V `

The electric field is determined by the gradient of the potential^{1} \(E = -\nabla V\):

`<- -1 %prod% gradient(V, c("r","theta","phi"), coordinates = "spherical") E `

Then, by Gauss’s law^{2}, the total charge enclosed within a given
volume is equal to the surface integral of the electric field \(q=\int E\cdot dA\) where \(\cdot\) denotes the scalar product between
the two vectors. In spherical coordinates, this reduces to the surface
integral of the radial component of the electric field \(\int E_rdA\). The following code computes
this surface integral on a sphere with fixed radius \(r=1\):

```
<- integral(E[1],
i bounds = list(r = 1, theta = c(0,pi), phi = c(0,2*pi)),
coordinates = "spherical")
$value
i#> [1] 1.000002
```

As expected \(q=\int E\cdot dA=\int E_rdA=1\), the unitary charge generating the electric potential.

Guidotti E (2022). “calculus: High-Dimensional Numerical and Symbolic
Calculus in R.” *Journal of Statistical Software*,
*104*(5), 1-37. doi:10.18637/jss.v104.i05

A BibTeX entry for LaTeX users is

```
@Article{calculus,
title = {{calculus}: High-Dimensional Numerical and Symbolic Calculus in {R}},
author = {Emanuele Guidotti},
journal = {Journal of Statistical Software},
year = {2022},
volume = {104},
number = {5},
pages = {1--37},
doi = {10.18637/jss.v104.i05},
}
```