[QuestionじゃなくてProblemの方の問題] 佐藤超函数を用いた数値積分法について

前置き

先日、私は「 求積アルゴリズム(+Excelのワンセルコード) 」という記事を投稿したのですが、なるべく色々な求積アルゴリズムを集めたいという理由で、緒方秀教氏、平山弘氏の発表した「 数値積分に対する超函数法 」などを参考にExcelのワンセルコードを作ったのですが、一部の関数を除いて結果が著しく合わないので、Excelのコードとか自分の理解や解釈が間違っているのではないかというのを確認するためにこの記事を作成しました。値や関数、パラメータなどを色々変えて実験したのですが、結果が合わず、自分一人で考えても埒があかないと感じたためです。

本題

詳しくは「 数値積分に対する超函数法(日本応用数理学会論文誌 Vol. 26, No. 1, 2016, pp. 33～43) 」や「 佐藤超函数論に基づく数値積分法(日本数学会 2016 年年会) 」などを参照してほしいのですが、佐藤超函数を用いて、
$[a,b]$を含む領域$D\subset \mathbb{C}$において解析的である関数$f(x)$と重み関数$w(x)$に対して
$$\int_a^bf(x)w(x)\: dx$$
を考えるとき、この被積分関数に指示関数$\chi_{(a,b)}(x)$を掛けたものを超函数とみなして、その定義関数$[F(x)]$を求めると、
$f(x)w(x)\chi_{(a,b)}(x) = \left[-\dfrac{1}{2\pi i}f(z)\Psi (x)\right]$
$$\Psi (z) = \int_a^b\dfrac{w(x)}{z-x}\: dx$$
となる。あとはこれを積分して、複素領域$D$内に正の向きの閉積分路$C$をとり、周期関数を$\varphi(u)$、周期を$u_{\text{period}}$と置いて台形則で近似すると、
\begin{eqnarray} \int_a^bf(x)w(x)\: dx &=& \dfrac{1}{2\pi i}\oint_Cf(z)\Psi (z)\: dz\\ &\approx& \dfrac{h}{2\pi i}\sum_{l=0}^{n-1}f(\varphi(lh))\Psi (\varphi (lh))\varphi'(lh) \end{eqnarray}
$h = \dfrac{u_{\text{period}}}{n}$
となります。
重み関数が$w(x)=1$、積分範囲が$(a,b)$のとき、
$$\Psi (z) = \log\left(\dfrac{z-a}{z-b}\right)$$
($-\pi \leq \Im (\log z) < \pi$)
となることと、
参考文献による周期関数の取り方から
$$\varphi (u) = \dfrac{a+b}{2} + \dfrac{b-a}{4}\left(k+\dfrac{1}{k}\right)\cos u + i\dfrac{b-a}{4}\left(k-\dfrac{1}{k}\right)\sin u$$
$$\left(\varphi'(u) = -\dfrac{b-a}{4}\left(k+\dfrac{1}{k}\right)\sin u + i\dfrac{b-a}{4}\left(k-\dfrac{1}{k}\right)\cos u\right)$$
($u_{\text{period}}=2\pi$、$k$は正の定数[パラメータ])
として、Excelのワンセルコードを作成しました。
とりあえずテスト積分は
$$\int_0^\pi e^{\cos x}\: dx = \pi I_0(1)\approx 3.97746326050642$$

$k=2$、$n=80$に設定
結果:$3.97746326050642$
$k=2.5$、$n=80$に設定
結果:$3.97746326050642$
$k=3$、$n=150$に設定
結果:$3.97746326050643$
この関数は非常に正確な値を返します。(しかし$k\geq 5$とかにすると結構ずれ始めます。おそらく$\cos z$の値が大きくなって大幅な桁落ちが起こるのが理由なのかと思われますがこちらも原因は不明です。)

ではテスト関数を変えてみます。
$$\int_0^1\sqrt{x}\: dx = \dfrac{2}{3}$$
すると
$k=2$、$n=100$に設定
結果:$0.693034765599263$
$k=2.5$、$n=100$に設定
結果:$0.718322768711128$
$k=3$、$n=100$に設定
結果:$0.746029367111263$
となってどんどん真値と離れていきます。(勿論$k$を1に近づけると真値に近づいていきますが、$k$は$\varphi(u)$が先ほどの解析的である領域$D$内であれば問題ないのでその範囲内に収まる$k$によって値が大幅に変わることはありません。)

今度は参考文献に載っている数値積分例で行います。
$$\int_0^1x^{\alpha-1}(1-x)^{\beta-1}e^x\: dx \approx 3.71819 70362 84701\times 10^4$$
$(\alpha = 10^{-4}, \beta = 10^{-4})$
$k=10$、$n=50$に設定
結果:$-0.847653395793876$
合わないです。$k$の値を1に近づけると一応真値に近づきそうな気配はあります。(ただ正負すら合ってないっていうのは本当に謎です。)

終わりに

そもそも$\text{Excel}$で複素数をこういった感じで扱うということ自体かなり無謀という気がします。それに$\log z$や$z^\alpha$の演算が正しいか(正しい分枝を取っているか)というのも一応確認しましたが、それが全域で正しいと言えるかどうかはもっと調べる必要があると思います。
ただExcelで縛った以上このコードも完成させたいと思うので、どうかこの記述に初歩的なミスがあるかどうかでも確認して頂けると幸いです。
(DE公式が苦手な積分もしっかり積分したいので...。)