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

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);

Answer 1

[Shashi Kiran]

Hi @Jhryssa,

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:

• ode45: https://www.mathworks.com/help/matlab/ref/ode45.html

Hope this helps.

[Jhryssa]

I appreciate the help! Thank you.

Answer 2

[Sam Chak]

Hi @Jhryssa

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')

[Sam Chak]

Hi @Jhryssa

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')

[Jhryssa]

Thank you so much!