如何快速计算Gamma函数和BesselJ函数？

dianyancao

\[
\begin{align*}
&\text{Gamma function: } \\
&-\ln [\Gamma (z)]=\ln z+\gamma z+\sum_{n=1}^{\infty}[\ln(1+\frac{z}{n})-\frac{z}{n}]\\
\\
&\text{BesselJ function:}\\
&J_m(x)=\begin{cases}
\sum_{i=0}^{\infty}{\frac{(-1)^l}{2^{2l+|m|}l!(|m|+l)!}x^{2l+|m|}}, & \text{ for } |m|\neq \frac{1}{2}\\
\sqrt{\frac{2}{\pi x}}\cos x, & \text{ for } m=-\frac{1}{2}\\
\sqrt{\frac{2}{\pi x}}\sin x, & \text{ for } m=\frac{1}{2}
\end{cases}
\end{align*}
\]
http://mathworld.wolfram.com/GammaFunction.html
http://mathworld.wolfram.com/BesselFunctionoftheFirstKind.html
mathematica中有Gamma和BesselJ函数，可用上述两个级数计算
有没有比这两个级数收敛更快的？

kastin

Mahtematica当然不会傻到直接用定义式计算，一般来说，超越函数的定义式收敛速度都会很慢的。
事实上，Mathematica里面计算各种超越函数，以及数论函数，组合函数，方程求根……都是采用了目前最优的算法来实现的。

比如Gamma函数，根据输入的数据类型和大小的不同（表达式还是数字，数字的大小范围），动态地采用递归、函数方程和 Binet 渐进公式计算。

同样，贝塞尔系列函数（J,Y,K,I）使用级数和渐进展开式求值； 对于整数次阶，有些也使用稳定前向递归。这个算法流程就更加复杂了——因为根据阶涉及到整数还是小数，输入变量会是表达式还是实数、复数以及范围大小——计算策略也会相应不同，尤其是对于特殊值的情况，会用特殊的算法来计算（这样速度更快）。

举个例子。比如计算Gamma函数，输入的是纯数值。那么根据http://mathworld.wolfram.com/GammaFunction.html中的公式59和61可知，只要计算出区间(0,1/2]上的gamma函数的值，则可以得到任何实数作为自变量的gamma函数值。而(0,1/2]区间较小，因此可以通过在该区间收敛的无穷级数来计算。另外，如果还想增加加速效果，可以通过级数加速收敛技术(比如欧拉转换公式，艾特金`\Delta^2`加速法等等)来使得收敛速度加快。

详细细节参见Mathemaitca参考资料中心网页上的：关于内部实现的一些注释

dianyancao

kastin 发表于 2014-9-3 15:08
Mahtematica当然不会傻到直接用定义式计算，一般来说，超越函数的定义式收敛速度都会很慢的。
事实上，Mat ...

Gamma函数可以通过递归，转换成(0,1]之间的Gamma值计算
然后在通过余元公式能转换成[-1/2,0)或(0,1/2]区间的计算？

用这些级数加速收敛技术会不会增加额外的计算开销?

282842712474

还有一种技巧是外推，利用已经得到的数字再进行一次计算，可以得到更多位有效数字，缺点是一般进行了一次外推之后，就没有办法继续后续过程了，也就是说外推是一个最终步骤~
http://spaces.ac.cn/index.php/archives/1292/

kastin

去年在查阅高振荡数值积分资料的时候，阅读了一下序列加速技术相关论文，这里简要介绍一下吧。

在数值分析中，高振荡数值积分、级数求和、方程迭代求根等过程中经常会遇到收敛较慢的问题。
【注意1：为了方便起见，我们只讨论标量序列情形，向量序列（多应用于高振荡向量函数数值积分、矩阵或向量级数求和、方程组迭代求根）稍微复杂一点，但是大体跟标量情形类似。】
【注意2：除非特别说明，下面出现的符号 `n` 默认取自然数（包括零）。】

这里有必要解释一下在高振荡数值积分中的应用。在很多情况下，通常的数值积分方法对于高振荡函数（比如核函数是高频三角函数、贝塞尔函数等）的积分都会失效（误差大，而且结果不稳定），这是由于高振荡核需要更小的积分步长才能使积分值稍微精细，但若积分步长太小，舍入误差就过大，并且积分值提高的精度并不增加很多。目前，已经发展出很多优秀的方法都解决这一类问题，而其中有一种方法是：按照被积函数在积分区间内的零点，分成许多个小的子区间（相邻零点确定的区间），然后分别积分，再将积分值求和（如果是在零上下振荡，就是是交错级数求和）——这相当于一个级数求和过程。但通常高振荡函数的振幅部分衰减太慢，因此求和收敛过程变得极其困难，这就需要应用加速收敛方法了。

本质上来说，级数求和 `S=\sum a_n` 的收敛过程可以看成是级数的部分和序列 `{S_n}` 的收敛。

上面的问题都可以抽象为求一个收敛速度较慢的序列 `\{s_n\}` 的极限。首先的问题是，存在加快其收敛速度的方法吗？答案是肯定的。

继续介绍之前，我们先给出收敛速度快慢的定义。

直觉上来说，某个序列比另一个序列收敛得更快的意思是，在同样极限过程中，该序列更快地接近对应的极限值。

如果有两个序列 `\{s_n\}` 和 `\{s_n'\}` 有相同的极限 `l`，我们称 `s_n'` 比 `s_n` 收敛得更快，若$$\lim_{n\to\oo}\frac{s_n'-l}{s_n-l}=0\tag{1}$$换句话说，新的序列的残差项（各项减去极限值）是原序列残差项的高阶无穷小。

定义`~~~~`设序列有极限 `\D\lim_{n\to \oo} a_n=l` （这里的 `l` 可以是无穷大量）。称 `a` 是 `p(>0)` 阶收敛于 `l` 的，当$$\lim_{n\to\oo}\frac{|a_{n+1}-l|}{|a_n-l|^p}=c\tag{2}$$其中 `c>0`.

特别地，若 `p=1`，且存在 `\mu\in(0,1)` 满足 `\D\lim_{n\to \oo} \frac{|a_{n+1}-l|}{|a_n-l|}=\mu`，我们称 `a_n` 是线性收敛的。
1. 若 `\mu=0`，则称 `a_n` 是超线性（superlinear）收敛的；
2. 若 `\mu=1`，则称 `a_n` 是亚线性（sublinear）收敛的。若此时还有`\D\lim_{n\to \oo} \frac{|a_{n+2}-a_{n+1}|}{|a_{n+1}-a_n|}=1`，则称 `a_n` 是对数速率收敛的。

经过数学家和学者的研究，大量的加速收敛方法以及特殊情况的加速收敛手段广泛被提出。其中，两大经典的加速收敛方法分别是Euler级数变换法、Kummer级数变换法。剩下的方法中还有常用的（广义）Rechardson外推法极其衍生方法、Aitken `\Delta^2`加速法， Wynn `\varepsilon` 方法、Brezinski `\theta`方法，Levin u-变换、Lubkin W-变换、Shanks 变换等。如果求和部分有规律的话（即可以写成整数n的某个函数），Euler-Maclaurin公式也是一个不错的方法。

这些方法的特点是，使用某种手续或者步骤，将原来的序列变换为新的收敛速度很快的序列：`f\colon S \mapsto S'`. 这时，只需要少数几项就能获得很高精度的收敛值。

一、Euler级数变换
这是是一种线性变换方法，源于对交错级数求和的研究。
考虑交错级数，用移位算子`E` 以及前向差分算子 `\Delta` 可以写成$$\sum_{n=0}^{\oo}(-1)^na_n=(I-E+E^2-E^3+\cdots)a_0=(I+E)^{-1}a_0=(2I+\Delta)^{-1}a_0=\frac{1}{2}(I+\frac{1}{2}\Delta )^{-1}a_0=\sum_{n=0}^{\oo}(-1)^n\frac{\Delta^na_0}{2^{n+1}}\tag{3}$$(不熟悉差分算子的请看这里http://bbs.emath.ac.cn/forum.php?mod=redirect&goto=findpost&ptid=5736&pid=55095&fromuid=8865)
因此，等式 `(3)` 最左端与最右端连结起来就是欧拉级数变换公式，欧拉用它研究超几何级数的相关性质。

在某些特殊情况下，需要用到欧拉级数变换公式的变形形式。对于序列`s_n`，若其高阶差分较小，那么可以忽略相应高阶差分部分。Longman给出对于较大正整数 `p` 和 `m（p\leqslant m）`$$\sum_{n=0}^{\oo}(-1)^na_n \mapsto \sum_{n=0}^{p-1}(-1)^n\frac{\Delta^na_0}{2^{n+1}}+(-1)^m\sum_{n=0}^{p-1}(-1)^n\frac{\Delta^na_{m-n}}{2^{n+1}}\tag*{(*)}$$

很多时候，我们并不希望从 `a_0` 开始进行欧拉级数变换。因为通常前面几项的和占主导，而后面的项则确定小量部分。因此，假若从第 `p` 项开始变换，那么欧拉级数变换应该是：$$\sum_{n=0}^{\oo}(-1)^na_n = \sum_{n=0}^{p-1}(-1)^na_n+\sum_{n=p}^{\oo}\frac{\Delta^{n-p}a_p}{2^{n-p+1}}\tag{4}$$

当然，如果慢收敛级数不是交错级数，如何应用欧拉级数变换呢？我们可以将正项级数`\sum b_n`变换为交错级数：$$\sum_{n=0}^{\oo}b_n \mapsto \sum_{n=0}^{\oo}(-1)^na_n \tag{5}$$其中 `\D a_n=\sum_{k=0}^{\oo} 2^kb_{2^kn}`
注意，虽然每个`a_n`是无穷项`2^kb_{2^kn}`的和，但是只需要少数几项就能达到很好的精度，而不必全部加起来。

二、Kummer级数变换
这也是一种线性变换。但是思想一样很简单：设 `\sum a_n` 是收敛级数，而 `\sum c_n=c` 是一个和为 `c` 的收敛级数，且满足 `\D\lim_{n\to\oo}\frac{a_n}{c_n}=\lambda\neq 0`. 则$$\sum_{n=0}^{\oo} a_n=\lambda c+\sum_{n=0}^{\oo}(1-\lambda \frac{c_n}{a_n})a_n\tag{6}$$推导过程也很简单，由于两个级数都是收敛的，因此等式右边乘进去，然后可以分组求和，最后化简就得到左边的式子。

三、Aitkin delta平方加速
这个方法与下面将要提到的Richardson外推法本质上都算是外推类（Extrapolation）的方法。
假设已知序列 `s_n` 以几何速率收敛于极限 `s`$$s_n=s+(c+\varepsilon_n)q^n\quad(c\neq 0,\varepsilon\to 0)\tag{7}$$其中，`|q|<1`
因为不知道 `c,\varepsilon_n,q`，但我们知道 `s-s_{n+m}\approx cq^{n+m}\approx q^m(s-s_n)`，所以$$\frac{s-s_{n+m}}{s-s_n}\approx q^m \approx \frac{s-s_n}{s-s_{n-m}}$$整理得$$\begin{align*}s &=\frac{s_{n+m}s_{n-m}-s_n^2}{s_{n+m}-2s_n+s_{n-m}}\\
&=s_{n+m}-\frac{(s_{n+m}-s_n)^2}{s_{n+m}-2s_n+s_{n-m}}\end{align*}$$这说明，按照上式右端构造的新序列$$s_n'=s_{n+m}-\frac{(s_{n+m}-s_n)^2}{s_{n+m}-2s_n+s_{n-m}}\tag{8}$$比原序列能更快地接近于极限 `s`。
当取 `m=1`，上式用差分符号表示就是$$s_n'=s_{n+1}-\frac{(\Delta s_n)^2}{\Delta^2 s_{n-1}}\tag{9}$$这就是Aitkin `\Delta^2` 加速法。

四、Rechardson外推
假设序列 `s_n` 收敛的阶能被估计，即序列是按照下列渐进展开式收敛于极限 `s` $$s=s_n+\frac{c}{n^\alpha}+\frac{d}{n^\beta}+\frac{e}{n^\gamma}+\cdots\quad(0<\alpha<\beta<\gamma<\cdots)\tag{10}$$其中 `\alpha,\beta,\gamma,...` 都是已知的，而系数 `c,d,e` 是未知量。
因此上式可以写为$$s=s_n+\frac{c}{n^{\alpha}}+o(\frac{1}{n^{\beta}})\tag{11.1}$$类似地，给定正整数 `k`，通过不同的n之间外推至无穷远$$s=s_{kn}+\frac{c}{k^{\alpha}n^{\alpha}}+o(\frac{1}{k^{\beta}n^{\beta}})\tag{11.2}$$将 `k^{\alpha}\times(11.2)-(11.1)` 得$$s=\frac{k^{\alpha}s_{kn}-s_n}{k^{\alpha}-1}+o(\frac{1}{n^{\beta}})\tag{11.3}$$这就消去了较大的误差项，只剩下高阶小量 `o(\frac{1}{n^{\beta}})`. 因此新序列$$s_n'=\frac{k^{\alpha}s_{kn}-s_n}{k^{\alpha}-1}\tag{12}$$将比原序列更快地收敛到 `s`.
通常在 `(12)` 式中取 `k=2`，这就是我们常见到的Richardson外推法。

当然，还有其他形式的外推法，无非就是假设不同收敛速度的形式。比如把序列假设是按照`1/n`的多项式的阶收敛的，然后外推至无穷远，得到一簇线性公式，这便是Salzer加速法。其中一组典型公式是$$s_n'=-\frac{8}{6}s_n+\frac{81}{6}s_{n+1}-\frac{192}{6}s_{n+2}+\frac{125}{6}s_{n+3}\tag{13}$$又如，假设序列是按照对数收敛的 `s=s_n+\frac{c}{(\ln n)^{\alpha}}+\cdots`，混合对数形式的 `s=s_n+\frac{c}{n^{\alpha}\ln n}+\cdots` 或者 `s=s_n+\frac{c\ln n}{n^{\alpha}}+\cdots`，等等。不同的形式就能得到不同的外推加速收敛方法。

需要注意的是，序列如果是近似或者恰好是按照假设的形式收敛的，那么对应的方法效果会很好，否则效果不太理想。

五、其他方法

Wynn epsilon方法以及Brezinski theta方法类似于构造不同层次的表（Table），具体过程比较复杂，因此符号上下标很多这里就不介绍了，感兴趣的话可以自己查查相关论文。


下列资料可以提供更进一步的参考

加速收敛技术的回顾与总结
C.Brezinski. Convergence acceleration during the 20th century

总结非线性变换用于加速收敛以及发散级数求和
Ernst Joachim Weniger. Nonlinear sequence transformations for the acceleration of convergence and summation of divergent series

一本介绍实用外插技术以及加速收敛方法理论与应用的书
Avram Sidi. Practical Extrapolation Methods:Theorem and Applications

dianyancao

kastin 发表于 2014-9-4 13:13
去年在查阅高振荡数值积分资料的时候，阅读了一下序列加速技术相关论文，这里简要介绍一下吧。

在数值分 ...

\[\sum_{n=0}^{\oo}(-1)^na_n=\sum_{n=0}^{p-1}(-1)^n\frac{\Delta^na_0}{2^{n+1}}+(-1)^m\sum_{n=0}^{p-1}(-1)^n\frac{\Delta^na_{m-n}}{2^{n+1}}\tag*{(*)}
\]
上述这个式子看不懂，\(=\)应该是\(\to\),\(m\)应该等于\(p\)？求推导过程

\[\D a_n=\sum_{k=0}^{\oo} 2^kb_{2^kn}\]
上述这个级数好像不一定收敛？

kastin

（下面的差分算子只对其后的系数起作用）
构造形式幂级数（formal power series）$$S(x)=a_0-a_1x+a_2x^2+ax^3-\cdots \tag{a.1}$$等式两边同时乘以 `x`$$xS(x)=a_0x-a_1x^2+a_2x^3+\cdots\tag{a.2}$$然后两式相加$$(1+x)S(x)=a_0-x(\Delta a_1+\Delta a_2 x+\Delta a_3 x^2+\cdots)\tag{a.3}$$整理得$$S(x)=\frac{a_0}{1+x}-\frac{x}{1+x}(\Delta a_1+\Delta a_2 x+\Delta a_3 x^2+\cdots)\tag{a.4}$$对(a.4)式等号右端括号内的式子反复进行上述操作，最后令`x=1`即可得到（*）式，以及一项更高阶的差分表达式，而当m,p都是非常大的数且(p<=m)时，这个多出来的高阶差分表达式认为近似为0（形式幂级数函数收敛域有限）。所以（*）式成立。

或者利用差分算子进行不太严密的直观解释：因为为无穷项的交错级数`a_0-a_1+a_2-\cdots`可以用欧拉级数变换`\D(I+E)^{-1}a_0=\D\frac{1}{2}\small[I-(\frac{1}{2}\Delta)+I-(\frac{1}{2}\Delta)^2-\cdots]a_0=\D\sum_{n=0}^{\oo}(-1)^n\frac{\Delta^na_0}{2^{n+1}}`严格得到。因此若对于部分和，只要级数中的项足够多，则这种欧拉变换仍然近似成立（高阶差分视为小量）。给定足够大的正整数`m`，我们可以直接取前`m`项做欧拉级数变换。但是一般是取一个比`m`小得多的整数`p`，从1到`p`项应用欧拉级数变换，而剩下的`m-p+1`项可以视作另一个交错级数的部分和，由于`m`远大于`p`，所以仍然可以应用欧拉级数变换，为了弥补前面操作的误差，我们倒序排列剩下的项，然后做欧拉级数变换`\D(I+E)^{-1}a_p=(2E-E\nabla)^{-1}a_p=\frac{1}{2}E^{-1}(I+\frac{\nabla}{2})^{-1}a_p=\frac{1}{2}E^{-1}(I+\frac{1}{2}E\Delta)^{-1}a_p=\frac{(-1)^m}{2}E^{-m}(I+\frac{\Delta}{2})^{-1}a_m+高阶差分项 \approx (-1)^m\sum_{n=0}^{p-1}(-1)^n\frac{\Delta^na_{m-n}}{2^{n+1}}`
两部分相加即为(*)式——部分和的欧拉级数变换近似代替无穷项级数。

注意，后项差分算子`\nabla=I-E^{-1}`.

补充内容 (2014-9-6 23:56):
由于疏忽，根据差分算子的定义，(a.3)和(a.4)中等式右端括号中的项的下标都应该减小1.

dianyancao

kastin 发表于 2014-9-5 18:10
（下面的差分算子只对其后的系数起作用）
构造形式幂级数（formal power series）$$S(x)=a_0-a_1x+a_2x^2+ ...

\[S(x)=\frac{a_0}{1+x}-\frac{x}{1+x}(\Delta a_1+\Delta a_2 x+\Delta a_3 x^2+\cdots)\tag{a.4}\]
上式，\(a_0\)的权是\(\frac{1}{2}\)?和原来的\(S(x)\)不相等了啊？
我重复上式得到的结果是：
\[\begin{align*}
S(x)&=\sum_{i=0}^{n} ({\frac{x}{x-1}}{\Delta})^i \frac{1}{1-x} a_0\quad+\quad{(\frac{x}{x-1})}^n{\Delta}^n(\sum_{i=0}^{\infty}(-1)^i x^i a_i ) \\
&=\sum_{i=0}^{n}(-1)^i \frac{1}{1-x}{(\frac{x}{1-x})}^{i}{\Delta}^i a_0\quad+\quad{\sum_{i=0}^{\infty}(-1)^i \frac{x^i {\Delta}^n a_i}{(\frac{x-1}{x})^n}}
\end{align*}\]

kastin

dianyancao 发表于 2014-9-6 21:53
\[S(x)=\frac{a_0}{1+x}-\frac{x}{1+x}(\Delta a_1+\Delta a_2 x+\Delta a_3 x^2+\cdots)\tag{a.4}\]
...

取给定大数 `m` 和 `p` 满足 `p\leqslant m`，那么无穷级数可以用m项部分和代替，因此(a.4)可以写为$$S(x)=\frac{a_0+(-1)^m a_mx^{m+1}}{1+x}+\frac{x}{1+x}\small[(\Delta a_0)-(\Delta a_1)x+(\Delta a_2)x^2-\cdots+(-1)^{m-1}(\Delta a_{m-1})x^{m-1}]\tag{b}$$
对括号中反复进行 `p` 次上述操作，并令 `x=1` 得到$$S(x)=\frac{1}{2}a_0-\frac{1}{4}\Delta a_0+\frac{1}{8}\Delta^2 a_0-\cdots+(-1)^{p-1}\frac{\Delta^{p-1}}{2^p}a_0\\+(-1)^m\small[\frac{1}{2}a_m-\frac{1}{4}\Delta a_{m-1}+\frac{1}{8}\Delta^2 a_{m-2}-\cdots+(-1)^{p-1}\frac{\Delta^{p-1}}{2^p}a_{m-p+1}]\\
+(-\frac{1}{2})^p\small[\Delta^p a_0-\Delta^p a_1+\Delta^p a_2-\cdots+(-1)^{m-p}\Delta^p a_{m-p}]\tag{c}$$因为 `p` 很大，所以(c)中最后一行可以忽略。这就得到(*)式.