关于ln2

今天我在读刚借的《数值实验》,黄友谦编,高教版。
上面第一个例子就是计算ln2,采用方式是把ln(1+x)在0处泰勒展开,也就是麦克劳林展开,然后让x等于1。
书上的例子是为了阐述舍入误差,他说他采用的是10位10进制的计算机上算的,my god.我不懂。
首先,我用bc算ln2的精确值,在保留100位精度的情况下
ln2=0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875
然后我在我的机子上写了个很小的程序

#include <iostream>
#include <cmath>

int main(int argc,char* argv[]){
    typedef double T;
    long n=strtol(argv[1],NULL,10);
    n++;
    T sum(0);
    for(int i=1;i!=n;++i)
        if(i%2) sum+=static_cast<T>(1)/i;
        else sum-=static_cast<T>(1)/i;
    std::cout<<sum<<std::endl;
    std::cout<<"实际误差:"<<log(2)-sum<<std::endl; 
    std::cout<<"理论最大误差:"<<(1.0/n)<<std::endl; 
    return 0;
}

发现误差与理论值有惊人的相似。
n等于5000的时候,误差0.00998
$ ./ln2 5000
0.693247
n等于50000的时候,误差0.00098
$ ./ln2 50000
0.693157
我把计算过程中的double换成了float,结果依旧,不知道是不是因为发生了类型提升所以和原来的代码其实是一样的。


迭代次数

计算值

绝对误差

理论最大误差

15

0.725372

-0.0322247

0.0588235

24

0.672747

0.0203997

0.0384615

64

0.685396

0.00775147

0.0153846

56768

0.693138

8.8077e-06

1.76152e-05

81954

0.693141

6.10095e-06

1.22018e-05

99274

0.693142

5.03654e-06

1.0073e-05

书上的值为


迭代次数

计算值

绝对误差

15

0.7253718504

3.22246698e-2

24

0.6727474994

-2.03996812e-2

64

0.6853957079

-7.7514727e-3

56768

0.6931383909

8.7897-e6

81954

0.6931411057

6.0749e-6

99274

0.6931421806

-5e-6

基本完全一致。
按拉格朗日余项方式计算,|R(n)|应小于1/(n-1)
实际误差恰好为理论最大误差的一半,不仅仅是巧合,我觉得这个可能与此级数是交错级数有关。
然后说问题出在了哪里。
书上说
"若要求|Xn-ln2|<p,只需要取p满足1/(n+1)<p,当p=0.5e-4时,n约等于5000;当p=0.5e-5时,n约等于50000;当p=0.5e-6时,n约等于50万;……首先要注意到,计算结果与估计式严重不符,例如要使|Xn-ln2|<0.5e-5,由估计式知n约等于50000,而实际计算结果n等于99274……"
我真信了作者的鬼话了,其实自己拿笔算算1/(n+1)<p,当p=0.5e-4时,n应该等于多少???不是5000,而是20000.后面的依次类推,全部算错了。所以,计算到99274次满足了他的要求那不过是巧合,理论需要200000次,而上面我们已经看出来了,实际误差恰好是理论误差的一半(说明理论可以进一步修正),于是我们拿100000次来算,
计算值:0.693142
实际误差:4.99997e-06
理论最大误差:9.9999e-06
他计算了99274得到了实际误差为5e-06的结果,无怪乎……
后记:
1。这是一个很好的教训。数学一定要严谨,推导过程中半点马虎不得,否则会被自己错误的发现而洋洋得意。
2。数值计算的关键已经不在于如何去算,难点在于误差估计。通常我们不是要告诉用户,你去迭代吧,迭代100次不够就迭代1000次,迭代1000次不够就迭代10000次,而是要通过误差估计理论,只要用户告诉我们需要精确到多少位,我们就可以给出他需要的结果。

此博客中的热门博文

少写代码,多读别人写的代码

在windows下使用llvm+clang

tensorflow distributed runtime初窥