Boundary Value Problems

Solves boundary value problems of linear second order differential equations.

bvp(f, g, h, x, y, n = 50)

Arguments

f, g, h

functions on the right side of the differential equation. If f, g or h is a scalar instead of a function, it is assumed to be a constant coefficient in the differential equation.

x

x[1], x[2] are the interval borders where the solution shall be computed.

y

boundary conditions such that y(x[1]) = y[1], y(x[2]) = y[2].

n

number of intermediate grid points; default 50.

Details

Solves the two-point boundary value problem given as a linear differential equation of second order in the form: $$y'' = f(x) y' + g(x) y + h(x)$$ with the finite element method. The solution \(y(x)\) shall exist on the interval \([a, b]\) with boundary conditions \(y(a) = y_a\) and \(y(b) = y_b\).

Value

Returns a list list(xs, ys) with the grid points xs and the values ys of the solution at these points, including the boundary points.

References

Kutz, J. N. (2005). Practical Scientific Computing. Lecture Notes 98195-2420, University of Washington, Seattle.

Note

Uses a tridiagonal equation solver that may be faster then qr.solve for large values of n.

See also

shooting

Examples

##  Solve y'' = 2*x/(1+x^2)*y' - 2/(1+x^2) * y + 1
##  with y(0) = 1.25 and y(4) = -0.95 on the interval [0, 4]:
f1 <- function(x) 2*x / (1 + x^2)
f2 <- function(x)  -2 / (1 + x^2)
f3 <- function(x) rep(1, length(x))     # vectorized constant function 1
x <- c(0.0,   4.0)
y <- c(1.25, -0.95)
sol <- bvp(f1, f2, f3, x, y)
if (FALSE) { # \dontrun{
plot(sol$xs, sol$ys, ylim = c(-2, 2),
     xlab = "", ylab = "", main = "Boundary Value Problem")
# The analytic solution is
sfun <- function(x) 1.25 + 0.4860896526*x - 2.25*x^2 + 
                    2*x*atan(x) - 1/2 * log(1+x^2) + 1/2 * x^2 * log(1+x^2)
xx <- linspace(0, 4)
yy <- sfun(xx)
lines(xx, yy, col="red")
grid()} # }