Shanks变换与Richardson外推法加速收敛级数求和

给定一个收敛数列,近似估计这个数列的收敛值至所要求的的精度。例如:

113+1517+19...=π41 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9} - ... = \frac{\pi}{4}

假设我们还不知道π\pi是多少,要近似到10910^{-9}精度则需要计算接近10910^{9}数量级这么多项,这是不能接受的。

有什么办法可以快速求这个收敛值的近似值呢?


Shanks变换(Shanks Transformation)

sn(A)s_n(A)表示序列AA的前nn项和。对于一个序列AA,现在想要快速求出其收敛值S=s+(A)=AiS = s_{+\infty}(A) = \sum A_i

构造序列BB,使得sn(B)=S+αqns_n(B) = S + \alpha q^n。其中0<q<10 \lt \vert q \vert \lt 1。我们可以通过联立方程组来获得这样一个解法:

{sn+1(B)=S+αqn+1sn(B)=S+αqnsn1(B)=S+αqn1(sn1(B)S)(sn+1(B)S)=(sn(B)S)2\begin{cases} s_{n+1}(B) = S + \alpha q^{n+1} \\ s_n(B) = S + \alpha q^n \\ s_{n-1}(B) = S + \alpha q^{n-1} \\ \end{cases} \\ \Rightarrow (s_{n-1}(B)-S)(s_{n+1}(B)-S)= (s_n(B)-S)^2

化简就得到与α,q\alpha,q无关的式子:

S=sn2(B)sn1(B)sn+1(B)2sn(B)sn1(B)sn+1(B)S = \frac{s_n^2(B)-s_{n-1}(B)s_{n+1}(B)}{2s_n(B)-s_{n-1}(B)-s_{n+1}(B)}

注意目前为止这个等式对于任何nn都严格成立(在极限意义下)而非一个近似估计。

但是毕竟我们不知道序列BB的表达式,不能直接拿来计算(如果能的话不就做完了)。于是接下来,Shanks变换假设sn(B)s_n(B)可以近似地等于sn(A)s_n(A)

S^1,n=sn2(A)sn1(A)sn+1(A)2sn(A)sn1(A)sn+1(A)\hat S_{1,n} = \frac{s_n^2(A)-s_{n-1}(A)s_{n+1}(A)}{2s_n(A)-s_{n-1}(A)-s_{n+1}(A)}

只要代入sn1(A),sn(A)s_{n-1}(A),s_n(A)sn+1(A)s_{n+1}(A)的函数值,我们就得到了一个关于SS的估计值S^1,n\hat S_{1,n}。特别地,如果分母是00就认为S^1,n=sn(A)\hat S_{1,n} = s_n(A)

显然在nn越大的时候由于sn(A)s_n(A)sn(B)s_n(B)都越来越接近最后的答案SSSS就越接近S^1,n\hat S_{1,n},估计越准确。这里构造αqn\alpha q^n作为一个附加项,就是因为αqn\alpha q^n衰减得很快,可以让序列BB在尽可能小的nn就与序列AA接近,使得将sn(B)s_n(B)替换为sn(A)s_n(A)的误差较小,而且最终limn+sn(B)=S\lim_{n \rightarrow +\infty}s_{n}(B) = S

事实上,Shanks变换的这个近似假设往往是比较准确的。如果我们接受Shanks变换比普通迭代更有效的事实,那么就可以嵌套Shanks:

S^i+1,n=S^i,n2S^i,n1S^i,n+12S^i,nS^i,n1S^i,n+1\hat S_{i+1,n} = \frac{\hat S_{i,n}^2-\hat S_{i,n-1}\hat S_{i,n+1}}{2\hat S_{i,n}-\hat S_{i,n-1}-\hat S_{i,n+1}}

代码实测一下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include <bits/stdc++.h>
using namespace std;
#define EPS 1e-15
#define MAXN 100000
#define db long double
#define rint register int
inline int rf(){int r;int s=0,c;for(;!isdigit(c=getchar());s=c);for(r=c^48;isdigit(c=getchar());(r*=10)+=c^48);return s^45?r:-r;}
db A[16][MAXN+5], _; int n, K; inline db Shanks(db a, db b, db c){return _=2*b-a-c,fabs(_)>EPS?(b*b-a*c)/_:b;}
int main(){
n = rf(); K = rf(); A[K][0] = 1; printf("pi = %.10Lf\n",(db)3.1415926535897932384626433832795);
for(rint i = 1; i <= n+K; A[K][i] = A[K][i-1]+(db)(i&1?-1.:1.)/(i<<1|1), i++);
for(rint j = K-1, i; ~j; j--) for(i = 1; i <= n+j; i++)
A[j][i] = Shanks(A[j+1][i-1],A[j+1][i],A[j+1][i+1]);
printf(" "); for(rint j = 0; j <= K; printf("Shanks %2d: ",j), j++); puts("");
for(rint i = 1, j; i <= n; i++)
for(printf("n = %3d: ",i), j = K; ~j; printf("%.10Lf%c",4*A[j][i],"\n "[!!j]), j--);
return 0;
}

下面是实测结果:

可以发现精度是非常高的,S^6,25\hat S_{6,25}已经达到了精度要求,完全不需要计算10910^9项。


Richardson外推法(Richardson Extrapolation)

Richardson外推法是另外一种加速级数求和的方式。Shanks变换一般对于交错数列(AiA_i每项一正一负)效果更好,而Richardson外推法对于单调数列效果更好。

Richardson外推法思路和Shanks变换几乎完全一样。不过,Shanks变换构造sn(B)=S+αqns_n(B) = S + \alpha q^n,而Richardson外推法则构造:

sn(B)=S+i=1+Ainis_n(B) = S + \sum_{i=1}^{+\infty}\frac{A_i}{n^i}

先只计算第一项,联立方程组:

{sn+1(B)=S+A1n+1sn(B)=S+A1n\begin{cases} s_{n+1}(B) = S + \frac{A_1}{n+1} \\ s_n(B) = S + \frac{A_1}{n} \\ \end{cases} \\

n+1n+1倍的①式减去nn倍的②式,得到的结果:

S=(n+1)sn+1(B)nsn(B)S = (n+1)s_{n+1}(B) - ns_n(B)

同样地,目前为止严格精确。因为同样满足limn+sn(B)=S\lim_{n \rightarrow +\infty}s_n(B) = S,所以同样地,Richardson外推法假设sn(B)s_n(B)可以近似地等于sn(A)s_n(A)

R^1,n=(n+1)sn+1(A)nsn(A)\hat R_{1,n} = (n+1)s_{n+1}(A) - ns_n(A)

这称为11阶Richardson外推法。

与Shanks变换不同,高阶Richardson外推法不是通过迭代这一过程而是通过增加展开的项数。在之前的Richardson外推法假设公式中,我们计算两项:

{sn+2(B)=S+A1n+2+A2(n+2)2sn+1(B)=S+A1n+1+A2(n+1)2sn(B)=S+A1n+A2n2\begin{cases} s_{n+2}(B) = S + \frac{A_1}{n+2} + \frac{A_2}{(n+2)^2} \\ s_{n+1}(B) = S + \frac{A_1}{n+1} + \frac{A_2}{(n+1)^2} \\ s_n(B) = S + \frac{A_1}{n} + \frac{A_2}{n^2} \\ \end{cases} \\

计算(n+2)2(n+2)^22(n+1)2-2(n+1)^2+n2+n^2③,化简后近似可以得到22阶Richardson:

R^2,n=(n+2)2sn+2(A)2(n+1)2sn+1(A)+n2sn(A)2\hat R_{2,n} = \frac{(n+2)^2s_{n+2}(A) - 2(n+1)^2s_{n+1}(A) + n^2s_n(A)}{2}

然后类似地计算33阶Richardson:

R^3,n=(n+3)3sn+3(A)3(n+2)3sn+2(A)+3(n+1)3sn+1(A)n3sn(A)6\hat R_{3,n} = \frac{(n+3)^3s_{n+3}(A)-3(n+2)^3s_{n+2}(A)+3(n+1)^3s_{n+1}(A)-n^3s_n(A)}{6}

就可以发现规律,对于任意的kk展开kk阶以后模仿上述步骤,就可以得到kk阶Richardson公式:

R^k,n=1k!i=0k(1)kiCki(n+i)ksn+i(A)\hat R_{k,n} = \frac{1}{k!} \sum_{i=0}^k (-1)^{k-i}C_k^i(n+i)^ks_{n+i}(A)

由于Richardson不擅长处理交错数列,我们改用ζ(2)\zeta(2)来检验其有效性:

1+122+132+142+...=π261 + \frac{1}{2^2} + \frac{1}{3^2} + \frac{1}{4^2} + ... = \frac{\pi^2}{6}

代码实现:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include <bits/stdc++.h>
using namespace std;
#define MAXN 100000
#define db long double
#define PI 3.1415926535897932384626433832795
#define rint register int
inline int rf(){int r;int s=0,c;for(;!isdigit(c=getchar());s=c);for(r=c^48;isdigit(c=getchar());(r*=10)+=c^48);return s^45?r:-r;}
db A[MAXN+5], _; int n, K; long long C[16][16], fac[16];
inline db Richardson(int k, int n){_ = 0; for(rint j = 0; j <= k; _ += ((k-j)&1?-1:1)*C[k][j]*powl(n+j,k)*A[n+j], j++); return _/fac[k];}
int main(){
for(rint i = C[0][0] = 1, j; i <= 15; i++) for(j = C[i][0] = 1; j <= 15; C[i][j] = C[i-1][j]+C[i-1][j-1], j++);
for(rint i = fac[0] = 1; i <= 15; fac[i] = 1ll*fac[i-1]*i, i++);
n = rf(); K = rf(); printf("pi = %.14Lf\n",(db)PI*(db)PI/6);
for(rint i = 1; i <= n+K; A[i] = A[i-1]+(db)1/(1ll*i*i), i++);
printf(" "); for(rint j = 0; j <= K; printf("Richardson %2d: ",j), j++); puts("");
for(rint i = 1, j; i <= n; i++) for(printf("n = %3d: ",i), j = 0; j <= K; printf("%.14Lf %c",Richardson(j,i),"\n "[j<K]), j++);
return 0;
}

下面是实测结果:

可以看到高阶Richardson出奇的有效,10位小数的精度都已经不足以体现出区别。

当然,在高阶Richardson的基础上在进行层层嵌套,应该可以得到更优秀的结果……


总结

Shanks变换和Richardson外推法都是通过构造一个带有附加项的近似序列来解决收敛级数求和的问题。

Shanks变换更适用于交错数列,Richardson外推法更适用于单调数列

注意两者都只用到了数列的某几项的信息。这意味着这两种方法不依赖于具体的数列/函数。同样这也意味着两种方法存在不足: 无法预测存在"突变"的函数。


扫描二维码即可在手机上查看这篇文章,或者转发二维码来分享这篇文章:


文章作者: Magolor
文章链接: https://magolor.cn/2019/08/04/2019-08-04-blog-01/
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 Magolor
扫描二维码在手机上查看或转发二维码以分享Magolor的博客