# contract() and related functions in the stokes package In stokes: The Exterior Calculus

knitr::opts_chunk$set(echo = TRUE) options(rmarkdown.html_vignette.check_title = FALSE) library("stokes") library("spray") library("disordR") library("magrittr") set.seed(0) ![](r system.file("help/figures/stokes.png", package = "stokes")){width=10%} contract contract_elementary ## Contractions 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})$. 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)) Thus$k=5$and we have$\phi=dx^1\wedge dx^2\wedge dx^3\wedge dx^4\wedge dx^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) 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) $$\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)) contract(phi,1:5) 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) 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)) (psi <- rform(2,5)) a <- 7 b <- 13 v <- 1:7 contract(a*phi + b*psi,v) == a*contract(phi,v) + b*contract(psi,v) ## Contraction of products Weintraub 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) The theorem is verified. We note in passing that the object itself is quite complicated: summary(contract(phi^psi,v)) We may also switch$\phi$and$\psi$, remembering to change the sign: contract(psi^phi,v) == contract(psi,v) ^ phi + psi ^ contract(phi,v) ## 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)) u <- c(1,3,2,4,5,4,6) v <- c(8,6,5,3,4,3,2) contract(contract(phi,u),v) 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) 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)))) 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) max(jj) - min(jj) # zero to numerical precision 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) 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) 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) 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]]
)

This is performed more succinctly by contract_elementary():

contract_elementary(o,v)

## 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)))
v <- 1:7

Then the inside bit would be

apply(index(K), 1, contract_elementary, v)

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)))

And finally use Reduce() to sum the terms:

Reduce("+",Map("*", apply(index(K), 1, contract_elementary, v), elements(coeffs(K))))

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("+",.)

Well it might be clearer to Hadley but frankly YMMV.

## References

• M. Spivak 1968. Calculus on manifolds. Perseus Books
• S. T. Weintraub 2014. Differential forms: theory and practice. Elsevier.

## Try the stokes package in your browser

Any scripts or data that you put into this service are public.

stokes documentation built on Aug. 19, 2023, 1:07 a.m.