玩计算器的发现
大家都玩过计算器吧, 不知注意到没有.
输入任意数, 然后不断按最后总会输出.
什么, 你说明明记得是:? 哦, 因为你用了角度制.
这一系列操作等价于求解方程, 角度制下就是.
当然对于现在的你来说求数值解没啥意思了, 要求就求解析解是吧.
不过这两个方程其实是一样的, 我们先变个形:
也就是说:
于是我们现在只要解决这一个方程了.
最早研究这个问题的是天文学家, 毕竟那时候也没什么计算器给你玩, 一切要从实际出发...
开普勒方程
你可能听说过, 三体问题很困难, 直到一百多年前的庞加莱时代才被搞定.
而二体问题则简单的多, 400年前开普勒时代就研究的差不多了.
你至少知道这个成果, 两个天体以一个为交点, 另一个必定在圆锥曲线上运动.
一般天体遵循椭圆轨道, 如图椭圆是实际运行的轨道, 与椭圆相切的是一个以半长轴为半径的辅助圆.
在一定的时间内, 椭圆轨道上的质点运行到了点, 而辅助圆上的假想质点运行到了点.
- 椭圆轨道上所转过的角度被称为真近点角(True Anomaly)
- 辅助圆轨道上假想质点所转过的角度被称为平近点角(Mean Anomaly)
- 将椭圆上的质点向上作延长线,交辅助圆于点所形成的角被称为偏近点角(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`}
拉格朗日反演
不能分离但, 展开,然后直接用级数反演即可.
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)
我们从上面的贝塞尔函数解开始, 还原掉贝塞尔函数:
然后交换积分求和顺序.
里面的部分圈起来叫, 用欧拉公式展开.
其中:
可以发现其实都是的结构.
我们引入多对数函数:
也就是说:
用这个函数化简等式:
同样的整理一下:
可以合并成两组, 然后再次展开, 运算量有点大.
化简的时候注意恒等式:.
注意到第二部分:
最后代回去大功告成!
代入数据就得到了 Stack Exchange 一样的结果.
我对这种写法感到很不爽.
这个当然不能直接抵消, 由于, 我们作复展开.
严格来说这两者不是完全相等的, 因为这样一来消掉了奇点.
不过积分的时候完全可以划等号, 因为区间开闭完全不影响积分值.
综上所述, 最后代入值, 我们得到了:
(*真男人从不回头看数值验证*) (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
参考书目
- On Taylor series and Kapteyn series of the first and second type
- Kepler's equation, radiation problems and Meissel's expansion
- An exact analytical solution of Kepler's Equation
题图: https://www.pixiv.net/member_illust.php?mode=medium&illust_id=65616708
来源:知乎 www.zhihu.com
作者:酱紫君
【知乎日报】千万用户的选择,做朋友圈里的新鲜事分享大牛。 点击下载
没有评论:
发表评论