My code primarily depends on the diameter or "d" as it is defined, however changing it to a smaller value beyond 0.001, all the results just go to infinity. How do i get around the rounding issue (it seems) ?

Asked by Neel Singh

Neel Singh (view profile)

on 15 Jun 2015
Latest activity Commented on by Neel Singh

Neel Singh (view profile)

on 15 Jun 2015
clear all clc rho = 2000; Cd = 0.44; %d = input('input diameter'); d = 0.01; %d = 9.99*10^-4; A = (d/2)^2*pi; %m = input('Enter multiple of mass: ')*2.8*10^-10; %m = 2*10^-10; m = rho * (4/3)*pi* (d/2)^3; dt = 0.01; N = 100; t = 0:dt:N*dt; T = dt * N; Sy = 0 ;
Vx(1) = 0; Ux = 0; Uy = 1; Sx(1) = 0;
for i = 1:N
Vy(1) = 0;
dSy(1) = 0;
if Sy <= 0.5
%For y direction
Vy(i+1) = (Cd*A*rho*(Uy-Vy(i))^2*dt)/(2*m) + Vy(i);
Ay(i) = (Vy(i+1)-Vy(i))/dt;
dSy(i) = ((Vy(i+1))^2 - (Vy(i))^2 )/ (2*Ay(i));
%disp(sum(Sy))
Sy = cumsum(dSy);
%For x direction
Vx(i+1) = Cd*A*rho/(2*m)*(Ux-Vx(i))^2*dt+Vx(i);
Ax(i) = (Vx(i+1)-Vx(i))/dt;
dSx(i) = (Vx(i)*dt);
Sx = cumsum(dSx);
else
Uy = 0;
Ux = 1;
Deacc = Cd*A*rho/(2*m)*(Uy-Vy(i))^2*dt
%For y direction
Vy(i+1) = Vy(i) - Deacc ;
%Vy(i+1) = Vy(i)- Vy(i-1);
Ay(i) = (Vy(i+1)-Vy(i))/dt;
dSy(i) = (Vy(i)*dt);
%disp(sum(Sy))
Sy = cumsum(dSy);
%For x direction
Vx(i+1) = Cd*A*rho/(2*m)*(Ux-Vx(i))^2*dt+Vx(i);
Ax(i) = (Vx(i+1)-Vx(i))/dt;
dSx(i) = ((Vx(i+1))^2 - (Vx(i))^2 )/ (2*Ax(i));
Sx = cumsum(dSx);
end
end
x=Sx;
y=Sy;
figure (1)
plot(x,y);
hold on;
h=plot(x(1),y(1),'r*');
This "d" value gives answers well, this is how its supposed to work. however changing it to 0.001 will give infinity to Vy, Vy, Sy, Sx and i cant figure out why. The only thing i found online is that it might be a round off from matlab at some point, driving everything divided by 0 to infinity. How can I work around this, so I can use small values of d ?
Thanks for any help,
Im new to this forum, so excuse my mistakes.

Tags

No tags entered yet.

Answer by Jonathan Poirier

Jonathan Poirier (view profile)

on 15 Jun 2015

First, never use i in a loop. Matlab use i for a complex number (Try to type i in command window >> i ans =
0.0000 + 1.0000i
Second, I never try to understand yous matematical equation, but -1.32e239 look inf, no?
If you really want to figure there big values, divide (by exemple) 1 e100 your value and in your y axis :
xlabel('data 1e100');

Neel Singh

Neel Singh (view profile)

on 15 Jun 2015
Hi thanks for the suggestions, however those are not the cause of the problem.

Katalin (view profile)

on 15 Jun 2015

I think it's not a rounding issue. In Vy(i+1) = (Cd*A*rho*(Uy-Vy(i))^2*dt)/(2*m) + Vy(i); The fraction (Cd*A*rho*dt)/(2*m) is either larger or smaller than 1 depending on the variable “d”.
It is equal to 1, when d = Cd*dt/(4/3) = 0.0033. When you are trying numbers below d = 0.0033 the system has a complete different behavior and it goes to extreme very fast.

Neel Singh

Neel Singh (view profile)

on 15 Jun 2015
Hi thanks for your answer, what you are saying seems to be true. Have you any suggestions or work arounds? I need to be working with d in 10^-6 to 10^-9 range, which is impossible at the moment.