Orbit 9-12-2020

I solved the problem using a Lagrangian approach. The Newtonian approach would be less straightforward and possibly harder. With this, I just yielded a system of equations, and numerically solved.

In Hossam's solution, a first order equation was produced using a lot of mathematical manipulation. In the end, Hossam's equation of motion is far better behaved and easier to numerically solve than the coupled system below.

With the system, there are 2 generalised coordinates, hence two coupled differential equations will result.

Later, I might try using generalised polar coordinates $r$ and $\theta$ , and solve the Euler-Lagrange equation for them.

It may lead to a first order equation, actually.

So, the Euler-Lagrange equations for this system is given by:

$\displaystyle \frac{d}{dt} \left( \frac{\partial L}{\partial\dot{x}} \right)= \frac{\partial L}{\partial x}$

$\displaystyle \frac{d}{dt} \left( \frac{\partial L}{\partial\dot{y}} \right)= \frac{\partial L}{\partial y}$

The Lagrangian for the system can be found by finding the kinetic energy (the easy part) and the potential energy (the slightly harder part).

The gravitational potential energy can be found by:

$\displaystyle U_g = -\frac{GMm}{r}$

Where $r$ is the distance from the origin. This distance can be found with a bit of Pythagoras:

$r = \sqrt{x^2 + y^2}$

The Lagrangian:

$\displaystyle L = \frac{\dot{x}^2}{2} + \frac{\dot{y}^2}{2} -\frac{10}{\sqrt{x^2 + y^2}}$

Solving the Euler-Lagrange equation will give:

$\displaystyle \ddot{x} = -\frac{10x}{(x^2+y^2)^{\frac{3}{2}}}$

$\displaystyle \ddot{y} = -\frac{10y}{(x^2+y^2)^{\frac{3}{2}}}$

The system is kind of badly behaved, even when you solve with a small timestep. But a timestep of $10^{-6}$ is enough to give you a good result.

Here's code if you want to see how I solved the coupled equation @Lil Doug

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

import math

time = 0
deltaT = 10**-6

#initial values
x = 1
y = 0

xDot = -2
yDot = 2

while time <= 5:
  #the acceleration in x and y; coupled second order non-linear differential equations
  xDotDot = -10*x/(x**2+y**2)**(1.5)
  yDotDot = -10*y/(x**2+y**2)**(1.5)


  #Euler's method
  xDot += xDotDot*deltaT
  yDot += yDotDot*deltaT

  x += xDot*deltaT
  y += yDot*deltaT

  #updating time
  time += deltaT

print(math.hypot(x,y))

Using polar coordinates, we can write

$\vec{r}(t) = (x, y) = r ( \cos \theta , \sin \theta )$

Differentiating with respect to time $t$ ,

$\dot{\vec{r}} = \dot{r} (\cos \theta, \sin \theta) + r ( - \sin \theta, \cos \theta ) \dot{\theta}$

Differentiating again,

$\ddot{\vec{r}} = \ddot{r} (\cos \theta, \sin \theta )+ 2 \dot{r} (- \sin \theta , \cos \theta) \dot{\theta} + r ( - \cos \theta, - \sin \theta) (\dot{\theta})^2 + r (-\sin \theta , \cos \theta ) \ddot{\theta}$

The force acting on the particle is $\vec{F} = \dfrac{-10}{r^2} (\cos \theta , \sin \theta) = m \ddot{\vec{r}}$

Equating radial and traverse components of the vectors, on both sides of the equation, we deduce that,

$\ddot{r} - r (\dot{\theta})^2 = \dfrac{-10}{r^2} \hspace{20pt}(1)$

and

$2 \dot{r} \dot{\theta} + r \ddot{\theta} = 0 \hspace{20pt} (2)$

Multiplying the second equation by $r$ , we get,

$2 r \dot{r} \dot{\theta} + r^2 \ddot{\theta} = 0 \hspace{20pt} (3)$

The left hand side is just the derivative of $r^2 \dot{\theta}$ , therefore,

$r^2 \dot{\theta} = K \hspace{20pt}(4)$

To determine K, we use the initial conditions, we're given that $r(0) = 1 , \theta(0) = 0 , \dot{\vec{r} } = (-2,2)$ , Using the expression for velocity developed above, we have,

$\dot{\vec{r}} = \dot{r} (\cos \theta, \sin \theta) + r ( - \sin \theta, \cos \theta ) \dot{\theta}$

so that,

$(-2, 2) = \dot{r} (1, 0) + (1) (0, 1) \dot{\theta}$

Therefore, $\dot{r}(0) = -2$ , and $\dot{\theta}(0) = 2$ , and it follows from this that the constant $K = 2$ . Plugging in equation (4) into (1),

$\ddot{r} - \dfrac{4}{r^3} = -\dfrac{10}{r^2} \hspace{20pt}(5)$

By relating the time-derivative of $r$ to the $\theta$ -derivative, equation (5) becomes,

$\dfrac{4 r''}{ r^4} - \dfrac{8 r'^2}{ r^5} - \dfrac{4 }{ r^3} = -\dfrac{10}{r^2}\hspace{20pt}(6)$

where $r' = \dfrac{dr}{d\theta}$ and $r'' = \dfrac{d^2 r}{d\theta^2}$ .

Now, we will use the change of variables, $u = \dfrac{1}{r}$ . This will change equation $(6)$ into,

$4 u'' + 4 u - 10 = 0 \hspace{20pt}(7)$

where $u'' = \dfrac{d^2 u}{d\theta^2}$ . Equation $(7)$ is a linear differential equation with constant coefficients. Dividing by 4,

$u'' + u - \frac{5}{2} = 0 \hspace{20 pt} (8)$

The solution of this differential equation is:

$u = \frac{5}{2} + A \cos \theta + B \sin \theta$

Since $r(0) = 1$ , then , $u(0) = 1$ , and therefore, $A = - \dfrac{3}{2}$ , and one can show that $u(0) = 1$ , hence $B = 1$ . And therefore,

$r(\theta) = \dfrac{1}{\frac{5}{2} - \frac{3}{2} \cos \theta + \sin \theta } \hspace{20pt} (9)$

This is an equation of an ellipse. Now from equation (4), we have,

$\dfrac{d \theta}{dt} = \dot{\theta} = \dfrac{2}{r^2} = 2 (\frac{5}{2} - \frac{3}{2} \cos \theta + \sin \theta )^2$

And this is a first-order differential equation in $\theta$ , which we want to integrate over $[0, 5]$ . For this I used the Runge-Kutta 4th-order method, with a step-size of $h = 10^{-4}$ , I got $\theta(5) = 23.725$ radians. Plugging this value for $\theta$ into equation $(9)$ , gives $r = 0.7876$ .