| 查看: 277 | 回复: 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 ] |
回复此楼» 猜你喜欢
今年审到国自然15份，谈谈感受
已经有8人回复
-
上海大学实验技术岗位非升即走
已经有7人回复
-
2026博士还有哪些学校有名额
已经有7人回复
-
光量子物理方向 博士招生 1人（2026.09）
已经有4人回复
-
西南大学考核制博士
已经有8人回复
-
这年头没有找到涵评专家，还有中面上的可能吗
已经有16人回复
-
26/27博士推荐
已经有3人回复
-
找博士生导师
已经有4人回复
-
重磅！青年科学基金项目（C类）资助增幅预计超过50%
已经有9人回复
-
售SCI一区T0P文章，我:8.O.5.5.1.O.5.4,科目齐全,可+急
已经有4人回复
lovein
金虫 (小有名气)
- 应助: 1 (幼儿园)
- 金币: 1419.6
- 散金: 4
- 红花: 1
- 帖子: 254
- 在线: 38.3小时
- 虫号: 534931
- 注册: 2008-03-28
- 性别: MM
- 专业: 其他无机非金属材料
2楼2009-10-26 13:39:16
精华III:
