Python Float Precision

I'm just kind of messing around...something Nodbox lends itself to very well. I'm playing with this imaginary physics animation pasted below. If you run it you will see that it is stable for a while (20 seconds or so) and then falls apart. I suspect this is due to floating point precision errors creaping in, but I'm not 100% sure it's not a consequence of the system. I wonder if someone here knows a good way to increase te precision. Is there something in the standard library? (I've been digging around to no avail). Or if you think it is not a float issue I would be curious to hear about that as well

from decimal import *
speed(100)
size(800,800)
strokewidth(.05)
from Numeric import *
 
class System:
    def __init__(self):
        self.universe = []
        self.center = None
    def addBody(self, body):
        self.universe.append(body)
    def distance(self, body1, body2):
        dist = sqrt(sum((body1.position - body2.position)**2))
        if dist < 10: return 10
        else: return dist
    def gravity(self, body1, body2, dist):
        k = .00001 # gravitational constant
        dist = self.distance(body1, body2)
        return (body1.mass * body2.mass) * (dist**2) * k
    def repulsion(self, body1, body2, dist):
        k = -200 # gravitational constant
        dist = self.distance(body1, body2)
        return (body1.mass * body2.mass) / (dist**2) * k
 
    def makeGravity(self):
        uCopy = self.universe[:]
        while uCopy:
            body1 = uCopy.pop()
            for body2 in uCopy:
                dist = self.distance(body1, body2)
                if body1.charge == body2.charge:
                    gForce = self. repulsion(body1, body2, dist )
                else:
                    gForce = self.gravity(body1, body2, dist )
 
                #calculate  acceleration 
                
                body1.velocity -= gForce/body1.mass * (body1.position - body2.position) /dist
                body2.velocity -= gForce/body2.mass * (body2.position - body1.position) /dist
                #print body1.mass, body2.mass, " : "
        map(lambda b: b.updatePosition(1), self.universe)
          
 
class Particle:
    def __init__(self, mass, position, velocity):
        #velocity and position should be x,y,z vectors
        self.mass = mass
        self.velocity = array(velocity, Float)
        self.position = array(position, Float)
        self.x = self.position[0]
        self.y = self.position[1]
        self.z = self.position[2]
        self.isStationary = False
        self.charge = 0
    def draw(self):
        oval(self.position[0]-2, self.position[1]-2, 4, 4 )
    def updatePosition(self, t):
        self.draw()
        if self.isStationary:
            return
        self.position += self.velocity * t
        self.x = self.position[0]
        self.y = self.position[1]
        self.z = self.position[2]
       
       
def setup():        
    global p1_mass, p2_mass, p3_mass, a, p1, p2, p3, p4, system
    p1_mass=.04
    p2_mass=.01
    p3_mass=.025
    a = 1 
    system = System()
    n = 16
    for i in range(1, n+1):
        Center_x = 0
        Center_y = 0
        r = 155
        x = Center_x + cos((2*i*pi)/(n)) * r
        y = Center_y + sin((2*i*pi)/(n)) * r
        if i % 2 == 0:
            vx = sin((2*i*pi)/(n))
            vy = -cos((2*i*pi)/(n))
        else:
            vx = 0
            vy = 0
 
        system.addBody( Particle(15, (x, y, 0), (vx*2, vy*2, 0)))
    sun  = Particle(2, (-0, -0, 0), (0, 0, 0))
    sun.isStationary = True
    sun.charge = 1
    system.addBody(sun)
    
def draw():
    global p1_mass, p2_mass, p3_mass, a, p1, p2, p3, system
    translate(400, 400)
    nofill()
    stroke(.2, .2, .5, .5)
    system.makeGravity()

Hi Mark,

Wonderful! I actually like how they start getting out of sync :-)

But I'll need some time to review the code. I have a hunch there may be something else going on besides decimal precision. In the meantime, here's my version of a physics system I have been messing with. In the near future we should probably merge code and create one solid Physics library out of it.

from math import sqrt
 
class particle(Point):
    
    def __init__(self, x, y, mass, life=1000):
        
        Point.__init__(self, x, y)
        self.mass = mass
        self.velocity = Point(0, 0)
        self.force = Point(0, 0)
        self.life = life
        self.fixed = False
 
class force:
    
    def __init__(self, particle1, particle2, strength=1.0, threshold=100):
        
        self.particle1 = particle1
        self.particle2 = particle2
        self.strength  = strength
        self.threshold = threshold    
    
    def update(self):
        
        # Horizontal and vertical distance.
        dx = self.particle2.x - self.particle1.x
        dy = self.particle2.y - self.particle1.y        
 
        # Distance between two particles.
        # Distance has a minimum threshold to keep forces from growing too large:
        # e.g. distance 100 divides force by 10000, distance 5 only by 25!
        d = sqrt(dx*dx + dy*dy)
        d = max(d, self.threshold)
        
        # The force between particles increases according to their weight.
        # The force decreases as distance between them increases.
        f = 10.0*self.strength * self.particle1.mass * self.particle2.mass
        f = f / (d*d)
        fx = f * dx / d
        fy = f * dy / d
 
        if not self.particle1.fixed:
            self.particle1.force.x += fx
            self.particle1.force.y += fy
        
        if not self.particle2.fixed:
            self.particle2.force.x -= fx
            self.particle2.force.y -= fy
    
class spring:
    
    def __init__(self, particle1, particle2, length, strength=1.0, damping=0.1):
        
        self.particle1 = particle1
        self.particle2 = particle2
        self.length    = length
        self.strength  = strength
        self.damping   = damping
        
    def update(self):
        
        # Horizontal and vertical distance.
        dx = self.particle2.x - self.particle1.x
        dy = self.particle2.y - self.particle1.y
 
        # Distance between two particles.
        d = sqrt(dx*dx + dy*dy)
        if d == 0: return
        
        # The attractive strength decreases for heavy particles.
        # The attractive strength increases when the spring is stretched.
        f = 10.0 * self.strength / (self.particle1.mass * self.particle2.mass)
        f *= d - self.length
 
        # Horizontal and vertical forces take into account
        # individual particles' velocity (i.e. bearing), although dampened.
        fx = f + (self.particle1.velocity.x - self.particle2.velocity.x) * dx/d * self.damping
        fy = f + (self.particle1.velocity.y - self.particle2.velocity.y) * dy/d * self.damping
        fx *= dx / d
        fy *= dy / d
        
        if not self.particle1.fixed:
            self.particle1.force.x += fx
            self.particle1.force.y += fy
            
        if not self.particle2.fixed:
            self.particle2.force.x -= fx
            self.particle2.force.y -= fy

class system:
    
    def __init__(self, gx, gy):
        
        self.gravity = Point(gx, gy)
        self.particles = []
        self.springs = []
        self.forces = []
    
    def update(self, t=1.0):
        
        for p in self.particles:
            p.force.x += self.gravity.x * p.mass
            p.force.y += self.gravity.y * p.mass
 
        for f in self.forces:
            f.update()
 
        for s in self.springs:
            s.update()
        
        for p in self.particles:
            p.force.x += p.velocity.x
            p.force.y += p.velocity.y
            p.x += p.force.x
            p.y += p.force.y
            p.force = Point(0,0)
            p.life -= 1
 
s = system(0, 0.75)
for i in range(100):
    s.particles.append(particle(random(300.0,310),random(300.0,310.0),random(5,20)))
    s.particles[-1].velocity = Point(0, -5)
 
for i in range(len(s.particles)):
    for j in range(i+1, len(s.particles)): #adjust 1 for optimisation
        f = force(
            s.particles[i],
            s.particles[j],
            -5.0
        )
        s.forces.append(f)
 
size(700,700)
speed(50)
def draw():
    s.update()
    stroke(0)
    d = 1
    for p in s.particles:
        nostroke()
        fill(0,0,0, p.life/50.0)
        oval(p.x* d-p.mass/2, p.y* d-p.mass/2, p.mass, p.mass)
    for sp in s.springs:
        stroke(0,0,0, 0.1*sp.strength)
        line(sp.particle1.x* d, sp.particle1.y* d, 
             sp.particle2.x* d, sp.particle2.y* d)