计算化学公社

标题: 数值计算编程中务必要注意【数值不稳定】问题 [打印本页]

作者
Author: coolrainbow 时间: 2018-9-22 12:15
标题: 数值计算编程中务必要注意【数值不稳定】问题
本帖最后由 coolrainbow 于 2018-9-22 12:42 编辑

问题的提出

昨天帮人查一个程序的bug。是一个看起来无害的纯数值程序，给定一个x，算出 f(x)。然而在x接近0的时候，给出了十分离谱的结果。最后发现是使用的公式数值不稳定造成的。这种bug常常难以识别，故这里记录下来，供大家参考。

程序中真正的f(x)非常复杂，这里为了简明，我们将问题的精华提取出来。f(x)类似于下面这个模样，其中x>=0：
(, 下载次数 Times of downloads: 104)
这个公式在数学上是完全正确的，而且0只是个可去奇点，意思是f(0)=0，这可以通过求极限证明。但是在用到计算机上时却十分危险。如果写这样的代码

((1.-3./x+3./x/x)-(1.+3./x+3./x/x)*exp(-2*x))/x2;

在x为比较大的数时候计算，没有问题。但是若是x比较小，例如0.001，则

x = 0.001
精确值：0.000000066600038079
计算值： 0.000000000000000000

计算值和精确值间有10^(-8)的误差。这个误差从人的角度看起来很小，但是在计算机的角度上看是个十分显著的。实际上，这个代码在真正的程序中带来的误差高达10^(-3)以上。因此，这段代码在x较小的时候是绝对不能用的，是错误的。

那么问题究竟处在哪里？观察原始公式，我们注意到括号里是两个数字相减。而当x<0.001时，1/x^2已经高达1000000，也就是说，f(x)是两个1000000数量级的数字相减，这几乎必然会导致0.1以下精度大大丢失，所以才会有显著误差。而且x再小的时候，1/x和1/x^2可能会溢出，此时的计算结果已经完全错误，没有价值。

我们来仔细看一下括号里的两项。在x=0.001时，这两项是这样的：

                        (1-3./x+3/x/x)                   (1+3./x+3/x/x)*exp(-2.*x)
精确值：2997001.000000000000000000    2997000.999999999866799923
计算值：2997001.000000000000000000    2997001.000000000000000000

果然，计算(1+3./x+3/x/x)*exp(-2.*x) 时，2997000.99......精度丢失成了2997001.000000000000000000，自然也就产生了误差。这个误差的来源可以从下面的分析看出：

                              (1+3./x+3/x/x)             exp(-2.*x)             (1+3./x+3/x/x)*exp(-2.*x)
精确值： 3003001.000000000000000000    0.998001998667333067 2997000.999999999866799923
计算值： 3003001.000000000000000000    0.998001998667333079 2997001.000000000000000000

因此，这个算法的误差源于

的数值不稳定。其本质原因是，这个表达式从数学上是有界的，可以证明其数值在[0,1]中，但这是两个无界项(...1/x...) 相减抵消的结果。每一项在x趋于0是都是无穷大的，计算机在计算x趋于零的过程中，精度逐渐不足以保证两个无穷大项的精确抵消，所以误差逐渐增大，也即不稳定。

再总结两点：

一般来说，只要数学表达式中包含类似1/x相减的项，而x又可能很小，那么几乎一定有数值稳定性问题。必须谨慎测试程序的可靠性。
这种数值稳定性是针对“本来应该稳定”的表达式而言的。比如本文中提到的。又比如f(x)=(exp(x)-1)/x，数学上知道x=0是它的值为1，但是直接用计算机算特别小的x时，可能会溢出从而不等于1，这说明这个表达式数值不稳定。而“本来”就发散的表达式，如f(x) = 1/x，x趋于零是本来就是无穷大，计算机也是无穷大，这个表达式就谈无所谓“不稳定”。

解决方法

说了半天，如何正确计算

在x特别小的时候的值呢？很遗憾，这个表达式，也许用一些trick，如调整数序等可以稍微提高一下精度，但是只要稍有意外还是会给出荒唐的结果。根本解决方法就是不要用这个表达式。那么怎么办？

考虑简单的问题：f(x) = (exp(x)-1)/x。我们把它Taylor展开：
(, 下载次数 Times of downloads: 108)
Taylor展开是没有奇点的（0是可去奇点），因此只要取前几项就可以精确的计算f(x)在x很小时候的值。

因此，对这个问题可以用类似的技术。不过直接exp(-2x)的Taylor展开会发现并不好用。我们把这个表达式变得对称了一些
(, 下载次数 Times of downloads: 111)
现在，把括号里的两个exp做Taylor展开，此时你会发现所有的1/x都消失了（x是可去奇点）。取前几项，就可以计算数值稳定的f(x)的值了。

作者
Author: granvia 时间: 2018-9-22 14:03
本帖最后由 granvia 于 2018-9-22 14:09 编辑

验证了一下，用Taylor展开exp(-2*x)前4项得到的表达式是：x^2/15*(2*x^2+x+1) = 6.67334666666667e-08
而直接用MATLAB算式子exp(-x)/2/x*( (1-3/x+3/x^2)*exp(x) - (1+3/x+3/x^2)*exp(-x) )得到的是2.32597929386742e-07

作者
Author: granvia 时间: 2018-9-22 14:03
本帖最后由 granvia 于 2018-9-22 14:05 编辑

我机子上的MATLAB的eps是2.22044604925031e-16

作者
Author: ene 时间: 2018-9-22 16:00
多谢分享，确实是应该仔细注意的问题

作者
Author: coolrainbow 时间: 2018-9-22 21:25

granvia 发表于 2018-9-22 14:03
验证了一下，用Taylor展开exp(-2*x)前4项得到的表达式是：x^2/15*(2*x^2+x+1) = 6.67334666666667e-08
而 ...

exp(-x)/2/x*( (1-3/x+3/x^2)*exp(x) - (1+3/x+3/x^2)*exp(-x) )也是要进行Taylor展开计算的。直接算还是有误差。

之所以选择第二个形式，是因为它的Taylor展开是x^2+x^4+x^6....形式的，比直接展开x^2/15(..)收敛要快的多。

作者
Author: Warm_Cloud 时间: 2018-9-23 08:52
开方一个很小的数也会有问题。

欢迎光临计算化学公社 (http://bbs.keinsci.com/)