真男人从不回头看数值验证
只有娘们才喜欢用特殊函数
玩计算器的发现大家都玩过计算器吧，不知注意到没有。
输入任意数，然后不断按 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}]//InverseSeriesSeries[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
本文由超级数学建模编辑整理
来自酱紫君（知乎）
https://zhuanlan.zhihu.com/p/36297534
发表于 2018-6-13 01:39
收藏
分享