# R/ex_andrews.R In deTestSet: Test Set for Differential Equations

#### Documented in andrews

```### ============================================================================
### andrews -- non-stiff second order differential algebraic equation of index
###          3, consisting of 14 differential and 13 algebraic equations
###          The problem describes the motion of 7 rigid bodies connected
###          by joints without friction
### ============================================================================

andrews <- function(times = seq(0,0.03, by=0.03/100) ,
yini = NULL, dyini = NULL, parms = list(),
printmescd = TRUE, method = mebdfi,
atol = 1e-7, rtol = 1e-7, maxsteps = 1e5, ...) {

### check input
parameter <- c(m1=.04325, m2=.00365, m3=.02373 ,m4=.00706 ,
m5=.07050 ,m6=.00706 ,m7=.05498 ,
xa=-.06934 ,ya=-.00227 ,xb=-0.03635 ,yb=.03273 ,
xc=.014 ,yc=.072 ,c0=4530 ,
i1=2.194e-6,i2=4.410e-7,i3=5.255e-6,i4=5.667e-7,
i5=1.169e-5,i6=5.667e-7,i7=1.912e-5,
d=28e-3,da=115e-4,e=2e-2,ea=1421e-5,
rr=7e-3,ra=92e-5,l0=7785e-5,
ss=35e-3,sa=1874e-5,sb=1043e-5,sc=18e-3,sd=2e-2,
ta=2308e-5,tb=916e-5,u=4e-2,ua=1228e-5,ub=449e-5,
zf=2e-2,zt=4e-2,fa=1421e-5,mom=33e-3)

parameter <- overrulepar(parameter, parms, 42)

prob <- andrewsprob()

if (is.null(yini))
yini <- c(-0.0617138900142764496358948458001, 0,
0.455279819163070380255912382449, 0.222668390165885884674473185609,
0.487364979543842550225598953530,-0.222668390165885884674473185609,
1.23054744454982119249735015568 ,0,     0,0,     0,0,          0,0,
14222.4439199541138705911625887,-10666.8329399655854029433719415,
0,0,   0,0,  0,98.5668703962410896057654982170,
-6.12268834425566265503114393122,0,         0,0,   0)

if (is.null(dyini)) {
dyini <- rep(0,27)
dyini[1:14] <- yini[8:21]
}
checkini(27, yini, dyini)

if (is.null(names(yini)))
names(yini) <- c(paste("q",1:7,sep=""),paste("dq",1:7,sep=""),
paste("dq2",1:7,sep=""),paste("lambda",1:6,sep=""))

if (is.null(times)) times <- prob\$times

### solve

ind <- c(7,7,13)  # index of the system

useres <- FALSE
if (is.character(method)) {
if (method %in% c("mebdfi", "daspk"))
useres <- TRUE
} else  if("res" %in% names(formals(method)))
useres <- TRUE

if (useres) {
AndOut <- dae(y = yini, dy = dyini, times = times,
res = "andres", nind = ind,
dllname = "deTestSet", jacres = "andjacres", initfunc = "andpar",
parms = parameter, jactype = "fullusr", method = method,
maxsteps = maxsteps, atol=atol, rtol=rtol,...)
}
else
{
if (prob\$numjac)
AndOut <- dae(y = yini, times = times, nind = ind,
func = "andfunc", mass =  prob\$mass,
massup = 0, massdown = 0,
dllname = "deTestSet", initfunc = "andpar", parms = parameter,
method = method, maxsteps = maxsteps, atol=atol, rtol=rtol, ...)
else{
fulljac = (prob\$mujac == prob\$neqn & prob\$mljac == prob\$neqn)
if (fulljac)
jactype <- "fullusr"
else
jactype <- "bandusr"
AndOut <- dae(y = yini, times = times, nind = ind,
func = "andfunc", mass =  prob\$mass,
massup = 0, massdown = 0,
dllname = "deTestSet", initfunc = "andpar", parms = parameter,
method = method, maxsteps = maxsteps, atol=atol, rtol=rtol,
jacfunc ="andjac", jactype = jactype, ...)
}
}
if (printmescd) AndOut <- printpr (AndOut, prob, "andrews", rtol, atol)
return(AndOut)
}

andrewsprob <- function(){
fullnm <- "Andrews' squeezing mechanism"
problm <- 'andrews'
type <- 'DAE'
neqn <- 27
ndisc <- 500
t <- matrix(1,2)
t[1] <- 0
t[2] <- 3e-2
numjac <- FALSE
mljac <- neqn
mujac <- neqn
mlmas <- 0
mumas <- 0
ind <- c(7,7,13)
mass =  as.double(c(rep(1, 14), rep(0, 13)))
return(list(fullnm=fullnm, problm=problm,type=type,neqn=neqn,ndisc=ndisc,
t=t,numjac=numjac,mljac=mljac,mujac=mujac,
mlmas=mlmas,mumas=mumas,nind=ind, mass=mass))
}

andrewsinit <- function( ){
yini <- c(-0.0617138900142764496358948458001, 0,
0.455279819163070380255912382449, 0.222668390165885884674473185609,
0.487364979543842550225598953530,-0.222668390165885884674473185609,
1.23054744454982119249735015568 ,0,     0,0,     0,0,          0,0,
14222.4439199541138705911625887,-10666.8329399655854029433719415,
0,0,   0,0,  0,98.5668703962410896057654982170,
-6.12268834425566265503114393122,0,         0,0,   0)
dyini <- rep(0,27)
dyini[1:14] <- yini[8:21]
return(list(yini,dyini))
}
```

## Try the deTestSet package in your browser

Any scripts or data that you put into this service are public.

deTestSet documentation built on May 29, 2017, 2:56 p.m.