DAE

## =============================================================================
## R-Code to solve example 5.5.1 from the book
## K. Soetaert, J. Cash and F. Mazzia, 2012.  
## Solving differential equations in R. UseR, Springer, 248 pp. 
## http://www.springer.com/statistics/computational+statistics/book/978-3-642-28069-6.
## implemented by Karline Soetaert
## =============================================================================

## A differential algebraic equation model
library(deSolve)

Caraxis <- function(t, y, dy, parms) {
 with(as.list(y), {
   f <- rep(0, 10)
   yb <- r * sin(w * t)
   xb <- sqrt(L^2  - yb^2)
   Ll <- sqrt(xl^2 + yl^2)
   Lr <- sqrt((xr - xb)^2 + (yr - yb)^2)
   f[1:4] <- y[5:8]
   f[5] <- 1/k*((L0-Ll)*xl/Ll + lam1*xb + 2*lam2*(xl-xr))
   f[6] <- 1/k*((L0-Ll)*yl/Ll + lam1*yb + 2*lam2*(yl-yr)) -g
   f[7] <- 1/k*((L0-Lr)*(xr - xb)/Lr    - 2*lam2*(xl-xr))
   f[8] <- 1/k*((L0-Lr)*(yr - yb)/Lr    - 2*lam2*(yl-yr)) -g
   f[9] <- xb * xl + yb * yl
   f[10]<- (xl - xr)^2 + (yl - yr)^2 - L^2

   delt       <- dy - f
   delt[9:10] <- -f[9:10]

   list(delt)
 })
}

eps <- 0.01; M <- 10; k <- M * eps * eps/2 
L <- 1; L0 <- 0.5; r <- 0.1; w <- 10; g <- 9.8

yini <- c(xl = 0,     yl = L0, xr = L,     yr = L0, 
         ul = -L0/L, vl = 0,  ur = -L0/L, vr = 0, 
         lam1 = 0, lam2 = 0)

library(deTestSet)
rootfun <- function (dyi, y, t)
  unlist(Caraxis(t, y, dy = c(dyi, 0, 0), 
        parms = NULL)) [1:8]

dyini <- multiroot(f = rootfun, start = rep(0,8), 
                  y = yini, t = 0)$root  
(dyini <- c(dyini, 0, 0))
Caraxis(t = 0, yini, dyini, NULL)

index <- c(4, 4, 2)
times <- seq(from = 0, to = 3, by = 0.01)
out <- mebdfi(y = yini, dy = dyini, times = times, 
             res = Caraxis, parms = parameter, nind = index)

par(mar = c(4, 4, 3, 2))
plot(out, lwd = 2, mfrow = c(4,3))
plot(out[,c("xl", "yl")], xlab = "xleft", ylab = "yleft",
    type = "l", lwd = 2)
plot(out[,c("xr", "yr")], xlab = "xright", ylab = "yright",
    type = "l", lwd = 2)