swh:1:snp:81eadaa089e8253d8469bcef66aa332632c6c669
Tip revision: e4641665e964c87f969d6efbb037ce5df2f85813 authored by Hans W. Borchers on 16 October 2018, 19:00:04 UTC
version 2.1.8
version 2.1.8
Tip revision: e464166
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))
}