Here’s a Python implementation of RK4, hardcoded for double-integrating the second derivative (acceleration up to position). For a more generalized solution, see my other implementation. I tried to keep this as simple as I could, so people can easily see the relation between the ‘math form’ and ‘code form’ of the algorithm.
def rk4(x, v, a, dt): """Returns final (position, velocity) tuple after time dt has passed. x: initial position (number-like object) v: initial velocity (number-like object) a: acceleration function a(x,v,dt) (must be callable) dt: timestep (number)""" x1 = x v1 = v a1 = a(x1, v1, 0) x2 = x + 0.5*v1*dt v2 = v + 0.5*a1*dt a2 = a(x2, v2, dt/2.0) x3 = x + 0.5*v2*dt v3 = v + 0.5*a2*dt a3 = a(x3, v3, dt/2.0) x4 = x + v3*dt v4 = v + a3*dt a4 = a(x4, v4, dt) xf = x + (dt/6.0)*(v1 + 2*v2 + 2*v3 + v4) vf = v + (dt/6.0)*(a1 + 2*a2 + 2*a3 + a4) return xf, vf
Here is an example usage of the function and a comparison to Euler integration:
def accel(x, v, dt): """Determines acceleration from current position, velocity, and timestep. This particular acceleration function models a spring.""" stiffness = 1 damping = -0.005 return -stiffness*x - damping*v t = 0 dt = 1.0/40 # Timestep of 1/40 second state = 50, 5 # Position, velocity euler = 50, 5 # For comparison with Euler integration print "Initial -position: %6.2f, velocity: %6.2f"%state # Run for 100 seconds while t < 100: t += dt state = rk4(state[0], state[1], accel, dt) # Integrate using Euler's method euler = ( euler[0] + euler[1]*dt, euler[1] + accel(euler[0],euler[1],dt)*dt ) print "Final RK4 -position: %6.2f, velocity: %6.2f"%state print "Final Euler-position: %6.2f, velocity: %6.2f"%euler
The output of this really shows how much more accurate RK4 integration can be:
Initial -position: 50.00, velocity: 5.00 Final RK4 -position: 52.18, velocity: 38.05 Final Euler-position: 178.38, velocity: 137.62
As the timestep is decreased (meaning more computation), Euler approaches RK4 (shown at timestep of 1/400 seconds):
Initial -position: 50.00, velocity: 5.00 Final RK4 -position: 52.28, velocity: 37.92 Final Euler-position: 59.20, velocity: 43.02
Pingback: Improved RK4 Implementation « Doswa
Fantastic code! Just one thing your damping should be positive in order to be consistent with your other signs; currently it is adding energy to the system not absorbing energy. Also I added a plotting function and did a pendulum. So using the rk4 function as above this will make a nifty animated pendulum:
def accel(x, v, dt): """Determines acceleration from current position, velocity, and timestep. This particular acceleration function models a pendulum. And plots!""" g=-9.8 l=100.0 return (g/l)*x l=100.0 #pendulum_length p=10 #postion for x axis c=(np.pi/180.0) #conversion to radians t = 0 dt = 1.0/10 # Timestep of 1/40 second state =20*c, 0*c # start Position, velocity do_plot=1 #switch off of zero to turn off plotting if (do_plot==1): ion() line, =plot([p,p],[l,0],'o-') # Run for 100 seconds while t < 100: t += dt state = rk4(state[0], state[1], accel, dt) if (do_plot==1): line.set_xdata([p,p+np.sin(state[0])]) line.set_ydata([l,l*(1.0-np.cos(state[0]))]) draw()Thanks soooo much for the simple way you put that RK4 code together, really helped me see how it works, as I’ve had difficulty with some of the other ways I’ve seen it written.
Just a few questions if you’ve got the time/inclination to help out!
I’m putting together an orbital code, and need to use an Runge Kutta Fehlberg routine, but can’t see how to adapt it to work. I think I’m just not sure how the time fractions fit in with the function fractions.
Thanks again for the easy to read code, worked a treat!
Hi simmuskhan:
I had never heard of the RKF method until now, but after looking into it, I feel that it overcomplicates things. For many situations, RK4 is “good enough”.
However, here’s a basic implementation of RKF. You’ll probably want to add some code to put minimum/maximum values on dt. This one is hardcoded for 1st order differential equations, but it shouldn’t be too hard to modify for higher orders.
And here’s an example usage of that:
Pingback: Improved RK4 Implementation | Doswa
What a pretty little snippet of code. I was playing with the Runge-Kutta method today and stumbled upon your code. I am playing around with some stellar atmosphere equations and decided I might be able to produce a computer model. I always assumed numerically solving a differential equation was very difficult, however this proved me wrong. Now I am off to figure out how to solve simultaneous DE’s, what a treat!
Cheers,
Jerry