1
応用数学解説
文献あり

不完全ベータ関数の連分数って何項で打ち切ればいいんですか? 前編

146
0
$$$$

はじめに

まさかの何項打ち切りシリーズ第三弾です。
こんなに続くとは思いませんでした。
第三弾は不完全ベータ関数の連分数表示です。
$$B_x(a,b) = \dfrac{x^a (1-x)^b}{a}\dfrac{1}{1+\dfrac{d_1}{1+\dfrac{d_2}{1+\cdots}}}$$
で、
$$\begin{align*} d_{2m} &= \dfrac{m(b-m)x}{(a+2m-1)(a+2m)}\\ \quad\\ d_{2m+1} &= -\dfrac{(a+m)(a+b+m)x}{(a+2m)(a+2m+1)} \end{align*} $$
です。結構複雑ですね。
で、この条件において、不可避かつもっとも収束しづらい点は$x = \dfrac{a+1}{a+b+2}$なので、この点について主に考えていきます。
なお、収束項数については$n$を割り当てます。$n$が偶数の場合が$2m$、奇数の場合が$2m+1$ということですね。

まずは数値実験

まず愚直にこれで修正Lentz法(連分数を前向きで計算するアルゴリズム)を動かすとどうなるかを 数値計算bot さん(a.k.a. Y.K. さん)に計算してもらいました。
...ちょっと値が暴れすぎてグラフには出来なかったのですが、
有効精度$d$桁に対して、収束項数$n$$1.5\times 10^d$項程度になりました。
(収束条件は$|\Delta-1|\leq 10^{-d}$)
この結果が得られたため、とりあえず項数$n$$a,b$に非依存的に収束することを示そうとしました。(ネタバレすると失敗しますが、なぜ修正Lentz法が$1.5\times 10^{d}$項で収束したと判定したかが分かるので書きます。留意してください。)

数式でシミュレーション

ここでは、修正Lentz法を近似しながら数式で追っていくことでシミュレーションしていきます。
簡単のため
$$1 + \dfrac{d_1}{1+\dfrac{d_2}{1+\dfrac{d_3}{1+\dfrac{d_4}{1+\cdots}}}}$$を計算します。つまり、最初の1ステップを省いた形で考えるという事です。
ここで、$b\gg a\gg n,m$とします。
$$d_{2m+1} = -\dfrac{(a+m)(a+b+m)(a+1)}{(a+2m)(a+2m+1)(a+b+2)} \approx \dfrac{(a+m)(a+1)}{(a+2m)(a+2m+1)}\approx -1 + \dfrac{3m}{a}$$
と近似できます。
ただし$d_1 = -1 + \dfrac{2}{b}$です。
また、
$$d_{2m} = \dfrac{m(b-m)x}{(a+2m-1)(a+2m)}\approx \dfrac{mba}{a^2b} = \dfrac{m}{a}$$
と近似できます。
修正Lentz法の初期値$C_0 = 1, \: D_0 = 0$を定めて、いざLentz法を回しましょう。
$n=1$では
$$C_1 = 1+\dfrac{d_1}{C_0} = 1 - 1 + \dfrac{2}{b} =\dfrac{2}{b}$$
これを微小値$\varepsilon$と置きます。つまり$C_1 = \varepsilon$です。
$$D_1 = \dfrac{1}{1+d_1D_0} = \dfrac{1}{1+0} = 1$$
となります。
$\Delta_1 = C_1D_1 = \varepsilon$ですね。
$n=2$ときは
$$C_2 = 1 + \dfrac{d_2}{C_1} =1 + \dfrac{1}{a\varepsilon} \approx \dfrac{1}{a\varepsilon}$$
$$D_2 = \dfrac{1}{1+d_2D_1} = \dfrac{1}{1+1/a} \approx 1-\dfrac{1}{a}\approx 1$$
と近似できます。
$\Delta_2 = C_2D_2 \approx \dfrac{1}{a\varepsilon}$となります。つまり、$n=1$での$\varepsilon$倍は、$n=2$で相殺されます。...なんとなくやばい感じが伝わるでしょうか。

では$n=3$いきましょう。
$$C_3 = 1 + \dfrac{d_3}{C_2} \approx 1 + \dfrac{-1 + 3/a}{1/(a\varepsilon)} \approx 1 + a\varepsilon \left(-1 + \dfrac{3}{a}\right)\approx 1$$
$$D_3 = \dfrac{1}{1 + d_3D_2} \approx \dfrac{1}{1+(-1+3/a)(1-1/a)}\approx \dfrac{1}{1-1+4/a}=\dfrac{a}{4}$$
よって$\Delta_3 = C_3D_3\approx \dfrac{a}{4}$ですね。
$n=4$のとき、
$$C_4 = 1 + \dfrac{d_4}{C_3} \approx 1 + \dfrac{2/a}{1} = 1 + \dfrac{2}{a}\approx 1$$
$$D_4 = \dfrac{1}{1+d_4D_3} \approx \dfrac{1}{1 + (2/a)\cdot(a/4)} = \dfrac{2}{3}$$
よって、$\Delta_4 =C_4D_4\approx \dfrac{2}{3},\: |\Delta_4-1| = \dfrac13$となります。
...これはチュートリアルです。いわば数学的帰納法の最初の場合を確かめたにすぎません。

ここで、天下り的ですが、$k\geq 1$について、$C_{4k}\approx 1+\dfrac{\gamma_k}{a},\: D_{4k} \approx \delta_k\: (\delta_k <1)$

を仮定します。何でこんな仮定が出来るのかというと計算結果を見ているからです。計算結果の規則性から見通しが立ったわけです。
また$w_k = \dfrac{1}{1-\delta_k}$を定めます。見やすさのためです。ここで$n=4k+4 = 4(k+1)$まで修正Lentz法を進めて、そこでもその仮定が守られることを証明しましょう。

$n=4k+1$のとき、
$$C_{4k+1} = 1 + \dfrac{d_{4k+1}}{C_{4k}} \approx 1 + \dfrac{-1 +\frac{6k}{a}}{1+\frac{\gamma_k}{a}}\approx 1 + \left(-1 + \dfrac{6k}{a}\right)\left(1-\dfrac{\gamma_k}{a}\right)\approx \dfrac{6k+\gamma_k}{a}$$
$$D_{4k+1} = \dfrac{1}{1+d_{4k+1}D_{4k}} \approx \dfrac{1}{1+\left(-1+\frac{6k}{a}\right)\delta_k}\approx \dfrac{1}{1-\delta_k}=w_k$$
よって、$\Delta_{4k+1} = C_{4k+1}D_{4k+1} \approx w_k\left(\dfrac{6k+\gamma}{a}\right),\: |\Delta_{4k+1}-1|\approx 1$
となります。

$n=4k+2$のとき、
$$C_{4k+2} = 1+\dfrac{d_{4k+2}}{C_{4k+1}}\approx 1 + \dfrac{\frac{2k+1}{a}}{\frac{6k+\gamma_k}{a}} = 1+\dfrac{2k+1}{6k+\gamma_k} = \dfrac{8k+\gamma_k+1}{6k+\gamma_k}$$
$$D_{4k+2} = \dfrac{1}{1+d_{4k+2}D_{4k+1}}\approx \dfrac{1}{1+\frac{2k+1}{a}w_k}\approx 1 - \dfrac{2k+1}{a}w_k$$
よって、$$\Delta_{4k+2} = C_{4k+2}D_{4k+2}\approx \dfrac{8k+\gamma_k+1}{6k+\gamma_k}\left(1 - \dfrac{2k+1}{a}w_k\right)\approx \dfrac{8k+\gamma_k+1}{6k+\gamma_k}$$
$|\Delta_{4k+2}-1|\approx \dfrac{2k+1}{6k+\gamma_k}$
です。

$n=4k+3$のとき、
$$C_{4k+3}= 1+\dfrac{d_{4k+3}}{C_{4k+2}} \approx 1 + \dfrac{-1 + \frac{6k+3}{a}}{\frac{8k+\gamma_k+1}{6k+\gamma_k}}\approx 1-\dfrac{6k+\gamma_k}{8k+\gamma_k + 1}=\dfrac{2k+1}{8k+\gamma_k+1}$$
$$D_{4k+3} = \dfrac{1}{1+d_{4k+3}D_{4k+2}} \approx \dfrac{1}{1+\left(-1+\frac{6k+3}{a}\right)\left(1-\frac{2k+1}{a}w_k\right))}\approx \dfrac{1}{1-1+\frac{2k+1}{a}w_k+\frac{6k+3}{a}}=\dfrac{a}{(2k+1)(w_k+3)}$$
よって、
$$\Delta_{4k+3} = C_{4k+3}D_{4k+3} \approx \dfrac{2k+1}{8k+\gamma_k+1}\cdot \dfrac{a}{(2k+1)(w_k+3)} = \dfrac{a}{(8k+\gamma_k+1)(w_k+3)}$$
$$|\Delta_{4k+3}-1| \approx \dfrac{a}{(8k+\gamma_k+1)(w_k+3)}$$
です。
いよいよ$n=4k+4 = 4(k+1)$の場合です。
$$C_{4k+4} = 1 + \dfrac{d_{4k+4}}{C_{4k+3}} \approx 1 + \dfrac{\frac{2k+2}{a}}{\frac{2k+1}{8k+\gamma_k + 1}} = 1 + \dfrac{1}{a}\left[\dfrac{(2k+2)(8k+\gamma_k+1)}{2k+1}\right] = 1 + \dfrac{\gamma_{k+1}}{a}$$
$$D_{4k+4} = \dfrac{1}{1+d_{4k+4}D_{4k+3}} \approx \dfrac{1}{1 + \frac{2k+2}{a}\cdot\frac{a}{(2k+1)(w_k+3)}}= \dfrac{1}{1+\frac{2k+2}{(2k+1)(w_k+3)}} = \dfrac{(2k+1)(w_k+3)}{(2k+1)(w_k+3)+2k+2} = \delta_{k+1}$$
よって、$\Delta_{4k+4} \approx \dfrac{(2k+1)(w_k+3)}{(2k+1)(w_k+3)+2k+2},\: |\Delta_{4k+4}-1| \approx\dfrac{2k+2}{(2k+1)(w_k+3)+2k+2}$です。
ここで、$\gamma_k$$w_k$の一般項を求めましょう。
$\gamma_k$について、以下の漸化式が近似的に成り立つことが分かります。(一応等号で結ぶ)
$$(2k+1)\gamma_{k+1} = (2k+2)(\gamma_k + 8k+1)$$
$\gamma_k$は2次以下です。3次以上だと係数が合わないことが、3次の項を確かめることで分かります。
あとは$\gamma_k =ak^2 +bk+c$と置いて係数比較すれば解けます。
$$\begin{align*} (2k+1)[a(k+1)^2 + b(k+1)+c] &= (2k+2)[ak^2 +(b+8)k+c+1]\\ (2k+1)[ak^2 + (2a+b)k + a+b+c] &= (2k+2)[ak^2 +(b+8)k+c+1]\\ (3a-16)k^2+(4a-b-18)k+a+b+c-2 &=0 \end{align*}$$
これを解くと、$a=\dfrac{16}{3},\: b=-\dfrac{10}{3},\: c=0$となります。
これは$\gamma_1 = 2$を満たすので$\gamma_k = \dfrac{16k^2-10k}{3}$と分かります。

次に$w_k$について考えましょう。(これも近似的だが等号で結ぶ)
$$w_{k+1} = \dfrac{1}{1-\delta_{k+1}} = \dfrac{1}{1-\frac{(2k+1)(w_k+3)}{(2k+1)(w_k+3)+2k+2}} = \dfrac{(2k+1)(w_k+3)+2k+2}{2k+2} = \dfrac{2k+1}{2k+2}w_k + \dfrac{8k+5}{2k+2}$$
これも解きましょう。$w_k$が2次以下であることは$\gamma_k$と同様の方法で確かめられます。$w_k = ak +b$とすると、
$$\begin{align*} (2k+2)[a(k+1)+b] &= (2k+1)(ak+b)+8k+5\\ (3a-8)k + 2a+b-5&=0 \end{align*}$$
これを解くと、$a = \dfrac{8}{3},\: b = -\dfrac13$です。しかし、実際の$w_1 = \dfrac{1}{1-\delta_k} = 3$であって、$\dfrac83-\dfrac13=\dfrac73$ではありません。
ここで、斉次解を求めましょう。
$$\begin{align*} w_{k+1}^{\star} &= \dfrac{2k+1}{2k+2}w_k^{\star}\\ w_k^{\star} &= C\prod_{j=1}^{k-1}\dfrac{2j+1}{2j+2}\\ &= C\dfrac{(2k-1)!!}{(2k)!!} \end{align*}$$
よって一般解は$w_k = \dfrac{8k-1}{3}+C\dfrac{(2k-1)!!}{(2k)!!}$であって、初期値$w_1=3$から$C=\dfrac{4}{3}$が定まり、
$$w_k = \dfrac{8k-1}{3}+\dfrac43\dfrac{(2k-1)!!}{(2k)!!}$$
これを用いて$|\Delta_{4k+2}-1|$を計算すると、
$$|\Delta_{4k+2}-1| \approx \dfrac{2k+1}{6k+\gamma_k} = \dfrac{2k+1}{6k + \frac{16k^2-10k}{3}}= \dfrac{3(2k+1)}{16k^2+8k} = \dfrac{3(2k+1)}{8k(2k+1)}=\dfrac{3}{8k}$$
となります。修正Lentz法の収束条件は$|\Delta-1|<10^{-d}$であるので、
$$\begin{align*} \dfrac{3}{8k}\approx \dfrac{3}{2n} &< 10^{-d}\\ n &> \dfrac{3}{2}\cdot 10^{d} \end{align*}$$
と、確かに$1.5\times 10^{d}$項で収束した、と出ることが分かります。
では、これは真の収束を意味するのでしょうか。
そこで、積$P_k = \Delta_{4k}\Delta_{4k+1}\Delta_{4k+2}\Delta_{4k+3}$を考えます。
$$\begin{align*} P_k &= \left(\dfrac{w_k-1}{w_k}\right)\cdot\left(\dfrac{6k+\gamma_k}{a}w_k\right)\cdot\left(\dfrac{8k+\gamma_k+1}{6k+\gamma_k}\right)\cdot\left(\dfrac{a}{(8k+\gamma_k+1)(w_k+3)}\right)\\ &=\dfrac{w_k-1}{w_k+3}\\ &\approx \dfrac{2k-1}{2k+2} \end{align*}$$
これを$k=1,2,\cdots,K$まで掛けます。
$$\begin{align*} \prod_{k=1}^K P_k &= 2\dfrac{(2K-1)!!}{(2K+2)!!}\\ &\approx \dfrac{1}{\sqrt{\pi}K^{3/2}} \end{align*}$$
連分数では、この積は0ではない定数に近づくことを期待しますが、$a\to\infty, K\to\infty$でこの積は0に収束します。つまり、収束項数$n$$a$$b$に依存せずに決まることはない、という事です。残念...。

ではどうすればよいのか

実は、この連分数をそのままの形で運用し、修正Lentz法で計算させることは全く得策ではありません
修正Lentz法の収束判定は、連分数の全ての項が正であることを仮定しています。その時に、極めて安全に機能するのです。
今回の場合、すでに$d_{2m+1}$の項が負で固定されているため、値が大きく揺れるのです。
ではどうすればよいのか、2項を1項に纏めて処理しましょう。
そうした場合の連分数はこうなります。
$$\dfrac{1}{\beta_0 + \dfrac{\alpha_1}{\beta_1+\dfrac{\alpha_2}{\beta_2 + \dfrac{\alpha_3}{\beta_3+\cdots}}}}$$
$$\alpha_m = \dfrac{(a+m-1)(a+b+m-1)(b-m)m}{(a+2m-1)^2}x^2$$
$$\beta_m = a + 2m + \left(\dfrac{m(b-m)}{a+2m-1}-\dfrac{(a+m)(a+b+m)}{a+2m+1}\right)x$$
(ごちゃごちゃするので、今$m$のインデックスをずらすことはせず、最後の収束項数$m$に1を足します。)
これは、$m< b$の範囲で正項連分数となるので安全に扱うことが出来るし、なにより正項連分数を仮定した収束定理を用いることが出来ます。

正項連分数の収束性

$m$項目まで計算したときの連分数の分子を$p_m$, 分母を$q_m$とする。
$\beta_m,\alpha_m \geq 0$なる連分数が収束することと、
$$\lim_{m\to\infty}\dfrac{1}{q_mq_{m-1}}\prod_{k=0}^{m-1}\alpha_k = 0$$
が成り立つことは同値であり、収束値を$Z$とすると
$$\dfrac{\beta_{m+1}}{q_mq_{m+2}}\prod_{k=0}^{m}\alpha_k\leq \left|Z-\dfrac{p_m}{q_m}\right|\leq \dfrac{1}{q_mq_{m+1}}\prod_{k=0}^{m}\alpha_k$$
が成り立つ。

ですね。
今度こそ正確な収束項数推定を行いましょう。$a,b\gg m$を仮定します。
$$\dfrac{1}{\beta_0 + \dfrac{\alpha_1}{\beta_1+\dfrac{\alpha_2}{\beta_2 + \dfrac{\alpha_3}{\beta_3+\cdots}}}}$$
$$\alpha_m = \dfrac{(a+m-1)(a+b+m-1)(b-m)m}{(a+2m-1)^2}\left(\dfrac{a+1}{a+b+2}\right)^2 \approx am\dfrac{t}{t+1}$$
$$\beta_m = a + 2m + \left(\dfrac{m(b-m)}{a+2m-1}-\dfrac{(a+m)(a+b+m)}{a+2m+1}\right)x\approx 2m\dfrac{2t+1}{t+1}$$
つまり、これの分母漸化式はこうなります。
$$Y_{m+1} = 2m\dfrac{2t+1}{t+1}Y_m + am\dfrac{t}{t+1}Y_{m-1}$$
$Y_m = e^{S_m}$とし、$Y_{m\pm 1} = e^{S_m \pm \Delta S_m}$とおくと、
$$\begin{align*} e^{S_m+\Delta S_m} &= 2m\dfrac{2t+1}{t+1}e^{S_m}+am\dfrac{t}{t+1}e^{S_m-\Delta S_m}\\ (t+1)e^{2\Delta S_m} &= 2m(2t+1)e^{\Delta S_m} +amt \end{align*}$$
ここで、$e^{\Delta S_m} = \lambda_m$とおくと、
$$(t+1)\lambda_m^2 = 2m(2t+1)\lambda_m +amt$$
$$(t+1)\lambda_m^2 - 2m(2t+1)\lambda_m -amt = 0$$
という二次方程式が出てきます。
$\lambda_m ^{\pm} = \dfrac{m(2t+1)\pm\sqrt{m^2(2t+1)^2 + amt(t+1)}}{t+1}$
という解が得られました。具体的な解の表示はあとで使います。
この解は局所的な隣接項の比を表します。つまり、$\dfrac{Y_{m+1}}{Y_m}\approx \lambda_m^{\pm}$と表せるわけです。これの値が大きい方が優越解、小さい方が劣位解といいます。
一般的には優越解$(\lambda_m^{+})$が優先されるので、
$$Y_m = A_m\prod_{k=0}^{m-1}\lambda_m^{+}$$
と置いてみましょう。また、$A_{m\pm 1} = A_m \pm\Delta A_m$としておきます。
$$(A_m + \Delta A_m)\prod_{k=0}^{m}\lambda_k^{+}= 2m\dfrac{2t+1}{t+1}A_m\prod_{k=0}^{m-1}\lambda_k^{+} +am\dfrac{t}{t+1} (A_m - \Delta A_m)\prod_{k=0}^{m-2}\lambda_k^{+}$$
$$(A_m + \Delta A_m)\lambda_m^{+}\lambda_{m-1}^{+}= 2m\dfrac{2t+1}{t+1}A_m\lambda_{m-1}^{+} +am\dfrac{t}{t+1} (A_m - \Delta A_m)$$
ここで、解と係数の関係より、
$$(A_m + \Delta A_m)\lambda_m^{+}\lambda_{m-1}^{+}= (\lambda_m^{+}+\lambda_m^{-})A_m\lambda_{m-1}^{+} -\lambda_m^{+}\lambda_m^{-} (A_m - \Delta A_m)$$
となります。これを変形すると、
$$\begin{align*}\dfrac{\Delta A_m}{A_m} &= -\dfrac{\lambda_m^{-}}{\lambda_m^{+}(\lambda_{m-1}^{+}-\lambda_m^{-})}(\lambda_m^{+}-\lambda_{m-1}^{+})\\ &\approx -\dfrac{\lambda_m^{-}}{\lambda_m^{+}(\lambda_m^{+}-\lambda_m^{-})}(\lambda_m^{+}-\lambda_{m-1}^{+})\end{align*}$$
ここで、
$$\begin{align*} \lambda_m ^{\pm} &= \dfrac{m(2t+1)\pm\sqrt{m^2(2t+1)^2 + amt(t+1)}}{t+1}\\ &\approx \pm\sqrt{am}\\ & = \pm R_m \end{align*}$$
と近似します。$(R_m = \sqrt{am})$
すると、
$$\begin{align*}\dfrac{\Delta A_m}{A_m} &\approx -\dfrac{\lambda_m^{-}}{\lambda_m^{+}(\lambda_m^{+}-\lambda_m^{-})}(\lambda_m^{+}-\lambda_{m-1}^{+})\\ &\approx -\dfrac{-R_m}{2R_m^2}\Delta R_m\\&= \dfrac{\Delta R_m}{2R_m}\end{align*}$$
となります。これを連続的なものだととらえて積分すると、
$$\begin{align*}\ln A_m &\approx \dfrac12 \ln R_m + \ln C\end{align*}$$
となります。つまり、$A_m \approx C\sqrt{R_m}$です。
これを誤差上界に突っ込むと、
$$\begin{align*} \dfrac{1}{Y_mY_{m+1}}\prod_{k=0}^m \alpha_k &\approx \dfrac{\displaystyle\prod_{k=0}^m \lambda_k^{+}\left|\lambda_k^{-}\right|}{\displaystyle C^2\sqrt{R_mR_{m+1}} \prod_{k=0}^{m-1}\lambda_{k}^{+}\prod_{k=0}^{m}\lambda_k^{+}}\\ &= \dfrac{|\lambda_m^{-}|}{C^2\sqrt{R_mR_{m+1}}}\prod_{k=0}^{m-1}\left|\dfrac{\lambda_k^{-}}{\lambda_k^{+}}\right| \end{align*}$$
ここで、$\rho_k = \left|\dfrac{\lambda_k^{-}}{\lambda_k^{+}}\right| = -\dfrac{\lambda_k^{-}}{\lambda_k^{+}}$について、
$$\rho_k = \dfrac{\sqrt{\beta_k^2 -4\alpha_k}-\beta_k}{\sqrt{\beta_k^2 -4\alpha_k}+\beta_k} = \dfrac{\sqrt{1+u_k^2}-u_k}{\sqrt{1+u_k^2}+u_k}\approx 1-2u_k$$
ただし、
$u_k = \dfrac{\beta_k}{2\sqrt{\alpha_k}} = \dfrac{2k\dfrac{2t+1}{t+1}}{2\sqrt{ak\dfrac{t}{t+1}}} = \sqrt{\dfrac{k}{a}}\dfrac{2t+1}{\sqrt{t(t+1)}}$
です。
$$\ln \rho_k \approx -2u_k = -2\sqrt{\dfrac{k}{a}}\dfrac{2t+1}{\sqrt{t(t+1)}}$$
$$\prod_{k=0}^{m-1}\rho_k = \exp \left(\sum_{k=0}^{m-1}\ln \rho_k\right) \approx \exp \left(\int_0^{m}\ln \rho_k\: \text{d}k\right)$$
より、
$$\int_0^{m}\ln \rho_k\: \text{d}k \approx -\dfrac{4}{3}\dfrac{m^{3/2}}{\sqrt{a}}\dfrac{2t+1}{\sqrt{t(t+1)}}$$
よって、誤差上界は、
$$\dfrac{|\lambda_m^{-}|}{C^2\sqrt{R_mR_{m-1}}}\cdot \exp\left(-\dfrac{4}{3}\dfrac{m^{3/2}}{\sqrt{a}}\dfrac{2t+1}{\sqrt{t(t+1)}}\right)$$
と近似されます。更に、$|\lambda_m^{-}|\approx R_m\approx R_{m-1}$なので$\dfrac{|\lambda_m^{-}|}{\sqrt{R_mR_{m-1}}}\approx 1$
と近似できます。これが$10^{-d}$未満になればよいので、
$$\dfrac{1}{C^2}\exp\left(-\dfrac{4}{3}\dfrac{m^{3/2}}{\sqrt{a}}\dfrac{2t+1}{\sqrt{t(t+1)}}\right) < 10^{-d}$$
を解けばよいわけです。解くと、
$$m > \left[\dfrac{3\sqrt{a}}{4}\dfrac{\sqrt{t(t+1)}}{2t+1}(d\ln 10 - 2\ln C)\right]^{2/3}$$
となります。$t = \dfrac{b}{a}$を用いて戻すと、
$$m > \left[ \dfrac{3\sqrt{ab(a+b)}}{4(a+2b)} \left(d\ln10 - 2\ln C\right) \right]^{2/3}$$
$m$のインデックスとしてはこれでOKですが、総項数としては最初に$\dfrac{1}{b_0+\cdots}$とあるので、$m+1$項としておきます。(勿論これによってオーダーが変わることはありません。)
これによると、収束項数は$O (\sqrt[3]{\min(a,b)})$という挙動を取ります。

$a\gg b$では、$m > \left[ \dfrac{3\sqrt{b}}{4} \left(d\ln10 - 2\ln C\right) \right]^{2/3} \sim O (b^{1/3})$
$b\gg a$では、$m > \left[ \dfrac{3\sqrt{a}}{8} \left(d\ln10 - 2\ln C\right) \right]^{2/3} \sim O (a^{1/3})$
といった感じですね。
一方、 Numerical Recipes in C++ (第3版) には以下のように書かれています。

This continued fraction converges rapidly for $x < (a+1)/(a+b+2)$, except when $a$ and $b$ are both large, when it can take $O(\sqrt{\min(a,b)})$ iterations.

実は、収束項数はこれよりも良い上界だったということです。(このほかにも文献を当たりましたが、これ以外に収束のオーダーについて書かれているものを見つけられませんでした。)
...しかし待ってください。この項数にはある仮定をしていることを忘れていました。$m< b$という仮定です。そうでなければ正項連分数という仮定が崩れます。

しかしこれには回避方法があって、具体的には
$$b > \left[ \dfrac{3\sqrt{ab(a+b)}}{4(a+2b)} \left(d\ln10 - 2\ln C\right) \right]^{2/3}$$
を満たせばよいことになります。変形すると、
$$\begin{align*} \left(d\ln10 - 2\ln C\right) & < \dfrac{4b}{3}\dfrac{a+2b}{\sqrt{a(a+b)}}\\ b &> \dfrac34 (d\ln 10 -2\ln C) \end{align*}$$
を満たす$b$まで繰り上げればよいことになります。(漸化式があります。 https://dlmf.nist.gov/8.17#E21  を参照してください。)
で、$C$はまたまた$\dfrac12$であることが数値実験的に分かっているので、
$$m > \left[ \dfrac{3\sqrt{ab(a+b)}}{4(a+2b)} \left(d\ln10 + 2\ln 2\right) \right]^{2/3}$$
です。$a,b$が小さい時でも、繰り上げたり、収束項数に定数を足したりすることで対処できるので、これで$a,b$全域で解決したことになります。

$x = s(a+1)/(a+b+1), \: (0< s<1)$のとき

このパターンを忘れていました。
$\alpha_m \approx ams^2\dfrac{t}{t+1}$
$$\begin{align*} \beta_m &= a+2m +s\dfrac{a+1}{a+b+2}\left(\dfrac{m(b-m)}{a+2m-1}-\dfrac{(a+m)(a+b+m)}{a+2m+1}\right)\\ &\approx a+2m + \dfrac{s}{t+1}\left(2tm -a(t+1)+t+1\right)\\ &= a+ 2m +2m\dfrac{s}{t+1}-as+s\\ &= a(1-s)+2m\left(1+ \dfrac{st}{t+1}\right)+s\\ &\approx a(1-s) \end{align*} $$
$$\lambda_m^{\pm} = \dfrac{\beta_m}{2}\pm\underbrace{\sqrt{\dfrac{\beta_m^2}{4}+\alpha_m}}_ {R_m}$$
$$R_m \approx \dfrac{\beta_m}{2}+\dfrac{\alpha_m}{\beta_m}\approx \dfrac{a}{2}(1-s)$$
$$\lambda_m^{+} \approx \beta_m\approx a(1-s)$$
$$\lambda_m^{-}\approx -\dfrac{\alpha_m}{\beta_m}\approx-\dfrac{ms^2t}{(1-s)(t+1)} $$
$$\left|\dfrac{\lambda_m^{-}}{\lambda_m^{+}}\right| \approx \dfrac{ms^2t}{a(1-s)^2(t+1)} = \dfrac{W_s}{a}m$$
$$\left(W_s := \dfrac{s^2t}{(1-s)^2(t+1)}\right)$$
これを
$$\dfrac{|\lambda_m^{-}|}{C^2\sqrt{R_mR_{m+1}}}\prod_{k=1}^{m-1}\left|\dfrac{\lambda_k^{-}}{\lambda_k^{+}}\right| \approx \dfrac{|\lambda_m^{-}|}{C^2R_m}\prod_{k=1}^{m-1}\left|\dfrac{\lambda_k^{-}}{\lambda_k^{+}}\right|$$
に代入すると、
$$\begin{align*}\dfrac{2ms^2t}{C^2a(1-s)^2(t+1)}\prod_{k=1}^{m-1}\dfrac{W_s}{a}k&= \dfrac{2W_s^m}{C^2a^{m}}m!\\ &\approx\dfrac{2W_s^m}{C^2a^{m}}\sqrt{2\pi m}\left(\dfrac{m}{e}\right)^m\\ &= \dfrac{2}{C^2}\sqrt{2\pi m}\left(\dfrac{mW_s}{ae}\right)^m \end{align*}$$
これが$10^{-d}$未満であればよいので、
$$\begin{align*}\dfrac{2}{C^2}\sqrt{2\pi m}\left(\dfrac{mW_s}{ae}\right)^m &< 10^{-d}\\ \dfrac{C^2}{2}\left(\dfrac{ae}{mW_s}\right)^m&> \sqrt{2\pi m}10^d\\ m\ln\left(\dfrac{ae}{mW_s}\right)-\dfrac{1}{2}\ln (m) &> \dfrac{1}{2}\ln (2\pi)-2\ln (C)+\ln(2)+d\ln(10)\\ m\ln\left(\dfrac{ae}{mW_s}\right)-\dfrac{1}{2}\ln (m) &> \dfrac{1}{2}\ln (2\pi)+3\ln (2)+d\ln(10) \quad (\because C=\dfrac12)\end{align*} $$
右辺は定数ですので$B$と置きます。また、$\dfrac{ae}{W_s} = A$と置きます。すると、
$$m\ln \left(\dfrac{A}{m}\right)-\dfrac12 \ln (m) > B$$
これを解けばよいことになります。が、これはちょっと解きづらいので、$m$の条件は厳しくなりますが、
$$m\ln \left(\dfrac{A}{m}\right)-\dfrac{m}2 \ln (m) > B$$
を解きます。私は今十分性のある項数推定式が欲しいので、こうした方が都合が良いのです。この解はLambertのW関数$W_{-1}$を用いて、
$m > -\dfrac{2B}{3W_{-1}\left(-\dfrac{2B}{3A^{2/3}}\right)}$
と求められます。(保険用に$+\text{const}.$項しておくとよいでしょう。)
ここで、$a$がとてつもなく大きいとき、$W_{-1}$の定義外になってしまうのではないかと思う方もいるでしょう。その場合は、$B<0$ならば$m=2$、そうでないならば$x=\dfrac{a+1}{a+b+2}$の場合と同じ項数までには必ず収束するでしょう。

結論

$x = \dfrac{a+1}{a+b+2}$$a,b\gg 1$において、不完全ベータ関数を、連分数を用いて誤差$10^{-d}$未満になるまで計算するためには、偶数縮約すれば
$$\left[ \dfrac{3\sqrt{ab(a+b)}}{4(a+2b)} \left(d\ln10 + 2\ln 2\right) \right]^{2/3} + \text{const.}$$
項で計算できる。
ただし、$b$が小さい時には、$b > \dfrac34 (d\ln 10 -2\ln C)$まで$b$を繰り上げればよい。
また、$x = s\dfrac{a+1}{a+b+2},\: (s<1)$の場合は、
$$A = \dfrac{ae}{W_s}$$
$$B = \dfrac{1}{2}\ln (2\pi)+3\ln (2)+d\ln(10)$$
として、
$m > -\dfrac{2B}{3W_{-1}\left(-\dfrac{2B}{3A^{2/3}}\right)}+\text{const.}$
項で求められる。なお、定義外では$B<0$ならば$m=2$、そうでないならば高々$x=\dfrac{a+1}{a+b+2}$の場合と同じ項数までで収束する

いかがでしたか?

今回は不完全ベータ関数の連分数の計算・収束項数についてより詳細に解析しました。
間違いや疑問等あれば教えてください。
最後に、数値計算や確認等に協力してくれた 数値計算bot さん(a.k.a. Y.K. さん)には改めて感謝申し上げます。

参考文献

投稿日:4日前
更新日:8時間前
数学の力で現場を変える アルゴリズムエンジニア募集 - Mathlog served by OptHub

この記事を高評価した人

高評価したユーザはいません

この記事に送られたバッジ

バッジはありません。

投稿者

vunu
vunu
70
8190
「西」の「西東南」の部分

コメント

他の人のコメント

コメントはありません。
読み込み中...
読み込み中