Runge-Kutta N-Body 实施无效

Question

我已经在 Java 中实现了本文所述的 Runge-Kutta 4 算法。 http://spiff.rit.edu/richmond/nbody/OrbitRungeKutta4.pdf

但是，动作不稳定。即使只有两个身体，也没有周期性运动。

protected void integrateRK4(double time) {

    final double H = time;
    final double HO2 = H/2;
    final double HO6 = H/6;

    Vector2[] currentVelocities = new Vector2[objects.length];
    Vector2[] currentPositions = new Vector2[objects.length];
    Vector2[] vk1 = new Vector2[objects.length];
    Vector2[] vk2 = new Vector2[objects.length];
    Vector2[] vk3 = new Vector2[objects.length];
    Vector2[] vk4 = new Vector2[objects.length];
    Vector2[] rk1 = new Vector2[objects.length];
    Vector2[] rk2 = new Vector2[objects.length];
    Vector2[] rk3 = new Vector2[objects.length];
    Vector2[] rk4 = new Vector2[objects.length];


    for (int i=0; i<objects.length; i++) {
        currentVelocities[i] = objects[i].velocity().clone();
        currentPositions[i] = objects[i].position().clone();
    }

        vk1 = computeAccelerations(objects);
    for (int i=0; i<objects.length; i++) {
        rk1[i] = currentVelocities[i].clone();
    }

    for (int i=0; i<objects.length; i++) {
        objects[i].setPosition(Vectors2.add(currentPositions[i], Vectors2.prod(rk1[i], HO2)));
    }
        vk2 = computeAccelerations(objects);

    for (int i=0; i<objects.length; i++) {
        rk2[i] = Vectors2.prod(vk1[i], HO2);
    }


    for (int i=0; i<objects.length; i++) {
        objects[i].setPosition(Vectors2.add(currentPositions[i], Vectors2.prod(rk2[i], HO2)));
    }
        vk3 = computeAccelerations(objects);

    for (int i=0; i<objects.length; i++) {
        rk3[i] = Vectors2.prod(vk2[i], HO2);
    }


    for (int i=0; i<objects.length; i++) {
        objects[i].setPosition(Vectors2.add(currentPositions[i], Vectors2.prod(rk3[i], H)));
    }
        vk4 = computeAccelerations(objects);

    for (int i=0; i<objects.length; i++) {
        rk4[i] = Vectors2.prod(vk3[i], H);
    }


    for (int i=0; i<objects.length; i++) {
        objects[i].setVelocity(Vectors2.add(currentVelocities[i], Vectors2.prod(Vectors2.add(vk1[i], Vectors2.prod(vk2[i], 2), Vectors2.prod(vk3[i], 2), vk4[i]), HO6)));
        objects[i].setPosition(Vectors2.add(currentPositions[i], Vectors2.prod(Vectors2.add(rk1[i], Vectors2.prod(rk2[i], 2), Vectors2.prod(rk3[i], 2), rk4[i]), HO6)));
    }
}

我的实现不知何故不正确吗？

注：

Vectors2是我自己实现的Vectors，它是大小为2的一阶张量。
所有的静态方法Vectors2.* return vector的副本。
在 Vector2 实例上调用的所有非静态方法都会修改向量，对于 objects[i].addVelocity() 和 objects[i].addPosition()

也是如此

Vectors2.add(Vector2... vectors2) 做 element-wise 加法。
Vectors2.prod(Vector2... vectors2) 做 element-wise 乘法。
Vectors2.prod(Vector2 vector2, double c) 做 element-wise 乘法。

computeAccelerations(PointBody[] objects) return 一个加速度向量数组。
变量 objects 是 class NBodyIntegrator 的属性，这些函数是其中的一部分。它是一个 PointBody[] 的数组。

为了确保这不是其他错误，我通过删除 k2、k3 和 k4 将 RK4 算法简化为显式欧拉算法。这个按预期工作。

protected void integrateRK1(double time) {
    final double H = time;
    final double HO2 = H/2;

    Vector2[] currentVelocities = new Vector2[objects.length];
    Vector2[] currentPositions = new Vector2[objects.length];
    Vector2[] vk1;
    Vector2[] rk1 = new Vector2[objects.length];


    for (int i=0; i<objects.length; i++) {
        currentVelocities[i] = objects[i].velocity().clone();
        currentPositions[i] = objects[i].position().clone();
    }


        vk1 = computeAccelerations(objects);
    for (int i=0; i<objects.length; i++) {
        rk1[i] = currentVelocities[i].clone();
    }


    for (int i=0; i<objects.length; i++) {
        objects[i].setVelocity(Vectors2.add(currentVelocities[i], Vectors2.prod(Vectors2.add(vk1[i]), H)));
        objects[i].setPosition(Vectors2.add(currentPositions[i], Vectors2.prod(Vectors2.add(rk1[i]), H)));
    }
}

Answer 1

您正在设置

rk1 = v0
pos2 = pos0 + rk1*h/2
vk2 = acc(pos2)

这是正确的。但是然后你继续

rk2 = vk1*h/2

它应该在哪里

rk2 = v0 + vk1*h/2

等等。当然那么累积位置更新也是错误的。

Answer 2

这不是一个很好的 RK 集成实现。

您正在使用共享的可变数据。那不是线程安全的。

我见过的任何正确的数值积分实现都将函数抽象为要集成到传递到积分方案的方法中。您应该能够通过将函数传递给新例程来更改集成方案。

从一个以函数和初始条件作为参数的集成接口开始。

Runge-Kutta N-Body 实施无效

Runge-Kutta N-Body Implementation not working

java

algorithm

math

runge-kutta