#################################################################
# Example 1: simple standard problem
# solve the BVP ODE:
# d2y/dt^2=-3py/(p+t^2)^2
# y(t= -0.1)=-0.1/sqrt(p+0.01)
# y(t= 0.1)= 0.1/sqrt(p+0.01)
# where p = 1e-5
#
# analytical solution y(t) = t/sqrt(p + t^2).
#
# The problem is rewritten as a system of 2 ODEs:
# dy=y2
# dy2=-3p*y/(p+t^2)^2
################################################################################
#--------------------------------
# Derivative function
#--------------------------------
fun <- function(t,y,pars)
{ dy1 <-y[2]
dy2 <- - 3*p*y[1]/(p+t*t)^2
return(list(c(dy1,
dy2))) }
# parameter value
p <-1e-5
# initial and final condition; second conditions unknown
init <- c(-0.1/sqrt(p+0.01), NA)
end <- c(0.1/sqrt(p+0.01), NA)
# Solve bvp
sol <- as.data.frame(bvptwp(yini=init,x=seq(-0.1,0.1,by=0.001),
func=fun, yend=end, guess=1))
plot(sol$time,sol[,2],type="l")
# add analytical solution
curve(x/sqrt(p+x*x),add=TRUE,type="p")
#################################################################
# Example 1b: simple
# solve d2y/dx2 + 1/x*dy/dx + (1-1/(4x^2)y = sqrt(x)*cos(x),
# on the interval [1,6] and with boundary conditions:
# y(1)=1, y(6)=-0.5
#
# Write as set of 2 odes
# dy/dx = y2
# dy2/dx = - 1/x*dy/dx - (1-1/(4x^2)y + sqrt(x)*cos(x)
#################################################################
f2 <- function(x,y,parms)
{
dy <- y[2]
dy2 <- -1/x*y[2]-(1-1/(4*x^2))*y[1] + sqrt(x)*cos(x)
list(c(dy,dy2))
}
x <- seq(1,6,0.1)
sol <- bvptwp(yini=c(1,NA),yend=c(-0.5,NA),x=x,func=f2,guess=1)
plot(sol)
# add the analytic solution
curve(0.0588713*cos(x)/sqrt(x)+1/4*sqrt(x)*cos(x)+0.740071*sin(x)/sqrt(x)+
1/4*x^(3/2)*sin(x),add=TRUE,type="l")
#################################################################
# Problem 2 - solved with specification of boundary, and jacobians
# d4y/dx4 =R(dy/dx*d2y/dx2 -y*dy3/dx3)
# y(0)=y'(0)=0
# y(1)=1, y'(1)=0
#
# dy/dx = y2
# dy2/dx = y3 (=d2y/dx2)
# dy3/dx = y4 (=d3y/dx3)
# dy4/dx = R*(y2*y3 -y*y4)
#################################################################
f2<- function(x,y,parms,R)
{
list(c(y[2],y[3],y[4],R*(y[2]*y[3]-y[1]*y[4]) ))
}
df2 <- function(x,y,parms,R)
{
df <- matrix(nr=4,nc=4,byrow=TRUE,data=c(
0,1,0,0,
0,0,1,0,
0,0,0,1,
-1*R*y[4],R*y[3],R*y[2],-R*y[1]))
}
g2 <- function(i,y,parms,R)
{
if ( i ==1) return(y[1])
if ( i ==2) return(y[2])
if ( i ==3) return(y[1]-1)
if ( i ==4) return(y[2])
}
dg2 <- function(i,y,parms,R)
{
if ( i ==1) return(c(1,0,0,0))
if ( i ==2) return(c(0,1,0,0))
if ( i ==3) return(c(1,0,0,0))
if ( i ==4) return(c(0,1,0,0))
}
require(bvpSolve)
init <- c(1,NA)
R <- 0.01
sol <- as.data.frame(bvptwp(x=seq(0,1,by=0.01), leftbc=2,
func=f2, guess=1,R=R,
bound=g2, jacfunc=df2, jacbound=dg2))
plot(sol[,1:2])Run the code above in your browser using DataLab