| 查看: 308 | 回复: 1 | |||
| 当前主题已经存档。 | |||
lovein金虫 (小有名气)
|
[交流]
【求助】麻烦帮忙看下一个“随机游走”的程序，为什么运行时老出错！！！
|
||
|
整个程序如下： (*----------------------------------------------------------------------*) BeginPackage["RW2D`"] RW2D::usage = "RW2D[maxstep, ktry, range, startarea, datadirectory] reads the two - dimensional digital image file (text file) and exports the mean - square displacement of random walkers as a text file, 'msdlist2d.out'. 'maxstep' is the final value of the time step of the random walker. 'ktry' is the number of average (i.e., number of walkers). 'range' is the lower and upper limit of the voxel intensity of the pore space. 'startarea' is the dimension of a square as the percentage to the system size. The random walker starts to diffuse from inside the square which is located at the center of the image. 'datadirectory' is the directory or folder name where a single two-dimensional digital image file (text file) is located. " ShowTrajectory::usage = "ShowTrajectory[kk] is for the graphic display of the random walk trajectory of the 'kk' - th trial. " Begin["`Private`"] (* Reading the two - dimensional data in 'datadirectory' *) readdata2d[range_List, datadirectory_] := Module[{dataname, datasize, x, m2d, rmin = range[[1]], rmax = \ range[[2]]}, dataname = First[FileNames["*", {datadirectory}]]; datasize = Round[Sqrt[Length[ReadList[ToString[dataname], Real]]]]; m2d = Partition[ Map[If[rmin (* random walk *) (* 'xylist' is the trajectory of 'maxstep' - step from 'o2d' *) tra2d[maxstep_, o2d_List, m2d_List] := Module[{n2d = {{1, 0}, {-1, 0}, {0, 1}, {0, -1}}, now = o2d, xylist = {o2d}, outoffield = 0, xmax = ymax = Length[m2d], step, new, x, y}, ran2d := Random[Integer, {1, 4}]; For[i = 1, i (* random walk of 'ktry' - walkers *) (* 'tra2dktry' is list of the trajectory of 'ktry' - walkers *) rw2d[maxstep_, ktry_, startarea_, m2d_List] := Module[{xdatasize = ydatasize = Length[m2d], xylist, outoffield, tra2dktry = {}, failure = 0}, rano2d[lx_, ly_, sa_] := {Random[Integer, {(1 - sa/100)/2 lx, (1 + sa/100)/2lx}], Random[Integer, {(1 - sa/100)/2 ly, (1 + sa/100)/2ly}]}; o2d := (For[i = 1; {x0, y0} = rano2d[xdatasize, ydatasize, startarea], m2d[[x0, y0]] != 0, i++, {x0, y0} = rano2d[xdatasize, ydatasize, startarea]]; {x0, y0}); Do[{xylist, outoffield} = tra2d[maxstep, o2d, m2d]; AppendTo[tra2dktry, xylist]; AddTo[failure, outoffield], {i, 1, \ ktry}]; Return[{tra2dktry, failure}]; ]; (* calculation of the mean - square displacement *) msd2d[data_List] := Module[{ktry = Dimensions[data][[1]], maxstep = Dimensions[data][[2]] - \ 1, d1, d2, d3, dd, ii}, d1 = Transpose[data, {2, 1}]; d2 = Map[(# - First[d1]) &, d1]; d3 = Transpose[Map[(#^2) &, d2]]; dd = Apply[Plus, Transpose[Apply[Plus, d3]/ktry]]; ii = Range[0, maxstep]; msdlist2d = Transpose[{ii, dd}]; Return[msdlist2d]; ]; RW2D[maxstep_, ktry_, range_List, startarea_, datadirectory_] := Module[{failure, msdlist2d}, (* reading the two - dimensional data in 'datadirectory' *) m2d = readdata2d[range, datadirectory]; (* execution of the random walk *) {tra2dktry, failure} = rw2d[maxstep, ktry, startarea, m2d]; Print["Failure = ", failure, " / Success = ", ktry]; (* calculation of the mean - square displacement of 'ktry' walkers *) msdlist2d = msd2d[tra2dktry]; (* file output of the mean - square displacement *) (* The file name is 'msdlist2d.out' *) Export["msdlist2d.out", N[msdlist2d], "Table"]; (* graphic output of the mean - square displacement *) ListPlot[msdlist2d, PlotJoined -> True, AxesLabel -> {"time \n step", "mean-square \n displacement"}]; ]; (* graphic output of the trajectory *) (* Show the trajectory of a single walker as an example. *) (* RGBColor[] of 0, 1, 2 are colors of the pore, obstacle, and trajectory, respectively. You can change the color as you like. *) gtra2d[tra2d_List, m2d_List] := Module[{color0 = RGBColor[0.9, 0.9, 0.9], color1 = RGBColor[0.5, 0.5, 0.5], color2 = RGBColor[0.8, 0.8, 0.8], plmin = 1 - 2, plmax = Length[m2d] + 1}, tra2dshift = Map[(# - {0.5, 0.5}) &, tra2d]; tra2dm2d = (aa2 = m2d; Map[(aa2[[#[[1]], #[[2]]]] = 2) &, tra2d]; aa2); Show[ Graphics[ RasterArray[tra2dm2d /. {0 -> color0, 1 -> color1, 2 -> color2}]], Graphics[{Line[Map[Reverse, tra2dshift]], AbsolutePointSize[5.0], Point[Reverse[First[tra2dshift]]], Point[Reverse[Last[tra2dshift]]]}], PlotRange -> {{plmin, plmax}, {plmin, plmax}}, AspectRatio -> 1]; ]; (* Only for the 'kk' - th trial, the trajectory of the random walk is shown graphically *) ShowTrajectory[kk_] := gtra2d[tra2dktry[[kk]], m2d]; End[]; EndPackage[]; (*----------------------------------------------------------------------*) (* example of usage *) (*----------------------------------------------------------------------*) (* As noted below, first of all, you put the command '<< RW2D.m'. *) (* Then you can put the command set, 'RW2D[500, 200, {0, 34000}, 10, "DATAsingle"];' and 'ShowTrajectory[2];'. \ *) (* 'RW2D[maxstep, ktry, range, startarea, datadirectory]' is the command for the performance of the random walk and the calculation of mean - square displacement. *) (* The command 'ShowTrajectory[kk]' is for the graphic display of the trajectory of the 'kk' - th trial. *) (* In the example below, a trajectory of the second trial will be shown graphically. *) (* In[1] := << RW2D.m *) (* In[2] := RW2D[500, 200, {0, 34000}, 10, "DATAsingle"]; ShowTrajectory[2]; \ *) 不知道为啥，运行的时候老是出错，我又看不大懂这程序，请求帮忙 [ Last edited by lovein on 2009-10-26 at 13:35 ] |
回复此楼» 猜你喜欢
售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急
已经有5人回复
-
售SCI一区T0P文章，我:8O.55.1.O.5.4,科目齐全,可+急
已经有4人回复
-
上海工程技术大学 激光智能制造课题组 2027级博士研究生招生
已经有6人回复
-
课题组招2027级博士 上海工程技术大学 激光智能制造方向
已经有6人回复
-
现代”学阀”该如何界定
已经有14人回复
-
求合成方法
已经有8人回复
-
师弟论文见刊大半年才想起来申请专利，还能抢救一下吗？
已经有4人回复
-
有机合成以后会不会被AI改变？做科研的虫友怎么看
已经有8人回复
-
我的奶奶
已经有3人回复
lovein
金虫 (小有名气)
- 应助: 1 (幼儿园)
- 金币: 1419.6
- 散金: 4
- 红花: 1
- 帖子: 254
- 在线: 38.3小时
- 虫号: 534931
- 注册: 2008-03-28
- 性别: MM
- 专业: 其他无机非金属材料
2楼2009-10-26 13:39:16