rkf54: Runge-Kutta-Fehlberg

In pracma: Practical Numerical Math Functions

Runge-Kutta-Fehlberg

Description

Runge-Kutta-Fehlberg with adaptive step size.

Usage

rkf54(f, a, b, y0, tol = .Machine$double.eps^0.5,
                   control = list(), ...)

Arguments

f	function in the differential equation y' = f(x, y).
a, b	endpoints of the interval.
y0	starting values at a.
tol	relative tolerance, used for determining the step size.
control	list for influencing the step size with components hmin, hmax, the minimal, maximal step size jmax, the maximally allowed number of steps.
...	additional parameters to be passed to the function.

Details

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.

Value

List with components x for grid points between a and b and y the function values of the numerical solution.

Note

This function serves demonstration purposes only.

References

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.

See Also

rk4, ode23

Examples

# 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)

pracma documentation built on Nov. 10, 2023, 1:14 a.m.