第一種・第二種ベッセル関数の漸近展開って何項で打ち切ればいいんですか?

はじめに

BesselJ,BesselY関数は以下の二階常微分方程式の特殊解で... っていうのは省略でいいですか。そういう数値計算では苦労する関数があるのです。そんなBesselJ,BesselY関数には、以下に示す便利な漸近展開が存在します。
$$J_{\nu}(x) \sim \sqrt{\dfrac{2 }{\pi x}}\sum_{n=0}^{\infty}\cos \left(\omega +\dfrac{n\pi}{2}\right)\prod_{k=1}^n \dfrac{4\nu ^2 - (2k-1)^2}{8kx}$$
$$Y_{\nu}(x) \sim \sqrt{\dfrac{2 }{\pi x}}\sum_{n=0}^{\infty}\sin \left(\omega +\dfrac{n\pi}{2}\right)\prod_{k=1}^n \dfrac{4\nu ^2 - (2k-1)^2}{8kx}$$
ただし空積は1とし、$\omega = x-\dfrac{\nu \pi}{2}-\dfrac{\pi}{4}$です。
見れば分かりますが、愚直に無限項足したら($\nu$が半整数でない限り)発散します。
つまり必要な精度と$\nu$を予め決めておいて、そこからある$x$に対してある程度の項数を計算し、その項数が精度を満たせばそれを返し、そうでなければ別の手法で計算するという事になります。
ただDesmosなどはこういうWhile的処理ができない(万が一出来るとしても相当複雑な手段を採らなければならないと思う)こともあり、必要な精度と$\nu$から最適な項数とつかえる$x$の範囲を決められると、色々と有利なのでは、と思ったのです。

まず具体的な値を求める

まず、項数を$p$(便宜上0項目からのカウントとする)、$x$をある正の実数とします。
そして、
$$T_p = \prod_{k=1}^{p}\dfrac{4\nu ^2-(2k-1)^2}{8kx} = \dfrac{1}{(8x)^p p!}\prod_{k=1}^{p}\left\{4\nu ^2-(2k-1)^2\right\}$$
を定めます。
ある必要精度を$\varepsilon$として$|T_p|<\varepsilon$を$x$について解きましょう。
$$\begin{eqnarray} \dfrac{1}{(8x)^p p!}\prod_{k=1}^{p}\left|4\nu ^2-(2k-1)^2\right| &<&\varepsilon\\ \dfrac{1}{\varepsilon p!}\prod_{k=1}^{p}\left|4\nu ^2-(2k-1)^2\right| &<& (8x)^p\\ \dfrac{1}{8(\varepsilon p!)^{1/p}}\prod_{k=1}^{p}\left|4\nu ^2-(2k-1)^2\right|^{1/p} &<& x\\ \end{eqnarray}$$
この左辺を$X_\nu(p)$としましょう。
見た目はこんな感じです。(必要精度はとりあえず$10^{-16}$で。)
$X_\nu(p)$($p$固定)
$X_\nu(p)$($\nu$固定)　なんか極小値みたいなものは見える
なんかすごく扱いづらい見た目をしています。それが先ほどあった「$\nu$が半整数の場合に限り有限項で終結する」という性質です。

とりあえず$\nu$を整数に限定して、$X_\nu(p)$が最小になる$p$を(非厳密的ですが)探索してみましょう。
$n=1,2,\cdots,100$について
$p=1,2,\cdots,500$まで探索した結果
なんか後半は線形に増加しているような気がします。ただ、今のままでは関数が扱いづらく何ともなりません。

変形して様子を見よう

よって、$X_\nu(p)$を変形してみましょう。
まず総積について、
$$\begin{eqnarray}\prod_{k=1}^{p}\left|4\nu ^2-(2k-1)^2\right| &=& 4^p \left(\prod_{k=1}^{p}\left|\nu+k-\dfrac{1}{2}\right|\right)\left(\prod_{k=1}^{p}\left|\nu-k+\dfrac{1}{2}\right|\right)\\ &=& 4^p \left|\dfrac{\Gamma(\nu + p+1/2)}{\Gamma\left(\nu + 1/2\right)}\cdot \dfrac{\Gamma(\nu+1/2)}{\Gamma(\nu-p+1/2)}\right|\\ &=& 4^p \left|\dfrac{\Gamma(\nu+p+1/2)}{\Gamma(\nu-p+1/2)}\right|\end{eqnarray}$$
ここで、先ほどの数値実験から$p>\nu-1/2$を仮定して、ガンマ関数の相補公式
$$\Gamma(-z) = -\dfrac{\pi}{\Gamma(z+1)\sin(\pi z)}$$
を適用すると、
$$4^p \left|\dfrac{\Gamma(\nu+p+1/2)}{\Gamma(\nu-p+1/2)}\right| = \dfrac{1}{\pi}4^p \left|\Gamma(p+\nu+1/2)\Gamma(p-\nu+1/2)\sin\left[\pi(p-\nu+1/2)\right]\right|$$
となります。罪を犯しているのは明らかに$\sin$ですね。

ここは大胆に$\sin$を外して考えましょう。
$$\tilde{X}_\nu(p) = \dfrac{1}{2}\left(\dfrac{\Gamma(p+\nu+1/2)\cdot\Gamma(p-\nu+1/2)}{\pi \varepsilon\cdot p!}\right)^{1/p}$$
(因みに$\nu$が半整数でなければ$p\rightarrow \infty$で$\sin (\pi(p-\nu+1/2))^{1/p} \rightarrow 1$となるため、大域的な挙動を調べるときには無視できます。)
Before
After
とりあえずこれを$\nu>0$、$p$が$\nu-1/2$より大きい実数を動く場合のモデル関数として考えることにします。

解析に踏み込む

ある$\nu>0$に対して$\tilde{X}_\nu(p)$を最小化する$p>\nu-1/2$を$p(\nu)$と置きます。
まず確かめるは一意性ですね。
前提として、$p\leq \nu-1/2$は解の候補になりえません。
この範囲では$|T_p|>|T_{p+1}|$が式より明らかに成り立ちますので、
より小さい$x$でも成り立つからです。
では$p>\nu-1/2$で最小値が一意に定まることを示しましょう。
今後使うことも見据えて
$$\begin{eqnarray}\Phi(\nu,p) &=& \ln \tilde{X}_\nu(p)\\&=& -\ln 2 +\dfrac{1}{p}\left[\ln\Gamma(p+\nu+1/2)+\ln\Gamma(p-\nu+1/2)-\ln\Gamma(p+1)-\ln(\pi \varepsilon)\right]\end{eqnarray}$$
を定めておきましょう。

では、この$\Phi(\nu,p)$の$p$に関する停留点から、$p(\nu)$の$\nu\to\infty$における漸近挙動を探りましょう。
まず、$\ln \Gamma(z)\sim z\ln z-z+\dfrac{1}{2}\ln(2\pi )-\dfrac12\ln(z)=\left(z-\dfrac12\right)\ln z -z+\dfrac12\ln (2\pi)$を代入して、さらに定数項を無視しましょう。
$$\Phi(\nu,p) \sim f(\nu,p) = \dfrac{1}{p}\left[(p+\nu)\ln(p+\nu+1/2)-(p+\nu+1/2)+(p-\nu)\ln(p-\nu+1/2)-(p-\nu+1/2)-p\ln p +p -\dfrac12\ln(p)+\tilde{C}\right]$$
ここで$C+\dfrac12\ln(2\pi)=\tilde{C}$としました。
さらに、$\nu\to\infty$で$(p+\nu)(\ln(p+\nu+1/2)-\ln(p+\nu))\to 1/2$なので、
$$f(\nu,p) = \dfrac{1}{p}\left[(p+\nu)\ln(p+\nu)-(p+\nu)+(p-\nu)\ln(p-\nu)-(p-\nu)-p\ln p +p -\dfrac12\ln(p)+\tilde{C}\right]$$
$p =s\nu$を代入してみましょう。
$$\begin{eqnarray}f(\nu,p) &=& \dfrac{1}{s}\left[(s+1)\ln(\nu(s+1))-(s+1)+(s-1)\ln(\nu(s-1))-(s-1)-s \ln (s\nu)+s -\dfrac1{2\nu} \ln(s\nu)+\dfrac{\tilde{C}}{\nu}\right]\\ &=& \ln \nu + \dfrac{1}{s}\left[(s+1)\ln(s+1)-(s+1)+(s-1)\ln(s-1)-(s-1)-s\ln s+s\right]-\dfrac{\ln \nu}{2s\nu}+\dfrac{1}{2s\nu}\left(2\tilde{C}-\ln s\right)\\ &=&\ln \nu +r_1(s)+r_2(s)\dfrac{\ln \nu}{\nu}+r_3(s)\dfrac{1}{\nu}\end{eqnarray}$$
$p$で微分するとこうなります。
$$\dfrac{\partial }{\partial p}f(\nu,p) =\dfrac{1}{\nu}\left[ r_1'(s)+r_2'(s)\dfrac{\ln \nu}{\nu}+r_3'(s)\dfrac{1}{\nu}\right]$$
これが0になる点が$p(\nu)$となります。
この関数の形に合わせて$s$を$\nu$の関数とみなしましょう。つまり、
$s=\alpha + \beta\dfrac{\ln \nu}{\nu}+\dfrac{\gamma}{\nu}$
とします。さらに、$r_1'(s),r_2'(s),r_3'(s)$に対して、点$\alpha$周りのTaylor展開
$$r_m'(s)\approx r_m'(\alpha)+r_m''(\alpha)(s-\alpha)$$
をします。(これは$\nu\to\infty$で成り立ちます。)
すると、
$$\dfrac{\partial }{\partial p}f(\nu,p) \approx \dfrac{1}{\nu}\left[r_1'(\alpha)+(\beta r_1''(\alpha)+r_2'(\alpha))\dfrac{\ln \nu}{\nu}+(\gamma r_1''(\alpha)+r_3'(\alpha))\dfrac{1}{\nu}\right]$$
これが0なので、恒等式から、
$$1.\:\:r_1'(\alpha)=0$$
$$2.\:\:\beta r_1''(\alpha)+r_2'(\alpha)=0$$
$$3.\:\: \gamma r_1''(\alpha)+r_3'(\alpha)=0$$
です。

1の式

$$r_1'(\alpha) = \dfrac{\alpha+\ln(\alpha-1)-\ln(\alpha+1)}{\alpha^2}=0$$で、$\alpha>0より、$
$$\begin{eqnarray}\alpha &=& \ln\left(\dfrac{\alpha+1}{\alpha-1}\right)\\&=& 2\coth^{-1}(\alpha)\end{eqnarray}$$
$\alpha\approx 1.5434046384182084479587097400533155536$

2の式

$$\begin{eqnarray}r_1''(\alpha)&=&\dfrac{1}{\alpha^3}\left[\dfrac{2\alpha}{\alpha^2-1}-\alpha+4\coth^{-1}(\alpha)\right]=\dfrac{\alpha^2+1}{\alpha^2(\alpha^2-1)}\end{eqnarray}$$
$$r_2'(\alpha) = \dfrac{1}{2\alpha^2}$$
より、
$\beta = -\dfrac{r_2'(\alpha)}{r_1''(\alpha)} = -\dfrac{\alpha^2-1}{2(\alpha^2+1)}\approx -0.20432552927071863993986485368419686$

3の式

$$r_3'(\alpha) = \dfrac{\ln \alpha-2\tilde{C}-1}{2\alpha^2}$$より、
$$\begin{eqnarray}\gamma = -\dfrac{r_3'(\alpha)}{r_1''(\alpha)}&=&\dfrac{\alpha^2-1}{2(\alpha^2+1)}(2\tilde{C}+1-\ln \alpha )\\&=&\dfrac{\alpha^2-1}{2(\alpha^2+1)}(-2\ln\varepsilon-\ln \dfrac{\pi}{2}+1-\ln \alpha)\\&\approx& 0.20432552927071863993986485368419686(-2\ln\varepsilon +0.11442651434895242834504107379726465)\end{eqnarray}$$
特に$\varepsilon=10^{-16}$のとき、$\gamma\approx15.078641628246962240520586$
$s = \alpha + \beta\dfrac{\ln \nu}{\nu}+\gamma\dfrac{1}{\nu},\:\:p=s\nu$より、
$$p(\nu)\sim\tilde{p}(\nu)= \alpha \nu+\beta\ln \nu+\gamma$$であることが分かりました。

欲

もっと項、知りたいですよね。
なんとなく、このTaylor展開から$\dfrac{(\ln\nu)^2}{\nu},\:\dfrac{\ln \nu}{\nu}$の項がありそうなのは分かります。しかも、これらは$r_m(s)$のTaylor展開によってのみ現れる項であり、スターリング近似の項を増やしたり$1/2$のシフトを考慮しても出てこないので、このまま
$s=\alpha+\beta\dfrac{\ln \nu}{\nu}+\gamma\dfrac{1}{\nu}+\delta\dfrac{(\ln \nu)^2}{\nu^2}+\eta\dfrac{\ln \nu}{\nu^2}$で計算続行。
$$r_m'(s)\approx r_m'(\alpha)+r_m''(\alpha)(s-\alpha)+\dfrac12r_m'''(\alpha)(s-\alpha)^2$$
$$\dfrac{\partial }{\partial p}f(\nu,p) \approx \dfrac{1}{\nu}\left[r_1'(\alpha)+(\beta r_1''(\alpha)+r_2'(\alpha))\dfrac{\ln \nu}{\nu}+(\gamma r_1''(\alpha)+r_3'(\alpha))\dfrac{1}{\nu}+\left(\delta r_1''(\alpha)+\beta r_2''(\alpha)+\dfrac12\beta^2 r_1'''(\alpha)\right)\dfrac{(\ln \nu)^2}{\nu^2}+\left(\eta \:r_1''(\alpha)+\beta \gamma r_1'''(\alpha)+\gamma r_2''(\alpha)+\beta r_3''(\alpha)\right)\dfrac{\ln \nu}{\nu^2}\right]$$
$$4.\:\:\delta r_1''(\alpha)+\beta r_2''(\alpha)+\dfrac12\beta^2 r_1'''(\alpha)=0$$
$$5.\:\:\eta \:r_1''(\alpha)+\beta \gamma r_1'''(\alpha)+\gamma r_2''(\alpha)+\beta r_3''(\alpha)=0$$

4の式

$r_1''(\alpha)=\dfrac{\alpha^2+1}{\alpha^2(\alpha^2-1)},\: r_2''(\alpha)=-\dfrac{1}{\alpha^3},\:r_1'''(\alpha)=-\dfrac{4(\alpha^4+\alpha^2-1)}{(\alpha-1)^2(\alpha+1)^2\alpha^3}$
$\beta= -\dfrac{\alpha^2-1}{2(\alpha^2+1)}$より、
$$\begin{eqnarray} \delta &=&\dfrac{\alpha^2(\alpha^2-1)}{\alpha^2+1}\cdot\dfrac{\alpha^2-1}{2(\alpha^2+1)}\cdot\left(\dfrac{1}{\alpha^3}-\dfrac{1}{2}\dfrac{\alpha^2-1}{2(\alpha^2+1)}\dfrac{4(\alpha^4+\alpha^2-1)}{\alpha^3(\alpha^2-1)^2}\right)\\&=&\dfrac{\alpha^2-1}{\alpha^2+1}\cdot\dfrac{\alpha^2-1}{2(\alpha^2+1)}\cdot\dfrac{\alpha}{(\alpha^2-1)(\alpha^2+1)}\\&=&\dfrac{\alpha(\alpha^2-1)}{2(\alpha^2+1)^3}\\&\approx& 0.02756957617749993357807088841\end{eqnarray}$$

5の式

$$r_3''(\alpha)=\dfrac{2\tilde{C}+3/2-\ln \alpha}{\alpha^3},\:\:-\beta(2\tilde{C}+3/2-\ln \alpha) =\gamma-\dfrac{\beta}{2}$$より、
$$\begin{eqnarray}\eta &=&\dfrac{\alpha^2-1}{\alpha(\alpha^2+1)}\left[\dfrac{4\beta\gamma(\alpha^4+\alpha^2-1)}{(\alpha^2-1)^2}+\gamma-\beta(2\tilde{C}+3/2-\ln \alpha)\right]\\ &=& \dfrac{\alpha^2-1}{\alpha(\alpha^2+1)}\left[\dfrac{4\beta\gamma(\alpha^4+\alpha^2-1)}{(\alpha^2-1)^2}+2\gamma-\dfrac{\beta}{2}\right] \end{eqnarray}$$
特に$\varepsilon=10^{-16}$のとき、$\eta\approx-4.042062381553583$

欲2

ではここで具体的な$p(\nu)$と、
$\tilde{p}(\nu)=\alpha+\beta\ln \nu +\gamma +\delta\dfrac{\ln \nu^2}{\nu}+\eta\dfrac{\ln \nu}{\nu}$の差分を見てみましょう。
$p(\nu)-\tilde{p}(\nu)$
んー。$\nu$掛けましょう。
$\nu(p(\nu)-\tilde{p}(\nu))$
どうやら$147.9$近辺という結構大きい値に収束してそうなので、そこまでは出したいですね。
という事で、スターリング近似を一段階進めて、さらに$1/2$シフトの近似をしっかり考慮して、$\dfrac{1}{\nu}$の項を出していきましょう。
$$\Phi(\nu,p)\sim f(\nu,p) = \ln \nu + r_1(s)+r_2(s)\dfrac{\ln \nu}{\nu}+r_3(s)\dfrac{1}{\nu}+r_4(s)\dfrac{1}{\nu^2}$$
なる$r_4(s)$を考えます。
まずスターリング近似の$\dfrac{1}{z}$の項で
$$\dfrac{1}{12}\left(\dfrac{1}{p+\nu+1/2}+\dfrac{1}{p-\nu+1/2}-\dfrac{1}{p}\right)\sim \dfrac{1}{12}\left(\dfrac{1}{p+\nu}+\dfrac{1}{p-\nu}-\dfrac{1}{p}\right)+O\left(\dfrac{1}{\nu^2}\right)$$
$(p\pm\nu)(\ln(p\pm\nu+1/2)-\ln(p\pm\nu))\sim \dfrac{1}{2}-\dfrac{1}{8(p\pm\nu)}+O\left(\dfrac{1}{\nu^2}\right)$
より、(すでに$1/2$は相殺済み)
これを$p=s\nu$で除して整理すると、
$$\begin{eqnarray}r_4(s)&=&-\dfrac{1}{24}\left(\dfrac{1}{s+1}+\dfrac{1}{s-1}+\dfrac{2}{s}\right)\end{eqnarray}$$
よって追加するべき項は
$$s=\alpha+\beta\dfrac{\ln\nu}{\nu}+\gamma\dfrac{1}{\nu}+\delta\dfrac{(\ln\nu)^2}{\nu^2}+\eta\dfrac{\ln\nu}{\nu}+\zeta\dfrac{1}{\nu^2}$$として、
$$\cdots + \left(\zeta r_1''(\alpha)+\dfrac12\gamma^2r_1'''(\alpha)+\gamma r_3''(\alpha)+r_4'(\alpha)\right)\dfrac{1}{\nu^2}$$
これが0になるので、
$$\zeta r_1''(\alpha)+\dfrac12\gamma^2r_1'''(\alpha)+\gamma r_3''(\alpha)+r_4'(\alpha)=0$$
$$r_4'(s) = \dfrac16\left(\dfrac{s}{s^2-1}+\dfrac{1}{s^3}\right),\: \: r_3''(s) = -\dfrac{\gamma}{\beta}+\dfrac12$$
その他諸々を代入して頑張ってまとめると、
$$\begin{eqnarray}\zeta &=&\frac{2 \alpha \gamma^2}{\alpha^4 - 1} + \frac{\beta \gamma}{\alpha} - \frac{2\alpha^4 - 2\alpha^2 + 1}{6\alpha(\alpha^4 - 1)}\\&=&\dfrac{2(3\beta\gamma-1)\alpha^4+2(6\gamma^2+1)\alpha^2-6\beta\gamma-1}{6\alpha(\alpha^4-1)}\end{eqnarray}$$
特に$\varepsilon=10^{-16}$のとき、$\zeta\approx 147.973021793977965$

結果

ある$\nu$に対するBessel関数の適切な項数$p(\nu)$は、
$p(\nu)\sim \tilde{p}(\nu)=s\nu=\alpha\nu+\beta\ln\nu+\gamma+\delta\dfrac{(\ln\nu)^2}{\nu}+\eta\dfrac{\ln \nu}{\nu}+\zeta\dfrac{1}{\nu}$
で表される。$(\nu\to\infty)$
ここで、
$\alpha = 2\coth^{-1}(\alpha)\quad (\alpha>0)$
$\alpha\approx1.5434046384182084479587097400533155536978837647192626945972413624898877396294736$
$\beta=-\dfrac{\alpha^2-1}{2(\alpha^2+1)}$
$\beta\approx-0.2043255292707186399398648536841968636454330742213819839124962484407731537145223$
$$ \begin{eqnarray}\gamma&=&\dfrac{\alpha^2-1}{2(\alpha^2+1)}(-2\ln\varepsilon-\ln \dfrac{\pi}{2}+1-\ln \alpha)\\&=&\beta(2\ln\varepsilon+\ln \dfrac{\pi}{2}-1+\ln \alpha)\end{eqnarray}$$
$\gamma\approx -0.2043255292707186399398648536841968636454330742213819839124962484407731537145223(2\ln \varepsilon-0.1144265143489524283450410737972646595395780315982045081110556037883564361127241)$
$\varepsilon=10^{-16}$のとき、
$\gamma\approx15.078641628246962240520586421951392099957661752771859179657196$
$$\begin{eqnarray}\delta&=&\dfrac{\alpha(\alpha^2-1)}{2(\alpha^2+1)^3}\\&=&-\beta\dfrac{\alpha}{(\alpha^2+1)^2}\end{eqnarray}$$
$\delta\approx0.0275695761774999335780708884105410615619800990553915104253666694606119555152969$
$$\begin{eqnarray}\eta &=& \dfrac{\alpha^2-1}{\alpha(\alpha^2+1)}\left[\dfrac{4\beta\gamma(\alpha^4+\alpha^2-1)}{(\alpha^2-1)^2}+2\gamma-\dfrac{\beta}{2}\right]\end{eqnarray}$$
$\varepsilon=10^{-16}$のとき、
$$\eta \approx -4.04206238155358331398716673019554456553065184556689$$
$$\begin{eqnarray}\zeta &=&\frac{2 \alpha \gamma^2}{\alpha^4 - 1} + \frac{\beta \gamma}{\alpha} - \frac{2\alpha^4 - 2\alpha^2 + 1}{6\alpha(\alpha^4 - 1)}\\&=&\dfrac{2(3\beta\gamma-1)\alpha^4+2(6\gamma^2+1)\alpha^2-6\beta\gamma-1}{6\alpha(\alpha^4-1)}\end{eqnarray}$$
$\varepsilon=10^{-16}$のとき、
$$\zeta \approx 147.97302179397796500561970060639325980493168471551414$$
ここに 最終的な結果と、それを項数として適用したBessel関数の漸近展開を載せました。($\nu$が小さい時は適当な近似式で代用しています。)

先行研究は見つけられませんでした。心当たりのある方や、この記事の間違いを見つけた方は、コメントして頂けると幸いです。