| 查看: 1038 | 回复: 1 | |||
revisionscc铜虫 (初入文坛)
|
[求助]
小弟初学MATLAB，编了一个仿真麦克斯韦速率分布的程序，求优化识错
|
大致就是如题，我已经分析出了在分子间碰撞那个循环会对一次碰撞做两次运算，但不知道怎么解决，并且这个程序运行时间太长，希望优化一下程序，求大牛们指导。。。 以下是文件内容以及附带的.m文件。 %以完全弹性碰撞模型为基础模拟N个粒子在初始状态为同速率的情况下，经过一段时间的速度以及速率分布 % %设模拟的分子数为10000，二维方形盒子边长为10^3。初始速度为400 %a为X方向速度矩阵，b为y方向速度矩阵，c为x方向位移矩阵，d为y方向位移矩阵 % clear k=rand(1,10000); a=400*cos(2*pi*k); %Vx b=400*sin(2*pi*k); %Vy c=10^3*k; %Sx k=rand(1,10000); %重置随机变量 d=10^3*k; %Sy t=0.01; %t for n=1:1000; c=c+t*a; d=d+t*b; e=(c>=1000&c<=0); %1为碰壁分子下标 17---22 解决分子碰壁问题 f=(d>=1000&d<=0); %0为不碰壁分子 c=c-e.*a*t; %重置碰壁分子位移 d=d-f.*b*t; % a=a-2*a.*e; %碰壁分子碰撞方向速度反向 b=b-2*b.*f; % g=round(c); %开始处理分子间碰撞问题，约化碰撞半径为1 23---39 解决分子间碰撞问题 h=round(d); %四舍五入分子芯的位置 for i=1:10000 %按顺序寻找分子芯位置相同的粒子 if length(find(and(not((g-g(i))),not((h-h(i))))))~=2 %排除不碰撞分子与两分子以上的碰撞 else j=find(and(not((g-g(i))),not((h-h(i))))); %找出相互碰撞的两分子坐标 u=j(1); % v=j(2); xx=c(u)-c(v); %求分子非对心碰撞模式 yy=d(u)-d(v); % theta=atan(yy/xx); % a(u)=(a(v)*cos(theta)+b(v)*sin(theta))*cos(theta)+(a(u)*sin(theta)+b(u)*cos(theta))*sin(theta); %碰后分子速度 b(u)=(a(v)*cos(theta)+b(v)*sin(theta))*sin(theta)+(a(u)*sin(theta)+b(u)*cos(theta))*cos(theta); % a(v)=(a(u)*cos(theta)+b(u)*sin(theta))*cos(theta)+(a(v)*sin(theta)+b(v)*cos(theta))*sin(theta); % b(v)=(a(u)*cos(theta)+b(u)*sin(theta))*sin(theta)+(a(v)*sin(theta)+b(v)*cos(theta))*cos(theta); % end end end hist(sqrt(a.^2+b.^2),100) %画出速率分布直方图 |
回复此楼» 猜你喜欢
梦想
已经有7人回复
-
学科评审组评审是指会评吗？
已经有4人回复
-
面上合作单位盖章
已经有8人回复
-
面上函评意见出来了，像什么等级？
已经有15人回复
-
为什么资助数各大高校都创新高，自己申请怎么就这么难
已经有14人回复
-
小白求助 投论文要求的highlights应该如何写
已经有4人回复
-
麻烦专家们看看评委们的意见（F口面上）
已经有12人回复
-
怎么看青基中了没有啊
已经有6人回复
-
国社科又开始会评了，不知道这次命运如何
已经有11人回复
-
有多少人是今天查系统知道结果的？
已经有19人回复
» 本主题相关商家推荐: (我也要在这里推广)
2楼2018-01-06 17:38:35
投票:
10