Classical Runge-Kutta
Classical Runge-Kutta of order 4.
rk4(f, a, b, y0, n) rk4sys(f, a, b, y0, n)
f |
function in the differential equation y' = f(x, y); |
a, b |
endpoints of the interval. |
y0 |
starting values; for m equations |
n |
the number of steps from |
Classical Runge-Kutta of order 4 for (systems of) ordinary differential equations with fixed step size.
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.
This function serves demonstration purposes only.
Fausett, L. V. (2007). Applied Numerical Analysis Using Matlab. Second edition, Prentice Hall.
## 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)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.