24小时热门版块排行榜    

查看: 1705  |  回复: 6
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

astringent

铜虫 (著名写手)

[交流] 【求助】相关函数归一化问题 已有4人参与

最近一直在看自相关函数的程序,有个问题想不通。
自相关函数的定义是A=(a(t)*a(0))/(a(0)*a(0)),当对体系中所有粒子求和时,是先算出每个例子的(a(t)*a(0)) /(a(0)*a(0))然后求和?还是先分别求a(t)*a(0)以及a(0)*a(0)和然后再相除呢?两者有什么区别?请大家指点。

[ Last edited by ghcacj on 2010-12-8 at 09:27 ]
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

astringent

铜虫 (著名写手)

引用回帖:
Originally posted by 老虎大王 at 2010-12-07 10:59:29:
问得好。我初学时也在这个问题上纠结过。后来我写程序把两种算法同时都用了,输出两列结果。经比较,二者在数值上差别是很小的。另外,按Leach的书(molecular modelling...),是对每个粒子,上下分别对时间原点求 ...

为什么我用两种方法得到的结果有很大差别呢,劳烦指教。
下面是我的程序
这是我的程序:
do i=1,10000 (粒子数)
     do ntel=1,npoints (共有npoints个数据)
          if (确定timeorigin的条件) then
             t0=t0+1
             time0(t0)=ntel
             tcf0(t0)=tcf(ntel)         
             do t=1,t0
                  delt=ntel-time0(t)+1
                  if(delt.lt.npoints)
                  tcf(delt)=tcf(delt)+tcf(ntel)*tcf0(t)
                  tcf0(delt)=tcf0(delt)+tcf0(t)*tcf0(t)
                  acf(delt)=tcf(delt)+tcf(ntel)*tcf0(t)/ tcf0(t)*tcf0(t)
              endif
            enddo
           endif
    enddo
enddo
do ical=1,npoints
tcf(ical)=(tcf(ical)/(npoints*ntime(t))/(tcf0(ical)/(npoints*ntime(t))
acf(ical)=acf(ical)/(npoints*ntime(t))
write(*,*) ical, tcf(ical)
write(*,*) ical, acf(ical)
enddo
下面是我的输出结果(我只截取了一部分)
  0.0000000000000000        1.0000000000000000     
   0.0000000000000000       1.14207400639561443E-004
   1.0000000149011612       0.67655786350148372     
   1.0000000149011612       7.72679149727596715E-005
   2.0000000298023224       0.62686567164179108     
   2.0000000298023224       7.15926989083817998E-005
   3.0000000447034836       0.60360360360360366     
   3.0000000447034836       6.89359985842397961E-005
   4.0000000596046448       0.58610271903323263     
   4.0000000596046448       6.69372680485647080E-005
   5.0000000745058060       0.57012195121951215     
   5.0000000745058060       6.51121460963353292E-005
   6.0000000894069672       0.55828220858895705     
   6.0000000894069672       6.37599598662582241E-005
   7.0000001043081284       0.54629629629629639     
   7.0000001043081284       6.23910799790196834E-005
   8.0000001192092896       0.51552795031055898     
   8.0000001192092896       5.88771071620099369E-005
   9.0000001341104507       0.52187499999999998     
   9.0000001341104507       5.96019872087711250E-005
   10.000000149011612       0.52515723270440251     
   10.000000149011612       5.99768424742350986E-005
5楼2010-12-10 00:17:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 7 个回答

老虎大王

木虫 (著名写手)

★ ★ ★ ★ ★ ★
小木虫(金币+0.5):给个红包,谢谢回帖交流
ghcacj(金币+5):谢谢 2010-12-07 16:26:07
问得好。我初学时也在这个问题上纠结过。后来我写程序把两种算法同时都用了,输出两列结果。经比较,二者在数值上差别是很小的。另外,按Leach的书(molecular modelling...),是对每个粒子,上下分别对时间原点求平均后再除,然后再对粒子求和,然后再除以总粒子数,也就是用的你的第一个方案。

也就是:
Cvv=(1/N)*sigma(i 从1到N)(/
2楼2010-12-07 10:59:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

astringent

铜虫 (著名写手)

引用回帖:
Originally posted by 老虎大王 at 2010-12-07 10:59:29:
问得好。我初学时也在这个问题上纠结过。后来我写程序把两种算法同时都用了,输出两列结果。经比较,二者在数值上差别是很小的。另外,按Leach的书(molecular modelling...),是对每个粒子,上下分别对时间原点求 ...

请问老虎大王,我在写程序的时候直接把t=0的相关函数做分母了,也就是 tcf(delt)=tcf(delt)+tcf(ntel)*tcf0(t)/ tcf0(t)*tcf0(t),可是输出错误,请问,我这样写为什么不可以呢?应该怎样写?
do i=1,10000 (粒子数)
     do ntel=1,npoints (共有npoints个数据)
          if (确定timeorigin的条件) then
             t0=t0+1
             time0(t0)=ntel
             tcf0(t0)=tcf(ntel)         
             do t=1,t0
                  delt=ntel-time0(t)+1
                  tcf(delt)=tcf(delt)+tcf(ntel)*tcf0(t)/ tcf0(t)*tcf0(t)
             enddo
           endif
    enddo
enddo
do ical=1,npoints
tcf(ical)=tcf(ical)/(npoints*ntime(t))
write(*,*) ical, tcf(ical)
enddo
3楼2010-12-07 15:57:21
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

astringent

铜虫 (著名写手)

问题已解决,请结贴
4楼2010-12-08 06:34:00
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 今日不放榜?网传国自然预计 8 月 27 日可查结果 +17 医学老男孩 2026-08-20 22/1100 2026-08-25 15:36 by 医学老男孩
[基金申请] 某些机构,以效率低为荣,以效率低作为存在感 +8 yuleib84 2026-08-25 9/450 2026-08-25 15:23 by Weed114114
[基金申请] 没有任何消息-是不是就凉了 +9 图啦图啦 2026-08-24 10/500 2026-08-25 11:59 by 南海小哥
[基金申请] 人气不行了 +11 fansofjerry 2026-08-21 11/550 2026-08-25 11:04 by 孤独的英雄6
[教师之家] 导师吐槽:我怎么摊上了这么个极品研究生! +3 苏东坡二世 2026-08-23 3/150 2026-08-25 10:35 by shisan1313
[基金申请] 2026年的国家社科基金项目通讯评审的新规则与新动向、新挑战 +5 process2012 2026-08-23 7/350 2026-08-25 09:42 by huixian257
[基金申请] 今天基金会出结果吗?20260819 +17 kkkl_v 2026-08-19 18/900 2026-08-25 09:41 by windflowerwy
[基金申请] 我面上完蛋了 +13 且听虎啸 2026-08-20 14/700 2026-08-25 09:10 by mrkang
[基金申请] 只有每年这种时候来逛逛小木虫 +25 yaoyewhu2008 2026-08-20 27/1350 2026-08-25 08:01 by Equinoxhua
[基金申请] 明天应该可查了!? +5 chengyan1220 2026-08-23 5/250 2026-08-24 23:28 by 我4大白菜
[基金申请] 能否退出参与的面上项目解除限项 +21 koalala 2026-08-24 24/1200 2026-08-24 19:25 by 家与远方
[基金申请] filecode,4个jtjc了 +14 ziyangfang 2026-08-19 17/850 2026-08-24 18:37 by 哈哈蛤?
[基金申请] 建议基金发布提前给出明确的时间点 +13 kulium 2026-08-21 16/800 2026-08-24 16:27 by superceng
[基金申请] 让我中一个面上吧! +13 大萍1987 2026-08-20 16/800 2026-08-24 10:23 by 太傻了
[基金申请] 时间戳今天,20号变了 +5 archvillain 2026-08-20 5/250 2026-08-22 06:12 by hui_daxiao
[基金申请] 看来今天不会放榜了? +8 chengyan1220 2026-08-21 11/550 2026-08-21 17:52 by dcqxinyang
[基金申请] 时间戳又变了 +13 wuchongjun 2026-08-20 19/950 2026-08-21 17:21 by 紫杉醇
[基金申请] 应该是下周三26日公布了吧? +4 哈哈蛤? 2026-08-21 4/200 2026-08-21 10:58 by Vivilian
[基金申请] 今天放榜没戏了吧 +9 yuleib84 2026-08-19 11/550 2026-08-21 10:06 by gltch
[基金申请] 基金啊基金 +4 longfie172 2026-08-20 4/200 2026-08-21 08:58 by mark mao
信息提示
请填处理意见