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

Question

Neel Singh 2015-6-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/223875-my-code-primarily-depends-on-the-diameter-or-d-as-it-is-defined-however-changing-it-to-a-smaller

评论： Neel Singh 2015-6-15

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Jonathan Poirier 2015-6-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/223875-my-code-primarily-depends-on-the-diameter-or-d-as-it-is-defined-however-changing-it-to-a-smaller#answer_182725

在 MATLAB Online 中打开

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');

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Neel Singh 2015-6-15

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

请先登录，再进行评论。

Answer 2

Katalin 2015-6-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/223875-my-code-primarily-depends-on-the-diameter-or-d-as-it-is-defined-however-changing-it-to-a-smaller#answer_182762

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

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.

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论