无存储Buffer的ShaderToy中如何用时间函数计算行星引力运动位置

核心结论 #

二体引力场景完全可以构造以时间为唯一输入的显式函数，直接计算任意时刻行星位置，不需要任何持久化存储逐帧递推状态。

递推依赖的本质说明 #

你提到的"位置决定加速度、加速度改变速度、速度更新位置"的循环关系，本质是牛顿运动定律对应的常微分方程初值问题，不是不可解的逻辑死锁：

• 逐帧递推是对连续物理规律的离散数值近似，只有用欧拉法、Verlet法这类离散数值方法时，才需要存储上一帧的状态。
• 这类初值问题的解本身就是时间t的一元函数：给定t=0时刻的初始位置、初始速度，任意t时刻的位置是唯一确定的，和计算路径无关，完全不需要依赖历史值。

非圆轨道二体场景的显式实现 #

你提到的非圆形、非稳定轨道场景，属于经典二体问题，存在成熟的解析解，轨道形态为圆锥曲线（椭圆、抛物线、双曲线），可以通过开普勒方程直接计算位置，不需要递推：

1. 基于初始位置、初始速度、中心天体引力参数GM，预计算轨道半长轴、偏心率、轨道平面朝向、初始平近点角等固定轨道参数，这些参数不随时间变化。
2. 对任意输入时间t，先计算当前时刻的平近点角，再通过3-5次牛顿迭代求解开普勒方程得到偏近点角，迭代精度足够满足渲染需求，全程不需要历史状态。
3. 把偏近点角转换为轨道平面坐标，再通过预计算的旋转参数转换到世界坐标系，直接得到t时刻的行星位置。

核心计算逻辑的无状态GLSL函数示例：

#define PI 3.1415926535
// 无状态纯函数：输入时间t，返回对应时刻行星位置
// 所有轨道常量可在shader初始化时预计算，不需要逐帧存储
vec3 calcPlanetPos(float t, vec3 r0, vec3 v0, float GM) {
    // 预计算固定轨道参数
    float r0Len = length(r0);
    float v0Sq = dot(v0, v0);
    float semiMajor = 1.0 / (2.0 / r0Len - v0Sq / GM);
    vec3 angMom = cross(r0, v0);
    float ecc = length(cross(v0, angMom)/GM - normalize(r0));
    float meanMotion = sqrt(GM / (semiMajor*semiMajor*semiMajor));
    
    // 计算初始平近点角，修正象限
    float M0 = acos(clamp((1.0 - r0Len/semiMajor)/ecc, -1.0, 1.0));
    if (dot(r0, v0) < 0.0) M0 = 2.0*PI - M0;

    // 当前时刻平近点角
    float M = mod(M0 + meanMotion * t, 2.0*PI);

    // 牛顿迭代解开普勒方程，4次迭代即可达到渲染足够精度
    float E = M;
    for (int i = 0; i < 4; i++) {
        E -= (E - ecc*sin(E) - M) / (1.0 - ecc*cos(E));
    }

    // 轨道平面坐标转换到世界坐标（省略轨道倾角旋转部分，平面场景可直接使用）
    float orbitX = semiMajor * (cos(E) - ecc);
    float orbitY = semiMajor * sqrt(1.0 - ecc*ecc) * sin(E);
    return vec3(orbitX, orbitY, 0.0);
}

多体场景的限制 #

如果是3个及以上天体存在互相引力作用的多体场景，不存在通用的初等函数闭式解析解，这时候无法写出简单的显式时间函数：

• 若坚持无buffer实现，只能对每个像素从t=0开始做数值积分到当前时刻，时间跨度稍大就会出现计算量爆炸的问题，仅适合短时间演示场景。
• 也可以采用级数展开、近似轨道拟合的方案，但会存在不同程度的精度损失。

通用问题解答 #

对于这类"取值依赖自身历史"的递推关系，只要对应的微分方程满足利普希茨连续条件，解就唯一存在，不存在本质上无法转为显式函数的障碍。只是这类显式解不一定是初等函数的简单组合，可能需要包含迭代、级数展开等运算，不一定符合"简单闭式表达式"的预期。

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