contract()
and
contract_elementary()
in the stokes
packagefunction (K, v, lose = TRUE)
{
if (is.vector(v)) {
out <- Reduce("+", Map("*", apply(index(K), 1, contract_elementary,
v), elements(coeffs(K))))
}
else {
stopifnot(is.matrix(v))
out <- K
for (i in seq_len(ncol(v))) {
out <- contract(out, v[, i, drop = TRUE], lose = FALSE)
}
}
if (lose) {
out <- lose(out)
}
return(disordR::drop(out))
}
function (o, v)
{
out <- zeroform(length(o) - 1)
for (i in seq_along(o)) {
out <- out + (-1)^(i + 1) * v[o[i]] * as.kform(rbind(o[-i]),
lose = FALSE)
}
return(out)
}
To cite the stokes
package in publications, please use
Hankin (2022). Given a k-form ϕ: Vk → ℝ
and a vector v ∈ V, the
contraction ϕv of
ϕ and v is a k − 1-form with
ϕv(v1, …, vk − 1) = ϕ(v, v1, …, vk − 1)
provided k > 1; if k = 1 we specify ϕv = ϕ(v).
If Spivak (1965) does discuss this, I have
forgotten it. Function contract_elementary()
is a low-level
helper function that translates elementary k-forms with coefficient 1 (in the
form of an integer vector corresponding to one row of an index matrix)
into its contraction with v; function
contract()
is the user-friendly front end. We will start
with some simple examples. I will use phi
and ϕ to represent the same object.
## An alternating linear map from V^5 to R with V=R^5:
## val
## 1 2 3 4 5 = 1
Thus k = 5 and we have ϕ = dx1 ∧ dx2 ∧ dx3 ∧ dx4 ∧ dx5. We have that ϕ is a linear alternating map with
$$\phi\left(\begin{bmatrix}1\\0\\0\\0\\0\end{bmatrix}, \begin{bmatrix}0\\1\\0\\0\\0\end{bmatrix}, \begin{bmatrix}0\\0\\1\\0\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\1\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\0\\1\end{bmatrix} \right)=1$$.
The contraction of ϕ with any vector v is thus a 4-form mapping V4 to the reals with ϕv(v1, v2, v3, v4) = ϕ(v, v1, v2, v3, v4). Taking the simplest case first, if v = (1, 0, 0, 0, 0) then
## An alternating linear map from V^4 to R with V=R^5:
## val
## 2 3 4 5 = 1
that is, a linear alternating map from V4 to the reals with
$$\phi_\mathbf{v}\left( \begin{bmatrix}0\\1\\0\\0\\0\end{bmatrix}, \begin{bmatrix}0\\0\\1\\0\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\1\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\0\\1\end{bmatrix}\right)=1$$.
(the contraction has the first argument of ϕ understood to be v = (1, 0, 0, 0, 0)). Now consider w = (0, 1, 0, 0, 0):
## An alternating linear map from V^4 to R with V=R^5:
## val
## 1 3 4 5 = -1
$$\phi_\mathbf{w}\left( \begin{bmatrix}0\\0\\1\\0\\0\end{bmatrix}, \begin{bmatrix}1\\0\\0\\0\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\1\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\0\\1\end{bmatrix}\right)=1 \qquad\mbox{or}\qquad \phi_\mathbf{w}\left( \begin{bmatrix}1\\0\\0\\0\\0\end{bmatrix}, \begin{bmatrix}0\\0\\1\\0\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\1\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\0\\1\end{bmatrix}\right)=-1$$.
Contraction is linear, so we may use more complicated vectors:
## An alternating linear map from V^4 to R with V=R^5:
## val
## 2 3 4 5 = 1
## 1 3 4 5 = -3
## An alternating linear map from V^4 to R with V=R^5:
## val
## 1 2 3 4 = 5
## 1 2 4 5 = 3
## 2 3 4 5 = 1
## 1 3 4 5 = -2
## 1 2 3 5 = -4
We can check numerically that the contraction is linear in its vector argument: ϕav + bw = aϕv + bϕw.
a <- 1.23
b <- -0.435
v <- 1:5
w <- c(-3, 2.2, 1.1, 2.1, 1.8)
contract(phi,a*v + b*w) == a*contract(phi,v) + b*contract(phi,w)
## [1] TRUE
We also have linearity in the alternating form: (aϕ + bψ)v = aϕv + bψv.
## An alternating linear map from V^5 to R with V=R^7:
## val
## 2 3 4 5 7 = -2
## 1 3 4 6 7 = 1
## An alternating linear map from V^5 to R with V=R^7:
## val
## 1 2 3 6 7 = 2
## 2 3 5 6 7 = 1
## [1] TRUE
Weintraub (2014) gives us the following theorem, for any k-form ϕ and l-form ψ:
(ϕ ∧ ψ)v = ϕvψ + (−1)kϕ ∧ ψv.
We can verify this numerically with k = 4, l = 5:
phi <- rform(terms=5,k=3,n=9)
psi <- rform(terms=9,k=4,n=9)
v <- sample(1:100,9)
contract(phi^psi,v) == contract(phi,v) ^ psi - phi ^ contract(psi,v)
## [1] TRUE
The theorem is verified. We note in passing that the object itself is quite complicated:
## A kform object with 47 terms. Summary of coefficients:
##
## a disord object with hash d6cdb7213e60a9d847f1752a839f18b1de98bc57
##
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2943.00 -516.00 48.00 44.47 768.00 2625.00
##
##
## Representative selection of index and coefficients:
##
## An alternating linear map from V^6 to R with V=R^9:
## val
## 1 2 4 6 8 9 = 390
## 1 2 3 5 6 7 = 420
## 1 2 3 4 6 8 = -840
## 2 5 6 7 8 9 = 1605
## 1 2 3 6 7 9 = 355
## 1 2 3 6 8 9 = -1200
We may also switch ϕ and ψ, remembering to change the sign:
## [1] TRUE
It is of course possible to contract a contraction. If ϕ is a k-form, then (ϕv)w is a k − 2 form with
(ϕu)v(w1, …, wk − 2) = ϕ(u, v, w1, …, wk − 2)
And this is straightforward to realise in the package:
## An alternating linear map from V^5 to R with V=R^7:
## val
## 1 4 5 6 7 = -2
## 2 4 5 6 7 = -1
## An alternating linear map from V^3 to R with V=R^7:
## val
## 2 5 6 = 10
## 1 4 7 = 2
## 4 5 7 = 73
## 1 4 5 = 20
## 2 4 7 = 1
## 1 5 6 = 20
## 2 4 6 = -14
## 1 6 7 = -2
## 2 6 7 = -1
## 2 4 5 = 10
## 5 6 7 = 73
## 4 6 7 = -90
## 1 4 6 = -28
## 4 5 6 = -122
But contract()
allows us to perform both contractions in
one operation: if we pass a matrix M to contract()
then
this is interpreted as repeated contraction with the columns of M:
## [1] TRUE
We can verify directly that the system works as intended. The lines
below strip successively more columns from argument V
and
contract with them:
## An alternating linear map from V^4 to R with V=R^9:
## val
## 3 7 8 9 = -0.1482116
## 1 5 6 7 = 0.4314737
V <- matrix(rnorm(36),ncol=4)
jj <- c(
as.function(o)(V),
as.function(contract(o,V[,1,drop=TRUE]))(V[,-1]), # scalar
as.function(contract(o,V[,1:2]))(V[,-(1:2),drop=FALSE]),
as.function(contract(o,V[,1:3]))(V[,-(1:3),drop=FALSE]),
as.function(contract(o,V[,1:4],lose=FALSE))(V[,-(1:4),drop=FALSE])
)
print(jj)
## [1] -0.4992204 -0.4992204 -0.4992204 -0.4992204 -0.4992204
## [1] 2.775558e-16
and above we see agreement to within numerical precision. If we pass
three columns to contract()
the result is a 0-form:
## [1] -0.4992204
In the above, the result is coerced to a scalar which is returned in
the form of a disord
object; in order to work with a formal
0-form (which is represented in the
package as a spray
with a zero-column index matrix) we can
use the lost=FALSE
argument:
## An alternating linear map from V^0 to R with V=R^0:
## val
## = -0.4992204
thus returning a 0-form. If we iteratively contract a k-dimensional k-form, we return the determinant, and this may be verified as follows:
o <- as.kform(1:5)
V <- matrix(rnorm(25),5,5)
LHS <- det(V)
RHS <- contract(o,V)
c(LHS=LHS,RHS=RHS,diff=LHS-RHS)
## LHS RHS diff
## 6.355108e+00 6.355108e+00 1.776357e-15
Above we see agreement to within numerical error.
Suppose we wish to contract ϕ = dxi1 ∧ ⋯ ∧ dxik with vector v = (v1e1, …, vkek). Thus we seek ϕv with ϕv(v1, …, vk − 1) = dxi1 ∧ ⋯ ∧ dxik(v, v1, …, vk − 1). Writing v = v1e1 + ⋯ + ek, we have
where we have exploited linearity. To evaluate this it is easiest and most efficient to express ϕ as $\bigwedge_{j=1}^kdx^{i_j}$ and cycle through the index j, and use various properties of wedge products:
(above, a hat indicates a term’s being omitted). From this, we see that l ∉ L → ϕ = 0 (where L = {i1, …ik} is the index set of ϕ), for all the dxp terms kill el. On the other hand, if l ∈ L we have
contract_elementary()
Function contract_elementary()
is a bare-bones low-level
no-frills helper function that returns ϕv for
ϕ an elementary form of the
form dxi1 ∧ ⋯ ∧ dxik.
Suppose we wish to contract ϕ = dx1 ∧ dx2 ∧ dx5
with vector v = (1, 2, 10, 11, 71)T.
Thus we seek ϕv with ϕv(v1, v2) = dx1 ∧ dx2 ∧ dx5(v, v1, v2). Writing v = v1e1 + ⋯ + v5e5 we have
(above, the zero terms are because the vectors e3 and e4 are killed by dx1 ∧ dx2 ∧ dx5). We can see that the way to evaluate the contraction is to go through the terms of ϕ [that is, dx1, dx2, and dx5] in turn, and sum the resulting expressions:
o <- c(1,2,5)
v <- c(1,2,10,11,71)
(
(-1)^(1+1) * as.kform(o[-1])*v[o[1]] +
(-1)^(2+1) * as.kform(o[-2])*v[o[2]] +
(-1)^(3+1) * as.kform(o[-3])*v[o[3]]
)
## An alternating linear map from V^2 to R with V=R^5:
## val
## 1 5 = -2
## 2 5 = 1
## 1 2 = 71
This is performed more succinctly by
contract_elementary()
:
## An alternating linear map from V^2 to R with V=R^5:
## val
## 1 5 = -2
## 2 5 = 1
## 1 2 = 71
contract()
Given a vector v
, and kform
object
K
, the meat of contract()
would be
Reduce("+", Map("*", apply(index(K), 1, contract_elementary, v), elements(coeffs(K))))
I will show this in operation with simple but nontrivial arguments.
## An alternating linear map from V^4 to R with V=R^7:
## val
## 2 4 5 7 = 2
## 1 2 3 6 = 1
Then the inside bit would be
## [[1]]
## An alternating linear map from V^3 to R with V=R^7:
## val
## 2 4 7 = 5
## 4 5 7 = 2
## 2 5 7 = -4
## 2 4 5 = -7
##
## [[2]]
## An alternating linear map from V^3 to R with V=R^6:
## val
## 1 2 6 = 3
## 2 3 6 = 1
## 1 3 6 = -2
## 1 2 3 = -6
Above we see a two-element list, one for each elementary term of
K
. We then use base R’s Map()
function to
multiply each one by the appropriate coefficient:
## [[1]]
## An alternating linear map from V^3 to R with V=R^7:
## val
## 2 4 5 = -14
## 2 5 7 = -8
## 4 5 7 = 4
## 2 4 7 = 10
##
## [[2]]
## An alternating linear map from V^3 to R with V=R^6:
## val
## 1 2 3 = -6
## 1 3 6 = -2
## 2 3 6 = 1
## 1 2 6 = 3
And finally use Reduce()
to sum the terms:
## An alternating linear map from V^3 to R with V=R^7:
## val
## 2 4 7 = 10
## 4 5 7 = 4
## 2 5 7 = -8
## 1 2 3 = -6
## 2 4 5 = -14
## 1 3 6 = -2
## 2 3 6 = 1
## 1 2 6 = 3
However, it might be conceptually easier to use magrittr
pipes to give an equivalent definition:
K %>%
index %>%
apply(1,contract_elementary,v) %>%
Map("*", ., K %>% coeffs %>% elements) %>%
Reduce("+",.)
## An alternating linear map from V^3 to R with V=R^7:
## val
## 2 4 7 = 10
## 4 5 7 = 4
## 2 5 7 = -8
## 1 2 3 = -6
## 2 4 5 = -14
## 1 3 6 = -2
## 2 3 6 = 1
## 1 2 6 = 3
Well it might be clearer to Hadley but frankly YMMV.