How to fix timestep in this code?

Question

BioZ 2020-4-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/519926-how-to-fix-timestep-in-this-code

回答： Gifari Zulkarnaen 2020-4-22

Hello everyone, I am new to Matlab and working on a code...

I have written this code.

dt = 1; %time step
tmax = 200; %max time
t = 0:dt:tmax;
G = (6.67408*10^-11); %Gravitaionals constant m^3.kg^-1.s^-2
Mearth = (5.9722*10^24); %Mass of earth kg
R = 6371; %Radius of the earth 6371 km
A = 75; %Area of the rocket m^2
Cd = 0.4; %Drag coefficient for rocket 
h = 0; %Initial altitude
Mempty = 54000; %Mass of the rocket Empty kg
Minitial = 894000; %Mass of the rocket Full (With propellant) kg
Ve = 4500; %Exhaust gas speed  ms^-1
dm = 5000; %Rate of change of mass kg/s (constant)
V = zeros(1,length(t));
V(1) = 1;
 % using dx/dt ~= (x(i)-x(i-1))/dt Discretization example
 % speed(i) = speed(i-1) + time_step*force/mass; Remember
 
 m = zeros(1,length(t));
 m(1) = Minitial; %Mass of the rocket kg
 h = zeros(1,length(t));
 h(1) = 1 %Altitude
 p = zeros(1,length(t));
 p(1) = 1.225; %Density of air sea lvl kg/m^3
for n=2:length(t)
    if m(n-1) >= Mempty;
        m(n) = Minitial-(dt*dm*n);
    else m(n-1) = Mempty;
    end
    V(n) = V(n-1)+dt*(Ve*dm-G*((Mearth*m(n-1))/(h(n-1))^2)-0.5*p(n-1)*A*(V(n-1))^2*Cd)/m(n-1);
    h(n) = V(n)*n;
    p(n) = 1.225*10^(-3*h(n-1)/50000);
end     
plot(t,V)
    %V(n) = V(n-1)+dt*(Ve*dm-G*((Mearth*m(n))/r^2)-0.5*p*A*V(n-1)^2*Cd
    %V(n)=V(n-dt)*(1+t(n-1)*dt);
 
 

This code is supposed to calculate the velocity and altitude of a rocket flying straight up. However, for some reason as it got evident from the graph my V(n) values seem to be calculate for only 2 times steps and somehow ignores tmax variable? My assumption is that there is something wrong in my V(n) Equation, line 32.

For my V(n) I am using this

and the graph that I am getting is

which is clearly not what I am trying to get here.

I would really appreciate if someone could help me out to fix the V(n) in order to find the velocity.