| 查看: 259 | 回复: 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 ] |
回复此楼» 猜你喜欢
0703 化学 求调剂，一志愿山东大学 342 分
已经有3人回复
-
0703化学调剂，求导师收
已经有9人回复
-
0856求调剂
已经有3人回复
-
085701求调剂初试286分
已经有4人回复
-
求调剂
已经有8人回复
-
316求调剂
已经有6人回复
-
315求调剂
已经有5人回复
-
一志愿北京化工大学材料与化工（085600）296求调剂
已经有11人回复
-
332求调剂
已经有4人回复
-
本科新能源科学与工程，一志愿华理能动285求调剂
已经有11人回复
lovein
金虫 (小有名气)
- 应助: 1 (幼儿园)
- 金币: 1419.6
- 散金: 4
- 红花: 1
- 帖子: 254
- 在线: 38.3小时
- 虫号: 534931
- 注册: 2008-03-28
- 性别: MM
- 专业: 其他无机非金属材料
2楼2009-10-26 13:39:16
