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 methodenter image description here

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 ?

  • 3
    $\begingroup$ 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. $\endgroup$ – spektr Feb 28 '17 at 16:18
  • $\begingroup$ 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 $\endgroup$ – Eray Xx Feb 28 '17 at 16:29
  • $\begingroup$ Im a beginner actually $\endgroup$ – Eray Xx Feb 28 '17 at 16:35
  • 2
    $\begingroup$ I assume your variable t is the timestep $\Delta t$? If so, I may add your equations in differential form to the question. $\endgroup$ – spektr Feb 28 '17 at 17:02
  • $\begingroup$ yes youre right it would be really good $\endgroup$ – Eray Xx 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.

| cite | improve this answer | |
  • 2
    $\begingroup$ 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. $\endgroup$ – spektr Feb 28 '17 at 19:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.