真男人从不回头看数值验证

只有娘们才喜欢用特殊函数

玩计算器的发现

大家都玩过计算器吧，不知注意到没有。

输入任意数，然后不断按cos ANS最后总会输出0.739085。

什么，你说明明记得是：0.999847哦，因为你用了角度制。

这一系列操作等价于求解方程x=cosx，角度制下就是

当然对于现在的你来说求数值解没啥意思了，要求就求解析解是吧。

不过这两个方程其实是一样的，我们先变个形：

也就是说：

于是我们现在只要解决Ax-B=sin(x)这一个方程了。

最早研究这个问题的是天文学家，毕竟那时候也没什么计算器给你玩，一切要从实际出发...

开普勒方程

你可能听说过，三体问题很困难，直到一百多年前的庞加莱时代才被搞定。

而二体问题则简单的多，400年前开普勒时代就研究的差不多了。

你至少知道这个成果，两个天体以一个为焦点，另一个必定在圆锥曲线上运动。

一般天体遵循椭圆轨道，如图椭圆是实际运行的轨道，与椭圆相切的是一个以半长轴为半径的辅助圆。

在一定的时间t内，椭圆轨道上的质点运行到了p点，而辅助圆上的假想质点运行到了y点。

椭圆轨道上所转过的角度∠T被称为真近点角(True Anomaly)

辅助圆轨道上假想质点所转过的角度∠M被称为平近点角(Mean Anomaly)

将椭圆上的质点向上作延长线，交辅助圆于x点所形成的角∠E被称为偏近点角(Eccentric Anomaly)

天文学家发现，它们满足如下关系式：

Kepler Equation：

抛物线就是

的特殊情况，双曲线有所不同。

Hyperbolic Kepler Equation：

但从数学上讲，这个式子其实就是：

也就是说不考虑物理意义其实是一样的。

开普勒方程的解析解

有了方程当然接下来就是求解了咯，古代计算力比较值钱，毕竟没有计算机，所以大家对解析解都有一种病态的追求。

怎么着推一天公式要比算一整天的牛顿迭代有趣吧?

作一下等价性检验：

In [] = FindRoot[x==Cos@x,{x,0}] x-Pi/2/.FindRoot[Pi/2==x-Sin@x,{x,1}] FindRoot[x==Cos[Pi x/180],{x,0}] 180x/Pi-90/.FindRoot[Pi/2==x-Pi Sin@x/180,{x,1}] Out[] = 0.7390851332151605` {x -> 0.7390851332151607`} 0.9998477415310987` {x -> 0.9998477415310881`}

拉格朗日反演

E不能分离但M，展开M(E)，然后直接用级数反演即可。

Mathematica 可以很方便的执行级数反演。

Series[M- Sin[M], {M, 0, 10}]//InverseSeries Series[M-e Sin[M], {M, 0, 10}]//InverseSeries

早期解这个方程使用了关于离心率

的麦克劳林展开。

这不是个整函数，所以引入了所谓的拉普拉斯极限。

超出收敛域的部分级数失效，级数反演则很好的解决了这个问题。

贝塞尔函数解

当然无穷级数不利于计算，能否使用微积分表达是我们接下来的探索重点。

我们来考虑函数方程：

我们假设它可以展开为傅里叶级数，分析原函数方程性态可以期望这是个正弦级数。

那么系数可以表达为：

我们来尝试计算，嗯？没思路怎么办...

无脑分部积分展开到能搞定为止呗。

而这正好是贝塞尔函数的定义式之一：

Bessel Function of the First Kind：

于是原式可以写成

赫维茨-勒奇超越函数解

Stack Exchange上有个用反三角函数和三角函数表示的解析解，这个解比较有难度。

特殊函数论中将以下级数称为赫维茨-勒奇超越函数(Lerch Transcendent Function)

我们从上面的贝塞尔函数解开始，还原掉贝塞尔函数：

然后交换积分求和顺序。

里面的部分圈起来叫F(M)，用欧拉公式展开。

其中：

可以发现其实都是

的结构。

我们引入多对数函数：

也就是说：

用这个函数化简等式：

同样的整理一下：

可以合并成两组，然后再次展开，运算量有点大。

化简的时候注意恒等式：

注意到第二部分：

最后代回去大功告成！

代入数据就得到了 Stack Exchange 一样的结果。

我对arctan(tan(x))这种写法感到很不爽。

这个当然不能直接抵消，由于arctan(tan(x))≠x，我们作复展开。

严格来说这两者不是完全相等的，因为这样一来消掉了奇点。

不过积分的时候完全可以划等号，因为区间开闭完全不影响积分值。

综上所述，最后代入值，我们得到了：

(*真男人从不回头看数值验证*) (2 + I Integrate[Log[-I/E^(I*(t - Sin[t]))], {t, 0, Pi}])/(2*Pi)//N (Pi + 90I Integrate[Log[(-I)*E^((-I)*t + (1/180)*I*Pi*Sin[t])], {t, 0, Pi}])/Pi^2//N > 0.7390851332151609` > 0.9998477415310951`

最后一个是百度贴吧上的，这个答案就非常魔幻了，它和上面两个方法不是一个系列的，和第一个方法有关。

暴力求解拉格朗日反演的解析形式，场面非常的少儿不宜...

我一时半会儿也没看懂，详情看参考书目(3)。

从这个结果上也能看出这个方法有多残暴...

(*怎么可以这么暴力的说*) \[Pi]/2 Exp[NIntegrate[1/(\[Pi] x) ArcTan[((\[Pi] x+2)Log[(Sqrt[1-x^2]+1)/x]x)/(x^2Log[(Sqrt[1-x^2]+1)/x]^2-\[Pi] x-1)],{x,0,1},WorkingPrecision->50]] ArcCot[1+1/(2\[Pi] ) NIntegrate[Log[((1-x^2)Pi^2+4(Sqrt[1-x^2]ArcTanh[x]+x)^2)/((1-x^2)Pi^2+4(Sqrt[1-x^2]ArcTanh[x]-x)^2)],{x,0,1},WorkingPrecision->50]] > 0.73908513321516064165531208767387340401341175890075746496567242428255184768807`50.12267193056545 > 0.73908513321516064165531208767387340401341175890075746496567993239614795659229`51.22422170141253

参考书目

1.On Taylor series and Kapteyn series of the first and second type

2.Kepler's equation, radiation problems and Meissel's expansion

3.An exact analytical solution of Kepler's Equation