8种机械键盘轴体对比
本人程序员，要买一个写代码的键盘，请问红轴和茶轴怎么选？
开普勒方程是航天动力学基础方程，是开普勒定律的数学描述。由于开普勒方程属于超越方程，因而无法通过解析法精确求解，这一问题在历史上困扰全世界数学家们达 200 年之久，直至牛顿迭代方法的出现。
本文将介绍一种实用的开普勒方程求解方法，并采用 C 语言实现其算法。该方法出自美国海军天文台 Marc A. Murison 名为 A Practical Method for Solving the Kepler Equation 的文章，由于其算法简单、高效且利于编码实现，因而具有较好的实用性。
1. 开普勒方程描述
开普勒方程可表达为：
它描述了线性时间 下偏近角 和平近角 之间的非线性转换关系。式中， 为轨道偏心率， 为平近角。其中， 表示二体模型中的平运动， 为天体(或卫星)飞越近地点后累积的时长， 为轨道半长轴。
实际计算中，有两种情况需要求解开普勒方程，一是已知真近角 ，求平近角 ；另外则是已知平近角 ，求真近角 。可见偏近角 只是过程量，真近角 和 才是最终想要的结果量。在航天动力学中偏近角 本就是不真实的虚拟量，但它与真近角 存在直观的几何关系，如图 1 所示。当然平近角 也是不真实的虚拟量，它必须通过偏近角和开普勒方程与真近角实现数值上的转换。
图 1: 真近角和偏近角之间的几何关系
由图 1 内部椭圆(真实运动轨迹)几何关系可得：
由图 1 外部圆(假想的平运动)几何关系可得：
于是可得真近角和偏近角之间的相互转换关系为：
开普勒方程属于超越方程，故需采用数值法才能求得精确值。理想情况下，一个数值方法是否实用，应该看它是否同时具备以下两点：计算 的 CPU 耗时是否最小化;
程序的复杂度是否最小化。
这两个需求往往是互相制约的，但幸运的是通过分析研究可以找到一个折中的方法来解决这个矛盾。开普勒方程只能通过迭代法求解，任何一个迭代过程的设计都有两个任务：第一个任务是设计迭代循环中的逼近算法。这个逼近算法需要重复执行直到结果达到令人满意的精度，一般说来迭代公式所用的阶数越高，迭代所需的次数就越少。然而，更高阶的迭代公式将使得公式的表达式变得十分复杂，这将极大地耗费 CPU 的运行时间。因此，无论我们选择何种算法，一个恰当(通常是比较低)的阶数会使得 CPU 的耗时最短；
第二个任务是选择一个迭代循环的初始值。初始值选的越精确，循环收敛的越快。选择初始值的方法不需要和迭代方法一样，哪怕两者极其相近。类似于迭代逼近算法，对于一个初值确定的算法，将有一个理想的阶数能够最大限度地减少 CPU 的耗时。
接下来将给出一个非常简单的初值确定方法，紧接着是一个快速迭代算法。最后，直接给出算法性能分析结果及其伪代码，并用 C 语言加以实现之。
2. 初值确定方法
由于必须用迭代法求解开普勒方程，因而循环迭代的初始值越精确，迭代效果就越好，至少在迭代表达式复杂得令人反感之前应该是如此。将开普勒方程写成：
在 的极限情况下，有 ，这是最简单的初始值近似。因此上面的方程可改写成如下简单的迭代近似公式：
其中初始条件 。我们可以对上述递归表达式进行反复迭代，直至得到足够高的偏心率阶数。例如，三阶近似公式为：
程序化后的三阶伪代码如下：1
2
3
4
5
6
7KeplerStart3 := proc(e,M)
local t33, t35, t34;
t34 := eˆ2;
t35 := e*t34;
t33 := cos(M);
return M+(-1/2*t35+e+(t34+3/2*t33*t35)*t33)*sin(M);
end proc;
3. 循环迭代方法
循环迭代的最终目的是要得到收敛的结果。具体到开普勒方程，当给定一个带有误差的 时，必须找到一个迭代公式，使其每次返回一个误差更小的近似值，同时它也必须收敛。基于此，按如下方式改写开普勒方程：
式中， 的解是 。令 为 逼近 时存在的误差，将 在 处进行泰勒展开，于是得到：
式中，假设 充分小。若只考虑一阶展开，可得：
上式可作为一阶迭代方案的的核心算法。假设一开始猜测 ，那么 比 更接近于 。于是得到一阶迭代方程：
其中初始值 的取值将在后面讨论。由上式可以得到单步一阶迭代法来估计 ，对应的伪代码如下：1
2
3eps1 := proc(e,M,x)
return (x-e*sin(x)-M)/(1-e*cos(x));
end proc;
现在把二阶项考虑进去，方程展开后写成 Horner 形式：
同理，令 ，整理得到如下形式：
类比单步一阶形式，创建单步二阶迭代方程：
重点来了，此处可以将 用一阶迭代公式替换，创建一个两步迭代过程，于是得到新的迭代方程：
上述方程经优化后的伪代码为：1
2
3
4
5
6
7eps2 := proc(e,M,x)
localt1,t2,t3;
t1 := -1+e*cos(x);
t2 := e*sin(x);
t3 := -x+t2+M;
return t3/(1/2*t3*t2/t1+t1);
end proc;
更进一步， 三阶展开方程的 Horner 形式为：
同理，令 ，整理得到如下形式：
重点又来了，此处可以先将 用二阶迭代公式替换，再将 用一阶迭代公式替换，创建一个三步迭代过程，于是得到新的迭代方程，由于方程过于复杂，直接给出优化后的伪代码：1
2
3
4
5
6
7
8
9
10eps3 := proc(e,M,x)
local t1, t2, t3, t4, t5, t6;
t1 := cos(x);
t2 := -1+e*t1;
t3 := sin(x);
t4 := e*t3;
t5 := -x+t4+M;
t6 := t5/(1/2*t5*t4/t2+t2);
return t5/((1/2*t3 - 1/6*t1*t6)*e*t6+t2);
end proc;
可以继续利用这种方式到更高阶的形式，本文的推导止于此，后文我们会看到三阶迭代已经处于 CPU 时耗最优区间。
4. 实用的开普勒方程求解方法
数值测试，令算法执行的 CPU 耗时为迭代阶数 (Niter) 和 初始值逼近阶数 (Nstart) 的二元函数，绘制求解 16000 个开普勒方程的 CPU 耗时等高线，其中 和 在 等间隔 网格域中选取。结果表明，不论是对于初值算法，还是迭代算法，三阶形式都接近最优，如图 2 所示。
图 2: 真近角和偏近角之间的几何关系
由此，给出优化后的三阶迭代和初始值方法求解开普勒方程的伪代码：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19KeplerSolve := proc( e, M, tol:=1.0e-14 )
local dE, E, E0, Mnorm, count;
global Estart3, eps3;
Mnorm := fmod(M,2*Pi);
E0 := KeplerStart3(e,Mnorm);
dE := tol + 1;
count := 0;
while dE > tol do
E := E0 - eps3(e,Mnorm,E0);
dE := abs(E-E0);
E0 := E;
count := count + 1;
if count=100 then
print “太令人惊讶了，KeplerSolve 竟然不收敛！”;
break;
end if;
end do;
return E;
end proc;
5. C 语言实现
根据优化后的实用开普勒方程求解伪代码，可编写其计算机语言实现代码。下面给出三个主要函数的 C 语言实现，这些代码已经过充分测试，可放心使用。函数中的 fmod 函数为取模函数，用于处理平近角周期问题，请读者自行补脑。至此，码完收工，拜了个拜！1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115**
** 函数名称: KeplerStart3
** 函数功能: 开普勒方程三阶初值估计
** 输入参数: ecc -> 偏心率 [N/A]
** M -> 平近角 [rad]
** 输出参数: 无
** 返回值 : E -> 偏近角 [rad]
**
** ---------------------------------------------------------------------*/
double (double ecc, double M){
double t34, t35, t33, E;
if(ecc < 0 || ecc >= 1){
printf("!!!错误的偏心率，该程序只适合椭圆轨道(0 =< ecc < 1)!!!n");
return -1;
}
if(M >= D2PI || M < 0){
M = fmod(M, 0, D2PI);
}
t34 = ecc * ecc;
t35 = ecc * t34;
t33 = cos(M);
E = M + (-0.5*t35 + ecc + (t34 + 3.0/2.0*t33*t35) * t33) * sin(M);
return E;
}
**
** 函数名称: eps3
** 函数功能: 开普勒方程三阶误差估计
** 输入参数: ecc -> 偏心率 [N/A]
** M -> 平近角 [rad]
** x -> 偏近角过程值 [rad]
** 输出参数: 无
** 返回值 : E_error -> 偏近角估计误差 [rad]
**
** ---------------------------------------------------------------------*/
double eps3(double ecc, double M, double x){
double t1, t2, t3, t4, t5, t6, E_error;
if(ecc < 0 || ecc >= 1){
printf("!!!错误的偏心率，该程序只适合椭圆轨道(0 =< ecc < 1)!!!n");
return -1;
}
if(M >= D2PI || M < 0){
M = fmod(M, 0, D2PI);
}
t1 = cos(x);
t2 = -1 + ecc * t1;
t3 = sin(x);
t4 = ecc * t3;
t5 = -x + t4 + M;
t6 = t5 / (0.5 * t5 * t4 / t2 + t2);
E_error = t5 / ((0.5*t3 - 1.0/6.0*t1*t6) * ecc * t6 + t2);
return E_error;
}
**
** 函数名称: KeplerSolve
** 函数功能: 求解开普勒方程
** 输入参数: ecc -> 偏心率 [N/A]
** M -> 平近角 [rad]
** 输出参数: 无
** 返回值 : E -> 偏近角 [rad]
**
** ---------------------------------------------------------------------*/
double KeplerSolve(double ecc, double M){
double E0, dE, E, Mnorm;
double tol = 1.0e-20;
int iter_count = 0;
if(ecc < 0 || ecc >= 1){
printf("!!!错误的偏心率，该程序只适合椭圆轨道(0 =< ecc < 1)!!!n");
return -1;
}
if(M >= D2PI || M < 0){
M = fmod(M, 0, D2PI);
}
Mnorm = fmod(M, 0, D2PI);
E0 = KeplerStart3(ecc, Mnorm);
dE = tol + 1;
while(dE > tol){
E = E0 - eps3(ecc, Mnorm, E0);
dE = fabs(E - E0);
E0 = E;
iter_count = iter_count + 1;
if(iter_count > 10){
printf("太令人惊讶了，KeplerSolve 竟然不收敛！n");
return -1;
}
// printf("迭代次数：iter_count = %dn", iter_count);
}
return E;
}
扫一扫
2321
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
