Runge-Kutta 4th Order with Fixed Step Size
Deriv = function(State)
{
nDim = length(State) / 2
r = State[1:nDim]
v = State[(nDim+1):(2*nDim)]
a = -GM/(sqrt(sum(r*r)))^3 * r
return(c(v, a))
}
## function (State, tau, Deriv)
## {
## F1 = Deriv(State)
## F2 = Deriv(State + F1 * tau/2)
## F3 = Deriv(State + F2 * tau/2)
## F4 = Deriv(State + F3 * tau)
## return(State + (F1 + 2 * F2 + 2 * F3 + F4) * tau/6)
## }
## <bytecode: 0x000001e2b2fd41d8>
## <environment: namespace:math>
StartTime = 0
EndTime = 5
tau = 0.01
nStep = (EndTime - StartTime)/tau + 1
nDim = 2
Time = vector(length=nStep)
State1 = matrix(nrow=nStep, ncol=(2*nDim))
State1[1,] = c(1, 0, 0, 2*pi)
State2 = matrix(nrow=nStep, ncol=(2*nDim))
State2[1,] = c(1, 0, 0, 2*pi)
t1 = Sys.time()
for (i in 1:(nStep-1)) State1[i+1,] = Euler(State1[i,], tau)
t2 = Sys.time()
difftime(t2, t1)
## Time difference of 0.0167 secs
for (i in 1:(nStep-1)) State2[i+1,] = RK4(State2[i,], tau, Deriv)
t3 = Sys.time()
difftime(t3, t2)
## Time difference of 0.01865 secs
#windows(width=7, height=8)
par(mfrow=c(1,2))
amax = max(abs(State1[,1:2]))
plot(State1[,1:2], type="l", xlim=c(-amax,amax), ylim=c(-amax,amax))
abline(v=0,h=0)
amax = max(abs(State2[,1:2]))
plot(State2[,1:2], type="l", xlim=c(-amax,amax), ylim=c(-amax,amax))
abline(v=0,h=0)
