calculus/vignettes/integrals.Rmd at master · eguidotti/calculus · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
---
title: "Multiple integrals in arbitrary orthogonal coordinates systems"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Multiple integrals in arbitrary orthogonal coordinates systems}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(calculus)
```

The package integrates seamlessly with [cubature](https://cran.r-project.org/package=cubature) for efficient numerical integration in `C`. The function [`integral`](https://calculus.eguidotti.com/reference/integral.html) provides the interface for multidimensional integrals of `functions`, `expressions`, and `characters` in arbitrary [orthogonal coordinate systems](https://calculus.eguidotti.com/articles/differential-operators.html). If the package [cubature](https://cran.r-project.org/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$.

## Examples

Univariate integral $\int_0^1xdx$:

```{r}
i <- integral(f = "x", bounds = list(x = c(0,1)))
i$value
```

that is equivalent to:

```{r}
i <- integral(f = function(x) x, bounds = list(x = c(0,1)))
i$value
```

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

```{r}
i <- integral(f = "y*x", bounds = list(x = c(0,1), y = 2))
i$value
```

Multivariate integral $\int_0^1\int_o^1yxdxdy$:

```{r}
i <- integral(f = "y*x", bounds = list(x = c(0,1), y = c(0,1)))
i$value
```

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

```{r}
i <- integral(f = 1,
              bounds = list(r = c(0,1), theta = c(0,2*pi)),
              coordinates = "polar")
i$value
```

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

```{r}
i <- integral(f = 1,
              bounds = list(r = c(0,1), theta = c(0,pi), phi = c(0,2*pi)),
              coordinates = "spherical")
i$value
```

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

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

The electric field is determined by the gradient of the potential^[https://en.wikipedia.org/wiki/Electric_potential] $E = -\nabla V$:

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

Then, by Gauss's law^[https://en.wikipedia.org/wiki/Gauss%27s_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$:

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

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](https://doi.org/10.18637/jss.v104.i05)

A BibTeX entry for LaTeX users is

```bibtex
@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},
}
```