Runge-Kutta time step effects solution in matlab simulation

Question

I am simulating a rocket launch using Runge-Kutta. I am trying to compute the max altitude of the rocket. When I use a different timestep I get different results.

$dt = 0.1$ gives $altitude = 1.0190e6$m

$dt = 0.01$ gives $altitude = 360020$ m

I am simulating a rocket launch from the equator and have initial velocities in the y and x direction. When I set those to zero, the results change very little as expected.

How should I choose the velocities.

Apollo.initialRocketMass = 750e3;                           % rocket mass
Apollo.initialFuelMass = 1840e3; %2150e3;                     % fuel mass
Apollo.totalMass = Apollo.initialRocketMass...
                    + Apollo.initialFuelMass;               % total mass
Apollo.burnRate = 15e3;                                     % burn rate
Apollo.initialThrust = 34e6*1.25;                                % initial thrust


dt = 0.001;               % time step
tstart = 0;             % start time
tstop = 1000;           % upper limit on how long to simulate
A = []; V = []; H = []; T = []; M = []; F = []; % empty vectors
k = 1;

y0 = 6.371e6; 
% v = 0; vx =0; vy = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vx =466*dt; vy = 29722*dt; v = sqrt((vx)^2+(vy)^2); % I SUSPECT THE PROBLEM IS HERE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x = 0; y = y0;          % initial height

M(1) = Apollo.totalMass;    % enter the loop with total mass    
F(1) = Apollo.initialThrust;
Y(1) = y;                   % enter the loop with initial position
m = M(1);
ft = F(1);
w = [x,vx,y,vy];
for t = tstart:dt:tstop

    if sqrt(w(1)^2+w(3)^2) >= y0 %ground is y0

        % track variable mass
        if ( t <= 0 )
            m = Apollo.totalMass;
        elseif ( t > 0 && fuel == true )
            m = Apollo.totalMass - Apollo.burnRate*t;
        else
            m = Apollo.initialRocketMass;
        end

        % track fuel consumption
        if m > Apollo.initialRocketMass
            fuel = true;
        else
            fuel = false;
        end

        % track thrust
        if ( t >= 0 && Apollo.initialFuelMass >= Apollo.burnRate && fuel == true)
            ft = Apollo.initialThrust;
        else
            ft = 0;

        end

        [t, w] = GetAcceleration_8rk( t,dt,Apollo,fuel, w,m,ft,v); % get current acceleration

        v = sqrt((w(2))^2+(w(4))^2); % save speed


    end
end

Answer 1

That's not really hard, and your suspicion is correct.

For $dt = 0.01$, the starting $V_y$ is $297.22$, while for $dt = 0.1$ your starting $V_y$ is $2972.2$. With $10$ times greater initial velocity, not only do you get an intrinsic altitude increase, the rocket gets to higher altitude more quickly, and gravity falls off faster.

For uniformly accelerated motion, of course, $$s = vt + \frac{at^2}{2} $$

The fix is to simply specify the starting velocities without referring to dt.

I'm actually puzzled why your altitudes are so small. You should be seeing a difference on the order of $$\Delta s = \Delta v\times t = 2675\times 1000 = 2.675\times10^6 $$ which is nearly 4 times larger than your stated numbers.