Runge-Kutta-Fehlberg
Runge-Kutta-Fehlberg with adaptive step size.
rkf54(f, a, b, y0, tol = .Machine$double.eps^0.5,
control = list(), ...)f |
function in the differential equation y' = f(x, y). |
a, b |
endpoints of the interval. |
y0 |
starting values at |
tol |
relative tolerance, used for determining the step size. |
control |
list for influencing the step size with components |
... |
additional parameters to be passed to the function. |
Runge-Kutta-Fehlberg is a kind of Runge-Kutta method of solving ordinary differential equations of order (5, 4) with variable step size.
“At each step, two different approximations for the solution are made and compared. If the two answers are in close agreement, the approximation is accepted. If the two answers do not agree to a specified accuracy, the step size is reduced. If the answers agree to more significant digits than required, the step size is increased.”
Some textbooks promote the idea to use the five-order formula as the
accepted value instead of using it for error estimation. This approach
is taken here, that is why the function is called rkf54. The idea
is still debated as the accuracy determinations appears to be diminished.
List with components x for grid points between a and b
and y the function values of the numerical solution.
This function serves demonstration purposes only.
Stoer, J., and R. Bulirsch (2002). Introduction to Numerical Analysis. Third Edition, Springer-Verlag, New York.
Mathematica code associated with the book:
Mathews, J. H., and K. D. Fink (2004). Numerical Methods Using Matlab.
Fourth Edition, Prentice Hall.
# Example: y' = 1 + y^2
f1 <- function(x, y) 1 + y^2
sol11 <- rkf54(f1, 0, 1.1, 0.5, control = list(hmin = 0.01))
sol12 <- rkf54(f1, 0, 1.1, 0.5, control = list(jmax = 250))
# Riccati equation: y' = x^2 + y^2
f2 <- function(x, y) x^2 + y^2
sol21 <- rkf54(f2, 0, 1.5, 0.5, control = list(hmin = 0.01))
sol22 <- rkf54(f2, 0, 1.5, 0.5, control = list(jmax = 250))
## Not run:
plot(0, 0, type = "n", xlim = c(0, 1.5), ylim = c(0, 20),
main = "Riccati", xlab = "", ylab = "")
points(sol11$x, sol11$y, pch = "*", col = "darkgreen")
lines(sol12$x, sol12$y)
points(sol21$x, sol21$y, pch = "*", col = "blue")
lines(sol22$x, sol22$y)
grid()
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.