计算机求解开普勒方程

Question

我正在尝试解决开普勒方程，以此作为给定时间找到轨道物体真近点角度的步骤。但事实证明，开普勒方程很难求解，维基百科页面使用微积分描述了这个过程。好吧，我不懂微积分，但我知道解方程涉及无限数量的集合，这些集合会产生越来越接近正确答案的近似值。

通过观察数学公式，我看不出如何在计算机上解决这个问题，所以我希望有更好的数学背景的人能够帮助我。我该如何在计算机上解决这个难题？

顺便说一句，我正在使用F#，我可以计算出该方程式所需的其他元素，只是这一部分让我困扰。

我也愿意接受通过时间、近心点距离和离心率逼近真近点角度的方法。

- Carson Myers

---

有趣，但我不使用F#。如果是Java的话，我可以帮忙，但我还在学习中。 - S L

@experimentX 我添加了一个与语言无关的标签。我不介意从伪代码或其他语言进行转换。 - Carson Myers

好像开普勒方程和我学过的椭圆方程有很大不同。还需要继续努力，可能需要花费很长时间。 - S L

一些相关的问答：解决开普勒方程，C++实现，逼真的n体太阳系模拟。 - Spektre

Answer 1

这篇论文介绍了一种实用的方法来解决开普勒方程，使用迭代计算方法。您可以将其翻译成您选择的语言。论文链接：http://murison.alpheratz.net/dynamics/twobody/KeplerIterations_summary.pdf。

---

你可能也会发现这很有趣。它是一个ocaml程序，其中一部分声称包含一个Kepler方程求解器。由于F#属于ML语言家族（与ocaml类似），这可能是一个很好的起点。

- Robert Harvey

---

你在编辑方面是正确的，我找到了开普勒方程求解器。谢谢。 - Carson Myers

Answer 2

如果有人在寻找类似材料，我希望能够通过这个页面回复。

以下内容是在Adobe After Effects软件中编写的“表达式”，因此它类似于javascript，尽管我有一个Python版本用于另一个应用程序（Cinema 4d）。思路是相同的：迭代执行牛顿法直到达到某个任意精度。

请注意，我发布这段代码不是为了示范或者有效地实现 任何 方式，只是发布我们在期限内完成特定任务（即按照开普勒定律将行星围绕焦点移动，并且准确地执行）。 我们不以编写代码为生，因此我们也不会发布此代码进行批判。 快速并且脏就是遇到期限时的最佳选择。

在After Effects中，任何“表达式”代码都会对动画的每一帧执行一次。这限制了在实现许多算法时可以做什么，因为很难轻松地处理全局数据（用于开普勒运动的其他算法使用迭代更新的速度向量，这是我们无法使用的方法）。代码留下的结果是该时刻对象的[x，y]位置（在内部，这是帧编号），代码旨在附加到时间轴上的对象层的位置元素。

该代码来自于http://www.jgiesen.de/kepler/kepler.html，并在这里提供给下一个人使用。

pi = Math.PI;
function EccAnom(ec,am,dp,_maxiter) { 
// ec=eccentricity, am=mean anomaly,
// dp=number of decimal places
    pi=Math.PI;
    i=0; 
    delta=Math.pow(10,-dp); 
    var E, F; 

    // some attempt to optimize prediction
    if (ec<0.8) {
        E=am;
    } else {
       E= am + Math.sin(am);
    }
    F = E - ec*Math.sin(E) - am; 
    while ((Math.abs(F)>delta) && (i<_maxiter)) {
        E = E - F/(1.0-(ec* Math.cos(E) )); 
        F = E - ec * Math.sin(E) - am; 
        i = i + 1;
    } 
    return Math.round(E*Math.pow(10,dp))/Math.pow(10,dp);
} 
function TrueAnom(ec,E,dp) { 
    S=Math.sin(E); 
    C=Math.cos(E); 
    fak=Math.sqrt(1.0-ec^2);

    phi = 2.0 * Math.atan(Math.sqrt((1.0+ec)/(1.0-ec))*Math.tan(E/2.0));
    return Math.round(phi*Math.pow(10,dp))/Math.pow(10,dp);
} 
function MeanAnom(time,_period) {
    curr_frame  = timeToFrames(time);
    if (curr_frame <= _period) {
        frames_done = curr_frame;
        if (frames_done < 1) frames_done = 1;
    } else {
        frames_done = curr_frame % _period;
    }
    _fractime = (frames_done * 1.0 ) / _period;
    mean_temp   = (2.0*Math.PI) * (-1.0 * _fractime);
    return mean_temp;
}
//==============================
// a=semimajor axis, ec=eccentricity, E=eccentric anomaly 
// delta = delta digits to exit, period = per., in frames
//----------------------------------------------------------
_eccen      = 0.9;              
_delta      = 14;
_maxiter    = 1000;                 
_period     = 300;          
_semi_a     = 70.0;
_semi_b     = _semi_a * Math.sqrt(1.0-_eccen^2); 
_meananom = MeanAnom(time,_period);
_eccentricanomaly = EccAnom(_eccen,_meananom,_delta,_maxiter);
_trueanomaly = TrueAnom(_eccen,_eccentricanomaly,_delta);
r = _semi_a * (1.0 - _eccen^2) / (1.0 + (_eccen*Math.cos(_trueanomaly)));
x = r * Math.cos(_trueanomaly);
y = r * Math.sin(_trueanomaly);
_foc=_semi_a*_eccen;
[1460+x+_foc,540+y];

- user1592212

Answer 3

你可以查看这个由Carl Johansen用C#实现

表示围绕一个质量巨大的中心天体的椭圆轨道上的物体

以下是代码的注释：

在这种情况下，真近点角（True Anomaly）是 物体与太阳之间的角度。 对于椭圆轨道，它有点复杂。 完成周期的百分比仍然是一个关键输入， 但我们还需要应用开普勒方程 （基于离心率）以确保我们在相等时间内扫过相等区域。 这个方程是超越方程（即不能代数解）， 所以我们要么使用一个近似方程， 要么使用数值方法进行求解。 我的实现使用牛顿-拉弗森迭代得到 一个非常好的近似答案（通常在2或3次迭代中）。

- Aivan Monceller

---

太棒了！我还可以再检查一个例子。 - Carson Myers

我必须接受这个答案，这个例子很容易理解并且运行得非常好。 - Carson Myers

循环 while(Math.Abs(E_old - E_new) > epsilon); 在第一次迭代时总是会结束，因为前一行代码说了 E_old = E_new。如果我有误，请纠正我。 - noisy cat