The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

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 \(\phi\colon V^k\longrightarrow\mathbb{R}\) and a vector \(\mathbf{v}\in V\), the contraction \(\phi_\mathbf{v}\) of \(\phi\) and \(\mathbf{v}\) is a \(k-1\)-form with

\[ \phi_\mathbf{v}\left(\mathbf{v}^1,\ldots,\mathbf{v}^{k-1}\right) = \phi\left(\mathbf{v},\mathbf{v}^1,\ldots,\mathbf{v}^{k-1}\right) \]

provided \(k>1\); if \(k=1\) we specify \(\phi_\mathbf{v}=\phi(\mathbf{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 \(\mathbf{v}\); function contract() is the user-friendly front end. We will start with some simple examples. I will use phi and \(\phi\) 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 \(\phi=\mathrm{d}x^1\wedge\mathrm{d}x^2\wedge \mathrm{d}x^3\wedge\mathrm{d}x^4\wedge\mathrm{d}x^5\). We have that \(\phi\) 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 \(\phi\) with any vector \(\mathbf{v}\) is thus a \(4\)-form mapping \(V^4\) to the reals with \(\phi_\mathbf{v}\left(\mathbf{v}^1,\mathbf{v}^2,\mathbf{v}^3,\mathbf{v}^4\right)=\phi\left(\mathbf{v},\mathbf{v}^1,\mathbf{v}^2,\mathbf{v}^3,\mathbf{v}^4\right)\). Taking the simplest case first, if \(\mathbf{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 \(V^4\) 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 \(\phi\) understood to be \(\mathbf{v}=(1,0,0,0,0)\)). Now consider \(\mathbf{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

\[\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:

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: \(\phi_{a\mathbf{v}+b\mathbf{w}}= a\phi_\mathbf{v}+b\phi_\mathbf{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\phi+b\psi)_\mathbf{v}=a\phi_\mathbf{v} + b\psi_\mathbf{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 \(\phi\) and \(l\)-form \(\psi\):

\[ \left(\phi\wedge\psi\right)_\mathbf{v} = \phi_\mathbf{v}\psi + (-1)^k\phi\wedge\psi_\mathbf{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 \(\phi\) and \(\psi\), 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 \(\phi\) is a \(k\)-form, then \(\left(\phi_\mathbf{v}\right)_\mathbf{w}\) is a \(k-2\) form with

\[ \left(\phi_\mathbf{u}\right)_\mathbf{v}\left(\mathbf{w}^1,\ldots,\mathbf{w}^{k-2}\right)=\phi\left(\mathbf{u},\mathbf{v},\mathbf{w}^1,\ldots,\mathbf{w}^{k-2}\right) \]

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.355108 6.355108 0.000000

Above we see agreement to within numerical error.

Contraction from first principles

Suppose we wish to contract \(\phi=dx^{i_1}\wedge\cdots\wedge dx^{i_k}\) with vector \(\mathbf{v}=(v_1\mathbf{e}_1,\ldots,v_k\mathbf{e}_k)\). Thus we seek \(\phi_\mathbf{v}\) with \(\phi_\mathbf{v}\left(\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right) = dx^{i_1}\wedge\cdots\wedge dx^{i_k}\left(\mathbf{v},\mathbf{v}_1,\ldots\mathbf{v}_{k-1}\right)\). Writing \(\mathbf{v}=v_1\mathbf{e}_1+\cdots+\mathbf{e}_k\), we have

\[\begin{eqnarray} \phi_\mathbf{v}\left(\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right) &=& dx^{i_1}\wedge\cdots\wedge dx^{i_k}\left(\mathbf{v},\mathbf{v}_1,\ldots\mathbf{v}_{k-1}\right)\\&=& dx^{i_1}\wedge\cdots\wedge dx^{i_k}\left(v_1\mathbf{e}_1+\cdots+v_k\mathbf{e}_k,\mathbf{v}_1,\ldots\mathbf{v}_{k-1}\right)\\&=& v_1 dx^{i_1}\wedge\cdots\wedge dx^{i_k}\left(\mathbf{e}_1,\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right)+\cdots+ v_k dx^{i_1}\wedge\cdots\wedge dx^{i_k}\left(\mathbf{e}_k,\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right). \end{eqnarray}\]

where we have exploited linearity. To evaluate this it is easiest and most efficient to express \(\phi\) as \(\bigwedge_{j=1}^kdx^{i_j}\) and cycle through the index \(j\), and use various properties of wedge products:

\[\begin{eqnarray} dx^{i_1}\wedge\cdots\wedge dx^{i_k} &=& (-1)^{j-1} dx^{i_j}\wedge\left(dx^{i_1}\wedge\cdots\wedge\widehat{dx^{i_j}}\wedge\cdots\wedge dx^{i-k}\right)\\ &=& (-1)^{j-1} k\operatorname{Alt}\left(dx^{i_j}\otimes\left(dx^{i_1}\wedge\cdots\wedge\widehat{dx^{i_j}}\wedge\cdots\wedge dx^{i-k}\right)\right) \end{eqnarray}\]

(above, a hat indicates a term’s being omitted). From this, we see that \(l\not\in L\longrightarrow\phi=0\) (where \(L=\left\lbrace i_1,\ldots i_k\right\rbrace\) is the index set of \(\phi\)), for all the \(dx^p\) terms kill \(\mathbf{e}_l\). On the other hand, if \(l\in L\) we have

\[\begin{eqnarray} \phi_{\mathbf{e}_l}\left(\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right) &=& \left(dx^{l}\wedge\left(dx^{i_1}\wedge\cdots\wedge\widehat{dx^{l}}\wedge\cdots\wedge dx^{i_k}\right)\right)\left(\mathbf{e}_l,\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right)\\ &=& (-1)^{l-1}k\left(dx^{l}\otimes\left(dx^{i_1}\wedge\cdots\wedge\widehat{dx^{l}}\wedge\cdots\wedge dx^{i_k}\right)\right)\left(\mathbf{e}_l,\left(\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right)\right)\\ &=& (-1)^{l-1}k\left(dx^{i_1}\wedge\cdots\wedge\widehat{dx^{l}}\wedge\cdots\wedge dx^{i_k}\right)\left(\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right) \end{eqnarray}\]

Worked example using contract_elementary()

Function contract_elementary() is a bare-bones low-level no-frills helper function that returns \(\phi_\mathbf{v}\) for \(\phi\) an elementary form of the form \(dx^{i_1}\wedge\cdots\wedge dx^{i_k}\). Suppose we wish to contract \(\phi=dx^1\wedge dx^2\wedge dx^5\) with vector \(\mathbf{v}=(1,2,10,11,71)^T\).

Thus we seek \(\phi_\mathbf{v}\) with \(\phi_\mathbf{v}\left(\mathbf{v}_1,\mathbf{v}_2 \right)=dx^1\wedge dx^2\wedge dx^5\left(\mathbf{v},\mathbf{v}_1,\mathbf{v}_2\right)\). Writing \(\mathbf{v}=v_1\mathbf{e}_1+\cdots+v_5\mathbf{e}_5\) we have

\[\begin{eqnarray} \phi_\mathbf{v}\left(\mathbf{v}_1,\mathbf{v}_2 \right) &=& dx^1\wedge dx^2\wedge dx^5\left(\mathbf{v},\mathbf{v}_1,\mathbf{v}_2\right)\\ &=& dx^1\wedge dx^2\wedge dx^5\left(v_1\mathbf{e}_1+\cdots+v_5\mathbf{e}_5,\mathbf{v}_1,\mathbf{v}_2\right)\\&=& v_1 dx^1\wedge dx^2\wedge dx^5\left(\mathbf{e}_1,\mathbf{v}_1,\mathbf{v}_2\right)+ v_2 dx^1\wedge dx^2\wedge dx^5\left(\mathbf{e}_2,\mathbf{v}_1,\mathbf{v}_2\right)\\ &{}&\qquad +v_3dx^1\wedge dx^2\wedge dx^5\left(\mathbf{e}_3,\mathbf{v}_1,\mathbf{v}_2\right)+ v_4 dx^1\wedge dx^2\wedge dx^5\left(\mathbf{e}_4,\mathbf{v}_1,\mathbf{v}_2\right)\\ &{}&\qquad\qquad +v_5dx^1\wedge dx^2\wedge dx^5\left(\mathbf{e}_5,\mathbf{v}_1,\mathbf{v}_2\right)\\&=& v_1 dx^2\wedge dx^5\left(\mathbf{v}_1,\mathbf{v}_2\right)- v_2 dx^1\wedge dx^5\left(\mathbf{v}_1,\mathbf{v}_2\right)+0+0+ v_5 dx^1\wedge dx^2\left(\mathbf{v}_1,\mathbf{v}_2\right) \end{eqnarray}\]

(above, the zero terms are because the vectors \(\mathbf{e}_3\) and \(\mathbf{e}_4\) are killed by \(dx^1\wedge dx^2\wedge dx^5\)). We can see that the way to evaluate the contraction is to go through the terms of \(\phi\) [that is, \(dx^1\), \(dx^2\), and \(dx^5\)] 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.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.