MATLAB: Plot graph using Tsiolkovsky’s rocket equations.

differential equationsdiscretizationgraph

Implement a program that can be used to calculate the speed and position of a rocket whose motion is expressed through the differential equation. A plot that shows a test of your implementation compared to the solution of Tsiolkovsky’s rocket equations. Use data from the tables below but set G and CD to zero.A plot of how the velocity changes in the first 1000 s of the rocket’s flight according to the solution of (1) using the parameters stated below.
m*dV/dt = ve*dm/dt -GMm/r^2 -0.5*rho*A*V^2*cd(1)
Initial values (At t=0)
V speed =0 ms-1
r radius =R
h altitude =0 km
this is my code so far.
clear;
clf;
%TIME CONDITION
End_time = 1000;
Time_step = 1;
%INTIAL CONDITION
Intial_speed = 0;
Intial_altitude = 0;
%FIXED PARAMETERS
G = 6.67408e-11;        % Universal Gravitational Constant
M = 5.9722e24;          % Mass of the earth
R = 6371e3;             % Radius of the earth(m)
A = 75;                 % Area of the rocket 
Cd = 0.4;               % Drag Coefficient for the rocket
m_e = 54000;            % Empty mass of the rocket
m_0 = 894000;           % Intial mass of the rocket
v_e = 4500;             % Exhaust gas speed
dmdt_f = 5000;          % Rate of change of mass with fuel left
dmdt_o = 0;             % Rate of change of mas with all fuel being used
%STTING ARRAYS
t = 0:Time_step:End_time;
Speed = zeros(1,length(t));
Altitude = zeros(1,length(t));
Speed(1) = Intial_speed;
Altitude(1) = Intial_altitude;
r(1) = Intial_altitude+R;
m(1)= m_0;
r(1) = Intial_altitude+R;
rho(1) = 1.225;
%FOR LOOP
for i=2:length(t)
  %CALCULATION FOR dmdt BEING NOT CONSTANT THROUGHT  
    if m(i-1)> m_e
        dmdt = dmdt_f;
    else
       dmdt = dmdt_o;
    end
    
    m(i) = m(i-1) - Time_step*dmdt;             %Calculate the current speed
    
    rho(i)=1.225*10.^((-3*Altitude(i-1)/50000));      %Calculate the current air density
  
    Speed(1) = v_e + Speed(i-1) - v_e*m(i)/m(i-1) - Time_step*G*M/(r(i-1))^2 - Time_step*0.5*rho(i)*A*Cd*Speed(i-1)^2/m(i-1);  %Calculate the current Velocity
    
    r(i) = r(i-1) + Time_step*0.5*(Speed(i)+Speed(i-1));    %Calculate the current radius
    
    Altitude(i) = Altitude(i-1) + Time_step*0.5*(Speed(i)+Speed(i-1));  %Calculate the current altitude
    
end
  
%DISPLAY THE RESULTS
subplot(2,1,1)
plot (t,Speed,'m')
xlabel('Time(s)')
ylabel('Speed(m/s)')
subplot(2,1,2)
plot(t,Altitude)
xlabel('Time(s)')
ylabel('Altitude(m)')
This are my graphs
As a result for this question i need to get the following graphs but i don't know how to change my code so that i can get these graphs
I need help to fix my code so that i can get the following graph as a result

Best Answer

This
Speed(1) = ...
needs to be this
Speed(i) = ...
and this
... v_e*m(i)/m(i-1) ...
should be this
... Time_step*v_e*dmdt/m(i-1) ...
Also, that V^2 term only works if the rocket is going up. If the rocket ever starts falling down that term will not work as is. You should change V^2 to abs(V)*V for this to work in that situation.
Related Question