Multiple integrals in arbitrary orthogonal coordinates systems

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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:

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\).

Examples

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

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

that is equivalent to:

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

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

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

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

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

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

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

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

i <- integral(f = 1, 
              bounds = list(r = c(0,1), theta = c(0,pi), phi = c(0,2*pi)), 
              coordinates = "spherical")
i$value
#> [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:

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

The electric field is determined by the gradient of the potential¹ \(E = -\nabla V\):

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

Then, by Gauss’s law², 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\):

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

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

Cite as

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},
}

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.