主题：Fortran精度问题
tianao87
[专家分：0] 发布于 2011-05-16 11:29:00
用以下程序计算tent map时间序列，结果大家可以看到迭代十几步之后精度的误差就这么大，完全偏离预期结果，不知道要怎么办，求解！！！！
program tent
implicit none
parameter N=1000
integer :: i
double precision :: x(N)
open(10,file='data.txt')
x(1)=0.2_8
do i=2,N
x(i)=1.0-2.0*abs(x(i-1)-0.5)
end do
do i=1,N
write(10,*)i,x(i)
enddo
end
结果：
1 0.200000000000000
2 0.400000000000000
3 0.800000000000000
4 0.400000000000000
5 0.800000000000000
6 0.400000000000000
7 0.800000000000001
8 0.399999999999999
9 0.799999999999997
10 0.400000000000006
11 0.800000000000011
12 0.399999999999977
13 0.799999999999955
14 0.400000000000091
15 0.800000000000182
PS：编译环境用的是CVF6.6
最后更新于：2011-05-16 11:34:00
回复列表 （共16个回复）
沙发
yeg001 [专家分：14390] 发布于 2011-05-16 12:03:00
x(i)=1.0_8-2.0_8*abs(x(i-1)-0.5_8)
这样会不会好一些?
不过如果遵守规范的话, 1.0, 2.0, 0.5转换成双精度应该没有问题.
看楼主的结果应该就是截断误差的累加. 每次循环都乘以2.
板凳
tianao87 [专家分：0] 发布于 2011-05-16 12:08:00
[quote]x(i)=1.0_8-2.0_8*abs(x(i-1)-0.5_8)
这样会不会好一些?[/quote]
试过了，结果还是那样！
3 楼
yeg001 [专家分：14390] 发布于 2011-05-16 12:36:00
先问楼主一个问题, 这个程序的目的是不是用来测量截断误差的?
4 楼
tianao87 [专家分：0] 发布于 2011-05-16 12:50:00
[quote]先问楼主一个问题, 这个程序的目的是不是用来测量截断误差的?[/quote]
当然不是，用来生成一个时间序列的，也就是x(N)
5 楼
yeg001 [专家分：14390] 发布于 2011-05-16 14:08:00
分析了一下, 截断误差是由初值x(1)=0.2_8引入的. 后面的算法怎么调整也没用.
把后面输出部分用格式化
write(10,*)i,x(i) => write(10,fmt="(I5,F23.18)")i,x(i)
你会发现二进制浮点数表示无法精确表示0.2, 当2.0*累计到一定程度就出现你所说的精度问题.
这个与语言无关, 不是fortran的问题.
6 楼
tianao87 [专家分：0] 发布于 2011-05-16 14:43:00
[quote]分析了一下, 截断误差是由初值x(1)=0.2_8引入的. 后面的算法怎么调整也没用.
把后面输出部分用格式化
write(10,*)i,x(i) => write(10,fmt="(I5,F23.18)")i,x(i)
你会发现二进制浮点数表示无法精确表示0.2, 当2.0*累计到一定程度就出现你所说的精度问题.
这个与语言无关, 不是fortran的问题.[/quote]
受教了，十分感谢，那么可有办法让我进行迭代之后使截断误差不影响我的结果，迭代十几次就有这么大的误差接受不了！
7 楼
yeg001 [专家分：14390] 发布于 2011-05-16 15:10:00
我不懂这个时间序列怎样才算是对的.
我是想算法可不可以反过来, 从你序列最终终止的1.0开始反过来算出每个时序的数? 而不是随手设定的一个例如0.2这样的常数?
8 楼
tianao87 [专家分：0] 发布于 2011-05-16 15:16:00
[quote]我不懂这个时间序列怎样才算是对的.
我是想算法可不可以反过来, 从你序列最终终止的1.0开始反过来算出每个时序的数? 而不是随手设定的一个例如0.2这样的常数?[/quote]
如果是精确的迭代的话，应该一直是0.8和0.4两个值一个周期序列，结果由于误差的影响使序列最终稳定在0.0了（60个点之后的值全为0.0），这个结果严重影响最终结论啊，搞得我很受伤啊。。。。
9 楼
dongyuanxun [专家分：7180] 发布于 2011-05-16 15:32:00
[quote]我不懂这个时间序列怎样才算是对的.
我是想算法可不可以反过来, 从你序列最终终止的1.0开始反过来算出每个时序的数? 而不是随手设定的一个例如0.2这样的常数?[/quote]
计算机只有计算整型是绝对精确的。
如果你需要精确计算，那么不要使用浮点数，请全部使用整型。
然后再把结果除以一个倍数就是你要的结果。
先明白理论解、数值解、误差限之间的关系。
10 楼
dongyuanxun [专家分：7180] 发布于 2011-05-16 15:38:00
这样误差是最小的
program tent
implicit none
integer,parameter :: N=1000
integer :: i
integer :: x(N)
open(10,file='data.txt')
x(1)=2
do i=2,N
x(i)=10-2*abs(x(i-1)-5)
end do
do i=1,N
write(10,*)i,x(i)/10._8
enddo
end
当然也可以输出精确解，以字符串之类的方式构造类型，可以参考各种高精度算法。
我来回复
您所在位置：