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 q1…qn the integral is computed
as:
∫J⋅f(q1…qn)dq1…dqn
where J=∏ihi is the Jacobian determinant of the transformation and is equal to the product of the scale factors h1…hn.
Univariate integral ∫10xdx:
<- 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 ∫10yxdx|y=2:
<- integral(f = "y*x", bounds = list(x = c(0,1), y = 2))
i $value
i#> [1] 1
Multivariate integral ∫10∫1oyxdxdy:
<- integral(f = "y*x", bounds = list(x = c(0,1), y = c(0,1)))
i $value
i#> [1] 0.25
Area of a circle ∫2π0∫10dA(r,θ)
<- 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 ∫π0∫2π0∫10dV(r,θ,ϕ)
<- 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=14πr arising from a unitary point charge:
<- "1/(4*pi*r)" V
The electric field is determined by the gradient of the potential1 E=−∇V:
<- -1 %prod% gradient(V, c("r","theta","phi"), coordinates = "spherical") E
Then, by Gauss’s law2, the total charge enclosed within a given volume is equal to the surface integral of the electric field q=∫E⋅dA where ⋅ 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 ∫ErdA. 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=∫E⋅dA=∫ErdA=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},
}