10.7 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))
}

RK4

## 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)