开普勒方程迭代求解C语言,一种实用的开普勒方程求解方法及其 C 语言实现

66b52468c121889b900d4956032f1009.png

8种机械键盘轴体对比

本人程序员,要买一个写代码的键盘,请问红轴和茶轴怎么选?

开普勒方程是航天动力学基础方程,是开普勒定律的数学描述。由于开普勒方程属于超越方程,因而无法通过解析法精确求解,这一问题在历史上困扰全世界数学家们达 200 年之久,直至牛顿迭代方法的出现。

本文将介绍一种实用的开普勒方程求解方法,并采用 C 语言实现其算法。该方法出自美国海军天文台 Marc A. Murison 名为 A Practical Method for Solving the Kepler Equation 的文章,由于其算法简单、高效且利于编码实现,因而具有较好的实用性。

1. 开普勒方程描述

开普勒方程可表达为:

它描述了线性时间 下偏近角 和平近角 之间的非线性转换关系。式中, 为轨道偏心率, 为平近角。其中, 表示二体模型中的平运动, 为天体(或卫星)飞越近地点后累积的时长, 为轨道半长轴。

实际计算中,有两种情况需要求解开普勒方程,一是已知真近角 ,求平近角 ;另外则是已知平近角 ,求真近角 。可见偏近角 只是过程量,真近角 和 才是最终想要的结果量。在航天动力学中偏近角 本就是不真实的虚拟量,但它与真近角 存在直观的几何关系,如图 1 所示。当然平近角 也是不真实的虚拟量,它必须通过偏近角和开普勒方程与真近角实现数值上的转换。

082d910756926d626c85e0534f97eb40.png图 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 所示。

12a91137745823d23b4f143882e139da.png图 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;

}

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值