在python中集成3体系统

Question

这是一个更数学的问题。

我们有 3 个身体系统，具有已知的初始参数，如位置、质量、速度。这使得系统像

i 和 j = 1,2,3
所以问题是如何处理这个问题，将该系统提供给 scipy odeint ？

更新
我写了代码

from scipy import linspace
from scipy.integrate import odeint

def rf(i, j, r):
    x1, y1, z1 = r[0]
    x2, y2, z2 = r[1]
    x3, y3, z3 = r[2]
    r12 = ((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)**2
    r13 = ((x1-x3)**2+(y1-y3)**2+(z1-z3)**2)**2
    r23 = ((x2-x3)**2+(y2-y3)**2+(z2-z3)**2)**2
    if i == 1:
        if j == 2:
            r = r12
        elif j == 3:
            r = r13
    if i == 2:
        if j == 1:
            r = r12
        elif j == 3:
            r = r23
    if i == 3:
        if j == 1:
            r = r13
        elif j == 2:
            r = r23
    return r


def dotf(r, t0):
    x1, y1, z1 = r[0]
    x2, y2, z2 = r[1]
    x3, y3, z3 = r[2]
    x = [x1, x2, x3]
    y = [y1, y2, y3]
    z = [z1, z2, z3]

    k2 = 6.67e-11
    m = [2e6, 2.2e7, 0.1e3]

    for i in range(1, 3, 1):
        xs, ys, zs = 0, 0, 0
        for j in range(1, 3, 1):
            if i != j:
                r = rf(i, j, r)
                xs = xs + m[j]*(x[i]-x[j])/r
                ys = ys + m[j]*(y[i]-y[j])/r
                zs = zs + m[j]*(z[i]-z[j])/r
        x[i] = -1 * k2 * xs
        y[i] = -1 * k2 * ys
        z[i] = -1 * k2 * zs
    return [[x1,y1,z1], [x2,y2,z2], [x3,y3,z3]]


t = linspace(0, 50, 1e4)
r1 = [1, 2, 1]
r2 = [0.5, 0.1, 3]
r3 = [0.6, 1, 1.5]
r = [r1, r2, r3] 

u = odeint(dotf, r, t)

并得到输出错误：

Traceback (most recent call last):
  File "3b.py", line 73, in <module>
    u = odeint(dotf, r, t)
  File "/usr/local/lib/python2.7/dist-packages/scipy/integrate/odepack.py", line 148, in odeint
    ixpr, mxstep, mxhnil, mxordn, mxords)
ValueError: object too deep for desired array

Answer 1

我已经更正了您代码中的两个明显错误；它运行，但我没有 确定它是正确的。代码在最后。

第一个错误是odeint 想要处理状态向量和 状态导数。向量是指 rank-1 数组；你正在提交 一个矩阵（rank-2 数组）作为您的初始条件。我改变了 r 赋值，可能令人困惑，如 r = r1 + r2 + r3 - + 列表上的运算符是连接，尽管在 numpy 数组上它是 逐元素加法。

此更改意味着我必须将分配更改为 x1 等，在 dotf 和 rf。我还将 dotf 的返回值更改为向量。

第二个错误是你使用r 来表示两个不同的东西： 首先是系统状态向量，其次是来自rf 的返回。我 将第二个更改为rr。

解决方案似乎不稳定；我不知道这是不是 似是而非。我想还有其他错误，但你至少有 现在运行的东西。

我的建议是从你知道的更简单的东西开始 解，例如，稳定的一阶线性系统，或稳定的 欠阻尼二阶系统，或者 Lorenz，如果你想要一些东西 混乱。

有关 Scipy 文档，请参阅 http://docs.scipy.org/doc/，例如，我使用 http://scipy-0.13.0/reference/generated/scipy.integrate.odeint.html#scipy.integrate.odeint 对于 odeint 对于 Scipy 0.13.0（这是我的 Ubuntu 附带的 14.04系统）。

from scipy import linspace
from scipy.integrate import odeint

def rf(i, j, r):
    x1, y1, z1 = r[:3]
    x2, y2, z2 = r[3:6]
    x3, y3, z3 = r[6:]
    r12 = ((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)**2
    r13 = ((x1-x3)**2+(y1-y3)**2+(z1-z3)**2)**2
    r23 = ((x2-x3)**2+(y2-y3)**2+(z2-z3)**2)**2
    if i == 1:
        if j == 2:
            r = r12
        elif j == 3:
            r = r13
    if i == 2:
        if j == 1:
            r = r12
        elif j == 3:
            r = r23
    if i == 3:
        if j == 1:
            r = r13
        elif j == 2:
            r = r23
    return r


def dotf(r, t0):
    x1, y1, z1 = r[:3]
    x2, y2, z2 = r[3:6]
    x3, y3, z3 = r[6:]
    x = [x1, x2, x3]
    y = [y1, y2, y3]
    z = [z1, z2, z3]

    k2 = 6.67e-11
    m = [2e6, 2.2e7, 0.1e3]

    for i in range(1, 3, 1):
        xs, ys, zs = 0, 0, 0
        for j in range(1, 3, 1):
            if i != j:
                rr = rf(i, j, r)
                xs = xs + m[j]*(x[i]-x[j])/rr
                ys = ys + m[j]*(y[i]-y[j])/rr
                zs = zs + m[j]*(z[i]-z[j])/rr
        x[i] = -1 * k2 * xs
        y[i] = -1 * k2 * ys
        z[i] = -1 * k2 * zs
    return [x1,y1,z1, x2,y2,z2, x3,y3,z3]


t = linspace(0, 50, 1e4)
r1 = [1, 2, 1]
r2 = [0.5, 0.1, 3]
r3 = [0.6, 1, 1.5]
r = r1 + r2 + r3

u = odeint(dotf, r, t)


# uncomment to plot
# from matplotlib import pyplot
#
# pyplot.plot(t,u)
# pyplot.show()