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

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)

Answer 1

[Walter Roberson]

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.

[bina]

i done this for for loop again it is runing but not ploting figure ......just paste this on your command window you will understand.........can you suggest something else.

[Walter Roberson]

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.