runge-kutta.R

library(pforeach)
library(ggplot2)

x <- 1 #初期値
h <- 0.01 #ステップ幅
maxit <- 100 #くり返しの最大値

iter <- seq(0, maxit, by = h) #初期時刻，くり返しの最大値，ステップ幅

dx.dt <- function(tt, xx) { #dx/dt の定義
    return(2*xx)
}

npforeach (t = iter, .c = rbind)({ #Runge-Kutta 法をくり返す
    k1 <- dx.dt(t, x)
    k2 <- dx.dt(t + h/2, x + h/2 * k1)
    k3 <- dx.dt(t + h/2, x + h/2 * k2)
    k4 <- dx.dt(t + h, x + h * k3)
    
    x <- x + h/6 * (k1 + 2*k2 + 2*k3 + k4)
    data.frame(x = x, t = t)
}) -> result

ggplot(result, aes(t, x)) + geom_line() + ylim(1, 10) + xlim(0, 1)