Classical Runge-Kutta

Classical Runge-Kutta of order 4.

rk4(f, a, b, y0, n)

rk4sys(f, a, b, y0, n)

Arguments

f

function in the differential equation \(y' = f(x, y)\);
defined as a function \(R \times R^m \rightarrow R^m\), where \(m\) is the number of equations.

a, b

endpoints of the interval.

y0

starting values; for \(m\) equations y0 needs to be a vector of length m.

n

the number of steps from a to b.

Details

Classical Runge-Kutta of order 4 for (systems of) ordinary differential equations with fixed step size.

Value

List with components x for grid points between a and b and y an n-by-m matrix with solutions for variables in columns, i.e. each row contains one time stamp.

References

Fausett, L. V. (2007). Applied Numerical Analysis Using Matlab. Second edition, Prentice Hall.

Note

This function serves demonstration purposes only.

See also

ode23, deval

Examples

##  Example1: ODE
# y' = y*(-2*x + 1/x) for x != 0, 1 if x = 0
# solution is x*exp(-x^2)
f <- function(x, y) {
  if (x != 0) dy <- y * (- 2*x + 1/x)
  else        dy <- rep(1, length(y))
  return(dy)
}
sol <- rk4(f, 0, 2, 0, 50)
if (FALSE) { # \dontrun{
x <- seq(0, 2, length.out = 51)
plot(x, x*exp(-x^2), type = "l", col = "red")
points(sol$x, sol$y, pch = "*")
grid()} # }

##  Example2: Chemical process
  f <- function(t, u) {
    u1 <- u[3] - 0.1 * (t+1) * u[1]
    u2 <- 0.1 * (t+1) * u[1] - 2 * u[2]
    u3 <- 2 * u[2] - u[3]
    return(c(u1, u2, u3))
  }
u0 <- c(0.8696, 0.0435, 0.0870)
a <- 0; b <- 40
n <- 40
sol <- rk4sys(f, a, b, u0, n)
if (FALSE) { # \dontrun{
matplot(sol$x, sol$y, type = "l", lty = 1, lwd = c(2, 1, 1),
        col = c("darkred", "darkblue", "darkgreen"),
        xlab = "Time [min]", ylab= "Concentration [Prozent]",
        main = "Chemical composition")
grid()} # }