天体动力学引擎的类或元类设计

Question

大师们：

用于模拟航天器运动的微分方程可以用一组加速度项来描述：

d2r/dt2 =  a0 + a1 + a2 + ... + an

通常a0是物体的点质量加速度（a0 = -mu * r/r^3）； “高阶”项可能与其他行星、太阳辐射压力、推力等有关。

我正在实现一组旨在用于此类系统的算法。我将从 Python 开始进行设计和原型设计，然后我将转向 C++ 或 Fortran 95。

我想设计一个类（或元类），它允许我为给定实例指定不同的加速项，类似于：

# please notice this is meant as "pseudo-code"
def some_acceleration(t):
    return (1*t, 2*t, 3*t)

def some_other_acceleration(t):
    return (4*t, 5*t, 6*t)

S = Spacecraft()
S.Acceleration += someacceleration + some_other_acceleration

在这种情况下，实例 S 将默认使用两个加速项，我将添加我想要的其他两个项：some acceleration 和 some_other_acceleration；他们返回一个向量（这里表示为一个三元组）。请注意，在我的“实现”中，我重载了 + 运算符。

这样，算法将为抽象的“航天器”设计，所有实际力场都将根据具体情况提供，使我能够使用简化模型、比较建模方法等。

您将如何实现一个类或元类来处理这个问题？

对于这个相当冗长且解释不清的问题，我深表歉意，但我的大脑有点模糊。

谢谢。

Answer 1

PyDSTool 让您可以将具有空间或物理意义的“组件”以及与之相关的数学表达式构建成更大的组件，这些组件知道如何总结等等。结果是一种指定微分方程的方法使用符号工具的模块化方式，然后 PyDSTool 将自动创建 C 代码以使用快速积分器模拟系统。在使用 C 和 Fortran 进行“繁重的工作”之前，无需将 python 仅视为一个缓慢的原型设计步骤。一旦您完全指定了问题，PyDSTool 已经将所有内容（包括您定义的结果向量场）移至 C 级别。

您的二阶 DE 示例与包含多种离子通道的生物细胞膜上的电位差的“电流平衡”一阶方程非常相似。根据基尔霍夫现行定律，p.d. 的变化率。是

dV/dt = 1/memb_capacitance *（电流总和）

PyDSTool/tests/ModelSpec_tutorial_HH.py 中的示例是捆绑在包中的几个示例之一，它从模块化规范组件（您继承以创建自己的 ModelSpec 类 - 例如“point_mass”中的“point_mass”）物理环境），并使用类似宏的规范将最终 DE 定义为添加到“膜”组件的“离子通道”组件中存在的任何电流的总和。求和是在makeSoma 函数中定义的，只需使用诸如“for(channel,current,+)/C”之类的语句，您可以直接在代码中使用。

希望这会有所帮助。如果是这样，请随时通过 sourceforge 上的帮助论坛询问我（PyDSTool 的作者）是否需要更多帮助来启动它。

Answer 2

对于那些想避免使用 numpy 并在纯 python 中执行此操作的人，这可能会给您一些好主意。我敢肯定这个小短剧也有缺点和缺陷。 "operator" 模块加快了您的数学计算，因为它们是使用 c 函数完成的：

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

打印以下内容：

加速度(2.1234, 3.1234, 4.1234)

物体的速度与加速度 一秒后的加速度(13.3574, 24.3574, 35.3574)： 速度(26.7148, 48.7148, 70.7148)

您可以享受更快的计算：

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

Answer 3

您是在问如何为航天器类存储任意数量的加速度源？

你不能只使用一个函数数组吗？ （进入c++时的函数指针）

即：

#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)

Answer 4

我会改用一些可以处理向量的库（在 python 中，尝试 numpy）并将加速度表示为向量。然后，您不会重新发明*，+ 运算符可以按照您的意愿工作。如果我误解了您的问题，请纠正我。