set.seed(0) library("stokes") library("spray") library("disordR") library("magrittr") knitr::opts_chunk$set(echo = TRUE) options(rmarkdown.html_vignette.check_title = FALSE) knit_print.function <- function(x, ...){dput(x)} registerS3method( "knit_print", "function", knit_print.function, envir = asNamespace("knitr") )
knitr::include_graphics(system.file("help/figures/stokes.png", package = "stokes"))
contract contract_elementary
To cite the stokes
package in publications, please use
@hankin2022_stokes. 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 @spivak1965 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))
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)
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)
@weintraub2014 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)
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.
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}
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)
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.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.