主题：光滑算子编程
wangxinling
[专家分：0] 发布于 2012-07-14 15:05:00
物理模型：上下两层，第一层速度为1600m/s,第二层速度为2500m/s。5*5网格。
用三点平均算子,即v(i,j)[n+1]=（v(i-1,j)[n]+v(i,j)[n]+v(i+1,j)[n]）/3,其中v(i,j)[n]是网格节点（xi,zj）处经过n次光滑以后的速度值。
用fortran编出此程序，我编过了，但和手算的就是不一样，这是我编的程序，求大侠们帮忙。
dimension v(1000,1000),vvv(1000,1000)
real v,vvv
do i=1,5
do j=1,2
v(i,j)=1600
enddo
enddo
do i=1,5
do j=3,5
v(i,j)=2500
enddo
enddo
do 10 k=1,50
do 20 i=2,4
do 30 j=2,4
vvv(i,j)=(v(i,j-1)+v(i,j)+v(i,j+1))/3
30 continue
20 continue
v(i,j)=vvv(i,j)
10continue
最后更新于：2012-07-14 15:10:00
回复列表 （共9个回复）
沙发
Maple42 [专家分：50] 发布于 2012-07-15 14:33:00
我想你倒数第二行是想将vvv的值全部传给v吧。你的逻辑有问题，自己仔细演算一遍就知道了。
板凳
wangxinling [专家分：0] 发布于 2012-07-15 16:30:00
是的，新算出的vvv（i,j）赋值给v(i,j)，因为第二次光滑处理时，所用的速度值为第一光滑算出来的，我确实找不到错误。大侠帮忙啊。
3 楼
xiaoyuan24 [专家分：100] 发布于 2012-07-17 03:40:00
那我觉得你这个写法有问题，v(i,j)=vvv(i,j)，你想把vvv的值全部传给v，但是vvv(i,j)只是一个数，你只传了一个点的值而已，并没有将v里面的9个点都更新。因为在你的第一，第二个循环结束以后i和j都是具体的定值，比如第一次i，j，循环结束后，i=4,j=4，（有可能i=5,j=5,总之是个具体数字）。所以你的v(i,j)=vvv(i,j)其实只相当于v(4,4)=vvv(4,4)，并没有达到全部传进去的目的。 所以你将v(i,j)=vvv(i,j)改成v=vvv,应该就可以了。你可以试试自己调调，我并没有调试，也许会有小错误。你试完再告诉我结果是不是这样。
4 楼
wangxinling [专家分：0] 发布于 2012-07-17 08:49:00
试了，出现好几个零，结果不对。结果应该是在1600-2500之间变化的数。
5 楼
wangxinling [专家分：0] 发布于 2012-07-17 08:50:00
试了，出现好几个零，结果不对。结果应该是在1600-2500之间变化的数
6 楼
xiaoyuan24 [专家分：100] 发布于 2012-07-18 01:24:00
当然会有零，你的程序这样写就是会有0，你自始自终只更新了vvv(2:4,2:4),那其他的点的速度不是0是什么？？你程序逻辑有问题。问题表述也不是很清楚，你的速度如何更新你需要描述清楚。
第一你的数组定义那么大。用到的却只有（1:5,1:5）。第二，你是4*4网格还是5*5网格，如果是后者，那会有36个节点，你只赋值1:5,6你根本没赋值。第三，你更新的时候只是2:4，然后6也没更新。第四，你说N次光滑，你的程序里面第一个循环比如先更新vvv（2,2），接着更新vvv(2,3)的时候你是用v(i,j-1)+v(i,j)+v(i,j+1)，还是用vvv(i,j-1)+v(i,j)+v(i,j+1)？如果是后者那你的v=vvv就放错地方。
所以以上几点问题在你的问题描述里要说清楚，尤其是第四点，具体是如何更新的。你自己一步步把你的算法思路理一下，看看每一步是不是达到你想要的效果。
7 楼
xiaoyuan24 [专家分：100] 发布于 2012-07-18 01:37:00
s
8 楼
wangxinling [专家分：0] 发布于 2012-07-18 12:03:00
其他点的速度就是最开始赋值语句给的值啊，如v(1,1)=1600
1，数组定义小了就出错，大点没事，所以编程我习惯定义大点。这样是不是有问题啊？
2，4*4网格，i=1,5；j=1,5，共16个网格，25个节点,
4，用的前者。把i=1,5：j=2,4之间所有节点算完后，作为旧的值，进行下一部光滑进行新值计算。
不知道这次说清楚了吗？
v(i,j)=vvv(i,j)，这里出错了，改为下面形式就行了，问题算是解决了。谢谢你的建议。高手就是高手。
[em2]
do i=1,5
do j=2,4
v(i,j)=vvv(i,j)
enddo
enddo
9 楼
xiaoyuan24 [专家分：100] 发布于 2012-07-18 13:00:00
恭喜恭喜，我也就小打小闹不是什么高手，问题解决了就行，呵呵。
我来回复
您所在位置：