2D asteroid problem

Question

I don't know how to solve the next problem, and if anyone could explain to me step by step how it is solved I would thank you a lot. I know it's not hard, but I'm not seeing how to do it. Thank you in advance!

An asteroid passes next to a planet with mass m, and position (0, 0). The asteroid has an initial velocity (Vx, Vy) and a starting position (Px, Py). Considering that gravity is 1, and ignoring the size of both objects, get the formula to calculate the asteroid position in the future.

Answer 1

Here is the coordinate system and the force diagram. At first, it's a good idea to forget about any initial velocities and just focus on the forces. There's only one force acting, which is the gravity on the asteroid. According to Newton's law of universal gravitation, it can be calculated as $$ \vec{F_g} = G \frac{m_{\text{asteroid}} m_{\text{planet}}}{r^2} \vec{r} $$ where $G=1$ is a constant (set to $1$ by the problem statement) and $\vec{r}$ is a unit vector pointing from the asteroid to the planet (overlapping with the blue vector in the picture). $r$ is just the distance between the two objects, so $r^2 = x ^2 + y^2$.

Next, we want to split this into $x$ and $y$ components. The angle that the blue vector makes (let's call it $\theta$) with the horizontal satisfies $$ \tan \theta = \frac{x}{y} $$ Now, the components of the force are $$ \begin{cases} x\text{-component}: F_x = -\sin \theta |\vec{F_g}| = -\sin \theta \frac{m_{\text{asteroid}} m_{\text{planet}}}{x ^2 + y^2} \\ y\text{-component}: F_x = -\cos \theta |\vec{F_g}| = -\cos \theta \frac{m_{\text{asteroid}} m_{\text{planet}}}{x ^2 + y^2} \\ \end{cases} $$ To wrap up the analysis with the forces, we apply Newton's second law of motion to the asteroid: $\sum \vec{F} = m_{\text{asteroid}} \vec{a}$. Since gravity is the only force, the previous expressions give the $x$ and $y$-components of the acceleration: $$ \begin{cases} x\text{-component}:-\sin \theta \frac{m_{\text{asteroid}} m_{\text{planet}}}{x ^2 + y^2} = m_{\text{asteroid}} a_x \\ y\text{-component}:-\cos \theta \frac{m_{\text{asteroid}} m_{\text{planet}}}{x ^2 + y^2} = m_{\text{asteroid}} a_y \\ \end{cases} $$ Here, $a_x$ and $a_y$ are the components of acceleration. The mass of the asteroid gets cancelled,and we set $m_{\text{planet}}=m$, and now we have the acceleration components $$ \begin{cases} a_x = \frac{d^2 x}{dt^2} = -\frac{m \sin \theta }{x ^2 + y^2} \\ a_y = \frac{d^2 y}{dt^2} = - \frac{m \cos \theta}{x ^2 + y^2} \\ \end{cases} $$ But we're not done yet. We have to form the equations of motion. So far, we only have commented on the second derivative of position. We have to incorporate the initial conditions as well (we know the initial values of $P_x$ and $P_y$, and also the initial values of the first derivative of these). So now we have the coupled second-order differential equations $$ \begin{cases} \frac{d^2 x}{dt^2} = -\frac{m \sin \theta }{x ^2 + y^2} \\ \frac{d^2 y}{dt^2} = - \frac{m \cos \theta}{x ^2 + y^2} \\ \end{cases} \qquad \text{with} \quad \begin{cases} \frac{d x}{dt} (0) = V_x \\ \frac{d y}{dt} (0) = V_y \\ \end{cases} \quad \text{and} \quad \begin{cases} x (0) = P_x \\ y(0) = P_y \\ \end{cases} $$ This can be solved numerically with a computer. When implementing the solution algorithm to a computer, one must be a bit careful with that, but there are several methods. The simplest is Euler's method, but there you have to choose a really small time step. You would have to make hundreds if not thousands of steps to have a decent accuracy. Another way to solve would be some version of Runge-Kutta.

---

Here's a quickly put together Python3 code that could serve as an example on how to numerically solve this problem:

import math
import sys
m = 10.0
delta_t = 0.001 # Don't put this to zero!
V_x = 10.0
V_y = -3.0
P_x = 2.0
P_y = 3.0
T_max = 500
total_iterations = int(T_max / delta_t) + 1
def calculateAcceleration(x,y):
    global m
    theta = math.atan2(x,y)
    acceleration_x = -mmath.sin(theta)/(x2 + y2)
    acceleration_y = -mmath.cos(theta)/(x2 + y2)
    return [acceleration_x,acceleration_y]
def iterate(x,y,velocity_x, velocity_y):
    global delta_t
    acceleration_vector = calculateAcceleration(x,y)
    acceleration_x = acceleration_vector[0]
    acceleration_y = acceleration_vector[1]
    # Using first-order Euler method
    velocity_x_new = acceleration_x + delta_tacceleration_x
    velocity_y_new = acceleration_y + delta_tacceleration_y
    x_new = x + velocity_x_newdelta_t
    y_new = y + velocity_y_newdelta_t
return [x_new,y_new,velocity_x_new,velocity_y_new]


def fly():
    global V_x, V_y
    global P_x, P_y
    global delta_t
    global total_iterations
iteration_number = 0
time = 0

# Initial values
x = P_x
y = P_y
velocity_x = V_x
velocity_y = V_y

while iteration_number < total_iterations:
    iteration_number += 1
    new_parameters = iterate(x,y,velocity_x, velocity_y)
    x       = new_parameters[0]
    y       = new_parameters[1]
    velocity_x  = new_parameters[2]
    velocity_y  = new_parameters[3]
    print(x,y)
    input() 

fly()