| 查看: 260 | 回复: 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 ] |
回复此楼» 猜你喜欢
求调剂 有机化学考研356分
已经有4人回复
-
0703 化学 求调剂，一志愿山东大学 342 分
已经有4人回复
-
086000生物与医药调剂
已经有7人回复
-
305求调剂
已经有7人回复
-
348求调剂
已经有3人回复
-
332求调剂
已经有6人回复
-
各位老师好，我的一志愿为北京科技大学085601材料专硕
已经有6人回复
-
085600 材料与化工 329分求调剂
已经有10人回复
-
356求调剂
已经有3人回复
-
一志愿武汉理工，总分321，英一数二，求老师收留。
已经有7人回复
lovein
金虫 (小有名气)
- 应助: 1 (幼儿园)
- 金币: 1419.6
- 散金: 4
- 红花: 1
- 帖子: 254
- 在线: 38.3小时
- 虫号: 534931
- 注册: 2008-03-28
- 性别: MM
- 专业: 其他无机非金属材料
2楼2009-10-26 13:39:16
