https://github.com/cran/pracma
Tip revision: 71455748623ef69836470c75c5f9384f6e872d45 authored by HwB on 28 June 2011, 00:00:00 UTC
version 0.6-3
version 0.6-3
Tip revision: 7145574
rk4.R
##
## r k 4 . R Classical Runge-Kutta
##
rk4 <- function(f, a, b, y0, n) {
h <- (b-a)/n
x <- seq(a+h, b, by = h)
y <- numeric(n)
k1 <- h * f(a, y0)
k2 <- h * f(a + h / 2 , y0 + k1 / 2 )
k3 <- h * f(a + h / 2 , y0 + k2 / 2 )
k4 <- h * f(a + h , y0 + k3)
y[1] <- y0 + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6
for (i in 1:(n-1)) {
k1 <- h * f(x[i], y[i])
k2 <- h * f(x[i] + h / 2, y[i] + k1 / 2)
k3 <- h * f(x[i] + h / 2 ,y[i] + k2 / 2 )
k4 <- h * f(x[i] + h , y[i] + k3)
y[i+1] = y[i] + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6
}
x <- c( a, x)
y <- c(y0, y)
return(list(x = x, y = y))
}
rk4sys <- function(f, a, b, y0, n){
m <- length(y0)
h <- (b-a)/n
x <- seq(a+h, b, by = h)
y <- matrix(0, nrow = n, ncol = m)
k1 <- h * f(a, y0)
k2 <- h * f(a + h / 2 , y0 + k1 / 2 )
k3 <- h * f(a + h / 2 , y0 + k2 / 2 )
k4 <- h * f(a + h , y0 + k3)
y[1, ] <- y0 + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6
for (i in 1:(n-1)) {
k1 <- h * f(x[i], y[i, ])
k2 <- h * f(x[i] + h / 2, y[i, ] + k1 / 2 )
k3 <- h * f(x[i] + h / 2, y[i, ] + k2 / 2)
k4 <- h * f(x[i] + h , y[i, ] + k3 )
y[i+1, ] <- y[i, ] + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6
}
x <- c(a, x)
y <- rbind(y0, y)
return(list(x = x, y = y))
}