i am unable to get a plot .....as my plot fig is coming empty .....please tel me if anyone notice my mistake i paste the code
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
x(1)=-4; %initial position of particle in x direction
y(1)=1; %initial position of particle in y direction
nsteps=100;
x(i-1)=zeros(1,nsteps);
y(i-1)=zeros(1,nsteps);
R(i-1)=zeros(1,nsteps);
R(1)=sqrt(x(1)^2+y(1)^2);
p=1780;
f=1.2;
dp=100e-6;
m = 15.11e-4; %viscosity of fluid
Ux(1)=-0.986; %initial velocity in horizontal direction
Uy(1)=-0.00503; %initial velocity in vertical direction.
Ux(i-1)=(-1)*(1-1./(R(i-1)).^3+3.*y(1).^2./(2.*(R(i-1)).^5));
Uy(i-1)=3.*x(i-1).*y(i-1)./(2.*R(i-1).^5);
Ux=zeros(1,nsteps+1);
Uy=zeros(1,nsteps+1);
c=p*dp/f;
i=2:1:101;
RX(i-1)=Ux(i-1)*dp/m; %reynolds number in x direction
RY(i-1)=Uy(i-1)*dp/m; %reynolds number in y direction
ts=p*dp^2/(18*m);
cdu(1)=24/RX(1);
cdv(1)=24/RY(1);
if(RX(i-1) >= 3.0)
cdu(i-1)=(24/RX(i-1))+(4/RX(i-1)^(1/3));
elseif(RX(i-1) >= 0.5)
cdu(i-1)=(24/RX(i-1))+(3.6/RX(i-1)^(0.313));
elseif(RX(i-1) >= 0.1)
cdu(i-1)=(24/RX(i-1))+4.5;
elseif(RX(i-1) >=0.0001)
cdu(i-1)=24/RX(i-1);
end
while i>2;
h1=(1/(cdu(i-1)*RX(i-1)))*(c)*(-1.33);
h2=(1/(cdu(i-1)*RX(i-1)))*c*(-1.33);
h3=(1/(cdu(i-1)*RX(i-1)))*c*(-1.33);
h4=(1/(cdu(i-1)*RX(i-1)))*c*(-1.33);
j=ts/10;
x(i-1)=j/3*(h1+(4*h2)+(2*h3)+h4);
end
if(RY(i-1) >= 3.0)
cdu(i-1)=(24/RY(i-1))+(4/RY(i-1)^(1/3));
elseif(RY(i-1) >= 0.5)
cdu(i-1)=(24/RY(i-1))+(3.6/RY(i-1)^(0.313));
elseif(RY(i-1) >= 0.1)
cdu(i-1)=(24/RY(i-1))+4.5;
elseif(RY(i-1) >=0.0001);
g=9.81; %acceleration due to gravity
while i>2;
k1=sqrt((dp*f*g)/(3*cdv(i-1)*p));
k2=sqrt((dp*f*g)/(3*cdv(i-1)*p));
k3=sqrt((dp*f*g)/(3*cdv(i-1)*p));
k4=sqrt((dp*f*g)/(3*cdv(i-1)*p));
j=ts/10;
y(i-1)=j/3*(k1+(4*k2)+(2*k3)+k4);
end
end
plot(x,y)
0 Commenti
Risposte (1)
Walter Roberson
il 24 Apr 2013
Modificato: Walter Roberson
il 24 Apr 2013
You have
i=2:1:101;
but in some places (such as the "if" statements) you treat "i" as if it was a scalar instead of a vector.
When you have a construct such as
if R(i) > 5
... do something ...
end
when "i" is a vector, it does not mean that the code should be done only for the values of "i" for which R(i) > 5. In MATLAB, it would instead being interpreted as doing the code only if all R(i) > 5 .
Use a "for" loop, or study logical indexing.
2 Commenti
Walter Roberson
il 25 Apr 2013
You initialize x(1) and y(1) first, and then assign zeros to x and y, so x and y will be all 0. Then you calculate R(1) based on x(1) and y(1), so R(1) will be 0. You divide by R (all 0) so you get NaN. Everything goes downhill from there.
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!