`contract()` and related functions in the `stokes` package

knitr::opts_chunk$set(echo = TRUE)
options(rmarkdown.html_vignette.check_title = FALSE)

![](`r system.file("help/figures/stokes.png", package = "stokes")`){width=10%}



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)

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)

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


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:


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)

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(contract(o,V[,1,drop=TRUE]))(V[,-1]), # scalar
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:


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:


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)

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


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

Well it might be clearer to Hadley but frankly YMMV.


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.