DAE mass

## =============================================================================
## R-Code to solve example 5.6.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.
##
## A DAE model, solved in linearly implicit form
## The transistor amplifier
## implemented by Karline Soetaert
## =============================================================================

library(deSolve)

Transistor <- function(t, u, du, pars) {
  delt <- vector(length = 8)
   uin  <- 0.1 * sin(200 * pi * t)
   g23  <- beta * (exp( (u[2] - u[3]) / uf) - 1)
   g56  <- beta * (exp( (u[5] - u[6]) / uf) - 1)

   delt[1] <- (u[1] - uin)/R0
   delt[2] <- u[2]/R1 + (u[2]-ub)/R2 + (1-alpha) * g23
   delt[3] <- u[3]/R3 - g23
   delt[4] <- (u[4] - ub) / R4 + alpha * g23
   delt[5] <- u[5]/R5 + (u[5]-ub)/R6 + (1-alpha) * g56
   delt[6] <- u[6]/R7 - g56
   delt[7] <- (u[7] - ub) / R8 + alpha * g56
   delt[8] <- u[8]/R9
   list(delt)
}

ub <- 6; uf <- 0.026; alpha <- 0.99; beta <- 1e-6; R0 <- 1000
R1 <- R2 <- R3 <- R4 <- R5 <- R6 <- R7 <- R8 <- R9 <- 9000
C1 <- 1e-6; C2 <- 2e-6; C3 <- 3e-6; C4 <- 4e-6; C5 <- 5e-6

mass <- matrix(nrow = 8, ncol = 8, byrow = TRUE, data = c(
      -C1,C1, 0,  0,  0,  0,  0,  0,
      C1,-C1, 0,  0,  0,  0,  0,  0,
      0,  0,-C2,  0,  0,  0,  0,  0,
      0,  0,  0,-C3, C3,  0,  0,  0,
      0,  0,  0, C3,-C3,  0,  0,  0,
      0,  0,  0,  0,  0,-C4,  0,  0,
      0,  0,  0,  0,  0,  0,-C5, C5,
      0,  0,  0,  0,  0,  0, C5,-C5
))

yini <- c(0, ub/(R2/R1+1), ub/(R2/R1+1),
         ub, ub/(R6/R5+1), ub/(R6/R5+1), ub, 0)
names(yini) <- paste("u", 1:8, sep = "")

ind   <- c(8, 0, 0)
times <- seq(from = 0, to = 0.2, by = 0.001)

out <- radau(func = Transistor, y = yini, parms = NULL, 
            times = times, mass = mass, nind = ind)

plot(out, lwd = 2, which = c("u1", "u4", "u5", "u8"))