主题：[讨论]基本功训练 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个回复）
沙发
jstzhurj [专家分：4680] 发布于 2010-10-13 10:44:00
我试一下，不知精度是否满足楼主的要求？
program main
implicit none
integer :: nd=365
real(kind=4) :: x=0.06,y,z
y=x/nd
z=x-0.5*x*y
write(*,*) 100.*(exp(z)-1.)/y
end program
[em18][em18]
板凳
yeg001 [专家分：14390] 发布于 2010-10-13 12:15:00
lim(1+(1/x))^x 当x趋于无穷极限趋于e.
于是
(1 + x / ND) ** ND 约等于 e^x 亦即exp(x)
(1 + x / ND) ** ND - 1 用泰勒展开得到: x+(x^2/2!)+(x^3/3!)+(x^4/4!)+... 截断到4或者5似乎就够了.
其他的应该计算误差比较小了.
不知道对不对, 没写代码.
3 楼
asymptotic [专家分：16630] 发布于 2010-10-13 12:21:00
1 楼网友用 幂指函数，Taylor 级数展开 Ln(1 + x) 在处理 x 比较小的情形很理想；但当 x 比较大时，不如 13.56，则不能保证至少 6 位精度了。 再好好想想，看有何好的方法。
提示：若是一味地增加 Taylor 级数的项数，效果未必明显，因为这个设计交错级数的求和，收敛很慢。
4 楼
jstzhurj [专家分：4680] 发布于 2010-10-13 12:24:00
2楼，写个代码试试！
5 楼
jstzhurj [专家分：4680] 发布于 2010-10-13 12:36:00
[quote]1 楼网友用 幂指函数，Taylor 级数展开 Ln(1 + x) 在处理 x 比较小的情形很理想；但当 x 比较大时，不如 13.56，则不能保证至少 6 位精度了。 再好好想想，看有何好的方法。
提示：若是一味地增加 Taylor 级数的项数，效果未必明显，因为这个设计交错级数的求和，收敛很慢。[/quote]
正因为x比较小才用这个方法，针对实际问题的，没有通用性。
6 楼
jstzhurj [专家分：4680] 发布于 2010-10-13 12:40:00
对于单精度，x小到没有必要用3次方项了。
7 楼
aliouying [专家分：1150] 发布于 2010-10-13 12:47:00
program main
implicit none
integer :: nd
real(kind=4) :: x,y
nd=365
x=0.06
y=x/nd
write(*,*) 100./y*( exp(ND*log(1+y)) -1. )
end program
37614.05
Press any key to continue
这个通用性如何？
8 楼
jstzhurj [专家分：4680] 发布于 2010-10-13 12:55:00
[quote]program main
implicit none
integer :: nd
real(kind=4) :: x,y
nd=365
x=0.06
y=x/nd
write(*,*) 100./y*( exp(ND*log(1+y)) -1. )
end program
37614.05
Press any key to continue
这个通用性如何？
[/quote]
看来我用泰勒展开是多此一举！
9 楼
jstzhurj [专家分：4680] 发布于 2010-10-13 13:15:00
[quote]program main
implicit none
integer :: nd
real(kind=4) :: x,y
nd=365
x=0.06
y=x/nd
write(*,*) 100./y*( exp(ND*log(1+y)) -1. )
end program
37614.05
Press any key to continue
这个通用性如何？
[/quote]
我这儿结果是：
37615.45
请按任意键继续. . .
[em18][em18][em18]
10 楼
hubangdie [专家分：0] 发布于 2010-10-13 13:44:00
http://www.soier.net/?45986.htm
我来回复
您所在位置：