[问题] 循环拆开计算后与原本值有微小误差

楼主: boa85391 ( )   2015-08-12 21:27:43
do i=1, 3999, 1
do j=i+1, 4000, 1
rx=bead(j)%x_p - bead(i)%x_p
ry=bead(j)%y_p - bead(i)%y_p
rz=bead(j)%z_p - bead(i)%z_p
rpq= dsqrt(rx**2+ry**2+rz**2)
if (rpq<=r_c) then
if (bead(i)%bead_number==1 .and. bead(j)%bead_number==1) then
rep_term=dabs(24.d0*eps_hh*(2.d0*sigma**12*rpq**(-14)&
-sigma**6*rpq**(-8)))
bead(i)%x_frep=bead(i)%x_frep - rep_term * rx
bead(i)%y_frep=bead(i)%y_frep - rep_term * ry
bead(i)%z_frep=bead(i)%z_frep - rep_term * rz
bead(j)%x_frep=bead(j)%x_frep + rep_term * rx
bead(j)%y_frep=bead(j)%y_frep + rep_term * ry
bead(j)%z_frep=bead(j)%z_frep + rep_term * rz
else if (bead(i)%bead_number/=1 .and. bead(j)%bead_number/=1) then
rep_term=dabs(24.d0*eps_tt*(2.d0*sigma**12*rpq**(-14)&
-sigma**6*rpq**(-8)))
bead(i)%x_frep=bead(i)%x_frep - rep_term * rx
bead(i)%y_frep=bead(i)%y_frep - rep_term * ry
bead(i)%z_frep=bead(i)%z_frep - rep_term * rz
bead(j)%x_frep=bead(j)%x_frep + rep_term * rx
bead(j)%y_frep=bead(j)%y_frep + rep_term * ry
bead(j)%z_frep=bead(j)%z_frep + rep_term * rz
if (bead(i)%Lipid_number/=bead(j)%Lipid_number) then
bead(i)%x_fatt=bead(i)%x_fatt - rep_term * rx
bead(i)%y_fatt=bead(i)%y_fatt - rep_term * ry
bead(i)%z_fatt=bead(i)%z_fatt - rep_term * rz
bead(j)%x_fatt=bead(j)%x_fatt + rep_term * rx
bead(j)%y_fatt=bead(j)%y_fatt + rep_term * ry
bead(j)%z_fatt=bead(j)%z_fatt + rep_term * rz
end if
else
rep_term=dabs(24.d0*eps_ht*(2.d0*sigma**12*rpq**(-14)&
-sigma**6*rpq**(-8)))
bead(i)%x_frep=bead(i)%x_frep - rep_term * rx
bead(i)%y_frep=bead(i)%y_frep - rep_term * ry
bead(i)%z_frep=bead(i)%z_frep - rep_term * rz
bead(j)%x_frep=bead(j)%x_frep + rep_term * rx
bead(j)%y_frep=bead(j)%y_frep + rep_term * ry
bead(j)%z_frep=bead(j)%z_frep + rep_term * rz
end if
else if (r_c<=rpq .and. rpq<=(r_c+w_tt) .and. &
bead(i)%Lipid_number/=bead(j)%Lipid_number .and. &
bead(i)%bead_number/=1 .and. bead(j)%bead_number/=1) then
att_term=dabs(eps_tt*pi/w_tt/rpq*dcos(pi*(rpq-r_c)/2.d0/w_tt)* &
dsin(pi*(rpq-r_c)/2.d0/w_tt))
bead(i)%x_fatt=bead(i)%x_fatt + att_term * rx
bead(i)%y_fatt=bead(i)%y_fatt + att_term * ry
bead(i)%z_fatt=bead(i)%z_fatt + att_term * rz
bead(j)%x_fatt=bead(j)%x_fatt - att_term * rx
bead(j)%y_fatt=bead(j)%y_fatt - att_term * ry
bead(j)%z_fatt=bead(j)%z_fatt - att_term * rz
end if
end do
end do
主要是有4000颗珠子,计算彼此间距离一条件算作用力,rep是排斥力,att是吸引力,
除了if判定式里面的number存成integer,其他都是real(8)
为了方便使用openMP,所以我把它写成两个与上式一样的双层循环,
一个只算i的部分,一个只算j的部分,可是最后得到的值有些会在小数点后十几位有误差
尝试输出每一步的rep_term去比较,存到小数点后15位都还是没什么问题
想请问一下这是合理的误差产生吗?还是我的写法会结构上有问题,
已经在这部分弄了快两天,因为是直观的同样循环算两次,还是无法去忽略这些误差。
作者: perceval (摸鱼中)   2015-08-12 23:45:00
浮点数误差
楼主: boa85391 ( )   2015-08-13 01:19:00
请问一下有办法修正这种误差吗?我把j部分的运算用另外的array储存最后再加总回来误差才终于不见这些小毛病真的让人闷闷的
作者: perceval (摸鱼中)   2015-08-13 02:25:00
arbitrary precision package 但是这个误差不用管他你把浮点数弄清楚就知道这个误差的来源
楼主: boa85391 ( )   2015-08-13 12:13:00
谢谢回复! 读了一下相关资料,猜测比较可能是大量相近位数的值相减造成的精确位数丧失。

Links booklink

Contact Us: admin [ a t ] ucptt.com