PyDSTool lets you build up "components" that have spatial or physical meaning, and have mathematical expressions associated with them, into bigger components that know how to sum things up etc. The result is a way to specify differential equations in a modular fashion using symbolic tools, and then PyDSTool will create C code automatically to simulate the system using fast integrators. There s no need to see python only as a slow prototyping step before doing the "heavy lifting" in C and Fortran. PyDSTool already moves everything, including the resulting vector field you defined, down to the C level once you ve fully specified the problem.

Your example for second order DEs is very similar to the "current balance" first-order equation for the potential difference across a biological cell membrane that contains multiple types of ion channel. By Kirchoff s current law, the rate of change of the p.d. is

dV/dt = 1/memb_capacitance * (sum of currents)

The example in PyDSTool/tests/ModelSpec_tutorial_HH.py is one of several bundled in the package that builds a model for a membrane from modular specification components (the ModelSpec class from which you inherit to create your own - such as "point_mass" in a physics environment), and uses a macro-like specification to define the final DE as the sum of whatever currents are present in the "ion channel" components added to the "membrane" component. The summing is defined in the makeSoma function, simply using a statement such as for(channel,current,+)/C , which you can just use directly in your code.

Hope this helps. If it does, feel free to ask me (the author of PyDSTool) through the help forum on sourceforge if you need more help getting it started.

For those who would like to avoid numpy and do this in pure python, this may give you a few good ideas. I m sure there are disadvantages and flaws to this little skit also. The "operator" module speeds up your math calculations as they are done with c functions:

from operator import sub, add, iadd, mul
import copy

class Acceleration(object):
   def __init__(self, x, y, z):
      super(Acceleration, self).__init__()
      self.accel = [x, y , z]
      self.dimensions = len(self.accel)

   @property
   def x(self):
      return self.accel[0]

   @x.setter
   def x(self, val):
      self.accel[0] = val


   @property
   def y(self):
      return self.accel[1]

   @y.setter
   def y(self, val):
      self.accel[1] = val

   @property
   def z(self):
      return self.accel[2]

   @z.setter
   def z(self, val):
      self.accel[2] = val

   def __iadd__(self, other):
      for x in xrange(self.dimensions):
         self.accel[x] = iadd(self.accel[x], other.accel[x])
      return self

   def __add__(self, other):
      newAccel = copy.deepcopy(self)
      newAccel += other
      return newAccel

   def __str__(self):
      return "Acceleration(%s, %s, %s)" % (self.accel[0], self.accel[1], self.accel[2])

   def getVelocity(self, deltaTime):
      return Velocity(mul(self.accel[0], deltaTime), mul(self.accel[1], deltaTime), mul(self.accel[2], deltaTime))

class Velocity(object):
   def __init__(self, x, y, z):
      super(Velocity, self).__init__()
      self.x = x
      self.y = y
      self.z = z

   def __str__(self):
      return "Velocity(%s, %s, %s)" % (self.x, self.y, self.z)

if __name__ == "__main__":
   accel = Acceleration(1.1234, 2.1234, 3.1234)
   accel += Acceleration(1, 1, 1)
   print accel

   accels = []
   for x in xrange(10):
      accel += Acceleration(1.1234, 2.1234, 3.1234)

   vel = accel.getVelocity(2)
   print "Velocity of object with acceleration %s after one second:" % (accel)
   print vel

prints the following:

Acceleration(2.1234, 3.1234, 4.1234)

Velocity of object with acceleration Acceleration(13.3574, 24.3574, 35.3574) after one second: Velocity(26.7148, 48.7148, 70.7148)

You can get fancy for faster calculations:

def getFancyVelocity(self, deltaTime):
   from itertools import repeat
   x, y, z = map(mul, self.accel, repeat(deltaTime, self.dimensions))
   return Velocity(x, y, z)

Are you asking how to store an arbitrary number of acceleration sources for a spacecraft class?

Can t you just use an array of functions? (Function pointers when you get to c++)

i.e:

#pseudo Python
class Spacecraft
    terms = []
    def accelerate(t):
       a = (0,0,0)
       for func in terms:
         a+= func(t)


s = Spacecraft
s.terms.append(some_acceleration)
s.terms.append(some_other_acceleration)
ac = s.accelerate(t)

I would use instead some library which can work with vectors (in python, try numpy) and represent the acceleration as vector. Then, you are not reinventing the wheel, the + operator works like you wanted. Please correct me if I misunderstood your problem.