untitled1.m

来自「基于流函数法求解流体颗粒二相绕流问题的源码」· M 代码 · 共 27 行

x0=[0.2;0;0;0;0];
Yr=x0';
t0=0,tmax=20000;
while 20000-t0>=0    
[t,YY]=ode45('DYDt2',[],[],[],[t0:0.01:tmax],x0)
Sz=size(YY);
lth=Sz(1)-1;
rt=YY(:,5);
a=find(rt-fix(rt)<0.01)
Yr=[Yr;YY(a,:)]
% Yr=[Yr;YY(1:lth,:)]
if fix(YY(end,1))>=1
    x0=[2-YY(end,1);YY(end,2);-YY(end,3);YY(end,4);YY(end,5)];
elseif YY(end,1)<0
    x0=[-YY(end,1);YY(end,2);-YY(end,3);YY(end,4);YY(end,5)];
else
    if fix(YY(end,2))>=1
    x0=[1-YY(end,1);YY(end,2)-1;-YY(end,3);YY(end,4);YY(end,5)];
    elseif YY(end,2)<0
    x0=[1-YY(end,1);YY(end,2)+1;-YY(end,3);YY(end,4);YY(end,5)];
    end
end
t0=YY(end,5);
end
Yr=Yr(2:end,:);
save('data1','Yr','t')