# Numerical Free Fall Analysis with RK4

I am trying to calculate real speed and time in free fall of a body. I wrote a code in Fortran and I am trying to improve it by using RK4 method

x=time y=total free fall

Purple line using:

DO WHILE ((K1-K2) > y)
y = y + t*((G*m*t/(2*((r-y)**2))) + sqrt(2*(G*m*y)/(r*(r-y))))
time = time + t
END DO


Green line using(RK4 method):

DO WHILE ((D1-D2) > y)
k1 = t*((G*m*(0)/(2*((r-y)**2))) + sqrt(2*(G*m*y)/(r*(r-y))))
k2 = (t/2)*((G*m*(t/2)/(2*((r-(y+k1))**2))) + sqrt(2*(G*m*(y+k1)/(r*(r-(y+k1))))))
k3 = (t/2)*((G*m*(t/2)/(2*((r-(y+k2))**2))) + sqrt(2*(G*m*(y+k2)/(r*(r-(y+k2))))))
k4 = t*((G*m*t/(2*((r-(y+k3))**2))) + sqrt(2*(G*m*(y+k3))/(r*(r-(y+k3)))))
y = y + (1.0/6.0)*(k1+2*k2+2*k3+k4)
time = time + t
END DO


All of other codes are the same (stepsize, for example) and I am getting better results without RK4. I guess I didn't understand RK4 right and I write it wrong but I cant see my mistake. What should I change to get better results with RK4 or did I do it wrong ?

• Can you write out the actual dynamics, in equation form, that you're solving? In terms of implementation, you might want to consider separating the RK4 code from the dynamics of your system. Separating the two will make it easier to read and debug. Feb 28 '17 at 16:18
• y(n+1) = y(n) + t((GMt/(2((r-y)^2))) + sqrt(2(GMy)/(r(r-y)) That was my normal equation. I tried to make it in runge kutta. I dont know if i did right Feb 28 '17 at 16:29
• Im a beginner actually Feb 28 '17 at 16:35
• I assume your variable t is the timestep $\Delta t$? If so, I may add your equations in differential form to the question. Feb 28 '17 at 17:02
• yes youre right it would be really good Feb 28 '17 at 17:18

1. You are using t for the time step. It would be clearer to use dt.
2. k1,... should be the derivatives, not the steps. So you should use

k1 = ((G*m*t/(2*((r-y)**2))) + sqrt(2*(G*m*y)/(r*(r-y))))

to calculate the first derivative. And also for k2,... don't multiply by time either.

3. Then when calculating k2 (and k3) don't use y+k1 but y + 0.5*dt*k1 inside the formula. The first two tentative steps in RK4 are only half step long.

To make the code look nicer, you should write a function that calculates the acceleration, instead of pasting the same formula 4 times. Less visual clutter also helps thinking.

Also your equations are wrong, but that's a separate issue. With these corrections at least you should be able to get the same result from your Euler and RK versions.

• I agree with what you said. I also noticed the equations seemed wrong. I started to write them out and they did not make sense the way they are coded. Feb 28 '17 at 19:06