########################################
# EXAMPLE 1
# 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 <- lsode(yini,times,f1,parms=0,jactype="fullint")
# stiff method, user-generated full jacobian
out2 <- lsode(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 <- lsode(yini,times,f1,parms=0,jactype="bandint",
bandup=1,banddown=1)
# stiff method, user-generated banded jacobian
out4 <- lsode(yini,times,f1,parms=0,jactype="bandusr",
jacfunc=bandjac,bandup=1,banddown=1)
# nono-stiff method
out5 <- lsode(yini,times,f1,parms=0,mf=10)
########################################
# example 2: diffusion on a 2-D grid
# partially specified jacobian
########################################
diffusion2D <- function(t,Y,par)
{
y <- matrix(nr=n,nc=n,data=Y)
dY <- r*y #production
#diffusion in X-direction; boundaries=0-concentration
Flux <- -Dx * rbind(y[1,],(y[2:n,]-y[1:(n-1),]),-y[n,])/dx
dY <- dY - (Flux[2:(n+1),]-Flux[1:n,])/dx
#diffusion in Y-direction
Flux <- -Dy * cbind(y[,1],(y[,2:n]-y[,1:(n-1)]),-y[,n])/dy
dY <- dY - (Flux[,2:(n+1)]-Flux[,1:n])/dy
return(list(as.vector(dY)))
}
# parameters
dy <- dx <- 1 # grid size
Dy <- Dx <- 1 # diffusion coeff, X- and Y-direction
r <- 0.025 # production rate
times <- c(0,1)
n <- 50
y <- matrix(nr=n,nc=n,0.)
pa <-par(ask=FALSE)
# initial condition
for (i in 1:n)
for (j in 1:n)
{dst <- (i-n/2)^2+(j-n/2)^2
y[i,j]<-max(0.,1.-1./(n*n)*(dst-n)^2)
}
filled.contour(y,color.palette=terrain.colors)
# jacfunc need not be estimated exactly;
# a crude approximation, with a smaller bandwidth will do.
# Here the half-bandwidth 1 is used, whereas the true
# half-bandwidths are equal to n.
# This corresponds to ignoring the y-direction coupling in the ODEs.
print(system.time(
for (i in 1:20)
{
out <- lsode(func=diffusion2D,y=as.vector(y),times=times,
parms=NULL,jactype="bandint", bandup=1,banddown=1)
filled.contour(matrix(nr=n,nc=n,out[2,-1]),zlim=c(0,1),
color.palette=terrain.colors,main=i)
y <- out[2,-1]
}
))
par(ask=pa)Run the code above in your browser using DataLab