bulirsch-stoer | R Documentation |
Bulirsch-Stoer algorithm for solving Ordinary Differential Equations (ODEs) very accurately.
bulirsch_stoer(f, t, y0, tol = 1e-08, ...)
midpoint(f, t0, tfinal, y0, tol = 1e-08, kmax = 101, ...)
f |
function describing the differential equation |
t |
vector of |
y0 |
starting values as (row or column) vector. |
... |
additional parameters to be passed to the function. |
tol |
relative tolerance (in the Ricardson extrapolation). |
t0 , tfinal |
start and end point of the interval. |
kmax |
number of grid pounts. |
The Bulirsch-Stoer algorithm is a well-known method to obtain high-accuracy solutions to ordinary differential equations with reasonable computational efforts. It exploits the midpoint method to get good accuracy in each step.
The (modified) midpoint method computes the values of the dependent
variable y(t)
from t0
to tfinal
by a sequence of
substeps, applying Richardson extrapolation in each step.
Bulirsch-Stoer and midpoint shall not be used with non-smooth functions or
singularities inside the interval. The best way to get intermediate points
t = (t[1], ..., t[n])
may be to call ode23
or ode23s
first and use the x
-values returned to start bulirsch_stoer
on.
bulirsch_stoer returns a list with x
the grid points input, and
y
a vector of function values at the se points.
Will be extended to become a full-blown Bulirsch-Stoer implementation.
Copyright (c) 2014 Hans W Borchers
J. Stoer and R. Bulirsch (2002). Introduction to Numerical Analysis. Third Edition, Texts in Applied Mathematics 12, Springer Science + Business, LCC, New York.
ode23
, ode23s
## Example: y'' = -y
f1 <- function(t, y) as.matrix(c(y[2], -y[1]))
y0 <- as.matrix(c(0.0, 1.0))
tt <- linspace(0, pi, 13)
yy <- bulirsch_stoer(f1, tt, c(0.0, 1.0)) # 13 equally-spaced grid points
yy[nrow(yy), 1] # 1.1e-11
## Not run:
S <- ode23(f1, 0, pi, c(0.0, 1.0))
yy <- bulirsch_stoer(f1, S$t, c(0.0, 1.0)) # S$x 13 irregular grid points
yy[nrow(yy), 1] # 2.5e-11
S$y[nrow(S$y), 1] # -7.1e-04
## Example: y' = -200 x y^2 # y(x) = 1 / (1 + 100 x^2)
f2 <- function(t, y) -200 * t * y^2
y0 <- 1
tic(); S <- ode23(f2, 0, 1, y0); toc() # 0.002 sec
tic(); yy <- bulirsch_stoer(f2, S$t, y0); toc() # 0.013 sec
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.