[Tutor] list of references to object properties

[Oscar Benjamin]

On 19 January 2013 15:47, Jose Amoreira <ljmamoreira at gmail.com> wrote:
> Thanks for the explanation, Alan
>
>> Is there a good reason to want to do such a thing?
>
> There is a reason, but maybe you won't consider it a good one...
>
> I was writing a small program to simulate the gravitational dynamics of a
> system of many planets,

I have many times written code for this kind of problem.

> using scipy's odeint to solve the equations of
> motion.

I haven't checked but I wouldn't automatically trust scipy's odeint to
solve these equations through a significant period of time. A simple
method to indicate (not prove) the accuracy of your program would be
to measure how well energy is conserved over time (assuming that your
equations are Hamiltonian).

> I defined a class, CelestialBody, that describes objects that
> represent planets in my simulation. These objects have three attributes:
> position, velocity and mass (the first two are 3D-vectors; as such, the
> number of attributes is actually 7).

It is useful in this kind of problem to draw a distinction between
parameters and variables. Mass is a parameter. Position and velocity
are variables (they change over time).

> The many-body system is represented in
> the simulation by a list of CelestialBody objects.
>
> The dynamical state of the system is represented by a 6N (N being the number
> of planets) component array storing the components of the position and
> linear momentum of each body, and the integration procedures (odeint in this
> case) usually take this array as argument.
>


I would make a class to represent the system as a whole.

class SystemAsAWhole(object):
    def __init__(self):
        self.numdimensions = 0
        self.subsystems = []

    def add_subsystem(self, subsys):
         # Add to the list of systems to consider when evaluating the derivative
        self.subsystems.append(subsys)
        # Adjust the variable indices in the subsystem
        self.subsystems.shift(self.numdimensions)
        self.numdimensions += subsys.numdimensions

    def f(self, x, t, dxdt):
        for sys in self.subsystems:
            sys.f(x, t, dxdt)

class SubSystem(object):
    def __init__(self):
        for n, var in self.VARIABLES:
            setattr(self, var, n)
    # Called when the subsystem is incorporated into a larger system
    # Need to adjust our indices to be part of a larger state vector
    def shift(self, N):
        for var in self.variables:
            setattr(self, var, getattr(self, var) + N)

    def f(self, x, t, dxdt):
        # Compute and return derivative (subclasses need to implement this)
        raise NotImplementedError

And then subclass SubSystem to implement the derivative and describe
the variables/parameters:

class SelestialBody(SubSystem):
    VARIABLES = ['x', 'y', 'z', 'vx', 'vy', 'vz']
    def __init__(self, mass):
        self.mass = mass
    def f(self, x, t, dxdt):
        dxdt[self.x] = x[self.vx]
        dxdt[self.y] = x[self.vy]
        dxdt[self.z] = x[self.vz]
        ax, ay, az = self.force(x)
        dxdt[self.vx] = ax
        dxdt[self.vy] = ay
        dxdt[self.vz] = az

>
>
> So, in my simulation code I have a list of planets (objects of class
> CelestialBody) because it is a natural way of representing the systems I
> want to simulate, and another list (an numpy.array, actually) storing the
> positions and momenta of each planet, needed for the numerical processing of
> the simulation. But the positions and momenta are already present in the
> planet list, since the objects themselves store that information. Then, the
> second list, even if necessary, is a duplication of things that already are
> defined in the code. I don't like it. It'd be better (no duplication) if it
> was just a list of references to the values stored in the planet objects, I
> think.

If you understand the code above you'll see that each planet instance
has attributes giving the indices of the particular variable into the
state vector. This is much better than storing the values as
attributes in the instances themselves. The position of the Earth is a
property of a state of the system which is in turn part of a
particular trajectory. It should be stored as part of a trajectory not
as a property of the object itself. To access the y velocity of the
Earth from the state vector x you can just use x[earth.vy].


Oscar