rk4, rk4sys | R Documentation |
Classical Runge-Kutta of order 4.
rk4(f, a, b, y0, n)
rk4sys(f, a, b, y0, n)
f |
function in the differential equation |
a , b |
endpoints of the interval. |
y0 |
starting values; for |
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.
ode23
, deval
## 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.