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

1 次查看(过去 30 天)
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.

回答(2 个)

Jonathan Poirier
Jonathan Poirier 2015-6-15
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');

Katalin
Katalin 2015-6-15
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.
  1 个评论
Neel Singh
Neel Singh 2015-6-15
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.

请先登录,再进行评论。

标签

尚未输入任何标签。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by