Matlab code does not run past a certain point: i.e. buffers at a certain line of code

3 visualizaciones (últimos 30 días)
So I'm trying to solve a problem (attached below) and I've written the code (below) to solve it. When I go to run the code it gets stuck at the following line:
% Solve the ODEs numerically
%[t, y] = ode45(odeFunc, [t0, 1.66*3600], [r0; v0], options);
It fully ran once, but continues to not run.
This is the problem:
Code:
G = 6.6742e-11; % universal gravitational constant (km^3/kg/s^2)
M = 5.974e24; % Earth mass estimate (kg)
R = 6378; % planet radius (km)
r0 = [3207 5459 2714]; % initial position vector (km)
v0 = [-6.532 0.7835 6.142]; % initial velocity vector (km/s)
t0 = 0; % initial time in seconds
% Define the ODEs
odeFunc = @(t, y) [y(4); y(5); y(6); -G*M*y(1)/(norm(y(1:3))^3); -G*M*y(2)/(norm(y(1:3))^3); -G*M*y(3)/(norm(y(1:3))^3)];
% Set options for ode45 solver
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
% Solve the ODEs numerically
[t, y] = ode45(odeFunc, [t0, 1.66*3600], [r0; v0], options);
% Calculate altitude above Earth's surface
altitudes = sqrt(sum(y(:, 1:3).^2, 2)) - R;
% Find the maximum altitude and the time at which it occurs
[maxAltitude, maxAltitudeIndex] = max(altitudes);
timeOfMaxAltitude = t(maxAltitudeIndex);
% Display the results
fprintf('Maximum Altitude: %.2f km\n', maxAltitude);
fprintf('Time of Maximum Altitude: %.2f hours\n', timeOfMaxAltitude / 3600);

Respuesta aceptada

Shashi Kiran
Shashi Kiran el 13 de Sept. de 2024
Editada: Shashi Kiran el 13 de Sept. de 2024
I understand that your code is getting stuck while solving the function using the ode45 solver.
To address the issue, please follow these suggestions:
  • Using the gravitational constant (G) and the Earth's mass (M) separately can lead to computational inefficiencies and precision errors. Instead, replace with the gravitational parameter for a more accurate and simplified calculation.
% Gravitational parameter and Earth's radius
mu = 398600; % Gravitational parameter in km^3/s^2
R = 6378; % Earth's radius in km
% Initial conditions
r0 = [3207, 5459, 2714]; % Initial position vector in km
v0 = [-6.532, 0.7835, 6.142]; % Initial velocity vector in km/s
t0 = 0; % Initial time in seconds
  • Modify the odefunc to use μ directly.
% Define the ODE function for orbital motion
odeFunc = @(t, y) [
y(4);
y(5);
y(6);
-mu * y(1) / (norm(y(1:3))^3);
-mu * y(2) / (norm(y(1:3))^3);
-mu * y(3) / (norm(y(1:3))^3)
];
% Set options for the ode45 solver
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
  • The initial conditions vectors were concatenated using a semicolon [r0; v0], causing the solver to misinterpret the input vector. Change [r0; v0] to [r0, v0] to ensure the initial state vector is correctly formatted as a single row, compatible with ode45.
% Solve the ODEs numerically over the specified time span
[t, y] = ode45(odeFunc, [t0, 1.66 * 3600], [r0, v0], options);
% Calculate altitude above Earth's surface
altitudes = sqrt(sum(y(:, 1:3).^2, 2)) - R;
% Find the maximum altitude and the time at which it occurs
[maxAltitude, maxAltitudeIndex] = max(altitudes);
timeOfMaxAltitude = t(maxAltitudeIndex);
Applying these changes will ensure that the implementation is accurate as demonstrated here:
% Display the results
fprintf('Maximum Altitude: %.2f km\n', maxAltitude);
Maximum Altitude: 9685.19 km
fprintf('Time of Maximum Altitude: %.2f hours\n', timeOfMaxAltitude / 3600);
Time of Maximum Altitude: 1.66 hours
Refer to the following documentation for more details about the function:
Hope this helps.

Más respuestas (1)

Sam Chak
Sam Chak el 13 de Sept. de 2024
If your professor is not a disciplinarian, then there is nothing wrong with using the standard gravitational parameter [m³/kg/s²] as a multiplicative term. However, your specified unit is incorrect because all other physical quantities of length in your code are specified in [km]. Since 1 km = 1000 m, you need to divide the product by .
Secondly, to simulate a system, it is recommended to formally describe the equations of motion within a function and then solve it using the ode45 solver. The defined output of the system should be the altitude, as you want to determine its maximum value.
Lastly, you should plot the altitude so that you can logically interpret the results visually. Your original simulation abruptly stops exactly at seconds, while the satellite's altitude is still increasing. You should allow the simulation to run for a longer period to observe whether the maximum has truly occurred.
%% Satellite's motion dynamics
function [dydt, alt] = satellite(t, y) % dydt = state derivatives, alt = altitude
dydt = zeros(6, 1);
% parameters
G = 6.6742e-11; % universal gravitational constant [m³/kg/s²]
M = 5.974e24; % Earth mass estimate (kg)
m = 1000; % satellite mass [kg]
R = 6378; % planet radius (km)
mu = G*(M + m)/(1000^3); % Standard gravitational parameter [km³/s²]
% Equations of motion
dydt(1) = y(4);
dydt(2) = y(5);
dydt(3) = y(6);
dydt(4) = - mu*y(1)/(norm(y(1:3))^3);
dydt(5) = - mu*y(2)/(norm(y(1:3))^3);
dydt(6) = - mu*y(3)/(norm(y(1:3))^3);
% altitude of satellite
alt = sqrt(sum(y(1:3).^2)) - R;
end
%% Call ode45 solver
r0 = [ 3207; 5459; 2714]; % initial position vector (km)
v0 = [-6.532; 0.7835; 6.142]; % initial velocity vector (km/s)
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);
[t, y] = ode45(@satellite, [0, 5*3600], [r0; v0], options);
%% Altitude
alt = zeros(1, numel(t)); % Pre-allocate for the altitude
for j = 1:numel(t)
[~, alt(j)] = satellite(t(j), y(j,:).');
end
%% Find the maximum altitude and the time at which it occurs
[maxAlt, maxAltIdx] = max(alt);
timeOfMaxAlt = t(maxAltIdx);
%% Display the results
fprintf('Maximum Altitude: %.2f km\n', maxAlt);
Maximum Altitude: 9675.90 km
fprintf('Time of Maximum Altitude: %.2f hours\n', timeOfMaxAlt/3600);
Time of Maximum Altitude: 1.70 hours
%% Plot the altitude
plot(t/3600, alt), grid on
xline(timeOfMaxAlt/3600, '--', sprintf('Max Time: %.4f hour', timeOfMaxAlt/3600), 'color', '#7F7F7F', 'LabelVerticalAlignment', 'bottom')
xlabel('Time [hour]'), ylabel('Altitude [km]')
title ('Altitude of Satellite')
  2 comentarios
Sam Chak
Sam Chak el 13 de Sept. de 2024
You can also visualize the satellite's orbit using the plot3() command.
%% Call ode45 solver
r0 = [ 3207; 5459; 2714]; % initial position vector (km)
v0 = [-6.532; 0.7835; 6.142]; % initial velocity vector (km/s)
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);
[t, y] = ode45(@satellite, [0, 5*3600], [r0; v0], options);
plot3(y(:,1), y(:,2), y(:,3)), grid on
xlabel('x'), ylabel('y'), zlabel('z')
title ('Satellite''s Orbit')

Iniciar sesión para comentar.

Categorías

Más información sobre Reference Applications en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by