如何高效计算对称力矩阵（牛顿第三定律）优化N体模拟加速度计算

在N体模拟任务中，需要计算每个天体受到其他所有天体作用的总受力Fi。该问题原生计算复杂度为O(n²)，需要逐对计算天体间的相互作用力。根据牛顿第三定律 f_i,j = -f_j,i，利用该特性可将作用力的计算总量减少一半。

如何在代码中落地该定律，实现加速度计算逻辑的性能优化？

目前已完成不包含该优化的天体总受力计算逻辑，对应C++实现代码如下：

std::vector<glm::vec3> SequentialAccelerationCalculationImpl::calcAccelerations(
        const std::vector<Body> &bodies,
        const float softening_factor
) {
    const float softening_factor_squared = softening_factor * softening_factor;
    const size_t num_bodies = bodies.size();
    std::vector<glm::vec3> accelerations(num_bodies);
    // O(n^2)
    for (size_t i = 0; i < num_bodies; ++i) {
        glm::vec3 received_total_force(0.0);
        for (size_t j = 0; j < num_bodies; ++j) { // TODO this can be reduced to half 
            if (i != j) {
                const glm::vec3 distance_vector = bodies[i].getCurrentPosition() - bodies[j].getCurrentPosition();
                float distance_squared =
                        (distance_vector.x * distance_vector.x) +
                        (distance_vector.y * distance_vector.y) +
                        (distance_vector.z * distance_vector.z);
                // to avoid zero in the following division
                distance_squared += softening_factor_squared;
                received_total_force += ((bodies[j].getMass() / distance_squared) * glm::normalize(distance_vector));
            }
        }
        accelerations[i] = GRAVITATIONAL_CONSTANT * received_total_force;
    }
    return accelerations;
}

对应的力矩阵结构示意图如下：

优化实现方案 #

核心思路是只遍历所有i < j的天体对，每对仅计算一次距离、方向等公共项，再根据牛顿第三定律同时更新两个天体的受力，直接把计算量从n²次两两运算降到n(n-1)/2次，同时省去原代码中i != j的判断逻辑。

优化后与原实现计算结果完全一致的代码如下：

std::vector<glm::vec3> SequentialAccelerationCalculationImpl::calcAccelerations(
        const std::vector<Body> &bodies,
        const float softening_factor
) {
    const float softening_factor_squared = softening_factor * softening_factor;
    const size_t num_bodies = bodies.size();
    // 初始化加速度数组全为0
    std::vector<glm::vec3> accelerations(num_bodies, glm::vec3(0.0f));
    
    // 仅遍历i<j的天体对，总计算量减半
    for (size_t i = 0; i < num_bodies; ++i) {
        for (size_t j = i + 1; j < num_bodies; ++j) {
            // 公共计算项仅算一次
            const glm::vec3 distance_vector = bodies[i].getCurrentPosition() - bodies[j].getCurrentPosition();
            float distance_squared =
                    (distance_vector.x * distance_vector.x) +
                    (distance_vector.y * distance_vector.y) +
                    (distance_vector.z * distance_vector.z);
            distance_squared += softening_factor_squared;
            const glm::vec3 force_direction = glm::normalize(distance_vector);
            
            // 累加i受到j的引力贡献
            accelerations[i] += (bodies[j].getMass() / distance_squared) * force_direction;
            // 反方向累加j受到i的引力贡献，注意质量项取天体i的质量
            accelerations[j] -= (bodies[i].getMass() / distance_squared) * force_direction;
        }
    }

    // 统一乘引力常数，减少循环内重复乘法
    const float G = GRAVITATIONAL_CONSTANT;
    for (auto& acc : accelerations) {
        acc *= G;
    }
    
    return accelerations;
}

优化点说明 #

• 内层循环从j = i+1开始，完全避免同一对天体的重复计算，也不需要额外判断i != j
• 距离向量、距离平方、力方向这些公共计算项每对仅计算一次，比原实现少做一半重复算术运算
• 把引力常数的乘法移到循环外统一执行，进一步减少循环内运算量
• 输出结果和原实现完全等价，无精度损失，串行场景下性能接近原实现的2倍

内容的提问来源于stack exchange，提问作者user13583700