x = cos x 的解析形式

05-08

玩計算器的發現

大家都玩過計算器吧, 不知注意到沒有.

輸入任意數, 然後不斷按 $mathtt{cos ANS}$ 最後總會輸出 $0.739085$ .

什麼, 你說明明記得是: $0.999847$ ? 哦, 因為你用了角度制.

這一系列操作等價於求解方程 $x=cos x$ , 角度制下就是 $x=cos x°=cosdfrac{pi x}{180}$ .

當然對於現在的你來說求數值解沒啥意思了, 要求就求解析解是吧.

不過這兩個方程其實是一樣的, 我們先變個形:

$egin{aligned} x=cos x;Longrightarrow& x-frac{pi}{2}=cos left(x-frac{pi}{2} ight)&=sin x\ x=cos x°Longrightarrow& frac{180 x}{pi }-90=cos left(frac{ pi}{180} left(frac{180 x}{pi }-90 ight) ight)&=sin x end{aligned}$

也就是說:

於是我們現在只要解決 $Ax-B=sin(x)$ 這一個方程了.

最早研究這個問題的是天文學家, 畢竟那時候也沒什麼計算器給你玩, 一切要從實際出發...

開普勒方程

你可能聽說過, 三體問題很困難, 直到一百多年前的龐加萊時代才被搞定.

而二體問題則簡單的多, 400年前開普勒時代就研究的差不多了.

你至少知道這個成果, 兩個天體以一個為交點, 另一個必定在圓錐曲線上運動.

一般天體遵循橢圓軌道, 如圖橢圓是實際運行的軌道, 與橢圓相切的是一個以半長軸為半徑的輔助圓.

在一定的時間 $t$ 內, 橢圓軌道上的質點運行到了 $p$ 點, 而輔助圓上的假想質點運行到了 $y$ 點.

橢圓軌道上所轉過的角度 $angle T$ 被稱為真近點角(True Anomaly)
輔助圓軌道上假想質點所轉過的角度 $angle M$ 被稱為平近點角(Mean Anomaly)
將橢圓上的質點向上作延長線,交輔助圓於 $x$ 點所形成的角 $angle E$ 被稱為偏近點角(Eccentric Anomaly)

天文學家發現, 它們滿足如下關係式:

Kepler Equation: $M= E-epsilon sin(E)$

拋物線就是 $epsilon=1$ 的特殊情況, 雙曲線有所不同.

Hyperbolic Kepler Equation: $M = epsilon sinh H - H quadmathrm{where}quad H=iE$

但從數學上講, 這個式子其實就是:

$M = i left( E - epsilon sin E ight)$

也就是說不考慮物理意義其實是一樣的.

開普勒方程的解析解

有了方程當然接下來就是求解了咯, 古代計算力比較值錢, 畢竟沒有計算機, 所以大家對解析解都有一種病態的追求.

怎麼著推一天公式要比算一整天的牛頓迭代有趣吧?

$left{egin{aligned} frac{pi}{2}&=x-sin x\ frac{pi}{2}&=x-frac{pi}{180}sin x end{aligned} ight.$

作一下等價性檢驗:

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)$ ,然後直接用級數反演即可.

$M(E) = (1-epsilon)E+epsilonsum _{n=1}^{infty } frac{(-1)^n}{(2 n+1)!}E^{2 n+1}$

$igstar E(M)= egin{cases} displaystyle sum _{n=1}^{infty }igg(lim_{ heta o 0^{+}}!{frac {mathrm {d} ^{,n-1}}{mathrm {d} heta ^{,n-1}}}{igg (}{frac { heta }{sqrt[{3}]{ heta -sin( heta )}}}{igg )}^{!!!n}{igg )}{frac {M^{frac{n}{3}}}{n!}}&epsilon=1\ displaystyle sum_{n=1}^{infty }igg(lim_{ heta o 0^{+}}{frac {mathrm {d} ^{,n-1}}{mathrm {d} heta ^{,n-1}}}{Big (}{frac { heta }{ heta -epsilonsin( heta )}}{Big )}^{!n}{igg )}{frac {M^{n}}{n!}}&epsilon eq 1 end{cases}$

Mathematica 可以很方便的執行級數反演.

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

早期解這個方程使用了關於離心率 $epsilon$ 的麥克勞林展開.

$E(M)=M+sum_{n=1}^infty a_n epsilon^n; epsilonleqmathrm{L}$

這不是個整函數, 所以引入了所謂的拉普拉斯極限.

$L=max_{xinmathbb{R}}frac{x}{cosh (x)}approx0.662743$

超出收斂域的部分級數失效, 級數反演則很好的解決了這個問題.

貝塞爾函數解

當然無窮級數不利於計算, 能否使用微積分表達是我們接下來的探索重點.

我們來考慮函數方程: $g (M) = E (M) - M$

我們假設它可以展開為傅里葉級數, 分析原函數方程性態可以期望這是個正弦級數.

$g (M) = sum_{n = 1}^{infty}a_nsin (n M)$

那麼係數可以表達為:

$a_n = frac {2}{pi}int_0^pi g(M) sin(nM),mathrm{d}M$

我們來嘗試計算, 嗯? 沒思路怎麼辦...

無腦分部積分展開到能搞定為止唄.

$egin{aligned} a_n&=frac{2}{pi n}int_0^pi cos(nM),mathrm{d}g(M) -frac{2}{pi}left[g(M)frac{cos(nM)}{n} ight]_0^pi\ &=frac{2}{pi n}int_0^pi cos(nM),mathrm{d}(E-M)\ &=frac{2}{pi n}int_0^pi cos[n(E-epsilonsin E)],mathrm{d}E -frac{2}{pi n}int_0^pi cos(nM),mathrm{d}M\ &=frac{2}{pi n}int_0^pi cos(nE-nepsilonsin E),mathrm{d}E end{aligned}$

而這正好是貝塞爾函數的定義式之一:

Bessel Function of the First Kind: $J_n(z)=frac{1}{pi}int_0^pi cos(n θ-z sin θ),mathrm{d}θ ;nin mathbb{Z}$

於是原式可以寫成

$igstar E(M)=M+sum _{n=1}^{infty } frac{2}{n}J_n(n epsilon )sin(nM)$

赫維茨-勒奇超越函數解

Stack Exchange上有個用反三角函數和三角函數表示的解析解, 這個解比較有難度.

$mathcal{D}=frac1pi int_0^{pi } arctanleft( an left(frac{t-sin t+frac{pi }{2}}2 ight) ight) ,mathrm{d}t+frac{1}{pi }$

特殊函數論中將以下級數稱為赫維茨-勒奇超越函數(Lerch Transcendent Function)

$Phi (z,t,h):=sum _{n=0}^{infty } frac{z^n}{(h+n)^t}$

我們從上面的貝塞爾函數解開始, 還原掉貝塞爾函數:

$E=M+sum _{n=1}^{infty } frac{2}{n}left[frac{1}{pi}int_0^pi cos(n θ-nepsilon sin θ),mathrm{d}θ ight]sin(nM)$

然後交換積分求和順序.

$E=M+frac{2}{pi}int_0^pileft[sum _{n=1}^{infty } frac{sin(nM)}{n} cos(n θ-nepsilon sin θ) ight],mathrm{d}θ$

裡面的部分圈起來叫 $F(M)$ , 用歐拉公式展開. $egin{aligned} F(M)=&frac{sin(nM)}{n} cos(n θ-nepsilon sin θ)\ =&frac{i}{4 n}left(e^{-i M n}-e^{i M n} ight) left(e^{-i ( heta n-n epsilon sin ( heta ))}+e^{i ( heta n-n epsilon sin ( heta ))} ight)\ =&frac{i}{4 n}left(e_1+e_2+e_3+e_4 ight) end{aligned}$

其中:

$egin{cases} e_1=+exp(-i M n+i heta n-i n epsilon sin ( heta ))\ e_2=-exp(i M n+i heta n-i n epsilon sin ( heta ))\ e_3=+exp(-i M n-i heta n+i n epsilon sin ( heta ))\ e_4=-exp(i M n-i heta n+i n epsilon sin ( heta ))\ end{cases}$

可以發現其實都是 $e^{nalpha}$ 的結構.

我們引入多對數函數:

$mathrm{Li}_s(z) := zPhi (z,s,1)=sum _{n=1}^{infty } frac{z^n}{n^s}$

也就是說:

$sum _{n=1}^{infty } frac{e^{a n}}{n}= ext{Li}_1left(e^a ight)=iarg (1-e^a)-ln |1-e^a|$

用這個函數化簡等式:

$egin{aligned} E=&M+frac{2}{pi}int_0^pileft[sum _{n=1}^{infty } frac{i}{4 n}left(e_1+e_2+e_3+e_4 ight) ight],mathrm{d}θ\ =&M+frac{i}{2pi}int_0^pileft[ ext{Li}_1(e^{a_1})+ ext{Li}_1(e^{a_2})+ ext{Li}_1(e^{a_3})+ ext{Li}_1(e^{a_4}) ight],mathrm{d}θ end{aligned}$

同樣的整理一下:

$egin{cases} a_1=+i ( heta -M-epsilon sin ( heta ))\ a_2=+i ( heta +M-epsilon sin ( heta ))\ a_3=-i ( heta +M-epsilon sin ( heta ))\ a_4=-i ( heta -M-epsilon sin ( heta ))\ end{cases}$

可以合併成兩組, 然後再次展開, 運算量有點大.

化簡的時候注意恆等式: $arg(e^{ix})=arctan( an (x))$ .

$egin{aligned} sum ext{Li}_1(e^{a}) =&frac{2}{i}arctan an left(frac{ heta -M-epsilon sin ( heta )+pi}{2} ight)\ +&frac{2}{i}arctan an left(frac{- heta -M+epsilon sin ( heta )+pi}{2} ight) end{aligned}$

注意到第二部分:

$int_0^{pi } arctancot left(frac{1}{2} ( heta +M-epsilon sin ( heta )) ight) , mathrm{d} heta =epsilon+frac{1}{4} left( +pi ^2-2 pi M ight)$

最後代回去大功告成!

$egin{aligned} igstar E=&M+frac{i}{2pi}int_0^pisum ext{Li}_1(e^{a}),mathrm{d}θ\ =&frac{1}{4} (2 M+pi )+frac{epsilon }{pi }+frac{1}{pi }int_0^piarctan an left(frac{ heta -M-epsilon sin ( heta )+pi}{2} ight)mathrm{d} heta end{aligned}$

代入數據就得到了 Stack Exchange 一樣的結果.

我對 $arctan( an (x))$ 這種寫法感到很不爽.

這個當然不能直接抵消, 由於 $arctan( an (x)) eq x$ , 我們作復展開.

$egin{aligned} arctan( an (x)) &=frac{1}{2} i log left(1+frac{e^{-i x}-e^{i x}}{e^{-i x}+e^{i x}} ight)-frac{1}{2} i log left(1-frac{e^{-i x}-e^{i x}}{e^{-i x}+e^{i x}} ight)\ &=frac{i}{2}log left(frac{2}{1+e^{2 i x}}iggl/frac{2 e^{2 i x}}{1+e^{2 i x}} ight)\ &=frac{i}{2}log (e^{-2 i x}) end{aligned}$

嚴格來說這兩者不是完全相等的, 因為這樣一來消掉了奇點.

不過積分的時候完全可以劃等號, 因為區間開閉完全不影響積分值.

$igstar E=frac{1}{4} (2 M+pi )+frac{epsilon }{pi }+frac{i}{2pi }int_0^pi log left(-e^{i (M+epsilon sin heta- heta)} ight)mathrm{d} heta$

綜上所述, 最後代入值, 我們得到了:

$egin{aligned} mathcal{D}_{;}=&frac{1}{pi }left[1+frac{i}{2}int_0^{pi}log left(-i e^{i (sin t-t)} ight);mathrm{d}t ight]\ mathcal{D}_°=&frac{1}{pi }left[1+frac{90i}{pi}int_0^{pi}log left(-i e^{i (pisin t/180-t)} ight);mathrm{d}t ight] end{aligned}$

(*真男人從不回頭看數值驗證*)(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).

$egin{aligned} mathcal{D}=&frac{1}{2} pi exp left(int_0^1 {frac{1}{pi x}arctanleft(frac{x log left(frac{sqrt{1-x^2}+1}{x} ight)(pi x+2) }{x^2 log ^2left(frac{sqrt{1-x^2}+1}{x} ight)-pi x-1} ight) , mathrm{d}x} ight)\ =&mathrm{arccot}left(1+frac{1}{2 pi }int_0^1{ log left(frac{pi ^2 left(1-x^2 ight)+4 left(sqrt{1-x^2} mathrm{arctanh}(x)+x ight)^2}{pi ^2 left(1-x^2 ight)+4 left(sqrt{1-x^2} mathrm{arctanh}(x)-x ight)^2} ight) , mathrm{d}x} ight) end{aligned}$

從這個結果上也能看出這個方法有多殘暴...

(*怎麼可以這麼暴力的說*)[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
Keplers equation, radiation problems and Meissels expansion
An exact analytical solution of Keplers Equation

題圖: https://www.pixiv.net/member_illust.php?mode=medium&illust_id=65616708