## =============================================================================
## 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(y=-0.1/sqrt(p+0.01), dy=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$x,sol$y,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(y=1,dy=NA),yend=c(-0.5,NA),x=x,func=f2,guess=1)
plot(sol,which="y")
# 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")
## =============================================================================
## Example 2 - using continuation
## Problem nr 24 from http://www.ma.ic.ac.uk/~jcash/BVP_software/PROBLEMS.PDF
## xi*A*y*y'' -((1+ga)/2-xi*A'*y*y'+y'/y+(A'/A*(1-(ga-1)/2*y^2)=0
## where A = 1 + x2, y(x=0) = 0:9129, y(x=1) = 0:375
## =============================================================================
Prob24 <- function(x, y, pars, xi) {
A <- 1+x*x
AA <- 2*x
ga <- 1.4
list(c(y[2],
(((1+ga)/2 -xi*AA)*y[1]*y[2]-y[2]/y[1]-
(AA/A)*(1-(ga-1)*y[1]^2/2))/(xi*A*y[1]) ))
}
ini <- c(0.9129,NA)
end <- c(0.375,NA)
# This problem needs an initial guess for y; (as it divides by y,
# the default, 0, does not work)
# guess of y; set equal to 1
xguess <- c(0,1)
yguess <- matrix(nrow=2,ncol=2,data=1)
# start with relatively large value of xi
xi <-0.01
print(system.time(
mod1 <- bvptwp(yini=ini,yend=end,x=seq(0,1,by=0.01),
func=Prob24,xi=xi, xguess=xguess,yguess=yguess)
))
plot(mod1[,1:2], ylim = c(0.4,1.4), type="l", lwd=2)
# use previous solution to initialise solver for smaller xi
xi <-0.005
mod1b <- bvptwp(yini=ini,yend=end,x=seq(0,1,by=0.01),
func=Prob24,xi=xi, xguess=mod1[,1],yguess=t(mod1[,2:3]))
lines(mod1b, lwd=2, col="red")
## =============================================================================
## Example 3 - 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 <- 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]) # columns do not have names
mf <- par ("mfrow")
sol <- bvptwp(x=seq(0,1,by=0.01), leftbc=2,
func=f2, guess=c(y=1,dy=1,dy2=1,dy3=1),R=R,
bound=g2, jacfunc=df2, jacbound=dg2)
plot(sol) # here they do
par(mfrow=mf)Run the code above in your browser using DataLab