vector_cross_product()
in the stokes
package
## 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>
Given two vectors u,v∈R3 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=u2v3−u3v2, and so on. Spivak (p83) gives a more rigorous definition and places it in a more general context. If v1,…,vn−1∈Rn and ϕ∈Λ1(Rn) is defined by
ϕ(w)=det(v1⋮vnw)=det(v11…vn1⋮⋱⋮v1n−1…vnn−1wn−1…wn−1)
then there is a unique z∈Rn 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,…,vn−1,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)
The R function takes a matrix with n rows and n−1 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
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:
## [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
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
H <- hodge(as.1form(M[,1]) ^ as.1form(M[,2]) ^ as.1form(M[,3]) ^ as.1form(M[,4]) ^ as.1form(M[,5]))
H - as.1form(vector_cross_product(M)) # zero to numerical precision.
## 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.