pracma (version 1.9.9)

rk4, rk4sys: Classical Runge-Kutta

Description

Classical Runge-Kutta of order 4.

Usage

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.

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.

Details

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

References

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

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)
## Not run: 
# x <- seq(0, 2, length.out = 51)
# plot(x, x*exp(-x^2), type = "l", col = "red")
# points(sol$x, sol$y, pch = "*")
# grid()## End(Not run)

##  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)
## Not run: 
# 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()## End(Not run)