# R/bvpshoot.R In bvpSolve: Solvers for Boundary Value Problems of Differential Equations

#### Documented in bvpshoot

```## EXTENDED FOR DAEs
##==============================================================================
## Solving boundary value problems of ordinary differential equations
## using the shooting method
##==============================================================================

bvpshoot<- function(yini = NULL, x, func, yend = NULL, parms = NULL,
order = NULL,  guess = NULL,
jacfunc = NULL, bound = NULL, jacbound = NULL,
leftbc = NULL, posbound = NULL, ncomp = NULL,
atol = 1e-8, rtol = 1e-8, extra = NULL, maxiter = 100,
positive = FALSE, method = "lsoda", ...)  {

## ---------------------
## check input
## ---------------------

if (is.null(yini)   && is.null(bound))
stop("either 'yini' and 'yend' or 'bound' should be inputted")
if (!is.null(yini)  && is.null(yend))
stop("if 'yini' is given, 'yend' should also be given")
if (!is.null(bound) && is.null(posbound)&& is.null(leftbc))
stop("if 'bound' is given, 'posbound' or 'leftbc' should also be given")
if (!is.null(yini)  && !is.null(bound))
stop("either 'yini' or bound should be given, not both")
if (!is.null(bound)  && !is.null(extra))
stop("cannot combine 'bound' with 'extra'; use 'yini' and 'yend' ")
if (is.character(func) | is.character(jacfunc) |
is.character(bound) | is.character(jacbound))
stop("cannot use 'bvpshoot' with compiled code")

lex      <- length(extra)

# dots for passing to function and dots for method
fdots <- list(...)
metargs <- names(formals(method))
fmet <- which(names(fdots) %in% metargs)
fdots [fmet] <- NULL

if (is.function(method))
Ode <- method  else
Ode <- function(...) ode (method = method, ...)

if (length(fdots) == 0)
Do.call <- function(x, ll) do.call(x, ll)
else
Do.call <- function(x, ll) do.call(x, c(ll,unlist(fdots)))

## ---------------------------
## The order of the equations
## ---------------------------
JacFunc <- jacfunc # save jacfunc for further passing it to ode()
Func <- func
testit <- FALSE
attrib <- rep(0,4) # number of steps, number of fn evaluations, # jacobians ,nr ivp

if (! is.null(order)) {
mstar <- sum(order)
neq   <- length(order)
if (! is.null(ncomp))
if (sum(order) != ncomp)
stop("'ncomp' and 'order' not compatible: ncomp should equal sum(order)")
if (max(order) > 1) { # wrapper over func and jacfunc
testit <- TRUE
stareq <- cumsum(order)            # from func to vector returned to c
higord <- (1:mstar)[-stareq]  # from state to vector returned to c

# expand func
Fret <- numeric(length = mstar)
Func    <- function(x, state, parms,...)  {
FF <- Do.call(func, list(x, state, parms))
Fret[stareq] <- unlist(FF[1])
Fret[higord] <- state[higord+1]
FF[1] <- list(Fret)
FF
}

if (! is.null(jacfunc))
stop ("can not combine analytical jacobian with higher-order equations - remove 'jacfunc'")
}
if (is.null(ncomp)) ncomp <- mstar
}

## ---------------------
## yini or bound
## ---------------------
if (! is.null(yini))  {
# initial value yini; a function or vector
inity <- function(X, Parms) {  # initialises yini and parms..
if (is.function(yini))
Y <- do.call(yini,list(X,Parms))
else Y <- yini
if (lini>0)
Y[inix] <- X[1:lini]
names(Y)<-Ynames
Y
}
# parameters; some may be estimated (lex)
initparms <- function (X) {
Initparms <- parms
if(lex>0)
Initparms[1:lex] <- X[(lini+1):(lini+lex)]
return(Initparms)
}

if (is.function(yini))
y <- Do.call(yini,list(extra, parms))
else
y <- yini

# root function to solve for  - note jacfunc = NULL to avoid error in devel R 2.12
rootfun <- function(X, jacfunc = NULL, ...)  {
times <- c(x[1], x[length(x)])
Parms <- initparms(X)
Y     <- inity(X,Parms)
out   <- Ode(y=Y, times=times, func=Func, jacfunc=JacFunc,
parms=Parms, #method=method,
atol=atol, rtol=rtol, ...)
attrib <<- attrib + c(attributes(out)\$istate[c(2,3,14)],1)
# deviation with yend should be =0
if (is.function(yend) )
Res   <- yend(out[nrow(out),2:(ly+1)], Y, Parms,...)
else {
Res <- yend - out[nrow(out),2:(ly+1)]
Res <- Res[! is.na(Res)]
}
return(Res)
}
# the jacobian of root function: not specified
JacBound <- NULL

} else {                  # bound is specified
posspecified <- FALSE
if (is.null(leftbc)& is.null(posbound))
stop("'leftbc' or 'posbound' should be inputted if 'bound' is given")
if (! is.null(posbound)) {
if (all(posbound %in% c(x[1],x[length(x)])))
leftbc <- sum(posbound == x[1])
else  {
posspecified <- TRUE
# check consistency
if (! is.null(ncomp)) {
if(length(posbound) != ncomp)
stop("'posbound' should have a length = number of variables ")
} else ncomp <- length(posbound)

iipos <- which (x %in% posbound)
if (length(iipos) != ncomp)
stop("all elements in 'posbound' should also be in 'x' ")

}
}
if (! is.null(guess) & ! is.null(ncomp))
if (length(guess) != ncomp) stop ("length of 'guess' should be = number of variables")

y <- guess

rootfun <- function(X, jacfunc = NULL, ...)  {
if (! posspecified)
times <- c(x[1], x[length(x)])
else
times <- x
out   <- Ode(y = X, times = times, func = Func, jacfunc = JacFunc,
parms = parms, #method = method,
atol = atol, rtol = rtol, ...)
attrib <<- attrib + c(attributes(out)\$istate[c(2,3,14)],1)
Res <- vector(length=ly)
if (! posspecified) {
Yend <- out[nrow(out),2:(ly+1)]
for (i in 1:leftbc)
Res[i] <- Do.call(bound,list(i,X,parms))
if (leftbc < ly)
for (i in  (leftbc+1):ly)
Res[i]<- Do.call(bound,list(i,Yend,parms))
} else
for (i in 1:length(posbound)) {
ii <- iipos[i]
Res[i] <- Do.call(bound,list(i,out[ii,-1],parms))
}
return(Res)
}
JacBound <- NULL
if (! is.null(jacbound)) {
JAC <- matrix(nrow=length(y),ncol=length(y),0)
JacBound <- function(x,...)  {
for (i in 1:ly)
JAC[i,]<- Do.call(jacbound,list(i,x,parms))
return(JAC)
}
}
}
## ---------------------
## names, guess of y
## ---------------------

Ynames <- attr(y,"names")
if (is.null(Ynames) & ! is.null(yend))
Ynames <- names(yend)
if (is.null(y)) {
if (is.null(ncomp ))
stop("don't know number of variables - provide 'ncomp'")
y <- rep(0,ncomp)
warning("estimates for initial conditions not given ('guess'); assuming 0's")
}
if (is.null(y))
stop ("either provide 'guess' for initial conditions or 'ncomp', number of compartments")

ly <-  length(y)

inix     <- which (is.na(y))
lini     <- length(inix)
if (lini > 0 & is.null(guess))  {
warning("estimates for unknown initial conditions not given ('guess'); assuming 0's")
guess <- rep(0,lini)
}

if (! is.null(yini) & lini != length(guess))  {
if (is.null(extra))
stop("length of guess should be equal to number of NAs in y") else
if (lex > length(parms))
stop("length of extra should be smaller than number of parameters")
}
if (lini > 0)
y[inix] <- guess

if (is.null(yini) & is.null(guess))
guess <- y
if (! is.null(yini) & lini+lex==0)
stop ("this is not a boundary value problem - use initial value problem solver instead")

if (testit) {
if (! is.null(yini))
Parms <- initparms(extra)
else
Parms <- parms
FF <- Do.call(func,list(x[1], y, Parms))
if (length(FF[[1]]) != neq)
stop("function 'func' should return as many elements as the length of 'order'")
}

## ---------------------
## root solver:
## ---------------------
# find unknown initial condition + extra parameter
sol <- multiroot(start = c(guess,extra), f = rootfun, atol=atol, rtol=rtol,
maxiter=maxiter, jacfunc = JacBound, positive =positive, ...)

## ---------------------
## Output
## ---------------------

if (! is.null(yini) ) {
Parms <- initparms(sol\$root)
Y     <- inity(sol\$root, Parms)
}  else {
Parms <- parms
Y     <-  sol\$root
}

out <- Ode (times=x, func=Func, y=Y, parms=Parms, jacfunc=jacfunc, #method=method,
atol=atol, rtol=rtol, ...)
attrib <- attrib + c(attributes(out)\$istate[c(2,3,14)],1)

attrib[2] <- attrib[2] +attrib[4] +1 # correct for one more per ivp solution

attr(out,"roots")  <- data.frame(root=sol\$root,
f.root=sol\$f.root, iter=sol\$iter)
class(out) <- c("bvpSolve","matrix","deSolve")  # a boundary value problem
colnames(out)[1] <- "x"
attr(out,"name") <- "bvpshoot"
attr(out,"istate2") <- attrib[c(2,3,1,4)]
names (attr(out,"istate2")) <- c("nfunc", "njac", "nstep", "nivp")

out
}
```

## Try the bvpSolve package in your browser

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

bvpSolve documentation built on June 30, 2020, 3:02 a.m.