Loading [MathJax]/jax/output/HTML-CSS/jax.js

Function vector_cross_product() in the stokes package

Robin K. S. Hankin

vector_cross_product
## function (M) 
## {
##     n <- nrow(M)
##     stopifnot(n == ncol(M) + 1)
##     (-1)^n * sapply(seq_len(n), function(i) {
##         (-1)^i * det(M[-i, ])
##     })
## }
## <bytecode: 0x55abbc5b0db0>
## <environment: namespace:stokes>

The vector cross product

Given two vectors u,vR3 with u=(u1,u2,u3) and v=(v1,v2,v3) then the standard vector cross product u×v is given by the mnemonic

u×v=det(ijku1u2u3v1v2v3)

so w=u×v has components w1=u2v3u3v2, and so on. Spivak (p83) gives a more rigorous definition and places it in a more general context. If v1,,vn1Rn and ϕΛ1(Rn) is defined by

ϕ(w)=det(v1vnw)=det(v11vn1v1n1vnn1wn1wn1)

then there is a unique zRn such that w,z=ϕ(w). The reason that w is at the bottom rather than the top is that it ensures that the the n-tuple (v1,,vn1,w) has positive orientation with respect to the standard basis vectors of Rn. In R3 we get the standard elementary mnemonic for u=(u1,u2,u3), v=(v1,v2,v3):

u×v=det(ijku1u2u3v1v2v3)

R implementation

The R function takes a matrix with n rows and n1 columns: the transpose of the work above. This is because stokes (and R) convention is to interpret columns of a matrix as vectors. If we wanted to take the cross product of u=(5,2,1) with v=(1,2,0):

##      [,1] [,2]
## [1,]    5    1
## [2,]   -2    2
## [3,]    1    0
## [1] -2  1 12

But of course we can work with higher dimensional spaces:

## [1] -7.892174  0.927769  1.070595 -5.334296  1.771056 -4.517277

Verification

We can demonstrate that the function has the correct orientation. We need to ensure that the vectors v1,vn,v1××vn constitute a right-handed basis:

det(cbind(M,vector_cross_product(M)))>0
## [1] TRUE

So it is right-handed in this case. Here is a more severe test:

f <- function(n){
  M <- matrix(rnorm(n^2+n),n+1,n)
  det(cbind(M,vector_cross_product(M)))>0
}

all(sapply(sample(3:10,100,replace=TRUE),f))
## [1] TRUE

Vector products and Hodge

The cross product has a coordinate-free definition as the Hodge conjugate of the wedge product of its arguments. This is not used in function vector_cross_product() because it is computationally inefficient and (I think) prone to numerical roundoff errors. We may verify that the definitions agree:

## [1] -21  -2   2  14
## An alternating linear map from V^1 to R with V=R^4:
##        val
##  2  =   -2
##  1  =  -21
##  4  =   14
##  3  =    2
## [1]   4.431826  -1.966102  -3.344998  -6.853352 -11.879641   7.170485
## An alternating linear map from V^1 to R with V=R^5:
##        val
##  4  =    0
##  2  =    0
##  1  =    0
##  5  =    0

Above, note that the output of vector_cross_product() is a 1-form (rather than an R vector), so has to be coerced to a 1-form.

Reference