不完全ガンマ関数の連分数って何項で打ち切ればいいんですか?　後編

はじめに

この記事は後編です。前編は こちら から。
今回絡む分だけ言うと、専ら$x\geq a> 0$で用いられる不完全ガンマ関数の連分数表示
$$\Gamma(a,x) = e^{-x}x^a\dfrac{1}{1+x-a+\dfrac{1\cdot(a-1)}{3+x-a + \dfrac{2\cdot(a-2)}{5+x-a + \dfrac{3\cdot (a-3)}{7+x-a +\cdots}}}}$$
を十分収束させるためには何項計算すればよいかを推定し、特に収束しづらい$x=a$について、高々の項数$n\sim O(a^{1/3}\ln(a)^{2/3})$という、通説の$n\sim O(\sqrt{a})$よりも厳しい評価を得ました。
今回はさらにこれを詰めようと思います。

なぜ詰めようと思ったのか?

問題は 数値計算bot さん(a.k.a Y.K. さん)に計算して頂いた項数のデータを解析していたときに起こりました。
$a$と$n/a^{1/3}$の関係(凡例は必要精度桁数$d$桁)
明らかに$a\to\infty $で$\dfrac{n}{a^{1/3}}$が収束しようとしています。
$n\sim O(a^{1/3}\ln(a)^{2/3})$ならば少しあり得ない挙動です。
つまり、オーダーを改善する余地があるという事です。

近似アプローチを変えてみる

前回は鞍点法で分母漸化式...というかラゲール陪多項式(ソニン多項式)$L_n^{-a}(-x)$を近似していましたが、
今度は、分母漸化式そのものを離散WKB近似という手法を用いて近似していこうと思います。

WKB近似とは(飛ばし可)

これは元々量子力学でシュレーディンガー方程式を近似で求める手法です。
$$\left[\dfrac{1}{2m}\left(\dfrac{\hbar}{i}\dfrac{\text{d}}{\text{d}x}\right)^2+V(x)-E\right]\psi(x)=0$$
について$\psi(x)=\exp\left(\dfrac{i}{\hbar}S(x)\right)$と決め打ちします。(これは微分方程式を解くときに普通に使います。)
$$\dfrac{1}{2m}\left[\dfrac{\hbar}{i}\dfrac{\text{d}^2 S(x)}{\text{d}x^2}+\left(\dfrac{\text{d}S(x)}{\text{d}x}\right)^2\right]+V(x)-E=0$$
代入するとこうなります。
で、$\hbar\to 0$で古典力学に対応していて、$\hbar$が十分小という仮定のもと、$S(x)$を$\dfrac{\hbar}{i}$に対して展開します。(これがWKB近似です。)
$$S(x)\approx S_0(x) + S_1(x)\dfrac{\hbar}{i}+\cdots$$
これの1次近似まで取り入れます。
\begin{align*} \dfrac{1}{2m}\left[\dfrac{\hbar}{i}\left(S_0''(x)+\dfrac{\hbar}{i}S_1''(x)\right) +\left(S_0'(x)^2 + 2\dfrac{\hbar}{i}S_0'(x)S_1'(x)\right)\right]+V(x)-E&=0\\ \dfrac{1}{2m}\left[S_0'(x)^2 +\dfrac{\hbar}{i}\left(S_1''(x) + 2S_0'(x)S_1'(x)\right)\right]+V(x)-E&=0\\ \end{align*}
$\left(\dfrac{\hbar}{i}\right)^2$の項は1次近似につき無視です。
0次について、
$$\dfrac{1}{2m}S_0'(x)^2=E-V(x)$$
$$S_0'(x) = \pm\sqrt{2m(E-V(x))}$$
$$S_0 (x) \pm\int_{x_0}^x\sqrt{2m(E-V(\xi))} \: \text{d}\xi$$
と求まります。今後のため$S_0'(x)=\pm p(x)$と置いておきます。
1次ついて、
\begin{align*} S_0''(x)+2S_0'(x)S_1'(x)=0\\ S_1'(x) = \dfrac{1}{2}\dfrac{S_0''(x)}{S_0'(x)}\\ S_1'(x) = \dfrac{1}{2}\dfrac{\pm p'(x)}{\pm p(x)} = \dfrac{1}{2}\dfrac{p'(x)}{p(x)}\\ S_1(x) = \dfrac{1}{2}\ln(p(x))+C \end{align*}
と求まりましたね。これを纏めれば、シュレーディンガー方程式を近似的に解けたという事になります。
$$\psi(x) = \exp\left(\dfrac{i}{\hbar}S(x)\right)\approx \dfrac{C}{\sqrt{p(x)}}\exp\left(\pm\dfrac{i}{\hbar}\int_{x_0}^x p(\xi)\: \text{d}\xi\right)$$
WKB近似のことを知るにはここまでで十分です。(古典的禁止領域とか$p(x)=0$の場合の話とかは省略)

離散WKB近似

これを差分方程式に利用したのが離散WKB近似です。
漸化式は$q_{n+1} = (2n+1)q_n +n(a-n)q_{n-1}$
まぁまぁ、まずは$q_n = \exp(S(n))$
と置いて、$S(n\pm 1) \approx S(n) \pm \Delta S(n)$と近似しましょう。
$$e^{S(n)+\Delta S(n)} = (2n+1)e^{S(n)}+n(a-n)e^{S(n)-\Delta S(n)}$$
$$e^{\Delta S(n)} = (2n+1) +n(a-n)e^{-\Delta S(n)}$$
$e^{\Delta S(n)}=\lambda_n$とおいて整理すると、
$$\lambda_n^2 -(2n+1)\lambda_n -n(a-n) = 0$$
これは三項間漸化式の特性方程式に対応します。この方程式は$n$によって変化するので、局所特性方程式とか呼ばれたりするかもしれません。
これの解を$\lambda_n^{\pm}$とします。大きい方の解が$\lambda_n^{+}$、小さい方が$\lambda_n^{-}$です。
一応解を並べておきます。
$$\lambda_n^{\pm} = \dfrac{(2n+1)\pm\sqrt{(2n+1)^2 +4n(a-n)}}{2} = \left(n + \dfrac{1}{2}\right)\pm R_n$$
$$R_n = \sqrt{\left(n+\dfrac12\right)^2+n(a-n)}\approx \sqrt{na}$$
$$\begin{align*}\\ \quad\\ \quad\\\end{align*}$$
つまり$\Delta S(n) = \ln (\lambda_n^{\pm})$、$\pm p(x)$が$\ln (\lambda_n^{\pm})$に対応するわけですね。

($\lambda_n^{-}$が負ですが、$\ln(\lambda_n^{-})$を直接扱うことは結果的に無くなるので一旦無視します。)

積分の離散ver.である累積和を取ります。
$$S(n) \approx \sum_{k=1}^{n-1}\ln(\lambda_k^{\pm}) = \ln\left(\prod_{k=1}^{n-1} \lambda_{k}^{\pm}\right)$$
つまり、$q_n \sim \displaystyle \prod_{k=1}^{n-1}\lambda_{k}^{\pm}$が成り立ちます。
($\lambda_k^{\pm}$は、$\lambda_k^{+}$の総積と$\lambda_k^{-}$の総積の線形結合と取るという意味で用いています。)
今回のような近似では$|\lambda_k^{+}|>|\lambda_k^{-}|$が、今回考える範囲では常に成り立つので、優越解である$\lambda_k^{+}$のみをとって、
$q_n\sim \displaystyle\prod_{k=1}^{n-1}\lambda_k^{+}$とします。
これはまだ0次近似です。しかし1次近似するための式がもうありません。
ここで、WKB近似の1次近似が$\dfrac{C}{\sqrt{p(x)}}$という振幅を規定していたことを思い出しましょう。つまり$q_n$に振幅$A(n)$を前因子として要請して、
$$q_n \sim A(n)\prod_{k=1}^{n-1}\lambda_k^{+}$$とすれば良いわけです。
では、またこれを漸化式に代入しましょう。
\begin{align*} A(n+1)\prod_{k=1}^{n}\lambda_k^{+} = (2n+1)A(n)\prod_{k=1}^{n-1}\lambda_k^{+}+n(a-n)A(n-1)\prod_{k=1}^{n-2}\lambda_k^{+}\\ A(n+1)\lambda_n^{+}\lambda_{n-1}^{+} = (2n+1)A(n)\lambda_{n-1}^{+}+n(a-n)A(n-1)\\ \end{align*}
解と係数の関係$2n+1 = \lambda_n^+ +\lambda_{n}^-$, $n(a-n) = -\lambda_n^+ \lambda_n^{-}$を用いると、
$$A(n+1)\lambda_n^{+}\lambda_{n-1}^{+} = A(n)\lambda_{n-1}^+ (\lambda_n^+ + \lambda_n^{-}) -A(n-1)\lambda_n^+\lambda_n^{-}$$
ここで、$A(n\pm1) \approx A(n)\pm \Delta A(n)$と近似して代入、整理すると、
$$A(n)\left[\lambda_n^+\lambda_{n-1}^{+}-\lambda_{n-1}^+(\lambda_n^++\lambda_n^{-})+\lambda_n^+\lambda_n^-\right]+\Delta A(n)\left[\lambda_n^{+}\lambda_{n-1}^{+}-\lambda_n^+\lambda_n^{-}\right]=0$$

$$A(n)\lambda_n^{-}\left(\lambda_n^{+}-\lambda_{n-1}^{+}\right)+\Delta A(n)\lambda_n^{+}(\lambda_{n-1}^{+}-\lambda_n^{-})=0$$
ここで、$\lambda_{n-1}^+-\lambda_n^{-}\approx\lambda_n^{+}-\lambda_n^{-}$として$\dfrac{\Delta A(n)}{A(n)}=$の形にします。
$$\dfrac{\Delta A(n)}{A(n)} = -\dfrac{\lambda_n^{-}}{\lambda_n^{+}(\lambda_n^{+}-\lambda_n^{-})}(\lambda_{n}^{+}-\lambda_{n-1}^{+})$$
ここで正当といえば正当ですがえぐい近似をします。
$$\lambda_n^{\pm}= \left(n+\dfrac{1}{2}\right)\pm R_n\approx\pm R_n\approx \pm \sqrt{an}$$
です。$n\sim a^{1/3},\:\:\sqrt{na}\sim a^{2/3}$なので、オーダー的には問題ありません。
(後々、矛盾するような演算を取りますが、それはオーダーが違うからであって、この$A(n)$の決定に関してはこれくらい大雑把でOKです)
また$\lambda_n^{+}-\lambda_n^{-} = 2R_n$より(これは正確)
$$\dfrac{\Delta A(n)}{A(n)}\approx \dfrac{R_n(R_n-R_{n-1})}{2R_n^2} \approx\dfrac{\Delta R_n}{2R_n}$$
これを積分すると、
$$\ln A(n) \approx \dfrac{1}{2}\ln R_n +C$$
となるので、$A(n)\approx C\sqrt{R_n}$
と表すことが出来ました。
よって、
$$q_n \approx C\sqrt{R_n}\prod_{k=1}^{n-1}\lambda_k^{+}$$
です。あえてこのまま誤差上界$\displaystyle\dfrac{1}{q_nq_{n+1}}\prod_{k=1}^{n}k(a-k)$につっこみます。
解と係数の関係から$\displaystyle \prod_{k=1}^{n}k(a-k) =-\prod_{k=1}^{n}\lambda_k^{+}\lambda_{k}^{-}$が成り立つので
\begin{align*} \dfrac{1}{q_nq_{n+1}}\prod_{k=1}^{n}k(a-k) &\approx \dfrac{\displaystyle\prod_{k=1}^{n}\lambda_k^{+}\left|\lambda_{k}^{-}\right|}{\displaystyle C^2\sqrt{R_nR_{n+1}}\prod_{k=1}^{n-1}\lambda_k^{+}\prod_{k=1}^{n}\lambda_k^{+}}\\ &= \dfrac{\lambda_n^{+}}{C^2\sqrt{R_nR_{n+1}}}\prod_{k=1}^{n}\dfrac{\left|\lambda_k^{-}\right|}{\lambda_k^{+}} \end{align*}
ここで、$\rho_k = \dfrac{|\lambda_k^{-}|}{\lambda_k^{+}}$とすると、
$$\rho_k = \dfrac{R_k-(k+1/2)}{R_k+(k+1/2)}\approx \dfrac{\sqrt{ak}-(k+1/2)}{\sqrt{ak}+(k+1/2)}\approx 1-\dfrac{2k+1}{\sqrt{ak}}+\dfrac{(2k+1)^2}{2ak}+\cdots$$
より、
$$\ln (\rho_k) \approx -\dfrac{2k+1}{\sqrt{ak}}-\dfrac{(2k+1)^3}{12(ak)^{3/2}}$$

\begin{align*} \ln\left(\prod_{k=1}^n\rho_k\right) &= \sum_{k=1}^n \ln(\rho_k)\\ &\approx \int_0^n \ln (\rho_x)\: \text{d}x \end{align*}
であるから、第一項目について、
\begin{align*} \int_0^n -\dfrac{2x+1}{\sqrt{ax}}\: dx &= \dfrac{1}{\sqrt{a}}\left[-\dfrac{4x^{3/2}}{3} -2\sqrt{x}\right]_0^n\\ &= -\dfrac{4n^{3/2}}{3\sqrt{a}}-2\sqrt{\dfrac{n}{a}}\\ &\approx -\dfrac{4n^{3/2}}{3\sqrt{a}} \qquad(\because n\sim O(a^{1/3})) \end{align*}
第二項目について、粗く見積もっても
\begin{align*} \int_0^n\dfrac{(2x+1)^3}{12(ax)^{3/2}}\: dx < n\cdot \dfrac{(2n+1)^3}{12(an)^{3/2}}=a^{-3/2}\cdot O(n^{5/2}) \end{align*}
$n\sim O(a^{1/3})$であることは分かっているので、これは$O(a^{-2/3})$となり、$a\to\infty$で$0$となります。
一回、今吟味している式をもう一度見せます。
$$\dfrac{\lambda_n^{+}}{C^2\sqrt{R_nR_{n+1}}}\prod_{k=1}^{n}\rho_k$$
この総積の部分はどうにかなりました。では、この$\dfrac{\lambda_n^{+}}{C^2\sqrt{R_nR_{n+1}}}$を見積もりましょう。
しかしながら、$\lambda_n^{+}\approx R_n,\:\sqrt{R_nR_{n+1}}\approx R_n$なので、
$$\dfrac{\lambda_n^{+}}{C^2\sqrt{R_nR_{n+1}}} \approx \dfrac{1}{C^2}$$
です。
よって
$$\dfrac{1}{q_nq_{n+1}}\prod_{k=1}^{n}k(a-k)\sim \dfrac{1}{C^2}\exp\left(-\dfrac{4n^{3/2}}{3\sqrt{a}}\right) + O(1)$$です。
$$\dfrac{1}{C^2}\exp\left(-\dfrac{4n^{3/2}}{3\sqrt{a}}\right) <10^{-d}$$
を$n$について解くと、
$$ n_h > a^{1/3}\left[\dfrac{3}{4}\biggl(d\ln(10)-2\ln (C)\biggr)\right]^{2/3}$$
よって、$n\sim O(a^{1/3})$

前編での$O(a^{1/3}\ln(a)^{2/3})$からさらに改善しました。

この$C$は数値的に決定することになります。
ここで、 前編 の数値計算の結果から、$C\approx\dfrac{1}{2}$であることが分かります。
こじつけながら理由を考えるならば、
$q_n$の一般項は本来以下のように与えられます。
$$ q_n \approx C_{+}\sqrt{R_n}\prod_{k=1}^{n-1}\lambda_k^{+}+C_{-}\sqrt{R_n}\prod_{k=1}^{n-1}\lambda_k^{-} $$
さきほど、$|\lambda_k^{+}|>|\lambda_k^{-}|$とおきましたが、大きさはどちらも$R_k$とほぼ同じです。
$\lambda_{k}^-$の符号は負なので時々打ち消しあいますが、大域的には
$$ q_n \approx \left(C_{+}+C_{-}\right)\sqrt{R_n}\prod_{k=1}^{n-1}R_k $$
と近似されて、$C_{+}\approx C_{-}$かつ初期因子が$1$と仮定すると、$C_{+}\approx\dfrac{1}{2}$が出てもおかしくはありません。
とはいえ、これは厳密性の欠片もない推定ですので、あくまで示唆程度にしておきましょう。

とはいえ、これで数値実験との整合が確かめられました。

結論

不完全ガンマ関数の連分数を、$x=a,\: a\to\infty$で、必要精度$10^{-d}$になるまで回した時の項数$n$は
$$a^{1/3}\left[\dfrac{3}{4}\biggl(d\ln(10)-2\ln (C)\biggr)\right]^{2/3}$$
特に数値実験的には
$$a^{1/3}\left[\dfrac{3}{4}\biggl(d\ln(10)+\ln(4)\biggr)\right]^{2/3}$$
で与えられる。

謝辞

この厳密なオーダーを明らかにした数値実験に協力して下さった 数値計算bot さん(a.k.a Y.K. さん)には改めて感謝申し上げます。