Value of variable not changing.

2 次查看（过去 30 天）
Sahil Kumar
Sahil Kumar 2021-8-8
评论： Sahil Kumar 2021-8-8
In the following code the value mdot is set to change once time t reaches a certain value but when i run it it does not change. Could anyone explain why this is happening?
clear all;
%Projectile motion with drag and thrust solution using RK4
%variation of density with height
%x'' = FDx + Tx, y'' = -g + FDy + Ty Equations to be solved
g=9.807;
vel=1000; th_deg=60;m=80; %input
cd=0.75; rho=1.225; S=2.176*10^(-4);
k=0.5*rho*cd*S; % constant k used in drag force F=kv^2
mdot = 0.2; %burn rate
Isp1 = 280; Isp2=480; %stage wise specific impulse
Isp=Isp1;
ms1=20; %structural mass of first stage
x0=0; y0=0; %initial condition
t0=0; tf=10000; %time span
vx=vel*cosd(th_deg); %velocity along x
vy=vel*sind(th_deg); %velocity along y
beta=deg2rad(th_deg);%angle made by drag force with horizontal
%transforming second order differential equation to first order
%x'=u represented by fx
%u'=-(k/m)*(u^2+v^2)*cos(beta)+mdot*g*Isp*cos(beta)/m represented by fu
%y'=v represented by fy
%v'=-g-(k/m)*v*(u^2+v^2)*sin(beta) +mdot*g*Isp*sin(beta)/m represented by fv
fx=@(t,x,u) u;
fu=@(t,x,u,v)(-(k/m)*(u^2+v^2)*cos(beta)+(mdot*g*Isp*cos(beta))/m);
fy=@(t,y,v) v;
fv=@(t,y,u,v)(-g-(k/m)*(u^2+v^2)*sin(beta)+(mdot*g*Isp*sin(beta))/m);
t(1)=0;
x(1)=0;y(1)=0;
u(1)=vx;v(1)=vy;
h=0.01;
N=ceil((tf-t(1))/h);
for j=1:N
t(j+1)=t(j)+h;
k1x=fx(t(j),x(j),u(j));
k1y=fy(t(j),y(j),v(j));
k1u=fu(t(j),x(j),u(j),v(j));
k1v=fv(t(j),y(j),u(j),v(j));
k2x=fx(t(j)+h/2,x(j)+h/2*k1x,u(j)+h/2*k1u);
k2y=fy(t(j)+h/2,y(j)+h/2*k1y,v(j)+h/2*k1v);
k2u=fu(t(j)+h/2,x(j)+h/2*k1x,u(j)+h/2*k1u,v(j)+h/2*k1v);
k2v=fv(t(j)+h/2,y(j)+h/2*k1y,u(j)+h/2*k1u,v(j)+h/2*k1v);
k3x=fx(t(j)+h/2,x(j)+h/2*k2x,u(j)+h/2*k2u);
k3y=fy(t(j)+h/2,y(j)+h/2*k2y,v(j)+h/2*k2v);
k3u=fu(t(j)+h/2,x(j)+h/2*k2x,u(j)+h/2*k2u,v(j)+h/2*k2v);
k3v=fv(t(j)+h/2,y(j)+h/2*k2y,u(j)+h/2*k2u,v(j)+h/2*k2v);
k4x=fx(t(j)+h,x(j)+h*k3x,u(j)+h*k3u);
k4y=fy(t(j)+h,y(j)+h*k3y,v(j)+h*k3v);
k4u=fu(t(j)+h,x(j)+h*k3x,u(j)+h*k3u,v(j)+h*k3v);
k4v=fv(t(j)+h,y(j)+h*k3y,u(j)+h*k3u,v(j)+h*k3v);
x(j+1)=x(j)+h/6*(k1x+2*k2x+2*k3x+k4x);
y(j+1)=y(j)+h/6*(k1y+2*k2y+2*k3y+k4y);
u(j+1)=u(j)+h/6*(k1u+2*k2u+2*k3u+k4u);
v(j+1)=v(j)+h/6*(k1v+2*k2v+2*k3v+k4v);
m=m-mdot*h;
%condition for burn time for first stage
switch(t(j+1))
case 120
mdot=0;
m=m-ms1;%first stage separation
%no thrust phase
%second stage ignition
case 240
mdot=0.8;
Isp=Isp2;
%condition for burn time for second stage
case 360
mdot=0;
end
%Density variation with height
if(y(j+1)<=11000)
T=15.04 - .00649*y(j+1);
p=101.29*((T+273.1)/288.05)^5.256;
end
if(y(j+1)>11000 && y(j+1)<=25000)
T=-56.46;
p=22.65*exp(1.73-0.000157*y(j+1));
end
if(y(j+1)>25000)
T=-131.21+0.00299*y(j+1);
p=2.488*((T+273.1)/216.6)^(-11.388);
end
rho=p/(0.2869*(T+273.1));
k=0.5*rho*cd*S;
beta = atan(v(j+1)/u(j+1));
%condition to stop simulation when particle touches earth
if(y(j+1)<0)
break;
end
end
hmax = max(y);
rmax = max(x);
plot(x,y);
grid on;
xlabel('X');
ylabel('Y');

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

采纳的回答

Walter Roberson
Walter Roberson 2021-8-8
Ran in:
Your code assumes that if you add up enough 0.01 that the result will be an exact integer such as 120 or 240. However,
format long
t=0; for K = 1 : 12000; t = t + 0.01; end
t
t =
1.200000000000245e+02
t - 120
ans =
2.448530267429305e-11
You might be thinking this is a bug, but it is normal.
MATLAB uses IEEE 754 double precision to store (most) numbers. IEEE 754 represents numbers as a 54 bit number, divided by 2 to a power. In order to represent 1/10th exactly, there would have to be some integers (P,Q) such that P/(2^Q) = 1/10 exactly. That would require that 10*P = 2^Q for integers P, Q -- there would have to be an integer multiple of 10 that exactly equalled 2 to a power. But other than 2^0 = 1, powers of 2 must end in a non-zero digit: 2, 4, 8, 6, 2, 4, 8, 6 and so no such P, Q exist.
Which is another way of saying that double precision cannot exactly represent 1/10, for the same reason that finite decimal cannot exactly represent 1/3 . 0.3333 + 0.3333 + 0.3333 = 0.9999 and 0.3334 + 0.3334 + 0.3334 = 1.0002 so you can see that if you try to represent 1/3 as 0. 3 repeated for any finite number of decimal places, or 0 . 3 repeated then final 4, then when you add three of them you cannot get exactly 1.0 back in decimal.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

产品


版本

R2015a

Community Treasure Hunt

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

Start Hunting!

Translated by