# How to fix timestep in this code?

5 views (last 30 days)
Ali Khabibulayev on 22 Apr 2020
Answered: Gifari Zulkarnaen on 22 Apr 2020
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.

Gifari Zulkarnaen on 22 Apr 2020
You can check your equation part by part to check whether your code or calculation is correct. For example, for n=2:
Ve*dm = 22500000 --> not weird (I don't know your field, so I don't know whether it is correct)
G*((Mearth*m(n-1))/(h(n-1))^2) = 3.5634e+20 --> Suspicious, is it correct? Maybe the r^2 in your formula means (radius_earth + h_rocket)^2...
And so on, check it part by part