[R-sig-dyn-mod] bvpshoot - different results?

[Radovan Omorjan]

Hello,

I just tried the help example about bvshoot()

Example 2 - initial condition is a function of the unknown x
## tubular reactor with axial dispersion

using both examples 1. by using yini() function and 2. by using bound() 
function.

I just used some slight changes, i.e. used -with(as.list() {})- 
construction. Here are four ways in two files. Only one differs from the 
other three (noted in the code2). I could not see what I was doing wrong 
- It must have been something trivial.

QUESTION: Could someone check these two codes, please, and tell me what 
I was doing wrong

# ###############  code1 ############################
require(bvpSolve)
##
Reactor<-function(x, y, parms)
{
   list(c(y[2],
          Pe * (y[2]+R*(y[1]^n))))
}

Pe <-  1
R <-  2
n <-  2

x <- seq(0, 1, by = 0.01)
# 1. yini is a function here
yini <- function (x, parms) c(x, Pe*(x-1))
system.time(
   sol <- bvpshoot(func = Reactor, yend = c(y = NA, dy = 0),
                    yini = yini, x = x, extra = 0.5)
)
plot(sol, which = "y", main = "Reactor", type = "l", lwd = 2)
attributes(sol)$roots
head(sol)


# 2. using boundary function rather than yini...
bound <- function(i, y, p) {
   if (i == 1) return(y[2] - Pe*(y[1]-1))
   if (i == 2) return(y[2])
}

# need to give number of left boundary conditions + guess of all initial
# conditions (+ names)
system.time(
   Sol2<- bvpshoot(func = Reactor, x = x, leftbc = 1,
                   bound = bound, guess = c(y = 1, dy = 1) )
)
plot(Sol2, which = "y", main = "Reactor", type = "l", lwd = 2)
attributes(Sol2)$roots
head(Sol2)
# ###############  code1 ############################


# ###############  code2 ############################
require(bvpSolve)
##
Reactor<-function(x, y, p){
   with(as.list(p),
          list(c(y[2],
                 Pe * (y[2]+R*(y[1]^n)))))
}
parms <- c(Pe = 1, R = 2, n = 2)

x <- seq(0, 1, by = 0.01)
# 1. yini is a function here
yini <- function (x, p) with(as.list(p),
                                    c(x, Pe*(x-1)))
system.time(
   sol <- bvpshoot(func = Reactor, yend = c(y = NA, dy = 0),
                   yini = yini, x = x, extra = 0.5, parms=parms)
)
plot(sol, which = "y", main = "Reactor", type = "l", lwd = 2)
attributes(sol)$roots
head(sol)###### THIS ONE IS DIFFERENT, WHY??? ########

# 2. using boundary function rather than yini...
bound <- function(i, y, p) {
   with(as.list(p),{
   if (i == 1) return(y[2] - Pe*(y[1]-1))
   if (i == 2) return(y[2])})
}

# need to give number of left boundary conditions + guess of all initial
# conditions (+ names)
system.time(
   Sol2<- bvpshoot(func = Reactor, x = x, leftbc = 1,
                   bound = bound, guess = c(y = 1, dy = 1), parms=parms )
)
plot(Sol2, which = "y", main = "Reactor", type = "l", lwd = 2)
attributes(Sol2)$roots
head(Sol2)
# ###############  code2 ############################