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) ?

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.

Answer 1

[Jonathan Poirier]

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]

Hi thanks for the suggestions, however those are not the cause of the problem.

Answer 2

[Katalin]

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]

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.