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 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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.
0 commentaires
Réponses (2)
Jonathan Poirier
le 15 Juin 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');
Katalin
le 15 Juin 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.
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!