## =======================================================================
## Example 1: ODE
## Various ways to solve the same model.
## =======================================================================
## the model, 5 state variables
f1 <- function (t, y, parms) {
ydot <- vector(len = 5)
ydot[1] <- 0.1*y[1] -0.2*y[2]
ydot[2] <- -0.3*y[1] +0.1*y[2] -0.2*y[3]
ydot[3] <- -0.3*y[2] +0.1*y[3] -0.2*y[4]
ydot[4] <- -0.3*y[3] +0.1*y[4] -0.2*y[5]
ydot[5] <- -0.3*y[4] +0.1*y[5]
return(list(ydot))
}
## the Jacobian, written as a full matrix
fulljac <- function (t, y, parms) {
jac <- matrix(nrow = 5, ncol = 5, byrow = TRUE,
data = c(0.1, -0.2, 0 , 0 , 0 ,
-0.3, 0.1, -0.2, 0 , 0 ,
0 , -0.3, 0.1, -0.2, 0 ,
0 , 0 , -0.3, 0.1, -0.2,
0 , 0 , 0 , -0.3, 0.1))
return(jac)
}
## the Jacobian, written in banded form
bandjac <- function (t, y, parms) {
jac <- matrix(nrow = 3, ncol = 5, byrow = TRUE,
data = c( 0 , -0.2, -0.2, -0.2, -0.2,
0.1, 0.1, 0.1, 0.1, 0.1,
-0.3, -0.3, -0.3, -0.3, 0))
return(jac)
}
## initial conditions and output times
yini <- 1:5
times <- 1:20
## default: stiff method, internally generated, full Jacobian
out <- radau(yini, times, f1, parms = 0)
plot(out)
## stiff method, user-generated full Jacobian
out2 <- radau(yini, times, f1, parms = 0, jactype = "fullusr",
jacfunc = fulljac)
## stiff method, internally-generated banded Jacobian
## one nonzero band above (up) and below(down) the diagonal
out3 <- radau(yini, times, f1, parms = 0, jactype = "bandint",
bandup = 1, banddown = 1)
## stiff method, user-generated banded Jacobian
out4 <- radau(yini, times, f1, parms = 0, jactype = "bandusr",
jacfunc = bandjac, bandup = 1, banddown = 1)
## =======================================================================
## Example 2: ODE
## stiff problem from chemical kinetics
## =======================================================================
Chemistry <- function (t, y, p) {
dy1 <- -.04*y[1] + 1.e4*y[2]*y[3]
dy2 <- .04*y[1] - 1.e4*y[2]*y[3] - 3.e7*y[2]^2
dy3 <- 3.e7*y[2]^2
list(c(dy1, dy2, dy3))
}
times <- 10^(seq(0, 10, by = 0.1))
yini <- c(y1 = 1.0, y2 = 0, y3 = 0)
out <- radau(func = Chemistry, times = times, y = yini, parms = NULL)
plot(out, log = "x", type = "l", lwd = 2)
## =============================================================================
## Example 3: DAE
## Car axis problem, index 3 DAE, 8 differential, 2 algebraic equations
## from
## F. Mazzia and C. Magherini. Test Set for Initial Value Problem Solvers,
## release 2.4. Department
## of Mathematics, University of Bari and INdAM, Research Unit of Bari,
## February 2008.
## Available at http://www.dm.uniba.it/~testset.
## =============================================================================
## Problem is written as M*y' = f(t,y,p).
## caraxisfun implements the right-hand side:
caraxisfun <- function(t, y, parms) {
with(as.list(y), {
yb <- r * sin(w * t)
xb <- sqrt(L * L - yb * yb)
Ll <- sqrt(xl^2 + yl^2)
Lr <- sqrt((xr - xb)^2 + (yr - yb)^2)
dxl <- ul; dyl <- vl; dxr <- ur; dyr <- vr
dul <- (L0-Ll) * xl/Ll + 2 * lam2 * (xl-xr) + lam1*xb
dvl <- (L0-Ll) * yl/Ll + 2 * lam2 * (yl-yr) + lam1*yb - k * g
dur <- (L0-Lr) * (xr-xb)/Lr - 2 * lam2 * (xl-xr)
dvr <- (L0-Lr) * (yr-yb)/Lr - 2 * lam2 * (yl-yr) - k * g
c1 <- xb * xl + yb * yl
c2 <- (xl - xr)^2 + (yl - yr)^2 - L * L
list(c(dxl, dyl, dxr, dyr, dul, dvl, dur, dvr, c1, c2))
})
}
eps <- 0.01; M <- 10; k <- M * eps^2/2;
L <- 1; L0 <- 0.5; r <- 0.1; w <- 10; g <- 1
yini <- c(xl = 0, yl = L0, xr = L, yr = L0,
ul = -L0/L, vl = 0,
ur = -L0/L, vr = 0,
lam1 = 0, lam2 = 0)
# the mass matrix
Mass <- diag(nrow = 10, 1)
Mass[5,5] <- Mass[6,6] <- Mass[7,7] <- Mass[8,8] <- M * eps * eps/2
Mass[9,9] <- Mass[10,10] <- 0
Mass
# index of the variables: 4 of index 1, 4 of index 2, 2 of index 3
index <- c(4, 4, 2)
times <- seq(0, 3, by = 0.01)
out <- radau(y = yini, mass = Mass, times = times, func = caraxisfun,
parms = NULL, nind = index)
plot(out, which = 1:4, type = "l", lwd = 2)Run the code above in your browser using DataLab