Runge-Kutta in rocketry - interconnected system of ODE's

Having a lot of trouble figuring out where my solver is going wrong here, solving for various metrics for a rocket's trajectory off the surface of the moon.

v is velocity, gamma is flight path angle (direction of rocket wrt surface - starts at 90 deg, straight up), and the dot versions are the derivatives wrt time.

According to the gammadot ODE, it should be zero while gamma is 90, but once I introduce the "kick" (at h = 100m), it should naturally increase, driving gamma towards horizontal, then decrease smoothly as gamma approaches 0 degrees at orbital velocity. The force of gravity, g, decreases as h increases, going by the gravconstant function. The mass of the rocket, m, decreases over time by its included function.

However, gammadot is decaying over time instead, back to zero from the instantaneous boost given by the kick. Just completely disobeying the ODE, and causing height to increase forever. I think it has to be the solver, though velocity curve looks right.

% Given Info

m0 = 4780; % kg, mass of rocket + fuel
mp = 2375; % kg, mass of fuel
ms = m0 - mp; % kg, mass of rocket
T_ascent = 15570; % N, thrust during ascent
Isp = 311; % sec, specific impulse (engine characteristic, only factors into rate of mass decrease)

g0_earth = 9.8067; % m/s^2, gravitational constant on surface of Earth (factors into Isp calculations)
g0_moon = 1.625; % m/s^2, gravitational constant on surface of Moon
R_moon = 1740000; % m, radius of Moon
v_p = 1712.313; % m/s, orbital velocity (thrust stops, point at which gamma SHOULD hit 0 deg)

mdot = T_ascent / (g0_earth * Isp); % constant rate of mass decrease, based on Earth gravity + Isp + thrust

% ODE's
vdot = @(time,gamma,g,m) (T_ascent / m) - g * sind(gamma);
gammadot = @(time,v,gamma,h,g) - ( (g / v) - (v / (R_moon + h)) ) * cosd(gamma);

% initialize time
step = 0.5;
time = 0 : step : 600;
heightcheck = 1;

% initialize variables
v = zeros( 1, length(time) ); % velocity, assumed in direction of gamma
v(1) = 0.000001;

gamma = zeros( 1, length(time) ); % flight path angle - angle b/w rocket and ground
gamma(1) = 90;

g = zeros( 1, length(time) ); % gravitational acceleration on the moon, decreases as h increases
g(1) = g0_moon;

vdottracker = zeros( 1, length(time) ); % acceleration tracker
gammadottracker = zeros( 1, length(time) ); % rate of change of gamma tracker

hdot = zeros( 1, length(time) ); % height, m
h = zeros( 1, length(time) );

xdot = zeros( 1, length(time) ); % horizontal distance, m
x = zeros( 1, length(time) );

m=zeros( 1, length(time) ); % mass of LM
for i = 1 : length(time)
    m(i) = mass( time(i) );
end

gamma_kick = -0.1; % degrees, to make gamma /= 90 so gammadot /= 0

% Runge Kutta

for i=1 : length(time) - 1
    
% solving for change in ODE variables

    k1vdot = vdot(   time(i),   gamma(i),   g(i),   mass( time(i) ));
    k1gammadot = gammadot(   time(i),   v(i),   gamma(i),   h(i),   gravconstant( h(i) ));
    
    k2vdot = vdot(   time(i) + step / 2,   gamma(i) + (step / 2) * k1gammadot,   g(i),   mass( time(i) + step / 2 ));
    k2gammadot = gammadot(   time(i) + step / 2,   v(i) + (step / 2) * k1vdot,   gamma(i) + (step / 2) * k1gammadot,   h(i),   gravconstant( h(i) ));

    k3vdot = vdot(   time(i) + step / 2,   gamma(i) + (step / 2) * k2gammadot,   g(i),   mass( time(i) + step / 2 ));
    k3gammadot = gammadot(   time(i) + step / 2, v(i) + (step / 2) * k2vdot,   gamma(i) + (step / 2) * k2gammadot,  h(i),  gravconstant( h(i) ));

    k4vdot = vdot(   time(i) + step,   gamma(i) + k3gammadot,   g(i),   mass( time(i) + step ));
    k4gammadot = gammadot(   time(i) + step,   v(i) + k3vdot,   gamma(i) + k3gammadot,   h(i),   gravconstant( h(i) ));

    k_vdot = (1/6) * (k1vdot + 2 * k2vdot + 2 * k3vdot + k4vdot);
    k_gammadot = (1/6) * (k1gammadot + 2 * k2gammadot + 2 * k3gammadot + k4gammadot);

% stepping forward in time

    v(i + 1) = v(i) + k_vdot;
    gamma(i + 1) = gamma(i) + k_gammadot;
    
    hdot(i + 1) = hdot(i) + v(i) * sind(gamma(i));
    h(i + 1) = h(i) + step * hdot(i);

    xdot(i + 1) = xdot(i) + v(i) * cosd(gamma(i));
    x(i + 1) = x(i) + step * xdot(i);

    g(i + 1) = gravconstant(h(i + 1));

% tracking ODE variables over time

    vdottracker(i) = (T_ascent / m(i)) - g(i) * sind(gamma(i));
    gammadottracker(i) = - (g(i) - ((v(i) ^ 2) / (R_moon + h(i)))) * cosd(gamma(i)) / v(i);
    
% checks

    if h(i) > 100  && heightcheck == 1 % checking for first instance of sufficient height 
        gamma(i+1) = gamma(i) + gamma_kick; % introducing kick to gamma to start decay to 0
        heightcheck = 0;
    end

    if m(i)==2405 || v(i)>v_p % checking if we're out of propellant or @ orbital velocity, either way the engine stops
        T_ascent = 0;
        vdot = @(time,gamma,g,m) (T_ascent / m) - g * sind(gamma); % eliminating thrust
    end
    
end


% plot to check curves
figure
title('gammadot vs time')
xlabel('time (s)')
ylabel('gammadot (degrees/s)')
plot(time,gammadottracker)

% The kick gives that instantaneous drop, and instead of continuing that curve down,
% it comes back up and tends towards zero!

% This is the equivalent of a pencil, balanced on its eraser, pushed a
% little with your finger - and then it just slows down and stabilizes in
% midair.

figure
plot(time, v)
title('velocity vs time')
xlabel('time (s)')
ylabel('velocity (m/s)')

% goes up to orbital velocity, thrust stops, velocity stays there. This looks fine.

% finding g as h increases
function g = gravconstant(h)
    g0_moon = 1.625; %m/s^2
    R_moon = 1740000; % m
    g = g0_moon / (1 + h / R_moon)^2;
end

% finding m as time increases
function m = mass(time)
    T_ascent = 15570; % N
    Isp = 311; % sec
    g0_earth = 9.8067; % m/s^2
    mdot = T_ascent / (g0_earth*Isp);

    m = 4780 - mdot*time;
    if m < 2405
        m = 2405; % out of propellant, left with mass of rocket
    end
end

Thanks very much in advance, I've been banging my head against this for hours.

[Torsten]

Didn't you forget to multiply the increment by the time step size ?

    v(i + 1) = v(i) + k_vdot;
    gamma(i + 1) = gamma(i) + k_gammadot;    

?

[Tanner]

You are absolutely right, good catch! I'd had the step at 1 so hadn't noticed. Unfortunately, this doesn't solve the overall issue with the solver.

[Torsten]

I'd suggest using a MATLAB integrator (e.g. ODE45) for your system of differential equations. If this gives the same results as your solver, you can at least be sure that the error is in your equations, not in the numerical method.

[Tanner]

After some finangling, it appears ode45 won't work for my problem. It would need h and x and g to be pre-defined arrays like m is, not dependent on the individual outcomes of the solver while it's working. Since they're all interconnected and not just dependent on the time at each point, it looks like I need the loop format, necessitating the more general Runge-Kutta. If I'm mistaken about that, let me know, I'll be glad to hear it.

[James Tursa]

@Tanner I don't uderstand this comment. ode45 does not require that the derivates of the states are dependent only on time. It is perfectly suited for derivatives to be dependent on current state and "interconnected". Can you explain your comment further? Maybe I am missing something.

[Torsten]

I agree with @James Tursa . Further I think that the expression for gammadot must be multiplied by 180/pi since it seems you want to write the equation for gamma in degrees instead of radians. To avoid these unit problems, I strongly recommend switching to radians for gamma instead.

[Tanner]

Maybe I just couldn't figure it out lol. But I switched to straight Euler's, then over to Excel using Euler's, and it's the same issue across all of them, including a change in syntax! Tells me there's something fundamental wrong with it, like equation-level.

I will try radians and get back to you.

[Tanner]

That is progress! Straight up changed all of the trig to radians and gamma0 to pi/2, and gammadot is still decaying instead of growing, but it's doing it slower. Unfortunately, I have no idea what that means.

[Sam Chak]

Hi @Tanner

Could you please ensure that all ODEs are placed inside the myRocketOnTheMoon function? Please verify that all constants are properly assigned. Additionally, check the mdot again as there should be a fixed rocket structure mass term, the remaining mass of the rocket itself, which must be significantly lower than its initial mass after the fuel is consumed.

%% Launching of a Rocket from the surface of the moon
function dstate = myRocketOnTheMoon(t, state)

    % definitions
    state(1) = m;
    state(2) = v;
    state(3) = gamma;
    state(4) = h;
    state(5) = x;

    % parameters
    m0       = 4780;        % kg, mass of rocket + fuel
    mp       = 2375;        % kg, mass of fuel
    ms       = m0 - mp;     % kg, mass of rocket
    T_ascent = 15570;       % N, thrust during ascent
    Isp      = 311;         % sec, specific impulse (engine characteristic, only factors into rate of mass decrease)
    g0_earth = 9.8067;      % m/s^2, gravitational constant on surface of Earth (factors into Isp calculations)
    g0_moon  = 1.625;       % m/s^2, gravitational constant on surface of Moon
    R_moon   = 1740000;     % m, radius of Moon
    v_p      = 1712.313;    % m/s, orbital velocity (thrust stops, point at which gamma SHOULD hit 0 deg)
    
    % ODEs
    mdot     = T_ascent/(g0_earth*Isp);                 % rate of change of rocket mass
    vdot     = T_ascent/m - g*sin(gamma);               % velocity
    gammadot =  - (g/v - v/(R_moon + h))*cos(gamma);    % rate of change of flight path angle
    hdot     = ...;                                     % rate of change of height
    xdot     = ...;                                     % rate of change of horizontal distance

    dstate   = [mdot
                vdot
                gammadot
                hdot
                xdot];
end

[Tanner]

Hi Sam,

That was my issue with ode45 - I couldn't figure out how to express those last two equations in anything other than a loop format, since they'd ideally be hdot(i+1) = hdot(i) + v(i) * sin(gamma(i)) and the equivalent for xdot. They build on the last value of themselves - how would I express that analytically in an ODE?

[Sam Chak]

Hi @Tanner

You may have misremembered the velocity components in projectile motion, which is a topic typically covered in introductory university physics. This is generally not addressed in rocket science or Tsiolkovsky rocket equation, as it is considered very basic. Always envision the physical laws of motion around you as a continuous flow of time, rather than as discrete snapshots of time.

https://www.vedantu.com/jee-main/physics-projectile-motion

[Torsten]

hdot = dh/dt = v*sin(gamma)
xdot = dx/dt = v*cos(gamma)

?

[Sam Chak]

Hi @Tanner

You already have the equations for hdot and xdot; however, in the code, they are expressed as snapshots of time. Please follow @Torsten.

Regarding the rocket mass, I believe you can remove the ODE function because it is directly solvable and independent of other states. Additionally, instead of creating a separate mass() function and using an if-else statement to represent the piecewise function, you can express the mass as a one-line mathematical function: , where is the dry mass.

T_ascent= 15570;    % N
Isp     = 311;      % sec
g0_earth= 9.8067;   % m/s^2
mdot    = T_ascent/(g0_earth*Isp);  % fuel mass consumption rate
m0      = 4780;     % initial mass (fuel + rocket structure)
m_dry   = 2405;     % dry mass (only rocket structure)
t       = linspace(0, 1000, 10001);

%% ----- put this into myRocketOnTheMoon() -----
m       = max(m0 - mdot*t, m_dry);
%% ---------------------------------------------

plot(t, m), grid on, title('Total rocket mass')
xlabel('Time / s')
ylabel('Mass / kg')

●

[Tanner]

Hey guys, thanks for all of the help - switching from degrees to radians for gamma did the trick. Problem wasn’t a mismatch between the units and the trig functions, or even the units of gamma. But GAMMADOT needed to be in radians for the equation to work correctly. Hate to see it, I’m sorry it ended up being so simple. You guys are awesome!