How to implement this coupled ODE

I want to solve this ODE but I don't know how to implement it. Basically I want du/dt to be some input function that I determine (e.g. heviside or constant function) and use its integral elsewhere, rho is obtained from some data fitting that is dependent on h = u * t.

Here is my try. I don't know how to define du/dt as an input, and at the same time use its solution to calcualte h and dm/dt

 

clc,clear,close all


global  Ve Cd  A_r g  

g = 9.8; % m/s^2
Cd = 0.47;
r = 3.7; % m
A_r = pi * r^2;
rho_fuel = 1000; % Kg/m^3
Ve = 3000; % m/s
Mbody = 60684; % Kg
Mfuel = 488370; % Kg
% Mfuel = Mfuel + 450000;
Mpayload = 31184; % Kg

t0 = 0;
Tf = 540;

m0 = Mbody + Mfuel + Mpayload;
a0 = 14;

% constant acceleration input
a_in = 14*ones(1,Tf);

tspan = [t0 Tf];

[t,m] = ode45(@(t,m) myode(t,m), tspan, m0);


function dydt = myode(t,y,a_in)
    global  Ve Cd  A_r g  
    %     dydt(2) = y(2);
    dydt(2) = interp1(a_in,t);
    H=[-1;0;1;2;3;4;5;6;7;8;9;10;11;13;15;17;20;25;30;32;35;40;45;47;50;51;60;70;71;80;84.9;89.7;100.4;105;110;120];
    H=H*1000;
    rho=[1.347;1.225;1.1116;1.0065;0.9091;0.8191;0.7361;0.6597;0.5895;0.5252;0.4664;0.4127;0.3639;0.2655;0.1937;0.1423;0.0880;0.0395;0.0180;0.0132;0.0082;0.0039;0.0019;0.0014;0.0010;0.00086;0.000288;0.000074;0.000064;0.000015;0.000007;0.000003;0.0000005;0.0000002;0.0000001;0.0000001];
    f_rho=fit(H, rho, 'linearinterp');
    h = y(2) * t; % what is y(2)??
    rho_air = f_rho(h);
    Fd = 0.5 * rho_air * (y(2))^2 * A_r * Cd; % what is y(2)??
    dydt(1) = (1/Ve) * (-m * g - Fd - m * dudt(2));

end

Answer 1

[Torsten]

clc,clear,close all


global  Ve Cd  A_r g  

g = 9.8; % m/s^2
Cd = 0.47;
r = 3.7; % m
A_r = pi * r^2;
rho_fuel = 1000; % Kg/m^3
Ve = 3000; % m/s
Mbody = 60684; % Kg
Mfuel = 488370; % Kg
% Mfuel = Mfuel + 450000;
Mpayload = 31184; % Kg

t0 = 0;
Tf = 540;

m0 = Mbody + Mfuel + Mpayload;
a0 = 14;

% constant acceleration input
a_in = 14*ones(1,Tf+1);

tspan = [t0 Tf];

H=[-1;0;1;2;3;4;5;6;7;8;9;10;11;13;15;17;20;25;30;32;35;40;45;47;50;51;60;70;71;80;84.9;89.7;100.4;105;110;120];
H=H*1000;
rho=[1.347;1.225;1.1116;1.0065;0.9091;0.8191;0.7361;0.6597;0.5895;0.5252;0.4664;0.4127;0.3639;0.2655;0.1937;0.1423;0.0880;0.0395;0.0180;0.0132;0.0082;0.0039;0.0019;0.0014;0.0010;0.00086;0.000288;0.000074;0.000064;0.000015;0.000007;0.000003;0.0000005;0.0000002;0.0000001;0.0000001];
f_rho=fit(H, rho, 'linearinterp');

[t,y] = ode45(@(t,m) myode(t,m,a_in,f_rho), tspan, [m0,a0]);
figure(1)
plot(t,y(:,1))

figure(2)
plot(t,y(:,2))

function dydt = myode(t,y,a_in,f_rho)
    global  Ve Cd  A_r g  
    m = y(1);
    u = y(2);
    dudt = interp1(0:540,a_in,t);
    h = u * t; 
    rho_air = f_rho(h);
    Fd = 0.5 * rho_air * u^2 * A_r * Cd; 
    dmdt = (1/Ve) * (-m * g - Fd - m * dudt);
    dydt = [dmdt;dudt];
end

[Sam Chak]

I think the final mass m should be Mbody + Mpayload, after the fuel is exhausted.

The velocity of the rocket when it runs out of fuel should be , not increasing continuously.

The altitude of the rocket is also unrealistically high.

I have not fixed the velocity. But I think it has something to do with the acceleration.

global  Ve Cd A_r g Tf Mbody Mpayload
g        = 9.8; % m/s^2
Cd       = 0.47;
r        = 3.7; % m
A_r      = pi * r^2;
rho_fuel = 1000; % Kg/m^3
Ve       = 3000; % m/s
Mbody    = 60684; % Kg
Mfuel    = 488370; % Kg
% Mfuel = Mfuel + 450000;
Mpayload = 31184; % Kg
t0       = 0;
Tf       = 1000;
m0       = Mbody + Mfuel + Mpayload;
a0       = 10;

% constant acceleration input
a_in   = a0*ones(1, Tf+1);

tspan  = [t0 Tf];
H      = [-1;0;1;2;3;4;5;6;7;8;9;10;11;13;15;17;20;25;30;32;35;40;45;47;50;51;60;70;71;80;84.9;89.7;100.4;105;110;120];
H      = H*1000;
rho    = [1.347;1.225;1.1116;1.0065;0.9091;0.8191;0.7361;0.6597;0.5895;0.5252;0.4664;0.4127;0.3639;0.2655;0.1937;0.1423;0.0880;0.0395;0.0180;0.0132;0.0082;0.0039;0.0019;0.0014;0.0010;0.00086;0.000288;0.000074;0.000064;0.000015;0.000007;0.000003;0.0000005;0.0000002;0.0000001;0.0000001];
f_rho  = fit(H, rho, 'linearinterp');
[t, y] = ode45(@(t, m) myode(t, m, a_in, f_rho), tspan, [m0, a0]);

figure(1)
plot(t, y(:,1)), grid on, xlabel('t'), title('Rocket Mass (kg)')

figure(2)
plot(t, y(:,2)/1e3), grid on, xlabel('t'), title('Vertical Velocity (km/s)')

figure(3)
plot(t, y(:,2).*t/1e3), grid on, xlabel('t'), title('Rocket Altitude (km)')

function dydt = myode(t, y, a_in, f_rho)
    global  Ve Cd A_r g Tf Mbody Mpayload
    m       = y(1) - Mbody - Mpayload;  % Final mass = Mbody + Mpayload
    u       = y(2);
    dudt    = interp1(0:Tf, a_in, t);
    h       = u * t; 
    rho_air = f_rho(h);
    Fd      = 0.5 * rho_air * u^2 * A_r * Cd; 
    dmdt    = (1/Ve) * (- m * g - Fd - m * dudt);
    dydt    = [dmdt; dudt];
end

Answer 2

[Jan]

[t,m] = ode45(@(t,m) myode(t,m), tspan, m0);

Here myode has 2 inputs only.

function dydt = myode(t,y,a_in)

Here it has 2 inputs. I assume you want to provide a_in also:

[t,m] = ode45(@(t,m) myode(t, m, a_in), tspan, m0);

Your equation has two variables. y(1) is related to m and y(2) to u, as far as I can see. But you do not integrate u with ODE45? Then it is not a part of the equation to be integrated.

Note, that ODE45 is designed to integrate smooth functions only. Linear interpolations are not smooth. If you are lucky the integration stops with an error message. In many cases ODE45 reduces the step size until the discontinuity is hidden by rounding errors, but the accumulated errors due to the huge number of steps can cause very random final values. Do not call this "result".