https://github.com/cran/pracma
Tip revision: c79a04b5074656b36e591191eb8137b70a349932 authored by Hans W. Borchers on 30 June 2014, 00:00:00 UTC
version 1.7.0
version 1.7.0
Tip revision: c79a04b
bvp.R
##
## b v p . R Boundary Value Problems
##
bvp <- function(f, g, h, x, y, n = 50) {
stopifnot(is.numeric(x), is.numeric(y), is.numeric(n))
if (length(x) != 2 || length(y) != 2)
stop("Arguments 'x' and 'y' must have length 2.")
if (length(n) != 1 || floor(n) != ceiling(n) || n < 2)
stop("Argument 'n' must be an integer greater or equal 2.")
if (is.numeric(f)) ffun <- function(x) rep(f[1], length(x))
else ffun <- match.fun(f)
if (is.numeric(g)) gfun <- function(x) rep(g[1], length(x))
else gfun <- match.fun(g)
if (is.numeric(h)) hfun <- function(x) rep(h[1], length(x))
else hfun <- match.fun(h)
xa <- x[1]; xb <- x[2]
ya <- y[1]; yb <- y[2]
xs <- linspace(xa, xb, n+2)[2:(n+1)]
dt <- (xb - xa) / (n+1)
a <- -2 - dt^2 * gfun(xs) # main diagonal
b <- 1 - dt/2 * ffun(xs[1:(n-1)]) # superdiagonal
d <- 1 + dt/2 * ffun(xs[2:n]) # subdiagonal
rhs <- dt^2 * hfun(xs) # right hand side
rhs[1] <- rhs[1] - ya * (1 + (dt/2) * ffun(xs[1]))
rhs[n] <- rhs[n] - yb * (1 - (dt/2) * ffun(xs[n]))
ys <- trisolve(a, b, d, rhs)
return(list(xs = c(xa, xs, xb), ys = c(ya, ys, yb)))
}
# bvp <- function(f, a, b, ya, yb, N, cc, ...) {
# stopifnot(is.numeric(a), length(a) == 1, is.numeric(b), length(b) == 1,
# is.numeric(ya), length(ya) == 1, is.numeric(yb), length(yb) == 1)
#
# if (!is.numeric(N) || length(N) != 1 || floor(N) != ceiling(N) || N < 0)
# stop("Argument 'N' must be an integer greater or equal 1.")
# if (!is.numeric(cc) || length(cc) != 3)
# stop("Argument 'cc' must be a real vector of length 3.")
#
# fun <- match.fun(f)
# f <- function(x) fun(x, ...)
#
# h <- (b-a)/(N+1)
# xh <- linspace(a, b, N+2)
# hm <- cc[1]/h^2
# hd <- cc[2]/(2*h)
#
# A <- diag((2*hm + cc[3]), N, N)
# A[col(A) == row(A)-1] <- -hm - hd
# A[col(A) == row(A)+1] <- -hm - hd
#
# F <- f(xh[2:(N+1)])
# F[1] <- F[1] + ya*(hm + hd)
# F[N] <- F[N] + yb*(hm - hd)
#
# yh <- qr.solve(A, F)
# yh <- c(ya, yh, yb)
#
# return(list(xh = xh, yh = yh))
# }