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

Run this code
##  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)

Run the code above in your browser using DataCamp Workspace