| 查看: 1493 | 回复: 4 | |||
[交流]
matlab的蒙特卡洛模拟提高运算速率 已有1人参与
|
|
%蒙特卡洛模拟光子在水雾中的 传播 kext=2.4056; ks=2.4053; ka=0.0003; r=250*10^(-6); n0=19646104; c=pi*(r^2)*n0*ks; for m=10:10:250 L=m; s1=0; s2=0; z=zeros(1,10000); x=zeros(1,10000); y=zeros(1,10000); n=1; while(n>=1&&n<=10000) if((z(1,n)>=L)&&((x(1,n))^2+(y(1,n))^2<2.88)) s1=s1+1; n=n+1; else if((z(1,n)>=L)&&((x(1,n))^2+(y(1,n))^2>=2.88)) n=n+1; else if(rand(1)<=(ka/kext)) n=n+1; s2=s2+1; else syms t si; % t=solve(int(4.40189*10^(-8)+1.39435*10^(-8)*sin((pi/1.58425)*(si-2.33743)),si,0,t)-rand(1)*1.3792*10^(-7)); y1=4.40189*10^(-8)+1.39435*10^(-8)*sin(pi/1.58425*(si-2.33743)); R=rand(1)*1.379*10^(-7); y2=int(y1,si,0,t)-R; t=vpasolve(y2); rand1=-log(rand(1))/c; rand2=rand(1)*2*pi; rand3=t-pi/2; z(1,n)=z(1,n)+rand1*cos(rand3); x(1,n)=x(1,n)+rand1*cos(rand2)*sin(rand3); y(1,n)=y(1,n)+rand1*sin(rand2)*sin(rand3); n=n+0; end end end end fprintf('Distance:%i\n',L); fprintf('The number of photons arriving at the probe plane:%i\n',s1); end 本程序需要对10000个光子进行模拟， 若先对1000个进行模拟，出来一个L对应的是s1耗时5个小时。 那改为10000个光子，出来所有L=10:10:250对应的s1就可能得几个月了  求问如何提高运行效率  |
回复此楼» 猜你喜欢
070200求调剂，一志愿某211，288分
已经有20人回复
-
法国博士后职位
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有235人回复
-
重庆交大26年硕士生招生拟调剂通知已出！欢迎加入机器视觉与3D光学成像课题组。
已经有0人回复
-
**
已经有1人回复
-
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有3人回复
-
基底STO,薄膜SRO，XRD里面的振荡，是laue震荡还是kiessig振荡? 怎么判断？
已经有2人回复
-
B区学生调剂-兰州交通大学材料科学与工程学院
已经有8人回复
-
山西大同大学物理学专业还有调剂名额，欢迎调剂！
已经有9人回复
persisthcw
新虫 (初入文坛)
- 应助: 0 (幼儿园)
- 金币: 195.4
- 帖子: 6
- 在线: 13.4小时
- 虫号: 1619804
- 注册: 2012-02-15
- 性别: GG
- 专业: 计算机应用技术
2楼2017-11-16 13:09:57
3楼2017-11-16 15:06:44
persisthcw
新虫 (初入文坛)
- 应助: 0 (幼儿园)
- 金币: 195.4
- 帖子: 6
- 在线: 13.4小时
- 虫号: 1619804
- 注册: 2012-02-15
- 性别: GG
- 专业: 计算机应用技术
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
小木虫: 金币+0.5, 给个红包，谢谢回帖
|
可以把程序对应的算法伪代码和公式贴出来，大家看到这个才能明白，方便提出修改建议。如果循环里面每次计算是互相独立的话，也可以用并行计算来加快速度。 发自小木虫Android客户端 |
4楼2017-11-16 17:38:16
|
对应的算法及原理如下，不太会写伪代码，麻烦大神帮忙看看怎么改  算法： 从某一平行于XOY面（XOY面为探测平面）的平面向Z轴负方向发射一定光子数量的平行光，光子的生命周期终结于被吸收或逃逸至探测平面（圆形，半径的平方2.88）。 （1）从距探测平面L处的光源向发射光线，设定光源性质为1万个随机生成的均匀分布的光子构成，光子发射方向均为平行于Z轴； （2）计算光子的传输步长d； （3）判断光子是否逃逸至探测平面，成功逃逸（z(1,n)>=L(1,n))&&((x(1,n))^2+(y(1,n))^2<2.88)）则记录储存为S1并返回（1）开始下一个光子； （4）否则根据传输步长d和传输方向和计算发生碰撞的位置，。m表示第m个光子，i表示该光子第i次发生碰撞； （5）判断该光子在碰撞时是被吸收还是被散射，若被吸收，则记录储存为S2并返回（1）开始下一个光子； （6）若被散射，则重复（2），直至光子到达探测平面（Z=0）或光子被吸收。 （7）返回到达探测平面的光子数量。 设光子由Z=L（L从10m开始取，到500）平面向下出发，发生第一次碰撞之前的传播方向平行于Z轴，第一次碰撞之后的是垂直于Z轴平面上新方向与X轴的夹角，是新方向与Z轴的夹角 （1）光子的传输步长d：传输步长d是光子发生连续两次碰撞之间所经历的距离。对于光子在介质中的传播，传输步长d由下式决定：=-lnr1/c，为一个取自[0,1]区间均匀分布的随机数，c=pi*(r^2)*n0*ks; kext=2.4056; ks=2.4053; ka=0.0003; r=250*10^(-6); n0=19646104; （2）在光子与水雾粒子发生碰撞后，光子可能被水雾粒子吸收或散射。当时，光子被水雾吸收，否则光子被水雾散射，为一个取自[0,1]区间均匀分布的随机数。 （3）光子遇到水雾粒子，若产生碰撞，且未被吸收而发生散射。该光子新的运动路径由介质的散射相位函数决定。 在垂直于光子传输路径的平面上，新传播方向的方向的取值为，为一个取自[0,1]区间均匀分布的随机数。 再求得新方向与原路径的夹角，即可获得新的传播方向，这个就是散射角。令为一个取自[0,1]区间均匀分布的随机数，满足，便可求得一个确切的散射角与相对应。 求出和两个方向便得到了光子的新的传播方向。 （4）至此，使用4个取自[0,1]区间的、均匀分布的、随机产生的数字，完整地描述了光子与水雾粒子的相互作用概率模型。这4个随机数在Monte Carlo概率模型中起到了不同作用：决定了光子两次碰撞之间的传播距离，决定了碰撞的结果，与决定了散射后光子新的运动方向。 |
5楼2017-11-20 09:42:48
