§AERE 351 Homework 2

Sep 25, 2025

§1.

Given:

\overset{...}{x} = -4 \ddot y + 3 \dot x - 5 x + 2 y + 8 \ddot x - 2 \dot y

\overset{...}{y} = 5 x - 7 \dot y + 2 \ddot x

Definitions:

x_1 = x, \quad x_2 = \dot x = \dot x_1, \quad x_3 = \ddot x = \dot x_2 \\ y_1 = y, \quad y_2 = \dot y = \dot y_1, \quad y_3 = \ddot y = \dot y_2

Substitution:

\dot x_3 = -4 y_3 + 3 x_2 - 5 x_1 + 2 y_1 + 8 x_3 - 2 y_2

\dot y_3 = 5 x_1 - 7 y_2 + 2 x_3

Definition of

z

involving

x

and

y

z = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ y_1 \\ y_2 \\ y_3 \end{bmatrix}, \quad \dot z = \begin{bmatrix} \dot x_1 \\ \dot x_2 \\ \dot x_3 \\ \dot y_1 \\ \dot y_2 \\ \dot y_3 \end{bmatrix}

Identities re-written:

\dot x_1 = x_2

\dot x_2 = x_3

\dot x_3 = - 5 x_1 + 3 x_2 + 8 x_3 + 2 y_1 - 2 y_2 - 4 y_3

\dot y_1 = y_2

\dot y_2 = y_3

\dot y_3 = 5 x_1 + 2 x_3 - 7 y_2

Everything put together in a matrix equation:

\begin{bmatrix} \dot x_1 \\ \dot x_2 \\ \dot x_3 \\ \dot y_1 \\ \dot y_2 \\ \dot y_3 \end{bmatrix} = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ -5 & 3 & 8 & 2 & -2 & -4 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 5 & 0 & 2 & 0 & -7 & 0 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ y_1 \\ y_2 \\ y_3 \end{bmatrix}

Hence:

\boxed{\dot z = A z}

Where:

\boxed{A = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ -5 & 3 & 8 & 2 & -2 & -4 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 5 & 0 & 2 & 0 & -7 & 0 \end{bmatrix}}

Given:

\ddot r = \frac{\mu}{r^3} r

The flattening:

x_1 = x, \quad x_2 = \dot x = \dot x_1, \quad x_3 = \ddot x = \dot x_2 \\ y_1 = y, \quad y_2 = \dot y = \dot y_1, \quad y_3 = \ddot y = \dot y_2 \\ z_1 = z, \quad z_2 = \dot z = \dot z_1, \quad z_3 = \ddot z = \dot z_2

The contents of

r

r = \begin{bmatrix} x \\ y \\ z \end{bmatrix}

Packing

r

into a column vector:

w = \begin{bmatrix} r \\ \dot r \end{bmatrix} \implies \begin{bmatrix} x \\ y \\ z \\ \dot x \\ \dot y \\ \dot z \end{bmatrix}

\dot w = \begin{bmatrix} \dot r \\ \ddot r \end{bmatrix} \implies \begin{bmatrix} \dot x \\ \dot y \\ \dot z \\ \ddot x \\ \ddot y \\ \ddot z \end{bmatrix} = \begin{bmatrix} x_2 \\ y_2 \\ z_2 \\ x_3 \\ y_3 \\ z_3 \end{bmatrix} = \begin{bmatrix} x_2 \\ y_2 \\ z_2 \\ \frac{\mu}{r^3} x_1 \\ \frac{\mu}{r^3} y_1 \\ \frac{\mu}{r^3} z_1 \end{bmatrix}

Resolving

r^3

r^3 = |r|^3 = \left( \sqrt{x^2 + y^2 + z^2} \right)^3 = (x^2 + y^2 + z^2)^{3/2}

MATLAB code:

r_0 = [5000; 10000; 2100];
r_dot_0 = [-5; 2; 1.5];
w_0 = [r_0; r_dot_0];

T = 4 * 60 * 60;

% the default settings for ode45 are too lax
options = odeset('RelTol', 1e-9, 'AbsTol', 1e-10);
[t, y] = ode45(@orbit, [0, T], w_0, options);

plot3(y(:, 1), y(:, 2), y(:, 3));
xlabel('x (km)');
ylabel('y (km)');
zlabel('z (km)');
grid on;

final_state = y(end, 1:3);
disp(final_state);

% ode45 gives the time in the first argument but my intellisense complains
% if I put a t there and don't use it so I have used ~ here
function dw_dt = orbit(~, w)
    % matlab complains if I declare this above like all other constants
    mu = 3.98 * 10 ^ 5;

    r = w(1:3);
    r_dot = w(4:6);
    % I am sure there's a built in function for this but I like component-wise
    r_cubed = (r(1) ^ 2 + r(2) ^ 2 + r(3) ^ 2) ^ (3/2);
    r_ddot = (-mu .* r) ./ r_cubed;

    dw_dt = [r_dot; r_ddot];
end

MATLAB • 881B

Plot:

Console output:

1.0e+03 *
-8.5505   -3.5230    0.4822

• 37B

§2.

m = 0.2kg

l = 3m

m_\text{beam} = 0

c = 0.06

g = 9.81m/s^2

Free body diagram:

Sum of forces:

\sum F_x = cv + mg \sin \theta = -ma

\frac{c}{m} v + g \sin \theta = -a

And assuming things stay in radians, dividing the velocity and acceleration of the tangential component by the length of the beam will result in the angular velocity and acceleration respectively:

\frac{c}{m} \dot \theta + \frac{g}{l} \sin \theta = -\ddot \theta

Rearranging stuff recovers the target equation:

\boxed{\ddot \theta + \frac{g}{l} \sin \theta + \frac{c}{m} \dot \theta = 0}

However, this will be more useful for the matrix form:

\ddot \theta = - \frac{g}{l} \sin \theta - \frac{c}{m} \dot \theta

Matrix form:

x_1 = \theta

x_2 = \dot \theta = \dot x_1

\dot x_2 = - \frac{g}{l} \sin x_1 - \frac{c}{m} x_2

\boxed{w = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}, \quad \dot w = \begin{bmatrix} x_2 \\ - \frac{g}{l} \sin x_1 - \frac{c}{m} x_2 \end{bmatrix}}

This can be plugged directly into MATLAB with ode45 and the following conditions:

t_0 = 0, \quad t_1 = 15s

\theta_0 = -60\degree

The code:

l = 3;

theta_0 = -60 * (pi / 180);
t_0 = 0;
t_1 = 15;

[t, w] = ode45(@pendulum, [t_0, t_1], [theta_0, 0]);

plot(t, w(:, 1) * l, 'b', 'DisplayName', 'Position (m)'); hold on;
plot(t, w(:, 2) * l, 'r', 'DisplayName', 'Velocity (m/s)');
legend('show');

xlabel('Time (s)');
ylabel('Position and Velocity (see respective units)');
grid on;

function dw_dt = pendulum(~, w)
    l = 3;
    g = 9.81;
    c = 0.06;
    m = 0.2;

    x_1_dot = w(2);
    x_2_dot = (-g / l) * sin(w(1)) - (c / m) * w(2);

    dw_dt = [x_1_dot; x_2_dot];
end

MATLAB • 534B

The plot:

Looks about right! Time to move onto Runge-Kutta of order 4. My updated script to handle Runge-Kutta of order 4 with the ode45 function both at once:

l = 3;

theta_0 = -60 * (pi / 180);
t_0 = 0;
t_1 = 15;
dt = 0.5;

[t, w] = ode45(@pendulum, [t_0, t_1], [theta_0, 0]);
[t_rk, w_rk] = rk4(@pendulum, [t_0, t_1], [theta_0, 0], dt);

plot(t, w(:, 1) * l, 'b', 'DisplayName', 'Position (m)'); hold on;
plot(t, w(:, 2) * l, 'r', 'DisplayName', 'Velocity (m/s)');

plot(t_rk, w_rk(:, 1) * l, '--b', 'DisplayName', 'RK4 Position (m)'); hold on;
plot(t_rk, w_rk(:, 2) * l, '--r', 'DisplayName', 'RK4 Velocity (m/s)');

legend('show');

xlabel('Time (s)');
ylabel('Position and Velocity (see respective units)');
grid on;

function [t, w] = rk4(f, t_span, w0, dt)
    tn = t_span(1);
    wn = w0(:)';

    % I initially did a for look until the dt's added up to t_span(2)
    % but that gave the w and t arrays mismatching sizes
    N = ceil((t_span(2) - t_span(1)) / dt);
    w = zeros(N + 1, length(wn));
    t = zeros(N + 1, 1);

    t(1) = tn;
    w(1, :) = wn;

    for i = 1:N
        % transposing every time here is really bad for performance
        % but it will suffice for a homework
        k1 = f(tn, wn).';
        k2 = f(tn + dt / 2, wn + (dt / 2) * k1).';
        k3 = f(tn + dt / 2, wn + (dt / 2) * k2).';
        k4 = f(tn + dt, wn + dt * k3).';

        wn = wn + (dt / 6) * (k1 + 2 * k2 + 2 * k3 + k4);
        tn = tn + dt;

        w(i + 1, :) = wn;
        t(i + 1) = tn;
    end

end

function dw_dt = pendulum(~, w)
    l = 3;
    g = 9.81;
    c = 0.06;
    m = 0.2;

    x_1_dot = w(2);
    x_2_dot = (-g / l) * sin(w(1)) - (c / m) * w(2);

    dw_dt = [x_1_dot; x_2_dot];
end

MATLAB • 1545B

The plot:

§3.

Given:

m_S = 2000kg

r_S = 42200km

m_E = 5.9636 * 10^{24} kg

r_M = 384400km

\frac{g_E}{g_M} = 6

\frac{R_M}{R_E} = 0.27

From the Internet:

G = 6.67430 * 10^{-11} \frac{N * m^2}{kg^2}

R_E = 6378 km

The force between the satellite and the Earth can be calculated right away:

F = G \frac{m_1 m_2}{r^2} = G \frac{m_S m_E}{(r_S + R_E)^2}

F_{ES} = 6.67430 * 10^{-11} \frac{N * m^2}{kg^2} * \frac{2000kg * 5.9636 * 10^{24} kg}{(6378km + 42200km)^2} = \boxed{337.34N}

For the force between the Moon and Earth, I will need the mass of the Moon first. The force on an object of mass

m

at the surface of a planet of mass

M

and radius

R

should be:

F = G \frac{m M}{R^2}

Newton's second law:

F = ma = \cancel{m}g = G \frac{\cancel{m} M}{R^2}

g = \frac{GM}{R^2}

Applying this to the given ratio:

\frac{g_E}{g_M} = 6 = \frac{\cancel{G} m_E}{R_E^2} \frac{R_M^2}{\cancel{G} m_M} = \frac{m_E}{m_M} \frac{R_M^2}{R_E^2}

The mass of the Moon can be solved for now:

m_M = m_E \frac{g_M}{g_E} \frac{R_M^2}{R_E^2}

m_M = 5.9636 * 10^{24} kg * \frac{1}{6} * 0.27^2 = 7.246 * 10^{22} kg

Now, the force between the two:

F = G \frac{m_E m_M}{r_M^2}

F_{EM} = 6.67430 * 10^{-11} \frac{N * m^2}{kg^2} * \frac{5.9636 * 10^{24} kg * 7.246 * 10^{22} kg}{(384400km)^2} = \boxed{1.9518 * 10^{20} N}

The force on the satellite from the Moon:

F_{MS} = G \frac{m_M m_S}{(r_M - r_S)^2}

F_{MS} = 6.67430 * 10^{-11} \frac{N * m^2}{kg^2} * \frac{7.246 * 10^{22} kg * 2000kg}{(384400km - 42200km)^2} = 0.0826N

And as for the percentage comparing the Moon's effect on the satellite vs. the Earth's:

\gamma_{EM} = \frac{F_{MS}}{F_{ES}} * 100\% = \frac{0.0826N}{337.34N} * 100\% = \boxed{0.0245\%}

Safe to say, the Moon can be ignored in this orbit for low fidelity calculations. Now it's time to evaluate the Sun's effect on the moon vs Earth's. Given:

m_{Su} = 1.9891 * 10^{30} kg

From the Internet:

r_{Su} = 149.6 * 10^6 km

Because of how far away the Sun is, the Earth's radius is practically the average radius of the moon around the Sun:

F_{SuM} = G \frac{m_{Su} m_M}{r_{Su}^2}

F_{SuM} = 6.67430 * 10^{-11} \frac{N * m^2}{kg^2} * \frac{1.9891 * 10^{30} kg * 7.246 * 10^{22} kg}{(149.6 * 10^6 km)^2} = \boxed{4.298 * 10^{20} N}

\gamma_{SuE} = \frac{F_{SuM}}{F_{EM}} = \frac{4.298 * 10^20 N}{1.9518 * 10^{20} N} = \boxed{2.2 = 220\%}

This is significant. I was not anticipating this! I was expecting it to be less than 1%.

§4.

The force between the two upon each other:

F_\text{gravity} = G \frac{m m}{d^2} = G \frac{m^2}{d^2}

There exists an imaginary point in the middle about which both the object rotate with a radius

r

r = \frac{1}{2}d

The centrifugal force on one of the objects orbiting this imaginary center:

F_\text{centrifugal} = m \omega^2 r = \frac{1}{2} m \omega^2 d

Both the forces are opposite and equal:

F_\text{gravity} - F_\text{centrifugal} = m a_\text{tangential}

Since the planets never leave the orbit, they will never have a tangential acceleration.

F_\text{gravity} - F_\text{centrifugal} = 0 \implies F_\text{centrifugal} = F_\text{gravity}

\frac{1}{2} m \omega^2 d = G \frac{m^2}{d^2}

The question asks for the angular velocity:

\omega^2 = 2 G \frac{m^2}{d^2 m d} = 2 G \frac{m}{d^3}

\boxed{\omega = \sqrt{\frac{2 G m}{d^3}}}

§5.

This problem seems to be identical to the one before, but just with trigonometry now. The contribution of object

1

on object

3

F_1 = G \frac{m^2}{d_0^2} \begin{bmatrix} - \sin 30 \degree \\ - \cos 30 \degree \end{bmatrix}

Object

2

applies a similar force on object

3

F_2 = G \frac{m^2}{d_0^2} \begin{bmatrix} \sin 30 \degree \\ - \cos 30 \degree \end{bmatrix}

The sum of the forces on object

3

F = F_1 + F_2 = G \frac{m^2}{d_0^2} \begin{bmatrix} 0 \\ - 2 \cos 30 \degree \end{bmatrix}

For simplicity sake, since object

3

is vertically aligned with the center of mass and the only component of the force is in the

y

axis, it's okay to omit the

x

component and just worry about the happenings of the

y

components now:

F_\text{gravity} = -2 G \frac{m^2}{d_0^2} \cos 30 \degree

Meanwhile, the centrifugal force is similar to what I used in the previous problem:

F_\text{centrifugal} = m \omega^2 r = \frac{2}{3} m \omega^2 h

The positive sign here means up. The sum of the forces is

0

again for the same reason as before, but this time the signs are built into the forces so I haven't added them again here:

F_\text{gravity} + F_\text{centrifugal} = 0

-2 G \frac{m^2}{d_0^2} \cos 30 \degree + \frac{2}{3} m \omega^2 h = 0

\frac{2}{3} m \omega^2 h = 2 G \frac{m^2}{d_0^2} \cos 30 \degree

\omega^2 = 3 G \frac{m}{d_0^2 h} \cos 30 \degree

\omega = \sqrt{3 G \frac{m}{d_0^2 h} \cos 30 \degree}

Oh, and I almost forgot, the

\cos

is analytically simplifiable:

\boxed{\omega = \sqrt{\frac{3 \sqrt{3}}{2} G \frac{m}{d_0^2 h}}}