0
$\begingroup$

I'm trying to simulate the orbits of Mars and Earth in MATLAB using the 4th order Runge-Kutta method for numerical integration. However, the orbits appear as straight lines rather than the expected ellipses.

I've used a variety of physical constants and the initial positions and velocities of the planets at perihelion. The orbits are calculated in a time loop, where the new positions and velocities are determined using the Runge-Kutta method.

However, the orbits still appear as straight lines on the plot. I have attempted to debug the code and have also tried to adjust the time step, but none of that helped. I believe the problem lies somewhere in the time loop or the Runge-Kutta method, but I can't find what I'm doing wrong.

Here is the MATLAB code:

% constants
G = 6.67408e-11; % m^3.kg^(-1).s^(-2)
M_s = 1.989e30;  % kg
m_e = 5.9722e24; % kg
m_m = 6.4171e23; % kg

% converting perihelion and semi-major axis from km to m
perihelion_e = 147100100e3;  % m
perihelion_m = 206650000e3;  % m
semi_major_e = 149597887e3;  % m
semi_major_m = 227939366e3;  % m

% initial velocity at perihelion
v_e = sqrt(G*M_s*(2/perihelion_e - 1/semi_major_e)); % m/s
v_m = sqrt(G*M_s*(2/perihelion_m - 1/semi_major_m)); % m/s

% martian year in seconds
martian_half_year = 1.88*365.25/2*24*60*60; % s

% time vector
dt = 500;  % Reduce the time step from 1000 seconds to 500 seconds
t = 0:dt:martian_half_year;

% initial conditions for x, y, vx and vy (martian perihelion has 45 deg angle)
initial_e = [perihelion_e, 0, 0, v_e];
initial_m = [perihelion_m*cosd(45), perihelion_m*sind(45), -v_m*sind(45), v_m*cosd(45)];

% pre-allocating position and velocity arrays
pos_e = zeros(length(t), 2);
pos_m = zeros(length(t), 2);
vel_e = zeros(length(t), 2);
vel_m = zeros(length(t), 2);

pos_e(1, :) = initial_e(1:2);
pos_m(1, :) = initial_m(1:2);
vel_e(1, :) = initial_e(3:4);
vel_m(1, :) = initial_m(3:4);

% the numerical integration method
runge_kutta = @(r, v, dt, m) deal(r + dt*v, v - dt*G*M_s*(r/norm(r)^3)/m);

% calculating the orbits
for i = 2:length(t)
    k1_r_e = dt * vel_e(i-1, :);
    k1_v_e = dt * (-G * M_s * (pos_e(i-1, :) / norm(pos_e(i-1, :))^3)) / m_e;
    k2_r_e = dt * (vel_e(i-1, :) + 0.5 * k1_v_e);
    k2_v_e = dt * (-G * M_s * ((pos_e(i-1, :) + 0.5 * k1_r_e) / norm(pos_e(i-1, :) + 0.5 * k1_r_e)^3)) / m_e;
    k3_r_e = dt * (vel_e(i-1, :) + 0.5 * k2_v_e);
    k3_v_e = dt * (-G * M_s * ((pos_e(i-1, :) + 0.5 * k2_r_e) / norm(pos_e(i-1, :) + 0.5 * k2_r_e)^3)) / m_e;
    k4_r_e = dt * (vel_e(i-1, :) + k3_v_e);
    k4_v_e = dt * (-G * M_s * ((pos_e(i-1, :) + k3_r_e) / norm(pos_e(i-1, :) + k3_r_e)^3)) / m_e;
    pos_e(i, :) = pos_e(i-1, :) + (k1_r_e + 2*k2_r_e + 2*k3_r_e + k4_r_e) / 6;
    vel_e(i, :) = vel_e(i-1, :) + (k1_v_e + 2*k2_v_e + 2*k3_v_e + k4_v_e) / 6;

    k1_r_m = dt * vel_m(i-1, :);
    k1_v_m = dt * (-G * M_s * (pos_m(i-1, :) / norm(pos_m(i-1, :))^3)) / m_m;
    k2_r_m = dt * (vel_m(i-1, :) + 0.5 * k1_v_m);
    k2_v_m = dt * (-G * M_s * ((pos_m(i-1, :) + 0.5 * k1_r_m) / norm(pos_m(i-1, :) + 0.5 * k1_r_m)^3)) / m_m;
    k3_r_m = dt * (vel_m(i-1, :) + 0.5 * k2_v_m);
    k3_v_m = dt * (-G * M_s * ((pos_m(i-1, :) + 0.5 * k2_r_m) / norm(pos_m(i-1, :) + 0.5 * k2_r_m)^3)) / m_m;
    k4_r_m = dt * (vel_m(i-1, :) + k3_v_m);
    k4_v_m = dt * (-G * M_s * ((pos_m(i-1, :) + k3_r_m) / norm(pos_m(i-1, :) + k3_r_m)^3)) / m_m;
    pos_m(i, :) = pos_m(i-1, :) + (k1_r_m + 2*k2_r_m + 2*k3_r_m + k4_r_m) / 6;
    vel_m(i, :) = vel_m(i-1, :) + (k1_v_m + 2*k2_v_m + 2*k3_v_m + k4_v_m) / 6;
end

% first 10 position vectors of Mars
disp(pos_m(1:10, :))

% plotting the orbits
hold on
plot(pos_e(:, 1), pos_e(:, 2), 'b', 'LineWidth', 1)  % Plot Earth's orbit
plot(pos_m(:, 1), pos_m(:, 2), 'r', 'LineWidth', 1)  % Plot Mars' orbit
title('Orbits of Mars and Earth')
xlabel('Distance (m)')
ylabel('Distance (m)')

% plotting symbols symbols
plot(pos_e(1, 1), pos_e(1, 2), 'bo')  % Start of Earth's orbit
plot(pos_e(end, 1), pos_e(end, 2), 'bs')  % End of Earth's orbit
plot(pos_m(1, 1), pos_m(1, 2), 'ro')  % Start of Mars' orbit
plot(pos_m(end, 1), pos_m(end, 2), 'rs')  % End of Mars' orbit

% range of the plot of both orbits
axis([-3e11, 3e11, -3e11, 3e11])

legend('Earth', 'Mars', 'Earth start', 'Earth end', 'Mars start', 'Mars end')
axis equal
grid on
hold off

Could someone help me understand why the orbits are not appearing as expected? Any insight into what might be going wrong in my calculations would be greatly appreciated.


The description of the problem:

enter image description here


The plot that I'm getting:

enter image description here

$\endgroup$
4
  • $\begingroup$ Without a deep look into your code, one thing that immediately comes to mind is that you are using SI units. Basically immediately you are going to fail because computers cannot handle that kind of large or small numbers properly. You are always supposed to compute with the units natural to a problem, so that every number you keep track of, is kept near 1 in size. $\endgroup$
    – naturallyInconsistent
    Commented Apr 29, 2023 at 16:33
  • $\begingroup$ In addition to the previous advice, it is possible to reduce the Kepler problem to a single first order differential equation for the angle. The radial coordinate simply follows from the equation for an ellipse. Info can be found in wiki. $\endgroup$
    – secavara
    Commented Apr 29, 2023 at 16:55
  • 2
    $\begingroup$ Would Computational Science be a better home for this question? $\endgroup$
    – Qmechanic
    Commented Apr 29, 2023 at 17:21
  • 1
    $\begingroup$ what might be going wrong in my calculations Implementation details of computational tasks are explicitly off-topic: “While computational physics is on topic, we are not a programming site. If your question is about implementing computational code - in particular, if it's about writing, compiling, debugging or optimizing code, or about a specific language or library - then it is off topic.” $\endgroup$
    – Ghoster
    Commented Apr 30, 2023 at 3:25

1 Answer 1

3
$\begingroup$

You accidentally divided by the planet's mass twice.

The force exerted by the sun on the planet is $$ \vec{F} = -\frac{G M_s m_p}{\|r\|^3} \vec{r} $$ Plugging that in to the equation of motion for the planet, $$ \vec{F} = m_p \vec{a}\\ \vec{a} = -\frac{G M_s}{\|r\|^3} \vec{r} $$ Notice that the mass of the planet cancels out and is not in the numerator or denominator.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.