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 views (last 30 days)
Show older comments
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 Comments
Answers (1)
Walter Roberson
on 24 Apr 2013
Edited: Walter Roberson
on 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 Comments
Walter Roberson
on 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.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!