请选择 进入手机版 | 继续访问电脑版

计算化学公社

 找回密码
 现在注册!
查看: 507|回复: 5

[数值算法] 数值计算编程中务必要注意【数值不稳定】问题

[复制链接]

75

帖子

5

威望

1082

eV
积分
1257

Level 4 (黑子)

发表于 2018-9-22 12:15:19 | 显示全部楼层 |阅读模式
本帖最后由 coolrainbow 于 2018-9-22 12:42 编辑

问题的提出

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

程序中真正的f(x)非常复杂,这里为了简明,我们将问题的精华提取出来。f(x)类似于下面这个模样,其中x>=0:
CodeCogsEqn (3).png
这个公式在数学上是完全正确的,而且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展开:
CodeCogsEqn (1).png
Taylor展开是没有奇点的(0是可去奇点),因此只要取前几项就可以精确的计算f(x)在x很小时候的值。

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






CodeCogsEqn (1).png
CodeCogsEqn (2).png
CodeCogsEqn.png
CodeCogsEqn (3).png

评分

参与人数 10威望 +1 eV +45 收起 理由
Jasminer + 5 谢谢提醒
highlight + 5 精品内容
wangxubo + 5 赞!
wbn + 5 谢谢
wangyj + 5 谢谢分享
asdf + 5 精品内容
一颗赛艇 + 5 精品内容
sobereva + 1
liyuanhe211 + 5 许久以前被坑过,印象深刻
ene + 5 谢谢

查看全部评分

151

帖子

0

威望

335

eV
积分
486

Level 3 能力者

发表于 2018-9-22 14:03:43 | 显示全部楼层
本帖最后由 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

151

帖子

0

威望

335

eV
积分
486

Level 3 能力者

发表于 2018-9-22 14:03:54 | 显示全部楼层
本帖最后由 granvia 于 2018-9-22 14:05 编辑

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

69

帖子

1

威望

890

eV
积分
979

Level 4 (黑子)

发表于 2018-9-22 16:00:24 | 显示全部楼层
多谢分享,确实是应该仔细注意的问题

75

帖子

5

威望

1082

eV
积分
1257

Level 4 (黑子)

 楼主| 发表于 2018-9-22 21:25:36 | 显示全部楼层
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(..)收敛要快的多。

183

帖子

3

威望

1524

eV
积分
1767

Level 5 (御坂)

发表于 2018-9-23 08:52:29 | 显示全部楼层
开方一个很小的数也会有问题。
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

手机版|北京科音自然科学研究中心|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949-1号 )

GMT+8, 2018-12-16 10:35 , Processed in 0.120124 second(s), 28 queries .

快速回复 返回顶部 返回列表