ode45 graph issue

1 次查看（过去 30 天）
Takura Nyatsuro
Takura Nyatsuro 2023-3-21
回答： Steven Lord 2023-3-21
Ran in:
I am trying to make a code which uses ode45 functions to integrate 2 equations which will allow me to produce a graphs of velocity and altuitude against time.
The code is meant to simulate a 2 stage rocket which launches with a mass of 700kg , a thrust of 6500N and burns for 101.87s with the mass reducing at a rate of 2.4540kg/s before the first stage is jettisoned leaving a mass of 100kg, afterards the rocket has a coast phase of 10 seconds after the first boost phase where the thrust is zero, after the coast phase the mass of the rocket is reduced to 100kg and the second boost stage begins for 75 seconds in which the thrust is 1422.5N and the mass is reduced by 0.5kg/s until the end of the second boost stage.
After attempting to implement ode45 into my code i keep getting negative values for velocity and altitude which i assume means theres an error in my code logic somewhere. Below i have attached an image the differential equations used.
Thank you
% Define the thrust as a function of time
T = @(t) 6500*(t<=101.81) + 0*(t>101.81 & t<=111.81) + 1422.5*(t>111.81 & t<=186.81) + 1422.5*(t>186.81 & t<=261.81) + 0*(t>261.81); % N
% Define the mass as a function of time
m = @(t) 700 - 2.4540*t - 0*(t>101.81 & t<=111.81) - 600*(t>111.81 & t<=186.81) - 0.5*(t>186.81 & t<=261.81) - 37.5*(t>261.81); % kg
% Define the initial conditions
h0 = 0; % m
v0 = 0; % m/s
% Define the time span for the integration
tspan1 = [0, 101.81]; % s (first boost phase)
tspan2 = [101.81, 111.81]; % s (coast phase)
tspan3 = [111.81, 186.81]; % s (second boost phase)
% Integrate the system of equations using ode45 for the first boost phase
f1 = @(t, y) [y(2); T(t)/m(t) - 9.81];
[t1, y1] = ode45(f1, tspan1, [h0; v0]);
% Integrate the system of equations using ode45 for the coast phase
f2 = @(t, y) [y(2); -9.81];
[t2, y2] = ode45(f2, tspan2, [y1(end,1); y1(end,2)]);
% Integrate the system of equations using ode45 for the second boost phase
f3 = @(t, y) [y(2); T(t)/m(t) - 9.81];
[t3, y3] = ode45(f3, tspan3, [y2(end,1); y2(end,2)]);
% Integrate the system of equations using ode45 for
% Concatenate the time and state vectors
t = [t1; t2(2:end); t3(2:end)];
y = [y1; y2(2:end,:); y3(2:end,:)];
% Plot the altitude and velocity as functions of time
figure;
subplot(2,1,1);
plot(t, y(:,1));
xlabel('Time (s)');
ylabel('Altitude (m)');
subplot(2,1,2);
plot(t, y(:,2));
xlabel('Time (s)');
ylabel('Velocity (m/s)');

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

回答（1 个）

Steven Lord
Steven Lord 2023-3-21
Ran in:
Instead of defining your integrand function piecewise like this:
% Define the thrust as a function of time
T = @(t) 6500*(t<=101.81) + 0*(t>101.81 & t<=111.81) + 1422.5*(t>111.81 & t<=186.81) + 1422.5*(t>186.81 & t<=261.81) + 0*(t>261.81); % N
% Define the mass as a function of time
m = @(t) 700 - 2.4540*t - 0*(t>101.81 & t<=111.81) - 600*(t>111.81 & t<=186.81) - 0.5*(t>186.81 & t<=261.81) - 37.5*(t>261.81); % kg
I would define one function for each of your time spans and solve using the appropriate function. It would help you avoid problems like the one shown by this plot of your m function:
figure
fplot(m, [0 200])
title("m")
I assume your "- 600*(t>111.81 & t<=186.81)" term was intended to simulate the drop of the first stage, but the first stage doesn't weigh 600 pounds.

类别

在 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!

Translated by