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⟶R
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
ϕ([10000],[01000],[00100],[00010],[00001])=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
ϕv([01000],[00100],[00010],[00001])=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
ϕw([00100],[10000],[00010],[00001])=1orϕw([10000],[00100],[00010],[00001])=−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
ϕv(v1,…,vk−1)=dxi1∧⋯∧dxik(v,v1,…,vk−1)=dxi1∧⋯∧dxik(v1e1+⋯+vkek,v1,…,vk−1)=v1dxi1∧⋯∧dxik(e1,v1,…,vk−1)+⋯+vkdxi1∧⋯∧dxik(ek,v1,…,vk−1).
where we have exploited linearity. To evaluate this it is easiest and most efficient to express ϕ as ⋀kj=1dxij and cycle through the index j, and use various properties of wedge products:
dxi1∧⋯∧dxik=(−1)j−1dxij∧(dxi1∧⋯∧^dxij∧⋯∧dxi−k)=(−1)j−1kAlt(dxij⊗(dxi1∧⋯∧^dxij∧⋯∧dxi−k))
(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
ϕel(v1,…,vk−1)=(dxl∧(dxi1∧⋯∧^dxl∧⋯∧dxik))(el,v1,…,vk−1)=(−1)l−1k(dxl⊗(dxi1∧⋯∧^dxl∧⋯∧dxik))(el,(v1,…,vk−1))=(−1)l−1k(dxi1∧⋯∧^dxl∧⋯∧dxik)(v1,…,vk−1)
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
ϕv(v1,v2)=dx1∧dx2∧dx5(v,v1,v2)=dx1∧dx2∧dx5(v1e1+⋯+v5e5,v1,v2)=v1dx1∧dx2∧dx5(e1,v1,v2)+v2dx1∧dx2∧dx5(e2,v1,v2)+v3dx1∧dx2∧dx5(e3,v1,v2)+v4dx1∧dx2∧dx5(e4,v1,v2)+v5dx1∧dx2∧dx5(e5,v1,v2)=v1dx2∧dx5(v1,v2)−v2dx1∧dx5(v1,v2)+0+0+v5dx1∧dx2(v1,v2)
(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.