inst/doc/wedge.R

In wedge: The Exterior Calculus

## ----message=FALSE-------------------------------------------------------
library("wedge")

## ------------------------------------------------------------------------
k <- 4
n <- 5
M <- matrix(c(5,1,1,1, 1,1,2,3, 1,3,4,2),3,4,byrow=TRUE)
M
S <- as.ktensor(M,coeffs= 0.5 + 1:3)
S

## ------------------------------------------------------------------------
E <- matrix(rnorm(n*k),n,k)   # A random point in V^k

## ------------------------------------------------------------------------
f <- as.function(S)
f(E)

## ------------------------------------------------------------------------

E1 <- E
E2 <- E
E3 <- E

x1 <- rnorm(n)
x2 <- rnorm(n)
r1 <- rnorm(1)
r2 <- rnorm(1)

E1[,2] <- x1
E2[,2] <- x2
E3[,2] <- r1*x1 + r2*x2

## ------------------------------------------------------------------------
f <- as.function(S)
way1 <- r1*f(E1) + r2*f(E2)
way2 <- f(E3)
c(way1,way2)

## ------------------------------------------------------------------------
E1 <- matrix(rnorm(n*k),n,k)
E2 <- matrix(rnorm(n*k),n,k)
way1 <- f(r1*E1+r2*E2)
way2 <- r1*f(E1)+r2*f(E2)
c(way1,way2)

## ------------------------------------------------------------------------
S1 <- ktensor(spray(cbind(1:3,2:4),1:3))
S2 <- as.ktensor(matrix(1:6,2,3))
S1
S2

## ------------------------------------------------------------------------
cross(S1,S2)

## ------------------------------------------------------------------------
as.function(cross(S1,S2))(matrix(rnorm(30),6,5))

## ------------------------------------------------------------------------
S1
Alt(S1)

## ------------------------------------------------------------------------
E <- matrix(rnorm(8),4,2)
Erev <- E[,2:1]
as.function(Alt(S1))(E) + as.function(Alt(S1))(Erev)  # should be zero

## ------------------------------------------------------------------------
M <- matrix(c(4,2,3,1,2,4),2,3,byrow=TRUE)
M
K <- as.kform(M,c(1,5))
K

## ------------------------------------------------------------------------
M1 <- matrix(c(3,4,5, 4,6,1),2,3,byrow=TRUE)
K1 <- as.kform(M1,c(2,7))
K1
M2 <- cbind(1:5,3:7)
K2 <- as.kform(M2,1:5)
K2
K1 %^% K2


## ------------------------------------------------------------------------
F1 <- as.kform(matrix(c(3,4,5, 4,6,1,3,2,1),3,3,byrow=TRUE))
F2 <- as.kform(cbind(1:6,3:8),1:6)
F3 <- kform_general(1:8,2)
(F1 %^% F2) %^% F3
F1 %^% (F2 %^% F3)

## ------------------------------------------------------------------------
(F1 %^% F2) %^% F3 - F1 %^% (F2 %^% F3)

## ------------------------------------------------------------------------
Krel <- kform_general(4,2,1:6)
Krel

## ------------------------------------------------------------------------
K1 <- as.kform(matrix(1:4,2,2),1:2)
K2 <- as.kform(matrix(c(1,3,7,8,2,4),ncol=2,byrow=TRUE),c(-1,5,4))
K1
K2
K1+K2

## ------------------------------------------------------------------------
U <- ktensor(spray(cbind(1:4,2:5),1:4))
U

## ------------------------------------------------------------------------
as.symbolic(U)

## ------------------------------------------------------------------------
K <- kform_general(3,2,1:3)
K
as.symbolic(K,d="d",symbols=letters[23:26])

## ------------------------------------------------------------------------
(o <- rform())  # a random 3-form
V <- matrix(runif(21),ncol=3)
LHS <- as.function(o)(V)
RHS <- as.function(contract(o,V[,1]))(V[,-1])
c(LHS=LHS,RHS=RHS,diff=LHS-RHS)

## ------------------------------------------------------------------------
as.function(contract(o,V[,1:2]))(V[,-(1:2),drop=FALSE])

## ------------------------------------------------------------------------
contract(o,V)

## ------------------------------------------------------------------------
contract(o,V,lose=FALSE)

## ------------------------------------------------------------------------
grad(c(0.4,0.1,-3.2,1.5))

## ------------------------------------------------------------------------
f <- function(x){
    n <- length(x)
    as.kform(t(apply(diag(n)<1,2,which)))
}

## ------------------------------------------------------------------------
f(1:5)

## ------------------------------------------------------------------------
df  <- function(x){
    n <- length(x)
    S <- sum(x^2)
    grad(rep(c(1,-1),length=n)*(S^(n/2) - n*x^2*S^(n/2-1))/S^n
    )
}

## ------------------------------------------------------------------------
df(1:5)

## ------------------------------------------------------------------------
x <- rnorm(9)
print(df(x) %^% f(x))  # should be zero

## ------------------------------------------------------------------------
f1 <- function(w,x,y,z){x + y^3 + x*y*w*z}
f2 <- function(w,x,y,z){w^2*x*y*z + sin(w) + w+z}
f3 <- function(w,x,y,z){w*x*y*z + sin(x) + cos(w)}

## ------------------------------------------------------------------------
dw <- as.kform(1)
dx <- as.kform(2)
dy <- as.kform(3)
dz <- as.kform(4)

## ------------------------------------------------------------------------
phi <-
  (
    +f1(1,2,3,4) %^% dw %^% dx
    +f2(1,2,3,4) %^% dw %^% dy
    +f3(1,2,3,4) %^% dy %^% dz
  )

## ------------------------------------------------------------------------
e1 <- dw %^% dx
e2 <- dw %^% dy
e3 <- dy %^% dz

phi <-
  (
    +f1(1,2,3,4) %^% e1
    +f2(1,2,3,4) %^% e2
    +f3(1,2,3,4) %^% e3
  )
phi

## ------------------------------------------------------------------------
library("Deriv")
Df1 <- Deriv(f1)(1,2,3,4)
Df2 <- Deriv(f2)(1,2,3,4)
Df3 <- Deriv(f3)(1,2,3,4)

## ------------------------------------------------------------------------
Df1

## ------------------------------------------------------------------------
Hf1 <- matrix(Deriv(f1,nderiv=2)(1,2,3,4),4,4)
Hf2 <- matrix(Deriv(f2,nderiv=2)(1,2,3,4),4,4)
Hf3 <- matrix(Deriv(f3,nderiv=2)(1,2,3,4),4,4)

## ---- echo=FALSE---------------------------------------------------------
rownames(Hf1) <- c("w","x","y","z")
colnames(Hf1) <- c("w","x","y","z")

## ------------------------------------------------------------------------
Hf1

## ------------------------------------------------------------------------
ddphi <- # should be zero
  (  
    +as.kform(which(!is.na(Hf1),arr.ind=TRUE),c(Hf1))
    +as.kform(which(!is.na(Hf2),arr.ind=TRUE),c(Hf2))
    +as.kform(which(!is.na(Hf3),arr.ind=TRUE),c(Hf3))
  )

ddphi

## ------------------------------------------------------------------------
phi <- function(x){
    n <- length(x)
    sum(x^seq_len(n)*rep_len(c(1,-1),n)) * as.kform(t(apply(diag(n)<1,2,which)))
}
phi(1:9)

## ------------------------------------------------------------------------
dphi <- function(x){
    nn <- seq_along(x)
    sum(nn*x^(nn-1)) * as.kform(seq_along(x))
}
dphi(1:9)