主题：[讨论]基本功训练 04】浮点数计算过程中的误差
asymptotic
[专家分：16630] 发布于 2010-10-13 09:10:00
用 kind = 4 的 real、integer 计算如下表达式：
100 * ( (1 + x / ND) ** ND - 1 ) / ( x / ND )
其中，ND = 365， x = 0.06； 精确答案 -- 37614.05。
编写程序，精确到有效数字第六位，第七位可以有区别。
提示：这个问题主要是设计算法，不能考虑更高精度的浮点数、整数，
以减少浮点数计算过程中的误差（当然，若是这样，则要求至少 14 位
有效数字正确）。
回复列表 （共57个回复）
11 楼
jstzhurj [专家分：4680] 发布于 2010-10-13 14:07:00
楼上捣什么乱呢，这样可以不？
program main
implicit none
integer :: nd=365,i
real(kind=4) :: x=0.06,y,z,term,esp=1.e-8
y=x/nd
term=1.
z=0.
i=1
do while(term.gt.esp)
term=term*nd*y/i
z=z+term
i=i+1
nd=nd-1
enddo
write(*,*) 100.*z/y
end program
37614.05
请按任意键继续. . .
12 楼
yeg001 [专家分：14390] 发布于 2010-10-13 14:30:00
我按自己的想法写了一下, 似乎误差是大了些.
program main
implicit none
integer, parameter :: SP=4
integer(kind=4) :: nd=365, i, factorial
real(kind=SP) :: x=0.06, y, z
y=0.0
i=1
factorial = 1
DO
factorial = i*factorial
z = x**i/REAL(factorial,kind=SP)
y = y + z
i = i +1
IF((z/x)<1e-6) EXIT
END DO
write(*, *) i-1
write(*, *) 100.*y/x*REAL(nd, kind=SP)
end program
5
37617.23
Press any key to continue
泰勒展开到5就稳定了, 但是结果倒数第三位不准. 看来直接用极限近似 x/ND 还不够小.
13 楼
yeg001 [专家分：14390] 发布于 2010-10-13 14:36:00
看来自然对数函数和e幂指数函数精度做得挺高的. 那近似展开似乎就真的多此一举了.
14 楼
jstzhurj [专家分：4680] 发布于 2010-10-13 14:37:00
[quote]我按自己的想法写了一下, 似乎误差是大了些.
program main
implicit none
integer, parameter :: SP=4
integer(kind=4) :: nd=365, i, factorial
real(kind=SP) :: x=0.06, y, z
y=0.0
i=1
factorial = 1
DO
factorial = i*factorial
z = x**i/REAL(factorial,kind=SP)
y = y + z
i = i +1
IF((z/x)<1e-6) EXIT
END DO
write(*, *) i-1
write(*, *) 100.*y/x*REAL(nd, kind=SP)
end program
5
37617.23
Press any key to continue
泰勒展开到5就稳定了, 但是结果倒数第三位不准. 看来直接用极限近似 x/ND 还不够小.[/quote]
我想 (1 + x / ND) ** ND - 1 用泰勒展开得到: x+(x^2/2!)+(x^3/3!)+(x^4/4!)+... 是有问题的！
15 楼
jstzhurj [专家分：4680] 发布于 2010-10-13 14:58:00
楼主，还有什么更精确的办法，拿出来分享一下，我是黔驴技穷了！！！[em20][em15][em19]
16 楼
yeg001 [专家分：14390] 发布于 2010-10-13 15:17:00
e^x 的收敛半径是整个实数区域, 只要项数足够多应该是没问题的啊.
17 楼
jstzhurj [专家分：4680] 发布于 2010-10-13 15:25:00
[quote]e^x 的收敛半径是整个实数区域, 只要项数足够多应该是没问题的啊.[/quote]
(1 + x / ND) ** ND 等于e^x 吗？ 尽管ND比较大，但不能等同于无穷！！这是问题所在！
18 楼
yeg001 [专家分：14390] 发布于 2010-10-13 15:48:00
[quote][quote]e^x 的收敛半径是整个实数区域, 只要项数足够多应该是没问题的啊.[/quote]
(1 + x / ND) ** ND 等于e^x 吗？ 尽管ND比较大，但不能等同于无穷！！这是问题所在！[/quote]
似乎兄弟你这个观点我已经在贴代码的帖子说了. 数值计算有不少类似的近似.
我不知道用exp(x)求解能保证精度到哪一位, 还有取一次ln对数也不知道能精确到那里, 所以我想法就是展开成多项式精度更好得到保证. 不过那个自然常数近似引入多少误差我没有取评估.
19 楼
jstzhurj [专家分：4680] 发布于 2010-10-13 15:58:00
[quote][quote][quote]e^x 的收敛半径是整个实数区域, 只要项数足够多应该是没问题的啊.[/quote]
(1 + x / ND) ** ND 等于e^x 吗？ 尽管ND比较大，但不能等同于无穷！！这是问题所在！[/quote]
似乎兄弟你这个观点我已经在贴代码的帖子说了. 数值计算有不少类似的近似.
我不知道用exp(x)求解能保证精度到哪一位, 还有取一次ln对数也不知道能精确到那里, 所以我想法就是展开成多项式精度更好得到保证. 不过那个自然常数近似引入多少误差我没有取评估.[/quote]
对于误差要求比较小的，你这种近似足够了，但是对于楼主提出的误差，这种近似是无法达到其要求的！
20 楼
yeg001 [专家分：14390] 发布于 2010-10-13 19:48:00
把幂运算用x=e^(ln(x))转换运算挺依赖lnx和e^x运算的精度. 我错觉一直认为他们运算的精度不高.
楼主, 你的准确值是什么算法得到的?
我来回复
您所在位置：