[½»Á÷] ¡¾ÌÖÂÛ¡¿¸ø¸öÄ£ÄâÆøÌåÀ©É¢µÄMATLAB´úÂë ÒÑÓÐ4È˲ÎÓë

function gas_diffusive
clc;
rand('state',0);set(gcf,'Doublebuffer','on');
N=300;M=120;h=2;xp=round(N/2);D=round(M/12);
B=ones(M,N);B(1:h,:)=0;B(:,1:h)=0;B(:,N-h+1:N)=0;B(M-h+1:M,:)=0;B(:,xp:...
    +h-1)=0;
B(round(M/2)-D+1:round(M/2)+D,xp:xp+h-1)=1;
S=zeros(M,N);S(h+1:M-h,h+1:xp-1)=rand(M-2*h,xp-h-1)>0.5;
Tc=700;C=min(1-S,B);subplot(211);asa=imshow(C,[]);
ti=title('time=0');
rL=sum(sum(S(:,1:xp)))/((xp-h)*(M-2*h));xR=sum(sum(S(:,xp+h:N)))/((N-xp...
-2*h+1)*(M-2*h));
subplot(212),p1=plot(1,rL,'r-');xlim([1,Tc]);hold on;p2=plot(1,xR,'k--');
xlabel('ʱ¼ä');ylabel('ʱ¼ä');legend('×ó²àÃܶÈ','ÓÒ²àÃܶÈ');
for k=1:Tc
    p=mod(k,2);x1=h+1+p:2:xp-3+p;y1=h+1+p:2:M-h-2+p;x2=xp-1+p:2:xp+h-1+p;
    y2=round(M/2)-D+2:2:round(M/2)+D-1;x3=xp+h+1+p:2:N-h-2+p;y3=h+1+p:2:...
M-h-2+p;
    St=S;St(y1,x1)=S(y1+1,x1+1);St(y1+1,x1+1)=S(y1,x1);St(y1,x1+1)=S(...
    y1+1,x1);
    St(y1+1,x1)=S(y1,x1+1);St(y2,x2)=S(y2+1,x2+1);St(y2+1,x2+1)=S(y2,x2);
    St(y2+1,x2)=S(y2,x2+1);St(y2,x2+1)=S(y2+1,x2+1);St(y3,x3)=S(y3+1,x3+1);
    St(y3+1,x3+1)=S(y3,x3);St(y3+1,x3)=S(y3,x3+1);St(y3,x3+1)=S(y3+1,x3);
    St=changep(St,S,x1,y1);St=changep(St,S,x2,y2);St=changep(St,S,x3,y3);
    S=St;
    set(ti,'string',['time=',num2str(k)]);set(asa,'CData',min(1-S,B));
    xL(k+1)=sum(sum(S(:,1:xp)))/((xp-h)*(M-2*h));
    xR(k+1)=sum(sum(S(:,xp+h:N)))/((N-xp-2*h+1)*(M-2*h));
    set(p1,'XData',1:k+1,'YData',xL);set(p2,'XData',1:k+1,'YData',xR);
    pause(0.2);
end

function St=changep(St,S,x,y)
Sr1=S(y+1,x)+S(y,x);Sc1=S(y,x+1)+S(y+1,x+1);Sr2=S(y,x+1)+S(y,x);
Sc2=S(y+1,x)+S(y+1,x+1);[p,q]=find(Sr1==1&Sc1==1&Sr2==1&Sc2==1);
for k=1:length(p)
    St(y(p(k)),x(q(k)))=S(y(p(k))+1,x(q(k)));St(y(p(k))+1,x(q(k)))...
        =S(y(p(k))+1,x(q(k)));
    St(y(p(k)),x(q(k))+1)=S(y(p(k))+1,x(q(k))+1);St(y(p(k))+1,x(q(k))+1)...
        =S(y(p(k))+1,x(q(k))+1);
end
   
   
   

Originally posted by zyj8119 at 2010-12-29 22:36:45:
[code]function gas_diffusive
clc;
rand('state',0);set(gcf,'Doublebuffer','on');
N=300;M=120;h=2;xp=round(N/2);D=round(M/12);
B=ones(M,N);B(1:h,=0;B(:,1:h)=0;B(:,N-h+1:N)=0;B(M-h+1:M,=0;B(:, ...

Originally posted by g1w2h31980 at 2010-12-31 21:31:40:
ÇëÎʺÍÐÅÏ¢À©É¢ÓÐʲô¹ØÁªÂð£¿

±¾ÎÄÀ´×Ô: Сľ³æÂÛ̳ http://muchong.com/bbs/viewthread.php?tid=2750503&fpage=1