How to implement this coupled ODE

1 次查看（过去 30 天）
Romio
Romio 2022-12-8
评论： Sam Chak 2022-12-8
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

请先登录，再回答此问题。

回答（2 个）

Torsten
Torsten 2022-12-8
编辑：Torsten 2022-12-8
Ran in:
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 个评论
Sam Chak
Sam Chak 2022-12-8
Ran in:
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

请先登录，再进行评论。

Jan
Jan 2022-12-8
[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".

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by