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;
}