How to implement this coupled ODE
1 view (last 30 days)
Show older comments
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1224097/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1224102/image.png)
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
0 Comments
Answers (2)
Torsten
on 8 Dec 2022
Edited: Torsten
on 8 Dec 2022
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
1 Comment
Sam Chak
on 8 Dec 2022
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1224577/image.png)
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
Jan
on 8 Dec 2022
[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".
0 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!