This is a Python implementation of the RK4 numerical integrator that works with differential functions of all orders. That is, any function in the form F(x, y, y’, y”, …, 
).
If you’re new to numerical integration or even RK4 integration, please read my other post first. It’s easier to understand because it’s a less generalized function.
def rk4(t0, h, s0, f):
"""RK4 implementation.
t = current value of the independent variable
h = amount to increase the independent variable (step size)
s0 = initial state as a list. ex.: [initial_position, initial_velocity]
f = function(state, t) to integrate"""
r = range(len(s0))
s1 = s0 + [f(t0, s0)]
s2 = [s0[i] + 0.5*s1[i+1]*h for i in r]
s2 += [f(t0+0.5*h, s2)]
s3 = [s0[i] + 0.5*s2[i+1]*h for i in r]
s3 += [f(t0+0.5*h, s3)]
s4 = [s0[i] + s3[i+1]*h for i in r]
s4 += [f(t0+h, s4)]
return t+h, [s0[i] + (s1[i+1] + 2*(s2[i+1]+s3[i+1]) + s4[i+1])*h/6.0 for i in r]
An example usage of this function is:
def damped_spring_accel(t, state):
x = state[0]
v = state[1]
stiffness = 1
damping = 0.005
return -stiffness*x - damping*v
t = 0
h = 1/40.0
s = [50, 5]
while t<100:
t, s = rk4(t, h, s, damped_spring_accel)
print "Final state: t=%.2f x=%.2f v=%.2f"%(t,s[0],s[1])
On a side note, I recently found out about Sage (sagemath.org and sagenb.org). Its plotting capabilities and convenient mathematical notation are especially useful. Plus, it uses Python!
I made a little Processing sketch based on the same physics as JavaScript physics.
(View the full post to try it out.)
Continue reading →
I made a very basic physics demo in JavaScript to practice using the MooTools library. It uses Verlet integration as outlined in Thomas Jakobsen’s paper.
To try out the demo, go to http://doswa.com/projects/physics_js/
Chain of particles attached by joints
In the previous post, I wrote about a way to generalize derivatives and integrals into one function. What happens if a number other than an integer is passed to that function?
Here is the generalization from the last post:
Set k=2 and 
:
Now there’s a problem. Since a! is only defined for integers 0 or greater, a different way of calculating factorials is needed. Luckily, there exists a Gamma function defined for all real and complex numbers such that:
Substitute the Gamma function into the other equation and simplify:
So, the half derivative of 
is approximately equal to 
.
This can be used for ‘fractional integration’ as well, if a negative number is used for n.
Start with a simple expression, 
, and take a few derivatives:
A pattern is emerging:
, where c is the coefficient.
Now the hard part is finding the pattern in the coefficient. This needs to be taken out of the ‘…’ form. Focus on that:
This is a series of numbers, each one larger than the next. This looks like a factorial, so divide k! by that:
Note how the top goes from 1 to k and the bottom goes from k-n+1 to k. That means that k-n+1 to k is a subset of 1 to k, so just divide that part out:
Now solve for c:
, so
Puts this back into the original equation to get:
So, the nth derivative of is equal to:
Verify this with a couple of derivatives:
And a couple of integrals:
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
