Processing math: 100%

Functions contract() and contract_elementary() in the stokes package

Robin K. S. Hankin

contract
function (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))
}
contract_elementary
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 ϕ:VkR and a vector vV, the contraction ϕv of ϕ and v is a k1-form with

ϕv(v1,,vk1)=ϕ(v,v1,,vk1)

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.

(phi <- as.kform(1:5))
## 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 ϕ=dx1dx2dx3dx4dx5. 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

v <- c(1,0,0,0,0)
contract(phi,v)
## 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):

w <- c(0,1,0,0,0)
contract(phi,w)
## 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:

contract(phi,c(1,3,0,0,0))
## An alternating linear map from V^4 to R with V=R^5:
##              val
##  2 3 4 5  =    1
##  1 3 4 5  =   -3
contract(phi,1:5)
## 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.

(phi <- rform(2,5))
## 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
(psi <- rform(2,5))
## 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
a <- 7
b <- 13
v <- 1:7
contract(a*phi + b*psi,v) == a*contract(phi,v) + b*contract(psi,v)
## [1] TRUE

Contraction of products

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:

summary(contract(phi^psi,v))
## 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:

contract(psi^phi,v) ==  contract(psi,v) ^ phi + psi ^ contract(phi,v)
## [1] TRUE

Repeated contraction

It is of course possible to contract a contraction. If ϕ is a k-form, then (ϕv)w is a k2 form with

(ϕu)v(w1,,wk2)=ϕ(u,v,w1,,wk2)

And this is straightforward to realise in the package:

(phi <- rform(2,5))
## 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
u <- c(1,3,2,4,5,4,6)
v <- c(8,6,5,3,4,3,2)
contract(contract(phi,u),v)
## 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:

M <- cbind(u,v)
contract(contract(phi,u),v) == contract(phi,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:

(o <- kform(spray(t(replicate(2, sample(9,4))), runif(2))))
## 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
max(jj) - min(jj) # zero to numerical precision
## [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:

contract(o,V)
## [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:

contract(o,V,lose=FALSE)
## 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.

Contraction from first principles

Suppose we wish to contract ϕ=dxi1dxik with vector v=(v1e1,,vkek). Thus we seek ϕv with ϕv(v1,,vk1)=dxi1dxik(v,v1,,vk1). Writing v=v1e1++ek, we have

ϕv(v1,,vk1)=dxi1dxik(v,v1,,vk1)=dxi1dxik(v1e1++vkek,v1,,vk1)=v1dxi1dxik(e1,v1,,vk1)++vkdxi1dxik(ek,v1,,vk1).

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:

dxi1dxik=(1)j1dxij(dxi1^dxijdxik)=(1)j1kAlt(dxij(dxi1^dxijdxik))

(above, a hat indicates a term’s being omitted). From this, we see that lLϕ=0 (where L={i1,ik} is the index set of ϕ), for all the dxp terms kill el. On the other hand, if lL we have

ϕel(v1,,vk1)=(dxl(dxi1^dxldxik))(el,v1,,vk1)=(1)l1k(dxl(dxi1^dxldxik))(el,(v1,,vk1))=(1)l1k(dxi1^dxldxik)(v1,,vk1)

Worked example using 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 dxi1dxik. Suppose we wish to contract ϕ=dx1dx2dx5 with vector v=(1,2,10,11,71)T.

Thus we seek ϕv with ϕv(v1,v2)=dx1dx2dx5(v,v1,v2). Writing v=v1e1++v5e5 we have

ϕv(v1,v2)=dx1dx2dx5(v,v1,v2)=dx1dx2dx5(v1e1++v5e5,v1,v2)=v1dx1dx2dx5(e1,v1,v2)+v2dx1dx2dx5(e2,v1,v2)+v3dx1dx2dx5(e3,v1,v2)+v4dx1dx2dx5(e4,v1,v2)+v5dx1dx2dx5(e5,v1,v2)=v1dx2dx5(v1,v2)v2dx1dx5(v1,v2)+0+0+v5dx1dx2(v1,v2)

(above, the zero terms are because the vectors e3 and e4 are killed by dx1dx2dx5). 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():

contract_elementary(o,v)
## An alternating linear map from V^2 to R with V=R^5:
##          val
##  1 5  =   -2
##  2 5  =    1
##  1 2  =   71

The “meat” of 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.

(K <- as.kform(spray(matrix(c(1,2,3,6,2,4,5,7),2,4,byrow=TRUE),1:2)))
## An alternating linear map from V^4 to R with V=R^7:
##              val
##  2 4 5 7  =    2
##  1 2 3 6  =    1
v <- 1:7

Then the inside bit would be

apply(index(K), 1, contract_elementary, v)
## [[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:

Map("*", apply(index(K), 1, contract_elementary, v), elements(coeffs(K)))
## [[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:

Reduce("+",Map("*", apply(index(K), 1, contract_elementary, v), elements(coeffs(K))))
## 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.

References

Hankin, R. K. S. 2022. “Stokes’s Theorem in R.” arXiv. https://doi.org/10.48550/ARXIV.2210.17008.
Spivak, M. 1965. Calculus on Manifolds. Addison-Wesley.
Weintraub, S. T. 2014. Differential Forms: Theory and Practice. Elsevier.