How to implement this coupled ODE
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(2 个)
Torsten
2022-12-8
编辑:Torsten
2022-12-8
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
2022-12-8
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
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".
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!